Add bin/sawsim_plot_histogram_file.py and Histogram.plot().
[sawsim.git] / pysawsim / parameter_scan.py
index b8e95c0df39d53fce92ccf725c3c52c3451b47ce..f986ce81aaed693989bb3aabc7deb0f7ab3fefc1 100644 (file)
@@ -25,21 +25,22 @@ import os.path
 import pickle
 from StringIO import StringIO
 
-import matplotlib
-matplotlib.use('Agg')  # select backend that doesn't require X Windows
 import numpy
-import pylab
 
 from . import log
-from .histogram import Histogram
+from . import PYSAWSIM_LOG_LEVEL_MSG as _PYSAWSIM_LOG_LEVEL_MSG
+from .histogram import Histogram as _Histogram
+from .histogram import FIGURE as _FIGURE
 from .sawsim_histogram import sawsim_histogram
 from .sawsim import SawsimRunner
 
+# import after .histogram, which selects an appropriate backend
+import pylab
 
-FIGURE = pylab.figure()  # avoid memory problems.
-"""`pylab` keeps internal references to all created figures, so share
-a single instance.
+_multiprocess_can_split_ = True
+"""Allow nosetests to split tests between processes.
 """
+
 EXAMPLE_HISTOGRAM_FILE_CONTENTS = """# Velocity histograms
 # Other general comments...
 
@@ -153,14 +154,14 @@ class HistogramMatcher (object):
     ...     '-s "unfolded,wlc,{0.39e-9,28e-9}" '
     ...     '-k "folded,unfolded,bell,{%g,%g}" -q folded')
     >>> m = ThreadManager()
-    >>> sr = SawsimRunner(sawsim='bin/sawsim', manager=m)
+    >>> sr = SawsimRunner(manager=m)
     >>> hm = HistogramMatcher(histogram_stream, param_format_string, sr, N=3)
     >>> hm.plot([[1e-5,1e-3,3],[0.1e-9,1e-9,3]], logx=True, logy=False)
     >>> m.teardown()
     """
     def __init__(self, histogram_stream, param_format_string,
                  sawsim_runner, N=400, residual_type='jensen-shannon',
-                 plot=True):
+                 plot=False):
         self.experiment_histograms = self._read_force_histograms(
             histogram_stream)
         self.param_format_string = param_format_string
@@ -169,7 +170,8 @@ class HistogramMatcher (object):
         self.residual_type = residual_type
         self._plot = plot
 
-    def _read_force_histograms(self, stream):
+    @staticmethod
+    def _read_force_histograms(stream):
         """
         File format:
 
@@ -226,7 +228,7 @@ class HistogramMatcher (object):
         for params,block in hist_blocks.iteritems():
             if params == None:
                 continue
-            h = Histogram()
+            h = _Histogram()
             h.from_stream(StringIO('\n'.join(block)))
             histograms[params] = h
         return histograms
@@ -256,7 +258,10 @@ class HistogramMatcher (object):
         log().debug('residual %s: %g' % (params, residual))
         return residual
 
-    def plot(self, param_ranges, logx=False, logy=False, contour=False):
+    def plot(self, param_ranges, logx=False, logy=False, contour=False,
+             csv=None):
+        if csv:
+            csv.write(','.join(('param 1', 'param 2', 'fit quality')) + '\n')
         xranges = param_ranges[0]
         yranges = param_ranges[1]
         if logx == False:
@@ -277,6 +282,8 @@ class HistogramMatcher (object):
                            % (i, j, i*(len(y)-1) + j, (len(x)-1)*(len(y)-1)))
                 params = (xi,yj)
                 r = self.residual(params)
+                if csv:
+                    csv.write(','.join([str(v) for v in (xi,yj,r)]) + '\n')
                 C[j,i] = numpy.log(r) # better resolution in valleys
                 if MEM_DEBUG == True:
                     log().debug('RSS: %d KB' % rss())
@@ -288,8 +295,8 @@ class HistogramMatcher (object):
         # import pickle
         # [X,Y,C] = pickle.load(file("histogram_matcher-XYC.pkl", "rb"))
         # ...
-        FIGURE.clear()
-        axes = FIGURE.add_subplot(111)
+        _FIGURE.clear()
+        axes = _FIGURE.add_subplot(1, 1, 1)
         if logx == True:
             axes.set_xscale('log')
         if logy == True:
@@ -300,18 +307,26 @@ class HistogramMatcher (object):
         else: # pseudocolor plot
             p = axes.pcolor(X, Y, C)
             axes.autoscale_view(tight=True)
-        FIGURE.colorbar(p)
-        FIGURE.savefig("figure.png")
+        _FIGURE.colorbar(p)
+        _FIGURE.savefig("figure.png")
 
     def _plot_residual_comparison(self, experiment_hist, theory_hist,
                                   residual, title, filename):
-        FIGURE.clear()
-        p = pylab.plot(experiment_hist.bin_edges[:-1],
-                       experiment_hist.probabilities, 'r-',
-                       theory_hist.bin_edges[:-1],
-                       theory_hist.probabilities, 'b-')
-        pylab.title(title)
-        FIGURE.savefig(filename)
+        _FIGURE.clear()
+        axes = _FIGURE.add_subplot(1, 1, 1)
+        axes.hold(True)
+        axes.hist(x=experiment_hist.bin_edges[:-1],  # one fake entry for each bin
+                  weights=experiment_hist.counts,    # weigh the fake entries by count
+                  bins=experiment_hist.bin_edges,
+                  align='mid', histtype='stepfilled',
+                  color='r')
+        axes.hist(x=theory_hist.bin_edges[:-1],  # one fake entry for each bin
+                  weights=theory_hist.counts,    # weigh the fake entries by count
+                  bins=theory_hist.bin_edges,
+                  align='mid', histtype='stepfilled',
+                  color='b', alpha=0.5)
+        axes.set_title(title)
+        _FIGURE.savefig(filename)
 
 
 def parse_param_ranges_string(string):
@@ -339,8 +354,7 @@ def main(argv=None):
     >>> f = tempfile.NamedTemporaryFile()
     >>> f.write(EXAMPLE_HISTOGRAM_FILE_CONTENTS)
     >>> f.flush()
-    >>> main(['-s', 'bin/sawsim',
-    ...       '-r', '[1e-5,1e-3,3],[0.1e-9,1e-9,3]',
+    >>> main(['-r', '[1e-5,1e-3,3],[0.1e-9,1e-9,3]',
     ...       '-N', '2',
     ...       f.name])
     >>> f.close()
@@ -373,6 +387,7 @@ def main(argv=None):
             '`<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).',
+            _PYSAWSIM_LOG_LEVEL_MSG,
             ])
     parser = OptionParser(usage, epilog=epilog)
     parser.format_epilog = lambda formatter: epilog+'\n'
@@ -380,41 +395,44 @@ def main(argv=None):
         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).",
+    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,%g}" -q folded'))
-    parser.add_option("-p","--initial-params", dest="initial_params",
-                      metavar="PARAMS",
-                      help="Initial params for fitting (%default).",
+    parser.add_option('-p','--initial-params', dest='initial_params',
+                      metavar='PARAMS',
+                      help='Initial params for fitting (%default).',
                       default='3.3e-4,0.25e-9')
-    parser.add_option("-r","--param-range", dest="param_range",
-                      metavar="PARAMS",
-                      help="Param range for plotting (%default).",
+    parser.add_option('-r','--param-range', dest='param_range',
+                      metavar='PARAMS',
+                      help='Param range for plotting (%default).',
                       default='[1e-5,1e-3,20],[0.1e-9,1e-9,20]')
-    parser.add_option("--logx", dest="logx",
-                      help="Use a log scale for the x range.",
-                      default=False, action="store_true")
-    parser.add_option("--logy", dest="logy",
-                      help="Use a log scale for the y 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()),
+    parser.add_option('--logx', dest='logx',
+                      help='Use a log scale for the x range.',
+                      default=False, action='store_true')
+    parser.add_option('--logy', dest='logy',
+                      help='Use a log scale for the y 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("-c","--contour-plot", dest="contour_plot",
-                      help="Select contour plot (vs. the default pseudocolor plot).",
-                      default=False, action="store_true")
+    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('-c','--contour-plot', dest='contour_plot',
+                      help='Select contour plot (vs. the default pseudocolor plot).',
+                      default=False, action='store_true')
+    parser.add_option('--csv', dest='csv', metavar='FILE',
+                      help='Save fit qualities to a comma-separated value file FILE.'),
 
     options,args = parser.parse_args(argv)
 
     initial_params = [float(p) for p in options.initial_params.split(",")]
     param_ranges = parse_param_ranges_string(options.param_range)
     histogram_file = args[0]
+    csv = None
     sr_call_params = sr.initialize_from_options(options)
 
     try:
@@ -424,7 +442,11 @@ def main(argv=None):
             sawsim_runner=sr, residual_type=options.residual,
             plot=options.plot_residuals, **sr_call_params)
         #hm.fit(initial_params)
+        if options.csv:
+            csv = open(options.csv, 'w')
         hm.plot(param_ranges, logx=options.logx, logy=options.logy,
-                contour=options.contour_plot)
+                contour=options.contour_plot, csv=csv)
     finally:
         sr.teardown()
+        if csv:
+            csv.close()