pysawsim/parameter_error.py

   1 # Copyright (C) 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 """Experiment vs. simulation comparison error bar calculation.
  21 """
  22
  23 import numpy
  24
  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
  32
  33
  34 _multiprocess_can_split_ = True
  35 """Allow nosetests to split tests between processes.
  36 """
  37
  38
  39 def find_error_bounds(histogram_matcher, param_range, log_scale, threshold):
  40     if log_scale == False:
  41         ps = numpy.linspace(*param_range)
  42     else:
  43         m,M,n = param_range
  44         ps = numpy.exp(numpy.linspace(numpy.log(m), numpy.log(M), n))
  45
  46     residual = {}
  47     for i,p in enumerate(ps):
  48         residual[p] = histogram_matcher.residual([p])
  49
  50     point = {}
  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))
  57     p_last = None
  58     for p in ps:
  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]
  64             break
  65         p_last = p
  66     p_last = None
  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]
  73             break
  74         p_last = p
  75     return point, residual
  76
  77
  78 def main(argv=None):
  79     """
  80     >>> import tempfile
  81     >>> f = tempfile.NamedTemporaryFile()
  82     >>> f.write(EXAMPLE_HISTOGRAM_FILE_CONTENTS)
  83     >>> f.flush()
  84     >>> main(['-r', '[1e-6,1e-4,3]', '--log',
  85     ...       '-N', '4', '-t', '0.8',
  86     ...       f.name])
  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)
  92     1e-06\t1...
  93     1e-05\t0...
  94     0.0001\t1...
  95     >>> f.close()
  96     """
  97     from optparse import OptionParser
  98     import sys
  99
 100     if argv == None:
 101         argv = sys.argv[1:]
 102
 103     sr = SawsimRunner()
 104
 105     usage = '%prog [options] histogram_file'
 106     epilog = '\n'.join([
 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:',
 112             '',
 113             EXAMPLE_HISTOGRAM_FILE_CONTENTS,
 114             ''
 115             '`#HISTOGRAM: <params>` lines start each histogram.  `params`',
 116             'lists the `sawsim` parameters that are unique to that',
 117             'experiment.',
 118             '',
 119             'Each histogram line is of the format:',
 120             '',
 121             '<bin_edge><whitespace><count>',
 122             '',
 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',
 125             'one ends).',
 126             _PYSAWSIM_LOG_LEVEL_MSG,
 127             ])
 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':
 132             continue
 133         parser.add_option(option)
 134     parser.add_option('-f','--param-format', dest='param_format',
 135                       metavar='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',
 139                       metavar='PARAMS',
 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',
 146                       metavar='STRING',
 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).')
 156
 157     options,args = parser.parse_args(argv)
 158
 159     param_ranges = parse_param_ranges_string(options.param_range)
 160     assert len(param_ranges) == 1, param_ranges
 161     param_range = param_ranges[0]
 162
 163     histogram_file = args[0]
 164     sr_call_params = sr.initialize_from_options(options)
 165
 166     try:
 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)
 173     finally:
 174         sr.teardown()
 175
 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]])