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 from __future__ import with_statement
  21
  22 import hashlib
  23 import os.path
  24 import shutil
  25
  26 import numpy
  27
  28 from . import __version__, log
  29 from .histogram import Histogram
  30 from .manager import InvokeJob, MANAGERS, get_manager
  31 from .sawsim import parse
  32
  33
  34 SAWSIM = 'sawsim'  # os.path.expand(os.path.join('~', 'bin', 'sawsim'))
  35 CACHE_DIR = os.path.expanduser(os.path.join('~', '.sawsim_histogram'))
  36 DEFAULT_PARAM_STRING = (
  37     '-s cantilever,hooke,0.05 -N1 '
  38     '-s folded,null -N8 '
  39     "-s 'unfolded,wlc,{0.39e-9,28e-9}' "
  40     "-k 'folded,unfolded,bell,{3.3e-4,0.25e-9}' "
  41     '-q folded -v 1e-6')
  42
  43
  44 class SawsimHistogram (object):
  45     def __init__(self, use_cache=False, clean_cache=False):
  46         self._use_cache = use_cache
  47         self._clean_cache = clean_cache
  48
  49     def _cache_dir(self, param_string):
  50         """
  51         >>> s = SawsimHistogram()
  52         >>> s._cache_dir(DEFAULT_PARAM_STRING)  # doctest: +ELLIPSIS
  53         '/.../.sawsim_histogram/...'
  54         """
  55         return os.path.join(
  56             CACHE_DIR, hashlib.sha256(param_string).hexdigest())
  57
  58     def _make_cache(self, param_string):
  59         cache_dir = self._cache_dir(param_string)
  60         os.makedirs(cache_dir)
  61         with open(os.path.join(cache_dir, 'param_string'), 'w') as f:
  62             f.write('# version: %s\n%s\n' % (__version__, param_string))
  63
  64     def _load_cached_data(self, cache_dir, N):
  65         data = {}
  66         for i in range(N):
  67             simulation_path = os.path.join(cache_dir, '%d.dat' % i)
  68             if os.path.exists(simulation_path):
  69                 with open(simulation_path, 'r') as f:
  70                     data[i] = parse(f.read())
  71             else:
  72                 break
  73         return data
  74
  75     def _cache_run(self, cache_dir, i, stdout):
  76         simulation_path = os.path.join(cache_dir, '%d.dat' % i)
  77         with open(simulation_path, 'w') as f:
  78             f.write(stdout)
  79
  80     def __call__(self, param_string=None, N=400, bin_edges=None, manager=None):
  81         """Run `N` simulations and return a histogram with `bin_edges`.
  82
  83         If `bin_edges == None`, return a numpy array of all unfolding
  84         forces.
  85
  86         Use the `JobManager` instance `manager` for asynchronous job
  87         execution.
  88
  89         If `_use_cache` is `True`, store an array of unfolding forces
  90         in `CACHE_DIR` for each simulated pull.  If the cached forces
  91         are already present for `param_string`, do not redo the
  92         simulation unless `_clean_cache` is `True`.
  93         """
  94         data = {}
  95         if self._use_cache == True:
  96             d = self._cache_dir(param_string)
  97             if os.path.exists(d):
  98                 if self._clean_cache == True:
  99                     shutil.rmtree(d)
 100                     self._make_cache(param_string)
 101                 else:
 102                     data = self._load_cached_data(d, N)
 103                     log().debug('using %d cached runs for %s'
 104                                 % (len(data), param_string))
 105             else:
 106                 self._make_cache(param_string)
 107
 108         jobs = {}
 109         for i in range(N):
 110             if i in data:
 111                 continue
 112             jobs[i] = manager.async_invoke(InvokeJob(
 113                     target='%s %s' % (SAWSIM, param_string)))
 114         complete_jobs = manager.wait([job.id for job in jobs.itervalues()])
 115         for i,job in jobs.iteritems():
 116             j = complete_jobs[job.id]
 117             assert j.status == 0, j.data
 118             if self._use_cache == True:
 119                 self._cache_run(d, i, j.data['stdout'])
 120             data[i] = parse(j.data['stdout'])
 121         del(jobs)
 122         del(complete_jobs)
 123
 124         # generate histogram
 125         events = []
 126         for d_i in data.values():
 127             events.extend([e.force for e in d_i
 128                            if e.initial_state == 'folded'])
 129         events = numpy.array(events)
 130         if bin_edges == None:
 131             return events
 132         h = Histogram()
 133         h.from_data(events, bin_edges)
 134         return h
 135
 136
 137 def main():
 138     import optparse
 139     import sys
 140
 141     global SAWSIM, CACHE_DIR
 142
 143     usage = "%prog [options] velocity_file"
 144
 145     parser = optparse.OptionParser(usage)
 146     parser.add_option("-s","--sawsim", dest="sawsim",
 147                       metavar="PATH",
 148                       help="Set sawsim binary (%default).",
 149                       default=SAWSIM)
 150     parser.add_option("-p","--params", dest="params",
 151                       metavar="PARAMS",
 152                       help="Initial params for fitting (%default).",
 153                       default=DEFAULT_PARAM_STRING)
 154     parser.add_option("-N", "--number-of-runs", dest="N",
 155                       metavar="INT", type='int',
 156                       help="Number of sawsim runs at each point in parameter space (%default).",
 157                       default=400)
 158     parser.add_option("-w", "--bin-width", dest="bin_width",
 159                       metavar="FLOAT", type='float',
 160                       help="Histogram bin width in newtons (%default).",
 161                       default=None)
 162     parser.add_option("-m", "--manager", dest="manager",
 163                       metavar="STRING",
 164                       help="Job manager name (one of %s) (%%default)."
 165                       % (', '.join(MANAGERS)),
 166                       default=MANAGERS[0])
 167     parser.add_option("-C","--use-cache", dest="use_cache",
 168                       help="Use cached simulations if they exist (vs. running new simulations) (%default)",
 169                       default=False, action="store_true")
 170     parser.add_option("--clean-cache", dest="clean_cache",
 171                       help="Remove previously cached simulations if they exist (%default)",
 172                       default=False, action="store_true")
 173     parser.add_option("-d","--cache-dir", dest="cache_dir",
 174                       metavar="STRING",
 175                       help="Cache directory for sawsim unfolding forces (%default).",
 176                       default=CACHE_DIR)
 177     options,args = parser.parse_args()
 178
 179     SAWSIM = options.sawsim
 180     CACHE_DIR = options.cache_dir
 181
 182     sh = SawsimHistogram(use_cache=options.use_cache,
 183                          clean_cache=options.clean_cache)
 184
 185     manager = get_manager(options.manager)()
 186     try:
 187         events = sh(param_string=options.params, N=options.N, manager=manager)
 188     finally:
 189         manager.teardown()
 190
 191     if options.bin_width == None:
 192         sys.stdout.write('# Unfolding force (N)\n')
 193         events.tofile(sys.stdout, sep='\n')
 194         sys.stdout.write('\n')
 195     else:
 196         h = Histogram()
 197         h.from_data(events, bin_edges=h.calculate_bin_edges(
 198                 events, options.bin_width))
 199         h.headings = ['Force (N)', 'Unfolding events']
 200         h.to_stream(sys.stdout)