pysawsim/parameter_scan.py

   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.
  19
  20 """Experiment vs. simulation comparison and scanning.
  21 """
  22
  23 from os import getpid  # for rss()
  24 import os.path
  25 import pickle
  26 from StringIO import StringIO
  27
  28 import matplotlib
  29 matplotlib.use('Agg')  # select backend that doesn't require X Windows
  30 import numpy
  31 import pylab
  32
  33 from . import log
  34 from .histogram import Histogram
  35 from .sawsim_histogram import sawsim_histogram
  36 from .sawsim import SawsimRunner
  37
  38
  39 _multiprocess_can_split_ = True
  40 """Allow nosetests to split tests between processes.
  41 """
  42
  43 FIGURE = pylab.figure()  # avoid memory problems.
  44 """`pylab` keeps internal references to all created figures, so share
  45 a single instance.
  46 """
  47 EXAMPLE_HISTOGRAM_FILE_CONTENTS = """# Velocity histograms
  48 # Other general comments...
  49
  50 #HISTOGRAM: -v 6e-7
  51 #Force (N)\tUnfolding events
  52 1.4e-10\t1
  53 1.5e-10\t0
  54 1.6e-10\t4
  55 1.7e-10\t6
  56 1.8e-10\t8
  57 1.9e-10\t20
  58 2e-10\t28
  59 2.1e-10\t38
  60 2.2e-10\t72
  61 2.3e-10\t110
  62 2.4e-10\t155
  63 2.5e-10\t247
  64 2.6e-10\t395
  65 2.7e-10\t451
  66 2.8e-10\t430
  67 2.9e-10\t300
  68 3e-10\t116
  69 3.1e-10\t18
  70 3.2e-10\t1
  71
  72 #HISTOGRAM: -v 8e-7
  73 #Force (N)\tUnfolding events
  74 8e-11\t1
  75 9e-11\t0
  76 1e-10\t0
  77 1.1e-10\t1
  78 1.2e-10\t0
  79 1.3e-10\t0
  80 1.4e-10\t0
  81 1.5e-10\t3
  82 1.6e-10\t3
  83 1.7e-10\t4
  84 1.8e-10\t4
  85 1.9e-10\t13
  86 2e-10\t29
  87 2.1e-10\t39
  88 2.2e-10\t60
  89 2.3e-10\t102
  90 2.4e-10\t154
  91 2.5e-10\t262
  92 2.6e-10\t402
  93 2.7e-10\t497
  94 2.8e-10\t541
  95 2.9e-10\t555
  96 3e-10\t325
  97 3.1e-10\t142
  98 3.2e-10\t50
  99 3.3e-10\t13
 100
 101 #HISTOGRAM: -v 1e-6
 102 #Force (N)\tUnfolding events
 103 1.5e-10\t2
 104 1.6e-10\t3
 105 1.7e-10\t7
 106 1.8e-10\t8
 107 1.9e-10\t7
 108 2e-10\t25
 109 2.1e-10\t30
 110 2.2e-10\t58
 111 2.3e-10\t76
 112 2.4e-10\t159
 113 2.5e-10\t216
 114 2.6e-10\t313
 115 2.7e-10\t451
 116 2.8e-10\t568
 117 2.9e-10\t533
 118 3e-10\t416
 119 3.1e-10\t222
 120 3.2e-10\t80
 121 3.3e-10\t24
 122 3.4e-10\t2
 123 """
 124
 125
 126 MEM_DEBUG = False
 127
 128
 129
 130 def rss():
 131     """
 132     For debugging memory usage.
 133
 134     resident set size, the non-swapped physical memory that a task has
 135     used (in kilo-bytes).
 136     """
 137     call = "ps -o rss= -p %d" % getpid()
 138     status,stdout,stderr = invoke(call)
 139     return int(stdout)
 140
 141
 142 class HistogramMatcher (object):
 143     """Compare experimental histograms to simulated data.
 144
 145     The main entry points are `fit()` and `plot()`.
 146
 147     The input `histogram_stream` should contain a series of
 148     experimental histograms with '#HISTOGRAM: <params>` lines starting
 149     each histogram.  `<params>` lists the `sawsim` parameters that are
 150     unique to that experiment.
 151
 152     >>> from .manager.thread import ThreadManager
 153     >>> histogram_stream = StringIO(EXAMPLE_HISTOGRAM_FILE_CONTENTS)
 154     >>> param_format_string = (
 155     ...     '-s cantilever,hooke,0.05 -N1 '
 156     ...     '-s folded,null -N8 '
 157     ...     '-s "unfolded,wlc,{0.39e-9,28e-9}" '
 158     ...     '-k "folded,unfolded,bell,{%g,%g}" -q folded')
 159     >>> m = ThreadManager()
 160     >>> sr = SawsimRunner(manager=m)
 161     >>> hm = HistogramMatcher(histogram_stream, param_format_string, sr, N=3)
 162     >>> hm.plot([[1e-5,1e-3,3],[0.1e-9,1e-9,3]], logx=True, logy=False)
 163     >>> m.teardown()
 164     """
 165     def __init__(self, histogram_stream, param_format_string,
 166                  sawsim_runner, N=400, residual_type='jensen-shannon',
 167                  plot=False):
 168         self.experiment_histograms = self._read_force_histograms(
 169             histogram_stream)
 170         self.param_format_string = param_format_string
 171         self.sawsim_runner = sawsim_runner
 172         self.N = N
 173         self.residual_type = residual_type
 174         self._plot = plot
 175
 176     def _read_force_histograms(self, stream):
 177         """
 178         File format:
 179
 180           # comment and blank lines ignored
 181           #HISTOGRAM: <histogram-specific params>
 182           <pysawsim.histogram.Histogram-compatible histogram>
 183           #HISTOGRAM: <other histogram-specific params>
 184           <another pysawsim.histogram.Histogram-compatible histogram>
 185           ...
 186
 187         >>> import sys
 188         >>> stream = StringIO(EXAMPLE_HISTOGRAM_FILE_CONTENTS)
 189         >>> hm = HistogramMatcher(StringIO(), None, None, None)
 190         >>> histograms = hm._read_force_histograms(stream)
 191         >>> sorted(histograms.iterkeys())
 192         ['-v 1e-6', '-v 6e-7', '-v 8e-7']
 193         >>> histograms['-v 1e-6'].to_stream(sys.stdout)
 194         ... # doctest: +NORMALIZE_WHITESPACE, +REPORT_UDIFF
 195         #Force (N)\tUnfolding events
 196         1.5e-10\t2
 197         1.6e-10\t3
 198         1.7e-10\t7
 199         1.8e-10\t8
 200         1.9e-10\t7
 201         2e-10\t25
 202         2.1e-10\t30
 203         2.2e-10\t58
 204         2.3e-10\t76
 205         2.4e-10\t159
 206         2.5e-10\t216
 207         2.6e-10\t313
 208         2.7e-10\t451
 209         2.8e-10\t568
 210         2.9e-10\t533
 211         3e-10\t416
 212         3.1e-10\t222
 213         3.2e-10\t80
 214         3.3e-10\t24
 215         3.4e-10\t2
 216         """
 217         token = '#HISTOGRAM:'
 218         hist_blocks = {None: []}
 219         params = None
 220         for line in stream.readlines():
 221             line = line.strip()
 222             if line.startswith(token):
 223                 params = line[len(token):].strip()
 224                 assert params not in hist_blocks, params
 225                 hist_blocks[params] = []
 226             else:
 227                 hist_blocks[params].append(line)
 228
 229         histograms = {}
 230         for params,block in hist_blocks.iteritems():
 231             if params == None:
 232                 continue
 233             h = Histogram()
 234             h.from_stream(StringIO('\n'.join(block)))
 235             histograms[params] = h
 236         return histograms
 237
 238     def param_string(self, params, hist_params):
 239         """Generate a string of options to pass to `sawsim`.
 240         """
 241         return '%s %s' % (
 242             self.param_format_string % tuple(params), hist_params)
 243
 244     def residual(self, params):
 245         residual = 0
 246         for hist_params,experiment_hist in self.experiment_histograms.iteritems():
 247             sawsim_hist = sawsim_histogram(
 248                 sawsim_runner=self.sawsim_runner,
 249                 param_string=self.param_string(params, hist_params),
 250                 N=self.N, bin_edges=experiment_hist.bin_edges)
 251             r = experiment_hist.residual(sawsim_hist, type=self.residual_type)
 252             residual += r
 253             if self._plot == True:
 254                 title = ", ".join(["%g" % p for p in params]+[hist_params])
 255                 filename = "residual-%s-%g.png" % (
 256                     title.replace(', ', '_').replace(' ', '_'), r)
 257                 self._plot_residual_comparison(
 258                     experiment_hist, sawsim_hist, residual=r,
 259                     title=title, filename=filename)
 260         log().debug('residual %s: %g' % (params, residual))
 261         return residual
 262
 263     def plot(self, param_ranges, logx=False, logy=False, contour=False,
 264              csv=None):
 265         if csv:
 266             csv.write(','.join(('param 1', 'param 2', 'fit quality')) + '\n')
 267         xranges = param_ranges[0]
 268         yranges = param_ranges[1]
 269         if logx == False:
 270             x = numpy.linspace(*xranges)
 271         else:
 272             m,M,n = xranges
 273             x = numpy.exp(numpy.linspace(numpy.log(m), numpy.log(M), n))
 274         if logy == False:
 275             y = numpy.linspace(*yranges)
 276         else:
 277             m,M,n = yranges
 278             y = numpy.exp(numpy.linspace(numpy.log(m), numpy.log(M), n))
 279         X, Y = pylab.meshgrid(x,y)
 280         C = numpy.zeros((len(y)-1, len(x)-1))
 281         for i,xi in enumerate(x[:-1]):
 282             for j,yj in enumerate(y[:-1]):
 283                 log().info('point %d %d (%d of %d)'
 284                            % (i, j, i*(len(y)-1) + j, (len(x)-1)*(len(y)-1)))
 285                 params = (xi,yj)
 286                 r = self.residual(params)
 287                 if csv:
 288                     csv.write(','.join([str(v) for v in (xi,yj,r)]) + '\n')
 289                 C[j,i] = numpy.log(r) # better resolution in valleys
 290                 if MEM_DEBUG == True:
 291                     log().debug('RSS: %d KB' % rss())
 292         C = numpy.nan_to_num(C) # NaN -> 0
 293         fid = file("histogram_matcher-XYC.pkl", "wb")
 294         pickle.dump([X,Y,C], fid)
 295         fid.close()
 296         # read in with
 297         # import pickle
 298         # [X,Y,C] = pickle.load(file("histogram_matcher-XYC.pkl", "rb"))
 299         # ...
 300         FIGURE.clear()
 301         axes = FIGURE.add_subplot(111)
 302         if logx == True:
 303             axes.set_xscale('log')
 304         if logy == True:
 305             axes.set_yscale('log')
 306         if contour == True:
 307             p = axes.contour(X[:-1,:-1], Y[:-1,:-1], C)
 308             # [:-1,:-1] to strip dummy last row & column from X&Y.
 309         else: # pseudocolor plot
 310             p = axes.pcolor(X, Y, C)
 311             axes.autoscale_view(tight=True)
 312         FIGURE.colorbar(p)
 313         FIGURE.savefig("figure.png")
 314
 315     def _plot_residual_comparison(self, experiment_hist, theory_hist,
 316                                   residual, title, filename):
 317         FIGURE.clear()
 318         p = pylab.plot(experiment_hist.bin_edges[:-1],
 319                        experiment_hist.probabilities, 'r-',
 320                        theory_hist.bin_edges[:-1],
 321                        theory_hist.probabilities, 'b-')
 322         pylab.title(title)
 323         FIGURE.savefig(filename)
 324
 325
 326 def parse_param_ranges_string(string):
 327     """Parse parameter range stings.
 328
 329     '[Amin,Amax,Asteps],[Bmin,Bmax,Bsteps],...'
 330       ->
 331     [[Amin,Amax,Asteps],[Bmin,Bmax,Bsteps],...]
 332
 333     >>> parse_param_ranges_string('[1,2,3],[4,5,6]')
 334     [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]
 335     >>> parse_param_ranges_string('[1,2,3]')
 336     [[1.0, 2.0, 3.0]]
 337     """
 338     ranges = []
 339     for range_string in string.split("],["):
 340         range_number_strings = range_string.strip("[]").split(",")
 341         ranges.append([float(x) for x in range_number_strings])
 342     return ranges
 343
 344
 345 def main(argv=None):
 346     """
 347     >>> import tempfile
 348     >>> f = tempfile.NamedTemporaryFile()
 349     >>> f.write(EXAMPLE_HISTOGRAM_FILE_CONTENTS)
 350     >>> f.flush()
 351     >>> main(['-r', '[1e-5,1e-3,3],[0.1e-9,1e-9,3]',
 352     ...       '-N', '2',
 353     ...       f.name])
 354     >>> f.close()
 355     """
 356     from optparse import OptionParser
 357     import sys
 358
 359     if argv == None:
 360         argv = sys.argv[1:]
 361
 362     sr = SawsimRunner()
 363
 364     usage = '%prog [options] histogram_file'
 365     epilog = '\n'.join([
 366             'Compare simulated results against experimental values over a',
 367             'range of parameters.  Generates a plot of fit quality over',
 368             'the parameter space.  The histogram file should look something',
 369             'like:',
 370             '',
 371             EXAMPLE_HISTOGRAM_FILE_CONTENTS,
 372             ''
 373             '`#HISTOGRAM: <params>` lines start each histogram.  `params`',
 374             'lists the `sawsim` parameters that are unique to that',
 375             'experiment.',
 376             '',
 377             'Each histogram line is of the format:',
 378             '',
 379             '<bin_edge><whitespace><count>',
 380             '',
 381             '`<bin_edge>` should mark the left-hand side of the bin, and',
 382             'all bins should be of equal width (so we know where the last',
 383             'one ends).',
 384             ])
 385     parser = OptionParser(usage, epilog=epilog)
 386     parser.format_epilog = lambda formatter: epilog+'\n'
 387     for option in sr.optparse_options:
 388         if option.dest == 'param_string':
 389             continue
 390         parser.add_option(option)
 391     parser.add_option('-f','--param-format', dest='param_format',
 392                       metavar='FORMAT',
 393                       help='Convert params to sawsim options (%default).',
 394                       default=('-s cantilever,hooke,0.05 -N1 -s folded,null -N8 -s "unfolded,wlc,{0.39e-9,28e-9}" -k "folded,unfolded,bell,{%g,%g}" -q folded'))
 395     parser.add_option('-p','--initial-params', dest='initial_params',
 396                       metavar='PARAMS',
 397                       help='Initial params for fitting (%default).',
 398                       default='3.3e-4,0.25e-9')
 399     parser.add_option('-r','--param-range', dest='param_range',
 400                       metavar='PARAMS',
 401                       help='Param range for plotting (%default).',
 402                       default='[1e-5,1e-3,20],[0.1e-9,1e-9,20]')
 403     parser.add_option('--logx', dest='logx',
 404                       help='Use a log scale for the x range.',
 405                       default=False, action='store_true')
 406     parser.add_option('--logy', dest='logy',
 407                       help='Use a log scale for the y range.',
 408                       default=False, action='store_true')
 409     parser.add_option('-R','--residual', dest='residual',
 410                       metavar='STRING',
 411                       help='Residual type (from %s; default: %%default).'
 412                       % ', '.join(Histogram().types()),
 413                       default='jensen-shannon')
 414     parser.add_option('-P','--plot-residuals', dest='plot_residuals',
 415                       help='Generate residual difference plots for each point in the plot range.',
 416                       default=False, action='store_true')
 417     parser.add_option('-c','--contour-plot', dest='contour_plot',
 418                       help='Select contour plot (vs. the default pseudocolor plot).',
 419                       default=False, action='store_true')
 420     parser.add_option('--csv', dest='csv', metavar='FILE',
 421                       help='Save fit qualities to a comma-separated value file FILE.'),
 422
 423     options,args = parser.parse_args(argv)
 424
 425     initial_params = [float(p) for p in options.initial_params.split(",")]
 426     param_ranges = parse_param_ranges_string(options.param_range)
 427     histogram_file = args[0]
 428     csv = None
 429     sr_call_params = sr.initialize_from_options(options)
 430
 431     try:
 432         hm = HistogramMatcher(
 433             file(histogram_file, 'r'),
 434             param_format_string=options.param_format,
 435             sawsim_runner=sr, residual_type=options.residual,
 436             plot=options.plot_residuals, **sr_call_params)
 437         #hm.fit(initial_params)
 438         if options.csv:
 439             csv = open(options.csv, 'w')
 440         hm.plot(param_ranges, logx=options.logx, logy=options.logy,
 441                 contour=options.contour_plot, csv=csv)
 442     finally:
 443         sr.teardown()
 444         if csv:
 445             csv.close()