pysawsim/histogram.py

   1 # Copyright (C) 2009-2010  W. Trevor King <wking@drexel.edu>
   2 #
   3 # This program is free software: you can redistribute it and/or modify
   4 # it under the terms of the GNU General Public License as published by
   5 # the Free Software Foundation, either version 3 of the License, or
   6 # (at your option) any later version.
   7 #
   8 # This program is distributed in the hope that it will be useful,
   9 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 # GNU General Public License for more details.
  12 #
  13 # You should have received a copy of the GNU General Public License
  14 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
  15 #
  16 # The author may be contacted at <wking@drexel.edu> on the Internet, or
  17 # write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
  18 # Philadelphia PA 19104, USA.
  19
  20 """Histogram generation and comparison.
  21 """
  22
  23 import numpy
  24
  25 from . import log
  26
  27
  28 class Histogram (object):
  29     """A histogram with a flexible comparison method, `residual()`.
  30
  31     >>> h = Histogram()
  32     """
  33     def calculate_bin_edges(self, data, bin_width):
  34         """
  35         >>> h = Histogram()
  36         >>> h.calculate_bin_edges(numpy.array([-7.5, 18.2, 4]), 10)
  37         array([-10.,   0.,  10.,  20.])
  38         >>> h.calculate_bin_edges(numpy.array([-7.5, 18.2, 4, 20]), 10)
  39         array([-10.,   0.,  10.,  20.])
  40         >>> h.calculate_bin_edges(numpy.array([0, 18.2, 4, 20]), 10)
  41         array([  0.,  10.,  20.])
  42         >>> h.calculate_bin_edges(numpy.array([18.2, 4, 20]), 10)
  43         array([  0.,  10.,  20.])
  44         >>> h.calculate_bin_edges(numpy.array([18.2, 20]), 10)
  45         array([ 10.,  20.])
  46         """
  47         m = data.min()
  48         M = data.max()
  49         bin_start = m - m % bin_width
  50         return numpy.arange(bin_start, M+bin_width, bin_width, dtype=data.dtype)
  51
  52     def from_data(self, data, bin_edges):
  53         """Initialize from `data`.
  54
  55         All bins should be of equal width (so we can calculate which
  56         bin a data point belongs to).
  57
  58         `data` should be a numpy array.
  59         """
  60         self.headings = None
  61         self.bin_edges = bin_edges
  62         bin_width = self.bin_edges[1] - self.bin_edges[0]
  63
  64         bin_is = numpy.floor((data - self.bin_edges[0])/bin_width)
  65         self.counts = []
  66         for i in range(len(self.bin_edges)-1):
  67             self.counts.append(sum(bin_is == i))
  68         self.total = float(len(data)) # some data might be outside the bins
  69         self.mean = data.mean()
  70         self.std_dev = data.std()
  71         self.normalize()
  72
  73     def from_stream(self, stream):
  74         """Initialize from `stream`.
  75
  76         File format:
  77
  78             # comment and blank lines ignored
  79             <bin_edge><whitespace><count>
  80             ...
  81
  82         For example:
  83
  84
  85         `<bin_edge>` should mark the left-hand side of the bin, and
  86         all bins should be of equal width (so we know where the last
  87         one ends).
  88
  89         >>> import StringIO
  90         >>> h = Histogram()
  91         >>> h.from_stream(StringIO.StringIO('''# Force (N)\\tUnfolding events
  92         ... 150e-12\\t10
  93         ... 200e-12\\t40
  94         ... 250e-12\\t5
  95         ... '''))
  96         >>> h.headings
  97         ['Force (N)', 'Unfolding events']
  98         >>> h.total
  99         55.0
 100         >>> h.counts
 101         [10.0, 40.0, 5.0]
 102         >>> h.bin_edges  # doctest: +ELLIPSIS
 103         [1.5e-10, 2.000...e-10, 2.500...e-10, 3e-10]
 104         >>> h.probabilities  # doctest: +ELLIPSIS
 105         [0.181..., 0.727..., 0.0909...]
 106         """
 107         self.headings = None
 108         self.bin_edges = []
 109         self.counts = []
 110         for line in stream.readlines():
 111             line = line.strip()
 112             if len(line) == 0 or line.startswith('#'):
 113                 if self.headings == None and line.startswith('#'):
 114                     line = line[len('#'):]
 115                     self.headings = [x.strip() for x in line.split('\t')]
 116                 continue # ignore blank lines and comments
 117             try:
 118                 bin_edge,count = line.split()
 119             except ValueError:
 120                 log().error('Unable to parse histogram line: "%s"' % line)
 121                 raise
 122             self.bin_edges.append(float(bin_edge))
 123             self.counts.append(float(count))
 124         bin_width = self.bin_edges[1] - self.bin_edges[0]
 125         self.bin_edges.append(self.bin_edges[-1]+bin_width)
 126         self.total = float(sum(self.counts))
 127         self.mean = 0
 128         for bin,count in zip(self.bin_edges, self.counts):
 129             bin += bin_width / 2.0
 130             self.mean += bin * count
 131         self.mean /=  float(self.total)
 132         variance = 0
 133         for bin,count in zip(self.bin_edges, self.counts):
 134             bin += bin_width / 2.0
 135             variance += (bin - self.mean)**2 * count
 136         self.std_dev = numpy.sqrt(variance)
 137         self.normalize()
 138
 139     def to_stream(self, stream):
 140         """Write to `stream` with the same format as `from_stream()`.
 141
 142         >>> import sys
 143         >>> h = Histogram()
 144         >>> h.headings = ['Force (N)', 'Unfolding events']
 145         >>> h.bin_edges = [1.5e-10, 2e-10, 2.5e-10, 3e-10]
 146         >>> h.counts = [10, 40, 5]
 147         >>> h.to_stream(sys.stdout)
 148         ... # doctest: +NORMALIZE_WHITESPACE, +REPORT_UDIFF
 149         #Force (N)\tUnfolding events
 150         1.5e-10\t10
 151         2e-10\t40
 152         2.5e-10\t5
 153         """
 154         if self.headings != None:
 155             stream.write('#%s\n' % '\t'.join(self.headings))
 156         for bin,count in zip(self.bin_edges, self.counts):
 157             stream.write('%g\t%g\n' % (bin, count))
 158
 159     def normalize(self):
 160         self.total = float(self.total)
 161         self.probabilities = [count/self.total for count in self.counts]
 162
 163     def mean_residual(self, other):
 164         return abs(other.mean - self.mean)
 165
 166     def std_dev_residual(self, other):
 167         return abs(other.std_dev - self.std_dev)
 168
 169     def chi_squared_residual(self, other):
 170         assert self.bin_edges == other.bin_edges
 171         residual = 0
 172         for probA,probB in zip(self.probabilities, other.probabilities):
 173             residual += (probA-probB)**2 / probB
 174         return residual
 175
 176     def jensen_shannon_residual(self, other):
 177         assert self.bin_edges == other.bin_edges
 178         def d_KL_term(p,q):
 179             """
 180             Kullback-Leibler divergence for a single bin, with the
 181             exception that we return 0 if q==0, rather than
 182             exploding like d_KL should.  We can do this because
 183             the way d_KL_term is used in the Jensen-Shannon
 184             divergence, if q (really m) == 0, then p also == 0,
 185             and the Jensen-Shannon divergence for the entire term
 186             should be zero.
 187             """
 188             if p==0 or q==0 or p==q:
 189                 return 0.0
 190             return p * numpy.log2(p/q)
 191         residual = 0
 192         for probA,probB in zip(self.probabilities, other.probabilities):
 193             m = (probA+probB) / 2.0
 194             residual += 0.5*(d_KL_term(probA,m) + d_KL_term(probB,m))
 195         return residual
 196
 197     def _method_to_type(self, name):
 198         return name[:-len('_residual')].replace('_', '-')
 199
 200     def _type_to_method(self, name):
 201         return '%s_residual' % name.replace('-', '_')
 202
 203     def types(self):
 204         """Return a list of supported residual types.
 205         """
 206         return sorted([self._method_to_type(x)
 207                        for x in dir(self) if x.endswith('_residual')])
 208
 209     def residual(self, other, type='jensen-shannon'):
 210         """Compare this histogram with `other`.
 211
 212         Supported comparison `type`\s may be found with `types()`:
 213
 214         >>> h = Histogram()
 215         >>> h.types()
 216         ['chi-squared', 'jensen-shannon', 'mean', 'std-dev']
 217
 218         Selecting an invalid `type` raises an exception.
 219
 220         >>> h.residual(other=None, type='invalid-type')
 221         Traceback (most recent call last):
 222           ...
 223         AttributeError: 'Histogram' object has no attribute 'invalid_type_residual'
 224
 225         For residual types where there is a convention, this histogram
 226         is treated as the experimental histogram and the other
 227         histogram is treated as the theoretical one.
 228         """
 229         r_method = getattr(self, self._type_to_method(type))
 230         return r_method(other)