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
  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(params_string)
  60         os.makedirs(cache_dir)
  61         with open(os.path.join(d, '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(d, '%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(d, '%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                     '%s %s' % (SAWSIM, param_string)))
 114         complete_jobs = manager.wait([job.id for job in jobs])
 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(d_i)
 130         if bin_edges == None:
 131             return events
 132         h = Histogram()
 133         h.from_data(events, bin_edges)
 134         return h