From: W. Trevor King Date: Wed, 27 Oct 2010 14:59:40 +0000 (-0700) Subject: Added bin/sawsim_param_error.py and pysawsim/parameter_error.py X-Git-Url: http://git.tremily.us/?p=sawsim.git;a=commitdiff_plain;h=a44755fbec6eeb54bd15d1b42c31381ef4c41cc7 Added bin/sawsim_param_error.py and pysawsim/parameter_error.py --- diff --git a/bin/sawsim_param_error.py b/bin/sawsim_param_error.py new file mode 100755 index 0000000..ed007e6 --- /dev/null +++ b/bin/sawsim_param_error.py @@ -0,0 +1,8 @@ +#!/usr/bin/env python + +from pysawsim.parameter_error import main +from pysawsim.manager.mpi import MPI_worker_death + + +MPI_worker_death() +main() diff --git a/pysawsim/parameter_error.py b/pysawsim/parameter_error.py new file mode 100755 index 0000000..8813e2f --- /dev/null +++ b/pysawsim/parameter_error.py @@ -0,0 +1,173 @@ +# Copyright (C) 2010 W. Trevor King +# +# 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 . +# +# The author may be contacted at 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 `(, )`', + '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: ` lines start each histogram. `params`', + 'lists the `sawsim` parameters that are unique to that', + 'experiment.', + '', + 'Each histogram line is of the format:', + '', + '', + '', + '`` 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]])