1 # Copyright (C) 2010 W. Trevor King <wking@drexel.edu>
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.
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.
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/>.
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 """Experiment vs. simulation comparison error bar calculation.
25 from . import PYSAWSIM_LOG_LEVEL_MSG as _PYSAWSIM_LOG_LEVEL_MSG
26 from .histogram import Histogram
27 from .parameter_scan import (
28 EXAMPLE_HISTOGRAM_FILE_CONTENTS, HistogramMatcher,
29 parse_param_ranges_string)
30 from .sawsim_histogram import sawsim_histogram
31 from .sawsim import SawsimRunner
34 _multiprocess_can_split_ = True
35 """Allow nosetests to split tests between processes.
39 def find_error_bounds(histogram_matcher, param_range, log_scale, threshold):
40 if log_scale == False:
41 ps = numpy.linspace(*param_range)
44 ps = numpy.exp(numpy.linspace(numpy.log(m), numpy.log(M), n))
47 for i,p in enumerate(ps):
48 residual[p] = histogram_matcher.residual([p])
51 r_best = min(residual.itervalues())
52 p_best = min([p for p,r in residual.iteritems() if r == r_best])
53 point['best fit'] = (p_best, r_best)
54 assert r_best <= threshold, (
55 'min residual %g (at %g) is over the threshold of %g'
56 % (r_best, p_best, threshold))
59 if residual[p] < threshold:
60 assert p_last != None, (
61 'range did not bracket lower error bound: residual(%g) = %g < %g'
62 % (p, residual[p], threshold))
63 point['error (low)'] = p_last, residual[p_last]
67 for p in reversed(ps):
68 if residual[p] < threshold:
69 assert p_last != None, (
70 'range did not bracket upper error bound: residual(%g) = %g < %g'
71 % (p, residual[p], threshold))
72 point['error (high)'] = p_last, residual[p_last]
75 return point, residual
81 >>> f = tempfile.NamedTemporaryFile()
82 >>> f.write(EXAMPLE_HISTOGRAM_FILE_CONTENTS)
84 >>> main(['-r', '[1e-6,1e-4,3]', '--log',
85 ... '-N', '4', '-t', '0.8',
87 ... # doctest: +ELLIPSIS, +NORMALIZE_WHITESPACE, +REPORT_UDIFF
88 #parameter value\tresidual\tdescription
89 #1e-05\t0...\tbest fit
90 #0.0001\t1...\terror (high)
91 #1e-06\t1...\terror (low)
97 from optparse import OptionParser
105 usage = '%prog [options] histogram_file'
107 'Compare simulated results against experimental over a sweep of',
108 'a single parameter. Prints `(<paramter value>, <residual>)`',
109 'pairs for every tested point, with comments listing the',
110 'best-fit, + error bound, and - error bound. The histogram file',
111 'should look something like:',
113 EXAMPLE_HISTOGRAM_FILE_CONTENTS,
115 '`#HISTOGRAM: <params>` lines start each histogram. `params`',
116 'lists the `sawsim` parameters that are unique to that',
119 'Each histogram line is of the format:',
121 '<bin_edge><whitespace><count>',
123 '`<bin_edge>` should mark the left-hand side of the bin, and',
124 'all bins should be of equal width (so we know where the last',
126 _PYSAWSIM_LOG_LEVEL_MSG,
128 parser = OptionParser(usage, epilog=epilog)
129 parser.format_epilog = lambda formatter: epilog+'\n'
130 for option in sr.optparse_options:
131 if option.dest == 'param_string':
133 parser.add_option(option)
134 parser.add_option('-f','--param-format', dest='param_format',
136 help='Convert params to sawsim options (%default).',
137 default=('-s cantilever,hooke,0.05 -N1 -s folded,null -N8 -s "unfolded,wlc,{0.39e-9,28e-9}" -k "folded,unfolded,bell,{%g,0.25e-9}" -q folded'))
138 parser.add_option('-r','--param-range', dest='param_range',
140 help='Param range for plotting (%default).',
141 default='[1e-5,1e-3,30]')
142 parser.add_option('--log', dest='log',
143 help='Use a log scale for the parameter range.',
144 default=False, action='store_true')
145 parser.add_option('-R','--residual', dest='residual',
147 help='Residual type (from %s; default: %%default).'
148 % ', '.join(Histogram().types()),
149 default='jensen-shannon')
150 parser.add_option('-P','--plot-residuals', dest='plot_residuals',
151 help='Generate residual difference plots for each point in the plot range.',
152 default=False, action='store_true')
153 parser.add_option('-t', '--threshold', dest='threshold',
154 metavar='FLOAT', type='float', default='0.2',
155 help='Residual value that defines the location of the error bar (%default).')
157 options,args = parser.parse_args(argv)
159 param_ranges = parse_param_ranges_string(options.param_range)
160 assert len(param_ranges) == 1, param_ranges
161 param_range = param_ranges[0]
163 histogram_file = args[0]
164 sr_call_params = sr.initialize_from_options(options)
167 hm = HistogramMatcher(
168 file(histogram_file, 'r'),
169 param_format_string=options.param_format,
170 sawsim_runner=sr, residual_type=options.residual,
171 plot=options.plot_residuals, **sr_call_params)
172 points,residuals = find_error_bounds(hm, param_range, options.log, options.threshold)
176 print '#%s' % '\t'.join(['parameter value', 'residual', 'description'])
177 for key,value in sorted(points.iteritems()):
178 print '#%s' % '\t'.join([str(x) for x in [value[0], value[1], key]])
179 for key,value in sorted(residuals.iteritems()):
180 print '%s' % '\t'.join([str(x) for x in [key, value]])