Add bin/sawsim_plot_histogram_file.py and Histogram.plot().

This moves some of the matplotlib initialization into the more basic
pysawsim.histogram module, and provides an easy way to display
histogram files.

diff --git a/bin/sawsim_plot_histogram_file.py b/bin/sawsim_plot_histogram_file.py
new file mode 100755 (executable)
index 0000000..25e2f96
--- /dev/null
+++ b/bin/sawsim_plot_histogram_file.py
@@ -0,0 +1,15 @@
+#!/usr/bin/env python
+
+from pysawsim.parameter_scan import HistogramMatcher
+
+
+if __name__ == '__main__':
+    import sys
+    
+    for histfile in sys.argv[1:]:
+        with open(histfile, 'r') as f:
+            histograms = HistogramMatcher._read_force_histograms(f)
+            for params,histogram in sorted(histograms.items()):
+                histogram.plot(
+                    title=params,
+                    filename='histogram{}.png'.format(params.replace(' ','_')))

diff --git a/pysawsim/histogram.py b/pysawsim/histogram.py
index acd135955e1223df0a0b0f591ceb8849cee54230..cdb0bbf086120789595bf4061969d3f1fd73b814 100644 (file)
--- a/pysawsim/histogram.py
+++ b/pysawsim/histogram.py
@@ -20,7 +20,10 @@
 """Histogram generation and comparison.
 """
 
 """Histogram generation and comparison.
 """
 

+import matplotlib
+matplotlib.use('Agg')  # select backend that doesn't require X Windows

 import numpy
 import numpy

+import pylab

 
 from . import log
 
 
 from . import log
 


@@ -29,6 +32,11 @@ _multiprocess_can_split_ = True
 """Allow nosetests to split tests between processes.
 """
 
 """Allow nosetests to split tests between processes.
 """
 

+FIGURE = pylab.figure()  # avoid memory problems.
+"""`pylab` keeps internal references to all created figures, so share
+a single instance.
+"""
+

 
 class Histogram (object):
     """A histogram with a flexible comparison method, `residual()`.
 
 class Histogram (object):
     """A histogram with a flexible comparison method, `residual()`.


@@ -243,3 +251,14 @@ class Histogram (object):
         """
         r_method = getattr(self, self._type_to_method(type))
         return r_method(other)
         """
         r_method = getattr(self, self._type_to_method(type))
         return r_method(other)

+
+    def plot(self, title=None, filename=None):
+        FIGURE.clear()
+        axes = FIGURE.add_subplot(1, 1, 1)
+        axes.hist(x=self.bin_edges[:-1],  # one fake entry for each bin
+                  weights=self.counts,    # weigh the fake entries by count
+                  bins=self.bin_edges,
+                  align='mid', histtype='stepfilled')
+        axes.set_title(title)
+        pylab.show()
+        FIGURE.savefig(filename)

diff --git a/pysawsim/parameter_scan.py b/pysawsim/parameter_scan.py
index 2dade8bd2624ac70ac8986d347442e6a83ce867c..f986ce81aaed693989bb3aabc7deb0f7ab3fefc1 100644 (file)
--- a/pysawsim/parameter_scan.py
+++ b/pysawsim/parameter_scan.py
@@ -25,26 +25,22 @@ import os.path
 import pickle
 from StringIO import StringIO
 
 import pickle
 from StringIO import StringIO
 

-import matplotlib
-matplotlib.use('Agg')  # select backend that doesn't require X Windows

 import numpy
 import numpy

-import pylab

 
 from . import log
 from . import PYSAWSIM_LOG_LEVEL_MSG as _PYSAWSIM_LOG_LEVEL_MSG
 
 from . import log
 from . import PYSAWSIM_LOG_LEVEL_MSG as _PYSAWSIM_LOG_LEVEL_MSG

-from .histogram import Histogram
+from .histogram import Histogram as _Histogram
+from .histogram import FIGURE as _FIGURE

 from .sawsim_histogram import sawsim_histogram
 from .sawsim import SawsimRunner
 
 from .sawsim_histogram import sawsim_histogram
 from .sawsim import SawsimRunner
 

+# import after .histogram, which selects an appropriate backend
+import pylab

 
 _multiprocess_can_split_ = True
 """Allow nosetests to split tests between processes.
 """
 
 
 _multiprocess_can_split_ = True
 """Allow nosetests to split tests between processes.
 """
 

-FIGURE = pylab.figure()  # avoid memory problems.
-"""`pylab` keeps internal references to all created figures, so share
-a single instance.
-"""

 EXAMPLE_HISTOGRAM_FILE_CONTENTS = """# Velocity histograms
 # Other general comments...
 
 EXAMPLE_HISTOGRAM_FILE_CONTENTS = """# Velocity histograms
 # Other general comments...
 


@@ -174,7 +170,8 @@ class HistogramMatcher (object):
         self.residual_type = residual_type
         self._plot = plot
 
         self.residual_type = residual_type
         self._plot = plot
 

-    def _read_force_histograms(self, stream):
+    @staticmethod
+    def _read_force_histograms(stream):

         """
         File format:
 
         """
         File format:
 


@@ -231,7 +228,7 @@ class HistogramMatcher (object):
         for params,block in hist_blocks.iteritems():
             if params == None:
                 continue
         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
             h.from_stream(StringIO('\n'.join(block)))
             histograms[params] = h
         return histograms


@@ -298,8 +295,8 @@ class HistogramMatcher (object):
         # import pickle
         # [X,Y,C] = pickle.load(file("histogram_matcher-XYC.pkl", "rb"))
         # ...
         # 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:
         if logx == True:
             axes.set_xscale('log')
         if logy == True:


@@ -310,18 +307,26 @@ class HistogramMatcher (object):
         else: # pseudocolor plot
             p = axes.pcolor(X, Y, C)
             axes.autoscale_view(tight=True)
         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):
 
     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):
 
 
 def parse_param_ranges_string(string):


@@ -411,7 +416,7 @@ def main(argv=None):
     parser.add_option('-R','--residual', dest='residual',
                       metavar='STRING',
                       help='Residual type (from %s; default: %%default).'
     parser.add_option('-R','--residual', dest='residual',
                       metavar='STRING',
                       help='Residual type (from %s; default: %%default).'

-                      % ', '.join(Histogram().types()),
+                      % ', '.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='jensen-shannon')
     parser.add_option('-P','--plot-residuals', dest='plot_residuals',
                       help='Generate residual difference plots for each point in the plot range.',