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