--- /dev/null
+# Copyright (C) 2010 W. Trevor King <wking@drexel.edu>
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+#
+# The author may be contacted at <wking@drexel.edu> on the Internet, or
+# write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
+# Philadelphia PA 19104, USA.
+
+"""Experiment vs. simulation comparison error bar calculation.
+"""
+
+import numpy
+
+from .histogram import Histogram
+from .parameter_scan import (
+ EXAMPLE_HISTOGRAM_FILE_CONTENTS, HistogramMatcher,
+ parse_param_ranges_string)
+from .sawsim_histogram import sawsim_histogram
+from .sawsim import SawsimRunner
+
+
+def find_error_bounds(histogram_matcher, param_range, log_scale, threshold):
+ if log_scale == False:
+ ps = numpy.linspace(*param_range)
+ else:
+ m,M,n = param_range
+ ps = numpy.exp(numpy.linspace(numpy.log(m), numpy.log(M), n))
+
+ residual = {}
+ for i,p in enumerate(ps):
+ residual[p] = histogram_matcher.residual([p])
+
+ point = {}
+ r_best = min(residual.itervalues())
+ p_best = min([p for p,r in residual.iteritems() if r == r_best])
+ point['best fit'] = (p_best, r_best)
+ assert r_best <= threshold, (
+ 'min residual %g (at %g) is over the threshold of %g'
+ % (r_best, p_best, threshold))
+ p_last = None
+ for p in ps:
+ if residual[p] < threshold:
+ assert p_last != None, (
+ 'range did not bracket lower error bound: residual(%g) = %g < %g'
+ % (p, residual[p], threshold))
+ point['error (low)'] = p_last, residual[p_last]
+ break
+ p_last = p
+ p_last = None
+ for p in reversed(ps):
+ if residual[p] < threshold:
+ assert p_last != None, (
+ 'range did not bracket upper error bound: residual(%g) = %g < %g'
+ % (p, residual[p], threshold))
+ point['error (high)'] = p_last, residual[p_last]
+ break
+ p_last = p
+ return point, residual
+
+
+def main(argv=None):
+ """
+ >>> import tempfile
+ >>> f = tempfile.NamedTemporaryFile()
+ >>> f.write(EXAMPLE_HISTOGRAM_FILE_CONTENTS)
+ >>> f.flush()
+ >>> main(['-s', 'bin/sawsim',
+ ... '-r', '[1e-6,1e-4,3]', '--log',
+ ... '-N', '4', '-t', '0.8',
+ ... f.name])
+ ... # doctest: +ELLIPSIS, +NORMALIZE_WHITESPACE, +REPORT_UDIFF
+ #parameter value\tresidual\tdescription
+ #1e-05\t0...\tbest fit
+ #0.0001\t1...\terror (high)
+ #1e-06\t1...\terror (low)
+ 1e-06\t1...
+ 1e-05\t0...
+ 0.0001\t1...
+ >>> f.close()
+ """
+ from optparse import OptionParser
+ import sys
+
+ if argv == None:
+ argv = sys.argv[1:]
+
+ sr = SawsimRunner()
+
+ usage = '%prog [options] histogram_file'
+ epilog = '\n'.join([
+ 'Compare simulated results against experimental over a sweep of',
+ 'a single parameter. Prints `(<paramter value>, <residual>)`',
+ 'pairs for every tested point, with comments listing the',
+ 'best-fit, + error bound, and - error bound. The histogram file',
+ 'should look something like:',
+ '',
+ EXAMPLE_HISTOGRAM_FILE_CONTENTS,
+ ''
+ '`#HISTOGRAM: <params>` lines start each histogram. `params`',
+ 'lists the `sawsim` parameters that are unique to that',
+ 'experiment.',
+ '',
+ 'Each histogram line is of the format:',
+ '',
+ '<bin_edge><whitespace><count>',
+ '',
+ '`<bin_edge>` should mark the left-hand side of the bin, and',
+ 'all bins should be of equal width (so we know where the last',
+ 'one ends).',
+ ])
+ parser = OptionParser(usage, epilog=epilog)
+ for option in sr.optparse_options:
+ if option.dest == 'param_string':
+ continue
+ parser.add_option(option)
+ parser.add_option('-f','--param-format', dest='param_format',
+ metavar='FORMAT',
+ help='Convert params to sawsim options (%default).',
+ 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'))
+ parser.add_option('-r','--param-range', dest='param_range',
+ metavar='PARAMS',
+ help='Param range for plotting (%default).',
+ default='[1e-5,1e-3,30]')
+ parser.add_option('--log', dest='log',
+ help='Use a log scale for the parameter range.',
+ default=False, action='store_true')
+ parser.add_option('-R','--residual', dest='residual',
+ metavar='STRING',
+ help='Residual type (from %s; default: %%default).'
+ % ', '.join(Histogram().types()),
+ default='jensen-shannon')
+ parser.add_option("-P","--plot-residuals", dest="plot_residuals",
+ help="Generate residual difference plots for each point in the plot range.",
+ default=False, action="store_true")
+ parser.add_option('-t', '--threshold', dest='threshold',
+ metavar='FLOAT', type='float', default='0.2',
+ help='Residual value that defines the location of the error bar (%default).')
+
+ options,args = parser.parse_args(argv)
+
+ param_ranges = parse_param_ranges_string(options.param_range)
+ assert len(param_ranges) == 1, param_ranges
+ param_range = param_ranges[0]
+
+ histogram_file = args[0]
+ sr_call_params = sr.initialize_from_options(options)
+
+ try:
+ hm = HistogramMatcher(
+ file(histogram_file, 'r'),
+ param_format_string=options.param_format,
+ sawsim_runner=sr, residual_type=options.residual,
+ plot=options.plot_residuals, **sr_call_params)
+ points,residuals = find_error_bounds(hm, param_range, options.log, options.threshold)
+ finally:
+ sr.teardown()
+
+ print '#%s' % '\t'.join(['parameter value', 'residual', 'description'])
+ for key,value in sorted(points.iteritems()):
+ print '#%s' % '\t'.join([str(x) for x in [value[0], value[1], key]])
+ for key,value in sorted(residuals.iteritems()):
+ print '%s' % '\t'.join([str(x) for x in [key, value]])
projects / sawsim.git / commitdiff