pysawsim/sawsim_histogram.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, Drudge's University, Physics Dept., 3141 Chestnut St.,
  18 # Philadelphia PA 19104, USA.
  19
  20 import numpy
  21
  22 from . import PYSAWSIM_LOG_LEVEL_MSG as _PYSAWSIM_LOG_LEVEL_MSG
  23 from .histogram import Histogram
  24 from .manager import MANAGERS, get_manager
  25 from .sawsim import SawsimRunner
  26
  27
  28 _multiprocess_can_split_ = True
  29 """Allow nosetests to split tests between processes.
  30 """
  31
  32
  33 def sawsim_histogram(sawsim_runner, param_string, N=400, bin_edges=None):
  34         """Run `N` simulations and return a histogram with `bin_edges`.
  35
  36         If `bin_edges == None`, return a numpy array of all unfolding
  37         forces.
  38         """
  39         events = []
  40         for run in sawsim_runner(param_string=param_string, N=N):
  41             events.extend([event.force for event in run
  42                            if (event.initial_state == 'folded'
  43                                and event.final_state == 'unfolded')])
  44         events = numpy.array(events)
  45         if bin_edges == None:
  46             return events
  47         h = Histogram()
  48         h.from_data(events, bin_edges)
  49         return h
  50
  51
  52 def main(argv=None):
  53     """
  54     >>> main(['-N', '2'])
  55     ... # doctest: +ELLIPSIS, +NORMALIZE_WHITESPACE
  56     #Force (N)  Unfolding events
  57     ...
  58     """
  59     from optparse import OptionParser
  60     import sys
  61
  62     if argv == None:
  63         argv = sys.argv[1:]
  64
  65     sr = SawsimRunner()
  66
  67     usage = '%prog [options]'
  68     epilog = '\n'.join([
  69             'Generate an unfolding force histogram from a series of `sawsim`',
  70             'runs.',
  71             _PYSAWSIM_LOG_LEVEL_MSG,
  72             ])
  73     parser = OptionParser(usage, epilog=epilog)
  74     parser.format_epilog = lambda formatter: epilog+'\n'
  75     for option in sr.optparse_options:
  76         parser.add_option(option)
  77     parser.add_option('-w', '--bin-width', dest='bin_width',
  78                       metavar='FLOAT', type='float',
  79                       help='Histogram bin width in newtons (%default).',
  80                       default=10e-12)
  81
  82     options,args = parser.parse_args(argv)
  83
  84     sr_call_params = sr.initialize_from_options(options)
  85     try:
  86         events = sawsim_histogram(
  87             sawsim_runner=sr, bin_edges=None, **sr_call_params)
  88     finally:
  89         sr.teardown()
  90
  91     if options.bin_width == None:
  92         sys.stdout.write('# Unfolding force (N)\n')
  93         events.tofile(sys.stdout, sep='\n')
  94         sys.stdout.write('\n')
  95     else:
  96         h = Histogram()
  97         h.from_data(events, bin_edges=h.calculate_bin_edges(
  98                 events, options.bin_width))
  99         h.headings = ['Force (N)', 'Unfolding events']
 100         h.to_stream(sys.stdout)