1 # Copyright (C) 2009-2010 W. Trevor King <wking@drexel.edu>
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.
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.
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/>.
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.
20 from __future__ import with_statement
28 from . import __version__, log
29 from .histogram import Histogram
30 from .manager import InvokeJob, MANAGERS, get_manager
31 from .sawsim import parse
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 '
39 "-s 'unfolded,wlc,{0.39e-9,28e-9}' "
40 "-k 'folded,unfolded,bell,{3.3e-4,0.25e-9}' "
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
49 def _cache_dir(self, param_string):
51 >>> s = SawsimHistogram()
52 >>> s._cache_dir(DEFAULT_PARAM_STRING) # doctest: +ELLIPSIS
53 '/.../.sawsim_histogram/...'
56 CACHE_DIR, hashlib.sha256(param_string).hexdigest())
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))
64 def _load_cached_data(self, cache_dir, 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())
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:
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`.
83 If `bin_edges == None`, return a numpy array of all unfolding
86 Use the `JobManager` instance `manager` for asynchronous job
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`.
95 if self._use_cache == True:
96 d = self._cache_dir(param_string)
98 if self._clean_cache == True:
100 self._make_cache(param_string)
102 data = self._load_cached_data(d, N)
103 log().debug('using %d cached runs for %s'
104 % (len(data), param_string))
106 self._make_cache(param_string)
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'])
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:
133 h.from_data(events, bin_edges)
141 global SAWSIM, CACHE_DIR
143 usage = "%prog [options] velocity_file"
145 parser = optparse.OptionParser(usage)
146 parser.add_option("-s","--sawsim", dest="sawsim",
148 help="Set sawsim binary (%default).",
150 parser.add_option("-p","--params", dest="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).",
158 parser.add_option("-w", "--bin-width", dest="bin_width",
159 metavar="FLOAT", type='float',
160 help="Histogram bin width in newtons (%default).",
162 parser.add_option("-m", "--manager", dest="manager",
164 help="Job manager name (one of %s) (%%default)."
165 % (', '.join(MANAGERS)),
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",
175 help="Cache directory for sawsim unfolding forces (%default).",
177 options,args = parser.parse_args()
179 SAWSIM = options.sawsim
180 CACHE_DIR = options.cache_dir
182 sh = SawsimHistogram(use_cache=options.use_cache,
183 clean_cache=options.clean_cache)
185 manager = get_manager(options.manager)()
187 events = sh(param_string=options.params, N=options.N, manager=manager)
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')
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)