c0d73b75c0fdd3a8810af5ceaf8fdcb8135d9e29
[sawsim.git] / 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 _multiprocess_can_split_ = True
29 """Allow nosetests to split tests between processes.
30 """
31
32
33 class Histogram (object):
34     """A histogram with a flexible comparison method, `residual()`.
35
36     >>> h = Histogram()
37     """
38     def calculate_bin_edges(self, data, bin_width):
39         """
40         >>> h = Histogram()
41         >>> h.calculate_bin_edges(numpy.array([-7.5, 18.2, 4]), 10)
42         array([-10.,   0.,  10.,  20.])
43         >>> h.calculate_bin_edges(numpy.array([-7.5, 18.2, 4, 20]), 10)
44         array([-10.,   0.,  10.,  20.])
45         >>> h.calculate_bin_edges(numpy.array([0, 18.2, 4, 20]), 10)
46         array([  0.,  10.,  20.])
47         >>> h.calculate_bin_edges(numpy.array([18.2, 4, 20]), 10)
48         array([  0.,  10.,  20.])
49         >>> h.calculate_bin_edges(numpy.array([18.2, 20]), 10)
50         array([ 10.,  20.])
51         """
52         m = data.min()
53         M = data.max()
54         bin_start = m - m % bin_width
55         return numpy.arange(bin_start, M+bin_width, bin_width, dtype=data.dtype)
56
57     def from_data(self, data, bin_edges):
58         """Initialize from `data`.
59
60         All bins should be of equal width (so we can calculate which
61         bin a data point belongs to).
62         """
63         self.headings = None
64         data = numpy.array(data)
65         self.bin_edges = numpy.array(bin_edges)
66         bin_width = self.bin_edges[1] - self.bin_edges[0]
67
68         bin_is = numpy.floor((data - self.bin_edges[0])/bin_width)
69         self.counts = numpy.zeros((len(self.bin_edges)-1,), dtype=numpy.int)
70         for i in range(len(self.counts)):
71             self.counts[i] = (bin_is == i).sum()
72         self.counts = numpy.array(self.counts)
73         self.total = float(len(data)) # some data might be outside the bins
74         self.mean = data.mean()
75         self.std_dev = data.std()
76         self.normalize()
77
78     def from_stream(self, stream):
79         """Initialize from `stream`.
80
81         File format:
82
83             # comment and blank lines ignored
84             <bin_edge><whitespace><count>
85             ...
86
87         `<bin_edge>` should mark the left-hand side of the bin, and
88         all bins should be of equal width (so we know where the last
89         one ends).
90
91         >>> import StringIO
92         >>> h = Histogram()
93         >>> h.from_stream(StringIO.StringIO('''# Force (N)\\tUnfolding events
94         ... 150e-12\\t10
95         ... 200e-12\\t40
96         ... 250e-12\\t5
97         ... '''))
98         >>> h.headings
99         ['Force (N)', 'Unfolding events']
100         >>> h.total
101         55.0
102         >>> h.counts
103         [10.0, 40.0, 5.0]
104         >>> h.bin_edges  # doctest: +ELLIPSIS
105         [1.5e-10, 2.000...e-10, 2.500...e-10, 3e-10]
106         >>> h.probabilities  # doctest: +ELLIPSIS
107         [0.181..., 0.727..., 0.0909...]
108         """
109         self.headings = None
110         self.bin_edges = []
111         self.counts = []
112         for line in stream.readlines():
113             line = line.strip()
114             if len(line) == 0 or line.startswith('#'):
115                 if self.headings == None and line.startswith('#'):
116                     line = line[len('#'):]
117                     self.headings = [x.strip() for x in line.split('\t')]
118                 continue # ignore blank lines and comments
119             try:
120                 bin_edge,count = line.split()
121             except ValueError:
122                 log().error('Unable to parse histogram line: "%s"' % line)
123                 raise
124             self.bin_edges.append(float(bin_edge))
125             self.counts.append(float(count))
126         bin_width = self.bin_edges[1] - self.bin_edges[0]
127         self.bin_edges.append(self.bin_edges[-1]+bin_width)
128         self.analyze()
129
130     def to_stream(self, stream):
131         """Write to `stream` with the same format as `from_stream()`.
132
133         >>> import sys
134         >>> h = Histogram()
135         >>> h.headings = ['Force (N)', 'Unfolding events']
136         >>> h.bin_edges = [1.5e-10, 2e-10, 2.5e-10, 3e-10]
137         >>> h.counts = [10, 40, 5]
138         >>> h.to_stream(sys.stdout)
139         ... # doctest: +NORMALIZE_WHITESPACE, +REPORT_UDIFF
140         #Force (N)\tUnfolding events
141         1.5e-10\t10
142         2e-10\t40
143         2.5e-10\t5
144         """
145         if self.headings != None:
146             stream.write('#%s\n' % '\t'.join(self.headings))
147         for bin,count in zip(self.bin_edges, self.counts):
148             stream.write('%g\t%g\n' % (bin, count))
149
150     def analyze(self):
151         """Analyze `.counts` and `.bin_edges` if the raw data is missing.
152
153         Generate `.total`, `.mean`, and `.std_dev`, and run
154         `.normalize()`.
155         """
156         bin_width = self.bin_edges[1] - self.bin_edges[0]
157         self.total = float(sum(self.counts))
158         self.mean = 0
159         for bin,count in zip(self.bin_edges, self.counts):
160             bin += bin_width / 2.0
161             self.mean += bin * count
162         self.mean /=  float(self.total)
163         variance = 0
164         for bin,count in zip(self.bin_edges, self.counts):
165             bin += bin_width / 2.0
166             variance += (bin - self.mean)**2 * count
167         self.std_dev = numpy.sqrt(variance)
168         self.normalize()
169
170     def normalize(self):
171         """Generate `.probabilities` from `.total` and `.counts`.
172         """
173         self.total = float(self.total)
174         self.probabilities = [count/self.total for count in self.counts]
175
176     def mean_residual(self, other):
177         return abs(other.mean - self.mean)
178
179     def std_dev_residual(self, other):
180         return abs(other.std_dev - self.std_dev)
181
182     def chi_squared_residual(self, other):
183         assert (self.bin_edges == other.bin_edges).all()
184         residual = 0
185         for probA,probB in zip(self.probabilities, other.probabilities):
186             residual += (probA-probB)**2 / probB
187         return residual
188
189     def jensen_shannon_residual(self, other):
190         assert (self.bin_edges == other.bin_edges).all()
191         def d_KL_term(p,q):
192             """
193             Kullback-Leibler divergence for a single bin, with the
194             exception that we return 0 if q==0, rather than
195             exploding like d_KL should.  We can do this because
196             the way d_KL_term is used in the Jensen-Shannon
197             divergence, if q (really m) == 0, then p also == 0,
198             and the Jensen-Shannon divergence for the entire term
199             should be zero.
200             """
201             if p==0 or q==0 or p==q:
202                 return 0.0
203             return p * numpy.log2(p/q)
204         residual = 0
205         for probA,probB in zip(self.probabilities, other.probabilities):
206             m = (probA+probB) / 2.0
207             residual += 0.5*(d_KL_term(probA,m) + d_KL_term(probB,m))
208         return residual
209
210     def _method_to_type(self, name):
211         return name[:-len('_residual')].replace('_', '-')
212
213     def _type_to_method(self, name):
214         return '%s_residual' % name.replace('-', '_')
215
216     def types(self):
217         """Return a list of supported residual types.
218         """
219         return sorted([self._method_to_type(x)
220                        for x in dir(self) if x.endswith('_residual')])
221
222     def residual(self, other, type='jensen-shannon'):
223         """Compare this histogram with `other`.
224
225         Supported comparison `type`\s may be found with `types()`:
226
227         >>> h = Histogram()
228         >>> h.types()
229         ['chi-squared', 'jensen-shannon', 'mean', 'std-dev']
230
231         Selecting an invalid `type` raises an exception.
232
233         >>> h.residual(other=None, type='invalid-type')
234         Traceback (most recent call last):
235           ...
236         AttributeError: 'Histogram' object has no attribute 'invalid_type_residual'
237
238         For residual types where there is a convention, this histogram
239         is treated as the experimental histogram and the other
240         histogram is treated as the theoretical one.
241         """
242         r_method = getattr(self, self._type_to_method(type))
243         return r_method(other)