--- /dev/null
+I've been spending some time comparing my [[force spectroscopy]] data
+with Marisa's, and the main problem has been normalizing the data
+collected using different systems. Marisa runs experiments using some
+home-grown [LabVIEW][] software, which saves the data in [[IGOR binary
+wave (IBW)|igor]] files. From my point of view, this is not the best
+approach, but it has been stable over the time we've been running
+experiments.
+
+I run my own experiments using my own code based on [[pyafm]], saving
+the data in [[HDF5]] files. I think this approach is much stronger,
+but it has been a bit of a moving target, and I've spent a good deal
+of my time working on the framework instead of running experiments.
+
+Both approaches save the digitized voltages that we read/write during
+an experiment, but other constants and calibration terms are not
+always recorded. My [[pyafm]]-based software has been gradually
+getting better in this regard, especially since I moved to
+[[h5config]]-based initialization. Anyhow, here's a quick runthough
+of all the important terms.
+
+For the [TL;DR][tldr] crowd, [[crunch.py]] is my [[Python]] module
+that crunches your input and prints out all the intermediate constants
+discussed below. To use it on your own data, you'll probably have to
+tweak the bits that read in the data to use your own format.
+
+<a id="piezo"></a>
+
+Calibrating the piezo
+=====================
+
+<a id="grid"></a>
+
+Calibration grid
+----------------
+
+We control the surface-tip distance by [[driving a piezo
+tube|pypiezo]]. We want to know the correspondence between the
+driving voltage and the piezo position, so we calibrate the piezo by
+imaging a sample with square pits of a known depth.
+
+[[!img piezo-good.png
+ alt="Piezo calibration trace"
+ title="Piezo calibration trace"]]
+
+In this trace, I swept the $x$ position across most of the range of my
+16-bit DAC. For each $x$ position, I adjusted the $z$ position until
+the measured cantilever deflection $d$ crossed a setpoint. This gives
+the $z$ voltage required to reach the surface as a function of $x$.
+You can see the 200 nm deep square pits, as well as a reasonable
+amount of $x$/$z$ crosstalk.
+
+Sometimes grid calibrations are even less convincing than the example
+shown above. For example, I have traces with odd shoulders on the
+sides of the pits:
+
+[[!img piezo.png
+ alt="Funky piezo calibration trace"
+ title="Funky piezo calibration trace"]]
+
+This is one of the better piezo calibration images I aquired for the
+piezo used in the pull below, so I'll use numbers from the funky
+imaging for the remainder of this analysis. This piezo has
+less $x$/$y$ range than the one used to measure the first calibration
+trace, which is why I was unable to aquire a full pit (or pits).
+
+The calibration constant $\alpha_z$ is the number of $z$ bits per
+depth meter:
+
+<a id="eq:alpha_z"></a>
+
+\[
+ \alpha_z = \frac{\text{grid piezo bits}}{\text{grid depth}}
+ = \frac{2779\pm87\text{ bits}}{200\text{ nm}}
+ = (1.39 \pm 0.04) \cdot 10^{10} \text{ bits/m} \;.
+\]
+
+Related conversion factors are
+
+<a id="eq:gamma_z"></a>
+<a id="eq:sigma_z"></a>
+
+\[
+\begin{aligned}
+ \gamma_z &= \frac{\text{piezo volts}}{\text{piezo bits}}
+ = \frac{20 \text{ piezo volts}}{\text{output volts}}
+ \cdot \frac{10 \text{ output volts}}{2^{15}\text{ piezo bits}}
+ = 6.10 \cdot 10^{-3} \text{ V/bit} \\
+ \sigma_z &= \frac{\text{grid piezo volts}}{\text{grid depth}}
+ = \gamma_z \alpha_z
+ = 8.48 \cdot 10^7 \text{ V/m} \;.
+\end{aligned}
+\]
+
+This is roughly in the ballpark of our piezo (serial number 2253E)
+which is spec'd at 8.96 nm/V along the $z$ axis, which comes out
+to $1.12\cdot10^8$ V/m.
+
+<a id="interference"></a>
+
+Laser interference
+------------------
+
+Another way to ballpark a piezo calibration that is much closer to a
+force spectroscopy pull is to use the laser interference pattern as a
+standard length. The laser for our [[microscope|MultiMode]] has a
+wavelength of 670 nm. We'll assume a geometric gain of
+
+\[
+ g_\lambda
+ = \frac{\text{increased laser path}}{\text{piezo displacement}}
+ \approx 2 \;.
+\]
+
+Measuring the length of an interference wave in bits then gives a
+standard equivalent to the known-depth pits of the calibration grid.
+
+[[!img pull-interference.png
+ alt="Laser interference during a velocity-clamp pull"
+ title="Laser interference during a velocity-clamp pull"]]
+
+\[
+\begin{aligned}
+ \alpha_z
+ &= \frac{\text{interference piezo bits}}{\text{interference depth}}
+ = \frac{g_\lambda}{\lambda} \cdot \text{interference piezo bits} \\
+ &= \frac{2}{670\cdot10^{-9}\text{ m}} \cdot (28935 - 23721)\text{ bits}
+ = 1.56 \cdot 10^{10}\text{ bits/m} \;.
+\end{aligned}
+\]
+
+The other piezo calibration parameters are found exactly as in the
+calibration grid case.
+
+\[
+ \sigma_z = \gamma_z \alpha_z
+ = 9.52 \cdot 10^7 \text{V/m} \;.
+\]
+
+which is fairly close to both the spec'd value and grid calibration
+values.
+
+<a id="photodiode"></a>
+
+Calibrating the photodiode
+==========================
+
+During experiments, we measure cantilever deflection via the top and
+bottom segments of a photodiode. We need to convert this deflection
+voltage into a deflection distance, so we'll use the
+already-calibrated piezo. When the tip is in contact with the
+surface, we can move the surface a known distance (using the piezo)
+and record the change in deflection.
+
+[[!img bump.png
+ alt="Surface bump for photodiode sensitivity"
+ title="Surface bump for photodiode sensitivity"]]
+
+The calibration constant $\alpha_d$ is the number of diode bits per
+piezo bit.
+
+<a id="eq:alpha_d"></a>
+
+\[
+ \alpha_d = \frac{\text{bump diode bits}}{\text{bump piezo bits}}
+ = 2.24 \text{ diode bits/piezo bits} \;.
+\]
+
+Related conversion factors are
+
+<a id="eq:gamma_d"></a>
+<a id="eq:sigma_d"></a>
+
+\[
+\begin{aligned}
+ \gamma_d &= \frac{\text{diode volts}}{\text{diode bits}}
+ = \frac{10 \text{ input volts}}{2^{15}\text{ diode bits}}
+ \cdot \frac{\text{diode volts}}{1 \text{input volts}}
+ = 3.05 \cdot 10^{-4} \text{ V/bit} \\
+ \sigma_d &= \frac{\text{bump diode volts}}{\text{bump tip position}}
+ = \gamma_d \alpha_d \alpha_z
+ = 9.52 \cdot 10^6 \text{ V/m} \;.
+\end{aligned}
+\]
+
+<a id="vibration"></a>
+
+Calibrating the cantilever spring constant
+==========================================
+
+To convert cantilever tip deflection to force, we need to know the
+spring constant of the cantilever. After bumping the surface, we move
+away from the surface and measure the cantilever's thermal vibration.
+We use the vibration data to calculate the spring constant [[using the
+equipartition theorem|calibcant]].
+
+[[!img vibration.png
+ alt="Thermal vibration measurement"
+ title="Thermal vibration measurement"]]
+
+The deflection variance $\langle d^2 \rangle$ is measured in frequency
+space, where the power spectral density (PSD) is fitted to the
+[[expected PSD of a damped harmonic oscillator|thesis]].
+
+\[
+\begin{aligned}
+ PSD_f &= \frac{G_{1f}}{(f_0^2 - f^2)^2 + \beta_f^2 f^2} \\
+ \langle d^2 \rangle &= \frac{\pi G_{1f}}{2\beta_f f_0^2} \\
+ \kappa &= \frac{2\beta_f f_0^2}{\pi G_{1f}}(\alpha_d\alpha_z)^2 k_B T
+ = \frac{2 \cdot (4.17 \cdot 10^3\text{ Hz})
+ \cdot (8.14 \cdot 10^3\text{ Hz})^2}
+ {\pi \cdot 3.10 \cdot 10^{13}\text{ bit}^2 \cdot\text{ Hz}^3}
+ \cdot (2.24 \cdot 1.39 \cdot 10^{10}\text{bit/m})^2
+ \cdot (1.38 \cdot 10^{-23}\text{ J/K}) \cdot (294 \text{K})
+ = 22.4 \cdot 10^{-3}\text{ N/m}
+\end{aligned}
+\]
+
+<a id="pull"></a>
+
+Analyzing a velocity-clamp pull
+===============================
+
+The raw data from a [[velocity-clamp pull|Force_spectroscopy]] is an
+array of output voltages used to sweep the piezo (moving the surface
+away from the cantilever tip) and an array of input voltages from the
+photodiode (measuring cantilever deflection). There are a number of
+interesting questions you can ask about such data, so I'll break this
+up into a few steps. Lets call the raw piezo data (in
+bits) $\alpha_v$, with the contact-kink located at $\alpha_{v0}$.
+Call the raw deflection data (also in bits) $\alpha_{F}$, with the
+contact-kink located at $\alpha_{F0}$.
+
+Piezo displacement
+------------------
+
+Using the [piezo calibration parameters](#piezo), we can calculate the
+raw piezo position using
+
+\[
+ z_\text{piezo} = \frac{\alpha_v - \alpha_{v0}}{\alpha_z} \;,
+\]
+
+measured in meters.
+
+Surface contact region
+----------------------
+
+During the initial portion of a velocity clamp pull, the cantilever
+tip is still in contact with the surface. This allows you to repeat
+the [photodiode calibration](#photodiode), avoiding problems due to
+drift in laser alignment or other geometric issues. This gives a new
+set of diode calibration parameters $\alpha_d'$ and $\sigma_d'$ (it is
+unlikely that $\gamma_d$ has changed, but it's easy to rework the
+following arguments to include $\gamma_d'$ if you feel that it might
+have changed).
+
+Tension
+-------
+
+We can use the new photodiode calibration and the [cantilever's spring
+constant](#vibration) to calculate the force from the Hookean
+cantilever:
+
+\[
+ F = \frac{\alpha_F - \alpha_{F0}}{\alpha_d' \alpha_z}\cdot\kappa
+ = (\alpha_F - \alpha_{F0})
+ \cdot\frac{2\beta_f f_0^2}{\pi G_{1f}}
+ \cdot\frac{\alpha_d^2\alpha_z}{\alpha_d'}
+ \cdot k_B T \;.
+\]
+
+[[!img pull-piezo.png
+ alt="Tension vs. piezo extension during a velocity-clamp pull"
+ title="Tension vs. piezo extension during a velocity-clamp pull"]]
+
+Protein extension
+-----------------
+
+As the piezo pulls on the cantilever/protein system, some of the
+increased extension is due to protein extension and the rest is due to
+cantilever extension. We can use Hooke's law to remove the cantilever
+extension, leaving only protein extension:
+
+\[
+ z_\text{protein} = z_\text{piezo} - \frac{F}{\kappa}
+ = \frac{1}{\alpha_z}
+ \left(\alpha_v - \alpha_{v0} - \frac{\alpha_F - \alpha_{F0}}{\alpha_d'}
+ \right)
+\]
+
+[[!img pull-protein.png
+ alt="Tension vs. protein extension during a velocity-clamp pull"
+ title="Tension vs. protein extension during a velocity-clamp pull"]]
+
+Contour-space
+-------------
+
+In order to confirm the piezo calibration, we look at changes in the
+unfolded contour length due to domain unfolding ($\Delta L$). There
+have been a number of studies of [[titin I27|I27-synthesis]], starting
+with [Carrion-Vazquez et al.][carrionvazquez99], that show an contour
+length increase of 28.1 ± 0.17 nm. Rather than fitting each loading
+region with a worm-like chain (WLC) or other polymer model, it is
+easier to calculate $\Delta L$ by converting the abscissa to
+contour-length space (following [Puchner et al.][puchner08]). While
+the WLC is commonly used, Puchner gets better fits using the [[freely
+rotating chain (FRC)|FRC]] model.
+
+In [[crunch.py]], I use either [Bustamante's formula][bustamante94]
+(WLC) or [Livadaru's equation 49][livadaru03] (FRC) to calculate the
+contour length of a particular force/distance pair.
+
+[[!img pull-contour-space.png
+ alt="Tension vs. unfolded contour length during a velocity-clamp pull"
+ title="Tension vs. unfolded contour length during a velocity-clamp pull"]]
+
+As you can see from the figure, my curves are spaced a bit too far
+appart. Because the contour length depends $F$ as well
+as $z_\text{protein}$, it depends on the cantilever calibration
+(via $\beta_f$, $f_0$, $G_{1f}$ and $\alpha_d$) as well as the piezo
+calibration (via $\alpha_z$). This makes adjusting calibration
+parameters in an attempt to match your $\Delta L$ with previously
+determined values a bit awkward. However, it is likely that my
+cantilever calibration is too blame. If we use the value of $\alpha_z$
+from the interference measurement we get
+
+[[!img pull-contour-space-interference.png
+ alt="Tension vs. unfolded contour length during a velocity-clamp pull,
+ using laser-interference to calibrate the piezo"
+ title="Tension vs. unfolded contour length during a velocity-clamp pull,
+ using laser-interference to calibrate the piezo"]]
+
+Which is still a bit too wide, but is much closer.
+
+[LabVIEW]: http://www.ni.com/labview/
+[tldr]: http://en.wiktionary.org/wiki/TL;DR
+[livadaru03]: http://dx.doi.org/10.1021/ma020751g
+[carrionvazquez99]: http://dx.doi.org/10.1073/pnas.96.20.11288
+[puchner08]: http://dx.doi.org/10.1529/biophysj.108.129999
+[bustamante94]: http://dx.doi.org/10.1126/science.8079175
+
+[[!tag tags/theory]]
--- /dev/null
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Copyright (C) 2012 W. Trevor King <wking@tremily.us>
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+"""Calculate all parameters associated with a velocity-clamp pull.
+
+This helps you track down where in the calibration process things may
+be going wrong.
+"""
+
+import collections as _collections
+import functools as _functools
+import os.path as _os_path
+
+try: # SciPy >= 0.10.0
+ from scipy.constants import Boltzmann as _kB
+except ImportError: # SciPy < 0.10.0
+ from scipy.constants import Bolzmann as _kB
+
+import h5py as _h5py
+
+# for some reason my scatterplot 'pixel' markers were taking up several pixels.
+import matplotlib.markers as _markers
+def _init(self, marker=None, fillstyle='full'):
+ self._fillstyle = fillstyle
+ self.set_marker(marker)
+ self.set_fillstyle(fillstyle)
+ if marker == ',':
+ self._transform = _markers.Affine2D().scale(0.1)
+_markers.MarkerStyle.__init__ = _init
+# end pixel hack
+
+import matplotlib.pyplot as _pyplot
+import matplotlib.widgets as _widgets
+import numpy as _numpy
+import scipy as _scipy
+
+import calibcant.calibrate as _calibcant_calibrate
+import calibcant.config as _calibcant_config
+import calibcant.bump_analyze as _calibcant_bump_analyze
+import calibcant.vibration_analyze as _calibcant_vibration_analyze
+import pypiezo.surface as _pypiezo_surface
+import FFT_tools as _FFT_tools
+
+## your data
+
+piezo_file = '~/tmp/rsrch/calibrate-piezo/data/calibrate-piezo-2012-02-06T16-23-30.dat'
+
+piezo_parameters = {
+ 'pits': [ # (z_bottom_bits, z_top_bits) pairs
+ (30044, 32736),
+ (30167, 33033),
+ ],
+ 'pit_depth': 200e-9, # meters
+ 'z_gain': 20, # volts applied at piezo / volts output from DAC
+ 'output_gamma_z': 10./2**15 # volts output from DAC / bit
+ }
+
+# cantilever calibration
+
+calibcant_file = '~/rsrch/data/calibrate_cantilever/2012-03-26T13-02-19-calibcant-data.h5'
+calibcant_bump = 3
+calibcant_vibration = 15
+deflection_parameters = {
+ 'd_gain': 1, # volts ADC input / volts photodiode differencer output
+ 'input_gamma_d': 10./2**15, # volts input at ADC / bit
+ }
+
+# unfolding pull
+
+pull_file = '~/rsrch/data/unfold-h5/good-2/2012-03-26T13-57-03.h5'
+_gamma = 22*_numpy.pi/180
+_l_over_b = _numpy.cos(_gamma/2.)/abs(_numpy.log(_numpy.cos(_gamma)))
+_b = 0.4e-9/_l_over_b
+pull_parameters = {
+ 'drift_slope': 0.06,
+ 'alpha_f_offset': -50,
+ 'smooth_steps': 600,
+ 'start_sawteeth': 50e-9, # meters from kink to start of sawteeth
+ 'stop_sawteeth': 300e-9, # meters from kink to end of sawteeth
+ 'dL': 29e-9, # meters, expected value for I27
+ 'contour_space_model': 'frc',
+ #'contour_space_model': 'wlc',
+ #'contour_space_parameters': None, # use defaults
+ 'contour_space_parameters': {'b': _b, 'gamma': _gamma},
+ 'contour_space_lines': _numpy.arange(75e-9, 350e-9, 29e-9),
+ }
+
+## readers for your data
+
+def load_piezo(filename=piezo_file):
+ data = _numpy.loadtxt(_os_path.expanduser(filename))
+ x = data[:,0]
+ z = data[:,1]
+ deflection = data[:,2]
+ return {'x': x, 'z': z, 'deflection': deflection}
+
+def load_bump(filename=calibcant_file, index=calibcant_bump):
+ calibrator,calibcant_data,calibcant_raw_data = (
+ _calibcant_calibrate.load_all(filename=_os_path.expanduser(filename)))
+ bump_data = calibcant_raw_data['bump'][index]
+ z = bump_data['raw']['z']
+ deflection = bump_data['raw']['deflection']
+ # HACK, 0 index not necessarily deflection
+ deflection_high_voltage_rail = (
+ calibrator.config['afm']['piezo']['inputs'][0]['maxdata'])
+ return {
+ 'z': z,
+ 'deflection': deflection,
+ 'deflection_high_voltage_rail': deflection_high_voltage_rail,
+ }
+
+def load_vibration(filename=calibcant_file, index=calibcant_vibration):
+ calibrator,calibcant_data,calibcant_raw_data = (
+ _calibcant_calibrate.load_all(filename=_os_path.expanduser(filename)))
+ vibration_data = calibcant_raw_data['vibration'][index]
+ temperature = calibcant_data['raw']['temperature'].mean()
+ deflection = vibration_data['raw']
+ config = vibration_data['config']['vibration']
+ psd_kwargs = {
+ 'chunk_size': config['chunk-size'],
+ 'overlap': config['overlap'],
+ 'window': config['window'],
+ }
+ fit_kwargs = {
+ 'min_frequency': config['minimum-fit-frequency'],
+ 'max_frequency': config['maximum-fit-frequency'],
+ 'offset': True,
+ }
+ #offset=config['model'] == _calibcant_config.OffsetBreitWigner)
+ return {
+ 'deflection': deflection,
+ 'temperature': temperature,
+ 'frequency': config['frequency'],
+ 'psd_kwargs': psd_kwargs,
+ 'fit_kwargs': fit_kwargs,
+ }
+
+def load_pull(filename=pull_file, with_temperature=True):
+ with _h5py.File(_os_path.expanduser(filename)) as f:
+ alpha_v = _numpy.array(f['unfold']['z'][...].flat, dtype=float)
+ alpha_f = _numpy.array(
+ f['unfold']['deflection'][...].flat, dtype=float)
+ if with_temperature:
+ temperature = float(f['environment']['temperature'][...])
+ else:
+ temperature = None
+ return {
+ 'alpha_v': alpha_v,
+ 'alpha_f': alpha_f,
+ 'temperature': temperature,
+ }
+
+
+## crunching
+
+class NamedObject (object):
+ def __init__(self, name=None):
+ self.name = name
+
+ def log(self, message):
+ print('{}: {}'.format(self.name, message))
+
+
+class SliderFigure (NamedObject):
+ def __init__(self, figure, axes=None, parameters=None, data=None,
+ callback=None, **kwargs):
+ super(SliderFigure, self).__init__(**kwargs)
+ self.figure = figure
+ if axes is None:
+ axes = list(figure.axes)
+ self.axes = axes
+ self.parameters = parameters
+ self.data = data
+ self.callback = callback
+
+ self._slider_height = 0.03
+ self._slider_margin = 0.02
+ self._total_slider_height = len(parameters)*(
+ self._slider_height+self._slider_margin) + self._slider_margin
+
+ self._shift_original_axes_up()
+ self._add_sliders()
+
+ @staticmethod
+ def _map_original_axes_y(y, new_bottom):
+ """map [0,1] -> [new_bottom,1]
+
+ >>> SliderFigure._map_original_axes_y(0, 0.5)
+ 0.5
+ >>> SliderFigure._map_original_axes_y(0.5, 0.4)
+ 0.7
+ """
+ return 1 - (1-y)*(1-new_bottom)
+
+ def _shift_original_axes_up(self):
+ for axes in self.figure.axes:
+ p = axes.get_position()
+ p.y0 = self._map_original_axes_y(p.y0, self._total_slider_height)
+ p.y1 = self._map_original_axes_y(p.y1, self._total_slider_height)
+ axes.set_position(p)
+
+ def _add_sliders(self):
+ self._sliders = {}
+ n = len(self.parameters)
+ h = self._slider_height
+ m = self._slider_margin
+ left = 0.4
+ w = 0.4
+ for i,(k,v) in enumerate(self.parameters.items()):
+ try:
+ x = float(v)
+ except (TypeError, ValueError):
+ continue # non-float parameter
+ axes = self.figure.add_axes([left, m+i*(h+m), w, h])
+ min,max = sorted([v/4., 4*v]) # v might be negative
+ slider = _widgets.Slider(
+ axes, k, min, max, valinit=v, valfmt='%1.5g',
+ closedmin=False, closedmax=False)
+ slider.on_changed(_functools.partial(
+ self.on_changed, parameter=k, slider=slider))
+ self._sliders[k] = slider
+
+ def on_changed(self, value, parameter, slider=None):
+ if slider is None: # external change (not direct slider movement)
+ slider = self._sliders[parameter]
+ slider.set_val(value)
+ return
+ self.log('changed {} to {}, recalculating...'.format(
+ parameter, value))
+ self.parameters[parameter] = value
+ if self.callback:
+ self.callback(
+ figure=self.figure, axes=self.axes, data=self.data,
+ parameters=self.parameters)
+ self.log('finished recalculating change of {} to {}'.format(
+ parameter, value))
+
+
+class Experiment (NamedObject):
+ def __init__(self, load_data=None, input_parameters=None,
+ output_parameters=None, output_parameter_units=None,
+ descendants=None, **kwargs):
+ super(Experiment, self).__init__(**kwargs)
+ self.load_data = load_data
+ if input_parameters is None:
+ input_parameters = tuple()
+ self.input_parameters = input_parameters
+ if output_parameters is None:
+ output_parameters = tuple()
+ self.output_parameters = output_parameters
+ if output_parameter_units is None:
+ output_parameter_units = {}
+ self.output_parameter_units = output_parameter_units
+ if descendants is None:
+ descendants = []
+ self.descendants = descendants
+ self._iparameters = self._oparameters = None
+
+ def crunch(self, parameters=None):
+ if self.load_data and (not hasattr(self, 'data') or self.data is None):
+ self.data = self.load_data()
+ self._iparameters = dict(
+ (k,parameters[k]) for k in self.input_parameters)
+ self._iparameters.update(self.data)
+ oparameters = self._crunch(**self._iparameters)
+ return self._crunch_callback(oparameters)
+
+ @staticmethod
+ def _crunch(**kwargs):
+ return {}
+
+ def _crunch_callback(self, oparameters):
+ if self._oparameters is not None:
+ for k in self.output_parameters:
+ if self._oparameters[k] != oparameters[k]:
+ self.log('{} changed to {} during recrunch'.format(
+ k, oparameters[k]))
+ for d in self.descendants:
+ if k in d.input_parameters:
+ self.log('{} depends on {}'.format(d.name, k))
+ if d.slider_figure:
+ #slider.set_val(value)
+ d.slider_figure.on_changed(
+ value=oparameters[k], parameter=k)
+ self._oparameters = oparameters
+ return dict((k,self._oparameters[k]) for k in self.output_parameters)
+
+ def plot(self, save=False, slider=False, filename=None, **kwargs):
+ """Wrapper around plotting functions for consistent handling."""
+ if save and not filename:
+ filename = '{}.png'.format(self.name)
+ assert self._oparameters is not None, 'crunch before plotting!'
+ kwargs.update(self._iparameters)
+ kwargs.update(self._oparameters)
+ kwargs['log'] = self.log
+ kwargs['slider'] = slider
+ figure = self._plot(**kwargs)
+ if filename:
+ figure.savefig(filename)
+ if slider:
+ parameters = dict(
+ (k,self._iparameters[k]) for k in self.input_parameters)
+ parameters['_crunch'] = self._crunch
+ parameters['_crunch_callback'] = self._crunch_callback
+ self.slider_figure = SliderFigure(
+ name='{}-slider'.format(self.name), figure=figure,
+ data=self.data, parameters=parameters, callback=self._replot)
+ return figure
+
+ @staticmethod
+ def _plot(**kwargs):
+ raise NotImplementedError()
+
+ @staticmethod
+ def _replot(figure, axes, data=None, parameters=None):
+ parameters = dict(parameters)
+ crunch = parameters.pop('_crunch')
+ crunch_callback = parameters.pop('_crunch_callback')
+ parameters.update(data)
+ oparameters = crunch(**parameters)
+ crunch_callback(oparameters)
+ return oparameters
+
+ # utilities
+
+ @staticmethod
+ def smooth(x, window_len=11):
+ """
+ >>> Experiment.smooth([1,2,3,4,5], window_len=2)
+ array([ 1. , 1.5, 2.5, 3.5, 4.5])
+
+ With the padding, the convolution was over `[1, 1, 2, 3, 4, 5]`.
+
+ >>> Experiment.smooth([1,2,3,4,5], window_len=3)
+ array([ 1.33333333, 2. , 3. , 4. , 4.66666667])
+
+ With the padding, the convolution was over `[1, 1, 2, 3, 4, 5, 5]`.
+ """
+ assert window_len > 1, window_len
+ window = _numpy.ones(window_len, 'd') # moving average
+ front_pad = int(window_len/2) # truncate
+ if int(window_len) % 2 == 0:
+ back_pad = front_pad - 1
+ else:
+ back_pad = front_pad
+ front_padding = [x[0]] * front_pad
+ back_padding = [x[-1]] * back_pad
+ x = _numpy.concatenate((front_padding, x, back_padding))
+ return _numpy.convolve(window/window.sum(), x, mode='valid')
+
+ @staticmethod
+ def true_set(dictionary):
+ """
+ >>> from collections import defaultdict
+ >>> Experiment.true_set({})
+ False
+ >>> Experiment.true_set({'a': True})
+ True
+ >>> Experiment.true_set(defaultdict(lambda: True))
+ True
+ >>> Experiment.true_set(defaultdict(lambda: 'abc'))
+ False
+ """
+ if True in dictionary.values():
+ return True
+ if isinstance(dictionary, _collections.defaultdict):
+ if dictionary.default_factory() is True:
+ return True
+ return False
+
+
+class Piezo (Experiment):
+ def __init__(self, name='piezo', *args, **kwargs):
+ kwargs['input_parameters'] = [
+ 'pits', 'pit_depth', 'z_gain', 'output_gamma_z']
+ kwargs['output_parameters'] = [
+ 'alpha_z', 'delta_alpha_z', 'gamma_z', 'sigma_z']
+ kwargs['output_parameter_units'] = {
+ 'alpha_z': 'bit/m',
+ 'delta_alpha_z': 'bit/m',
+ 'gamma_z': 'V/bit',
+ 'sigma_z': 'V/m',
+ }
+ super(Piezo, self).__init__(name=name, *args, **kwargs)
+
+ @staticmethod
+ def _crunch(pits, pit_depth, z_gain, output_gamma_z, **kwargs):
+ ret = {}
+ ret['dz'] = _numpy.array([(top-bot) for bot,top in pits], dtype=float)
+ ret['alpha_z'] = ret['dz'].mean()/pit_depth
+ ret['delta_alpha_z'] = ret['dz'].std()/pit_depth
+ ret['gamma_z'] = z_gain * output_gamma_z
+ ret['sigma_z'] = ret['gamma_z'] * ret['alpha_z']
+ return ret
+
+ @staticmethod
+ def _plot(x, z, deflection, **kwargs):
+ figure = _pyplot.figure()
+ axes = figure.add_subplot(1, 1, 1)
+ color_dimension = 'deflection'
+ color_dimension = 'trace'
+ #color_dimension = 'debug trace'
+ if color_dimension in ['trace', 'debug trace']:
+ trace = _numpy.zeros(x.shape, dtype=bool)
+ direction = None
+ for i,xi in enumerate(x):
+ last_x = x[max(0, i-1)]
+ new_direction = xi > last_x
+ if direction is None and xi != last_x: # first movement
+ trace[0:i] = direction = new_direction
+ else: # subsequent movement
+ if xi != last_x:
+ direction = new_direction
+ trace[i] = direction
+ if color_dimension == 'debug trace':
+ axes.plot(x, 'b.', 2**15*trace, 'r.-')
+ elif color_dimension == 'trace':
+ scatter = axes.scatter(
+ x, z, c=trace, marker=',', edgecolors='face')
+ else: # color_dimension == 'deflection'
+ scatter = axes.scatter(
+ x, z, c=deflection, marker=',',
+ edgecolors='face')
+ if color_dimension in ['trace', 'deflection']:
+ colorbar = figure.colorbar(scatter)
+ if color_dimension == 'trace':
+ colorbar.set_label('trace/retrace')
+ elif color_dimension == 'deflection':
+ colorbar.set_label('deflection (bits)')
+ axes.set_title('piezo calibration')
+ axes.set_xlabel('x piezo (bits)')
+ axes.set_ylabel('z piezo (bits)')
+ return figure
+
+
+class Bump (Experiment):
+ def __init__(self, name='bump', *args, **kwargs):
+ kwargs['input_parameters'] = [
+ 'input_gamma_d', 'd_gain', 'alpha_z']
+ kwargs['output_parameters'] = [
+ 'alpha_d', 'gamma_d', 'sigma_d']
+ kwargs['output_parameter_units'] = {
+ 'alpha_d': 'diode_bit/piezo_bit',
+ 'gamma_d': 'V/bit',
+ 'sigma_d': 'V/m',
+ }
+ super(Bump, self).__init__(name=name, *args, **kwargs)
+
+ @staticmethod
+ def _crunch(z, deflection, deflection_high_voltage_rail,
+ input_gamma_d, d_gain, alpha_z):
+ parameters = _calibcant_bump_analyze.fit(
+ z, deflection, high_voltage_rail=deflection_high_voltage_rail,
+ param_guesser=_calibcant_bump_analyze.limited_linear_param_guess,
+ model=_calibcant_bump_analyze.limited_linear,
+ sensitivity_from_fit_params=lambda parameters: parameters)
+ ret = {}
+ ret['deflection_fit'] = _calibcant_bump_analyze.limited_linear(
+ z, parameters, high_voltage_rail=deflection_high_voltage_rail)
+ ret['alpha_d'] = parameters[2] # bump diode bits / bump piezo bit
+ ret['gamma_d'] = input_gamma_d / d_gain
+ ret['sigma_d'] = ret['gamma_d'] * ret['alpha_d'] * alpha_z
+ return ret
+
+ @staticmethod
+ def _plot(alpha_d, gamma_d, sigma_d, z, deflection, deflection_fit,
+ **kwargs):
+ figure = _pyplot.figure()
+ axes = figure.add_subplot(1, 1, 1)
+ axes.plot(z, deflection, 'b.',
+ z, deflection_fit, 'r-')
+ axes.set_title('photodiode calibration')
+ axes.set_xlabel('z piezo (bits)')
+ axes.set_ylabel('deflection (bits)')
+ return figure
+
+ @staticmethod
+ def _replot(figure, axes, data=None, parameters=None):
+ oparameters = Experiment._replot(
+ figure=figure, axes=axes, data=data, parameters=parameters)
+ deflection = data['deflection']
+ z = data['z']
+ deflection_fit = oparameters['deflection_fit']
+ axes[0].lines[0].set_data(z, deflection)
+ axes[0].lines[1].set_data(z, deflection_fit)
+ _pyplot.draw()
+
+
+class Vibration (Experiment):
+ def __init__(self, name='vibration', *args, **kwargs):
+ kwargs['input_parameters'] = [
+ 'alpha_z', 'alpha_d']
+ kwargs['output_parameters'] = [
+ 'vibration_temperature', 'vibration_variance', 'f_0', 'beta_f',
+ 'G_1f', 'Q', 'kappa']
+ kwargs['output_parameter_units'] = {
+ 'vibration_temperature': 'K',
+ 'vibration_variance': 'bit^2',
+ 'f_0': 'Hz',
+ 'beta_f': 'Hz',
+ 'G_1f': 'bit^2 Hz^3',
+ 'Q': '',
+ 'kappa': 'N/m',
+ }
+ super(Vibration, self).__init__(name=name, *args, **kwargs)
+
+ @staticmethod
+ def _crunch(deflection, temperature, frequency, alpha_z, alpha_d,
+ psd_kwargs, fit_kwargs):
+ freq_axis,power = _FFT_tools.unitary_avg_power_spectrum(
+ (deflection - deflection.mean()), frequency, **psd_kwargs)
+ f,beta,G,offset = _calibcant_vibration_analyze.fit_psd(
+ freq_axis, power, **fit_kwargs)
+ power_fit = _calibcant_vibration_analyze.breit_wigner(
+ freq_axis, f, beta, G, offset)
+ variance = _calibcant_vibration_analyze.breit_wigner_area(f, beta, G)
+ kappa = 1/variance * (alpha_d*alpha_z)**2 * _kB*temperature
+ return {
+ 'vibration_temperature': temperature,
+ 'freq_axis': freq_axis,
+ 'power': power,
+ 'f_0': f,
+ 'beta_f': beta,
+ 'G_1f': G,
+ 'Q': f/beta,
+ 'offset': offset,
+ 'power_fit': power_fit,
+ 'vibration_variance': variance,
+ 'kappa': kappa,
+ }
+
+ @staticmethod
+ def _plot(deflection, freq_axis, power, power_fit, frequency, **kwargs):
+ figure = _pyplot.figure()
+ axes = figure.add_subplot(2, 1, 1)
+ max_time = len(deflection)/frequency
+ time = _numpy.linspace(0, max_time, len(deflection))
+ axes.plot(time, deflection, '.')
+ axes.set_xbound([0, max_time])
+ axes.set_title('cantilever calibration')
+ axes.set_xlabel('time (s)')
+ axes.set_ylabel('deflection (bit)')
+
+ axes = figure.add_subplot(2, 1, 2)
+ axes.hold(True)
+ axes.set_yscale('log')
+ axes.plot(freq_axis, power, 'b.')
+ xmin,xmax = axes.get_xbound()
+ ymin,ymax = axes.get_ybound()
+ axes.plot(freq_axis, power_fit, 'r-')
+ axes.axis([xmin, xmax, ymin, ymax])
+ axes.set_xlabel('frequency (Hz)')
+ axes.set_ylabel('deflection PSD (bit^2/Hz)')
+ return figure
+
+ @staticmethod
+ def _replot(figure, axes, data=None, parameters=None):
+ oparameters = Experiment._replot(
+ figure=figure, axes=axes, data=data, parameters=parameters)
+ deflection = data['deflection']
+ frequency = data['frequency']
+ freq_axis = oparameters['freq_axis']
+ power = oparameters['power']
+ power_fit = oparameters['power_fit']
+ max_time = len(deflection)/frequency
+ time = _numpy.linspace(0, max_time, len(deflection))
+ axes[0].lines[0].set_data(time, deflection)
+ axes[1].lines[0].set_data(freq_axis, power)
+ axes[1].lines[1].set_data(freq_axis, power_fit)
+ _pyplot.draw()
+
+
+class Pull (Experiment):
+ def __init__(self, name='pull', *args, **kwargs):
+ kwargs['input_parameters'] = [
+ 'alpha_z', 'alpha_d', 'vibration_variance',
+ 'vibration_temperature', 'drift_slope', 'alpha_f_offset',
+ 'smooth_steps', 'start_sawteeth', 'stop_sawteeth',
+ 'contour_space_model', 'contour_space_parameters',
+ 'contour_space_lines']
+ kwargs['output_parameters'] = [
+ 'pull_alpha_d', 'pull_temperature']
+ kwargs['output_parameter_units'] = {
+ 'pull_alpha_d': 'diode_bit/piezo_bit',
+ 'pull_temperature': 'K',
+ }
+ super(Pull, self).__init__(name=name, *args, **kwargs)
+
+ @staticmethod
+ def _crunch(alpha_v, alpha_f, temperature,
+ alpha_z, alpha_d, vibration_variance, vibration_temperature,
+ drift_slope, alpha_f_offset, start_sawteeth, stop_sawteeth,
+ contour_space_model, contour_space_parameters,
+ **kwargs):
+ alpha_f = alpha_f - (alpha_v-alpha_v.min())*drift_slope
+ ddict = {
+ 'approach': { # for _analyze_surface_position_data
+ 'z': alpha_v[::-1],
+ 'deflection': alpha_f[::-1],
+ }
+ }
+ (non_contact_offset,non_contact_slope,kink,alpha_dp
+ ) = _pypiezo_surface.analyze_surface_position_data(
+ ddict, return_all_parameters=True)
+ # The fitted kink offset, but thrown off by protein presense
+ #alpha_f0 = non_contact_offset + non_contact_slope*(kink-alpha_v.min())
+ alpha_f0 = alpha_f[-100:].mean()
+ alpha_v0 = alpha_v.max() - (alpha_f.max()-alpha_f0)/alpha_dp
+ d_alpha_f = -(alpha_f - alpha_f0 + alpha_f_offset)
+ d_alpha_v = -(alpha_v - alpha_v0)
+ kBT = _kB * vibration_temperature
+ f = d_alpha_f/vibration_variance*(alpha_d**2 * alpha_z)/alpha_dp * kBT
+ z_piezo = 1./alpha_z * d_alpha_v
+ z = 1./alpha_z * (d_alpha_v - d_alpha_f/alpha_dp)
+ z_contour = contour_values(
+ z, f, temperature=temperature, model=contour_space_model,
+ parameters=contour_space_parameters)
+ in_sawteeth = (z > start_sawteeth) * (z < stop_sawteeth)
+ z_sawteeth = z[in_sawteeth]
+ f_sawteeth = f[in_sawteeth]
+ return {
+ 'pull_alpha_d': alpha_d,
+ 'pull_temperature': temperature,
+ 'z': z,
+ 'f': f,
+ 'z_piezo': z_piezo,
+ 'z_contour': z_contour,
+ 'z_sawteeth': z[in_sawteeth],
+ 'f_sawteeth': f[in_sawteeth],
+ }
+
+ @staticmethod
+ def _plot(z_piezo, z, z_contour, deflection=None, f=None,
+ pull_temperature=None, smooth_steps=None,
+ xlabel='z piezo (m)', ylabel='deflection (N)',
+ bounds=None, log=None, slider=None,
+ contour_space_model=None, contour_space_parameters=None,
+ contour_space_lines=None, **kwargs):
+ """Generate pull-plots with three possible abscissa.
+
+ The contour-space plots are following Puchner et al., 2008.
+ http://dx.doi.org/10.1529/biophysj.108.129999
+ """
+ if deflection is None:
+ assert f is not None, "must set 'deflection' or 'f'"
+ deflection = f
+ else:
+ assert f is None, "cannot set both 'f' and 'deflection'"
+ figure = _pyplot.figure()
+ axes = figure.add_subplot(1, 1, 1)
+ axes._z_data = {
+ 'z piezo (m)': z_piezo,
+ 'protein extension (m)': z,
+ 'contour-space (m)': z_contour,
+ }
+ try:
+ zp = axes._z_data[xlabel]
+ except KeyError as e:
+ raise NotImplementedError(xlabel) # from e (for Python 3)
+ plots = []
+ plots.extend(
+ axes.plot(
+ zp, deflection,
+ color='b', linestyle='none', marker=',', markeredgecolor='b'))
+ if smooth_steps:
+ window_len = int(len(z)/smooth_steps)
+ axes.hold(True)
+ z_smooth = Experiment.smooth(z, window_len)
+ d_smooth = Experiment.smooth(deflection, window_len)
+ try:
+ c_smooth = contour_values(
+ z_smooth, d_smooth, temperature=pull_temperature,
+ model=contour_space_model,
+ parameters=contour_space_parameters)
+ except ValueError: # avoid crashing figure.py
+ c_smooth = None
+ axes._smooth_z_data = {
+ 'z piezo (m)': Experiment.smooth(z_piezo, window_len),
+ 'protein extension (m)': z_smooth,
+ 'contour-space (m)': c_smooth,
+ }
+ zp_smooth = axes._smooth_z_data[xlabel]
+ if xlabel in ['z piezo (m)', 'protein extension (m)']:
+ plots.extend(axes.plot(
+ zp_smooth, d_smooth, color='r', markeredgecolor='r'))
+ else: # 'contour-space (m)'
+ plots.extend(axes.plot(
+ zp_smooth, d_smooth,
+ color='r', linestyle='none', marker=',',
+ markeredgecolor='r'))
+ if xlabel == 'contour-space (m)' and contour_space_lines is not None:
+ for zi in contour_space_lines:
+ axes.axvline(zi, color='g')
+ if bounds:
+ axes.axis(bounds)
+ axes.set_title('velocity clamp unfolding')
+ axes.set_xlabel(xlabel)
+ axes.set_ylabel(ylabel)
+
+ if slider: # add radio buttons for plot-type selection
+ figure.subplots_adjust(left=0.5)
+ left = 0.05
+ bottom = 0.2
+ width = 0.3
+ height = 0.6
+ rax = figure.add_axes(
+ [left, bottom, width, height], axisbg='none')
+ labels = (
+ 'z piezo (m)',
+ 'protein extension (m)',
+ 'contour-space (m)',
+ )
+ radio = _widgets.RadioButtons(
+ rax, labels=labels, active=labels.index(xlabel))
+ def change_type(label): # closure; upvalues: log, plots, axes
+ log('change plot to {}'.format(label))
+ axes.set_xlabel(label)
+ for i,plot in enumerate(plots):
+ old_z,deflection = plot.get_data()
+ if i == 0: # regular plot
+ z = axes._z_data[label]
+ elif i == 1: # smooth plot
+ if label in ['z piezo (m)', 'protein extension (m)']:
+ plots[-1].set_linestyle('-')
+ plots[-1].set_marker(None)
+ else: # 'contour-space (m)'
+ plots[-1].set_linestyle('None')
+ plots[-1].set_marker(',')
+ z = axes._smooth_z_data[label]
+ else:
+ break
+ plot.set_data(z, deflection)
+ _pyplot.draw()
+ log('changed plot to {}'.format(label))
+ axes._callback = change_type # avoid garbage collection
+ axes._radio = radio # avoid garbage collection
+ radio.on_clicked(change_type)
+ return figure
+
+ @staticmethod
+ def _replot(figure, axes, data=None, parameters=None):
+ oparameters = Experiment._replot(
+ figure=figure, axes=axes, data=data, parameters=parameters)
+ axes[0]._z_data = {
+ 'z piezo (m)': oparameters['z_piezo'],
+ 'protein extension (m)': oparameters['z'],
+ 'contour-space (m)': oparameters['z_contour'],
+ }
+ if parameters.get('smooth_steps', None):
+ window_len = int(len(oparameters['z'])/parameters['smooth_steps'])
+ z_smooth = Experiment.smooth(oparameters['z'], window_len)
+ d_smooth = Experiment.smooth(oparameters['f'], window_len)
+ axes[0]._smooth_z_data = {
+ 'z piezo (m)': Experiment.smooth(oparameters['z_piezo']),
+ 'protein extension (m)': z_smooth,
+ 'contour-space (m)': contour_values(
+ z_smooth, d_smooth,
+ temperature=oparameters['pull_temperature'],
+ model=parameters['contour_space_model'],
+ parameters=parameters['contour_space_parameters']),
+ }
+ for i,plot in enumerate(axes[0].lines):
+ z,deflection = plot.get_data()
+ if i == 0: # regular plot
+ deflection = oparameters['f']
+ elif i == 1: # smooth plot
+ deflection = d_smooth
+ plot.set_data(z, deflection)
+ axes[0]._callback(axes[0].get_xlabel())
+
+
+def inverse_wlc(F):
+ """
+ f = kB T / p ( 1/4 * (1/(1-z/L)**2 - 1) + z/L )
+
+ With the replacements x = z/L and F = fp/(kB T), this gives
+
+ F = ( 1/4 * (1/(1-x)**2 - 1) + x )
+
+ After some math (see the sawsim manual for details), we get
+
+ 0 = x**3 - (F+9/4) x**2 + (2F+3/2) x - F
+
+ To get x, just solve this cubic equation.
+
+ >>> def wlc(x):
+ ... return (0.25*(1/(1-x)**2 - 1) + x)
+ >>> F = wlc(0.8)
+ >>> '{:g}'.format(F)
+ '6.8'
+ >>> '{:g}'.format(inverse_wlc(F))
+ '0.8'
+ """
+ if F < 0:
+ return 0
+ roots = _scipy.roots([1, -(F+9./4), (2*F+3./2), -F])
+ for r in roots:
+ if r.imag == 0 and r.real > 0 and r.real < 1:
+ return r.real
+ raise ValueError((F, roots))
+
+def inverse_frc(F, gamma, c=2, beta=2, g=float('inf')):
+ """
+ From Livadaru et al., 2003, equations 22, 48, and 49 [1].
+
+ l = b*cos(γ/2)/|log(cos(γ))| (eqn. 22)
+ F_WLC[x] = 3/4 - 1/x + x^2/4 (eqn. 48)
+ x = 1 - {(F_WLC^{-1}[Fl/b])^β + (cF)^β}^{-1/β} + F/g (eqn. 49)
+
+ where b is the segment length, γ is the fixed bond angle, x=R_z/L,
+ and F = fb/kBT. They use l for persistence length.
+
+ They also use f/g for chain extension instead of F/g, but it's
+ easy enough to make appropriate adjustments to g, and using F
+ keeps temperature out of the formula. Combining equations 22 and
+ 49 to minimize constants:
+
+ x = 1 - {(F_WLC^{-1}[F cos(γ/2)/|log(cos(γ))|])^β + (cF)^β}^{-1/β} + F/g
+
+ The F_WLC form (48) can be used to reconstruct Bustamante's more
+ familiar version [2] via Livadaru's equation 47:
+
+ fl/(kBT) = F_WLC[(1-R_z/L)^{-1}]
+
+ I found Livadaru via Puchner et al., 2008, equation 3 [3]. Note
+ that there are typos in Livadaru's equation 46 which were
+ propogated in Puchner's equation 3. It should have read
+
+ { fa/(3kBT) for fb/(kBT) < b/l
+ R_z/L ~= { 1-(4fl/(kBT))^{-1/2} for b/l < fb/(kBT) < l/b
+ { 1-(cfb/(kBT))^{-1} for l/b < fb/(kBT)
+
+ Also, the form's low-force regime is missing a prefactor to allow
+ it to match up with the mid-force form at the crossover point.
+ They avoid this issue by using Bustamante's WLC interpolation
+ formula [2] for both the low- and mid-force portions in equation
+ 49.
+
+ Check against values from Livadaru's figure 14.
+
+ >>> for gamma in [1, 20, 70]:
+ ... for F in [0.01, 0.1, 1, 10]:
+ ... x = inverse_frc(F=F, gamma=gamma*_numpy.pi/180)
+ ... print(('(1-Rz/L)^(-1) = {:g} '
+ ... '(for gamma={:g} deg, F={:g})').format(
+ ... 1./(1-x), gamma, F))
+ (1-Rz/L)^(-1) = 16.1199 (for gamma=1 deg, F=0.01)
+ (1-Rz/L)^(-1) = 51.2165 (for gamma=1 deg, F=0.1)
+ (1-Rz/L)^(-1) = 162.053 (for gamma=1 deg, F=1)
+ (1-Rz/L)^(-1) = 512.833 (for gamma=1 deg, F=10)
+ (1-Rz/L)^(-1) = 1.11106 (for gamma=20 deg, F=0.01)
+ (1-Rz/L)^(-1) = 2.26794 (for gamma=20 deg, F=0.1)
+ (1-Rz/L)^(-1) = 8.05245 (for gamma=20 deg, F=1)
+ (1-Rz/L)^(-1) = 32.1006 (for gamma=20 deg, F=10)
+ (1-Rz/L)^(-1) = 1.0053 (for gamma=70 deg, F=0.01)
+ (1-Rz/L)^(-1) = 1.07101 (for gamma=70 deg, F=0.1)
+ (1-Rz/L)^(-1) = 2.56046 (for gamma=70 deg, F=1)
+ (1-Rz/L)^(-1) = 20.6952 (for gamma=70 deg, F=10)
+
+ [1]: http://dx.doi.org/10.1021/ma020751g
+ [2]: http://dx.doi.org/10.1126/science.8079175
+ [3]: http://dx.doi.org/10.1529/biophysj.108.129999
+ """
+ if F < 0:
+ return 0
+ l_over_b = _numpy.cos(gamma/2.)/abs(_numpy.log(_numpy.cos(gamma)))
+ xWLC = inverse_wlc(F*l_over_b) # R_z/L for a WLC
+ invFWLC = 1./(1 - xWLC) # F_WLC^{-1}[Fl/b] (see eqn.s 47 and 48)
+ return 1 - (invFWLC**beta + (c*F)**beta)**(-1./beta) + F/g
+
+def contour_values(z, deflection, temperature=300., model='wlc',
+ parameters=None):
+ """Convert (z, f) points to (L, f) space using a polymer model.
+
+ Default parameters are from Puchner et al., 2008 [1]. The FJC
+ parameters are corrected to match the corrected version of
+ Livadaru et al., 2003 [2].
+
+ On Wed, Sep 12, 2012 at 01:37:23AM +0000, Puchner, Georg Elias Michael wrote:
+ > Hello Trevor,
+ >
+ > thanks for pointing this out, I wasn't aware of this typo and
+ > also think they should have published and erratum.
+ >
+ > I dug into my old data and code (I am not doing force
+ > spectroscopy andy more since some years) and implemented your
+ > corrected version. I attached a comparis on of my old
+ > transformation equation (black line with the parameters γ=22°
+ > and b=0.4 nm as published) and the fixed equation (dashed red
+ > line with γ=41° and b=0.11 nm)
+ >
+ > [snip extra details]
+
+ [1]: http://dx.doi.org/10.1529/biophysj.108.129999
+ [2]: http://dx.doi.org/10.1021/ma020751g
+ """
+ if parameters is None:
+ parameters = {}
+ if model == 'wlc':
+ p = parameters.get('p', 0.4e-9) # meters, persistence length
+ d_scale = p/(_kB*temperature)
+ elif model == 'frc':
+ b = parameters.get('b', 0.11e-9) # meters, segment length
+ gamma = parameters.get(
+ 'gamma', 41*_numpy.pi/180) # radians, fixed bond angle
+ d_scale = b/(_kB*temperature)
+ else:
+ raise NotImplementedError(model)
+ L_values = _numpy.zeros(z.shape, dtype=float)
+ for i,(z_i,d_i) in enumerate(zip(z, deflection)):
+ if model == 'wlc':
+ x = inverse_wlc(F=d_i*d_scale)
+ else:
+ x = inverse_frc(F=d_i*d_scale, gamma=gamma)
+ if x == 0:
+ L_values[i] = 0
+ else:
+ L_values[i] = z_i/x # x = z/L
+ return L_values
+
+def load_experiments(load_piezo=None, piezo_parameters=None,
+ load_bump=None, bump_parameters=None,
+ load_vibration=None, vibration_parameters=None,
+ load_pull=None, pull_parameters=None,
+ override_parameters={}):
+ if load_piezo is None:
+ load_piezo = globals()['load_piezo']
+ if piezo_parameters is None:
+ piezo_parameters = globals()['piezo_parameters']
+ if load_bump is None:
+ load_bump = globals()['load_bump']
+ if bump_parameters is None:
+ bump_parameters = globals()['deflection_parameters']
+ if load_vibration is None:
+ load_vibration = globals()['load_vibration']
+ if vibration_parameters is None:
+ vibration_parameters = globals()['deflection_parameters']
+ if load_pull is None:
+ load_pull = globals()['load_pull']
+ if pull_parameters is None:
+ pull_parameters = globals()['pull_parameters']
+
+ parameters = {}
+ experiments = []
+ for experiment,external_parameters in [
+ (Piezo(load_data=load_piezo), piezo_parameters),
+ (Bump(load_data=load_bump), deflection_parameters),
+ (Vibration(load_data=load_vibration), deflection_parameters),
+ (Pull(load_data=load_pull), pull_parameters),
+ ]:
+ external_parameters.update(parameters)
+ parameters.update(experiment.crunch(external_parameters))
+ parameters.update(override_parameters)
+ for exp in experiments:
+ exp.descendants.append(experiment)
+ experiments.append(experiment)
+ return (parameters, experiments)
+
+def print_parameters(parameters, experiments):
+ for experiment in experiments:
+ for k in experiment.output_parameters:
+ if k.startswith('delta_'):
+ continue
+ delta_k = 'delta_{}'.format(k)
+ if delta_k in experiment.output_parameters:
+ print('{} = {:g} +/- {:g} {}'.format(
+ k, parameters[k], parameters[delta_k],
+ experiment.output_parameter_units[k]))
+ else:
+ print('{} = {:g} {}'.format(
+ k, parameters[k],
+ experiment.output_parameter_units[k]))
+
+
+if __name__ == '__main__':
+ import argparse
+
+ parser = argparse.ArgumentParser(description='Play a pure tone')
+ parser.add_argument(
+ '-f', '--frequency', type=float, help='tone frequency in Hz')
+
+ args = parser.parse_args()
+
+ sigma_z = 1/9.2e-9 # V/m
+ gamma_z = piezo_parameters['z_gain'] * piezo_parameters['output_gamma_z']
+ alpha_z = sigma_z / gamma_z
+ print(('hacked alpha_z = {:g} to match {:g} m/V piezo calibration'
+ ).format(alpha_z, 1/sigma_z))
+ parameters,experiments = load_experiments(
+ override_parameters={ # decreasing sigma_z spreads the peaks
+ 'sigma_z': sigma_z,
+ 'alpha_z': alpha_z,
+ 'delta_alpha_z': 0,
+ })
+ print_parameters(parameters, experiments)
+ for experiment in experiments:
+ experiment.plot(slider=True)
+ _pyplot.show()
--- /dev/null
+#!/usr/bin/env python
+#
+# Freely rotating chain (FRC) simulation following Livadaru et al. [1].
+#
+# [1]: http://pubs.acs.org/doi/abs/10.1021/ma020751g
+
+import numpy as _n
+import matplotlib.pyplot as _pyplot
+from mpl_toolkits.mplot3d import Axes3D as _Axes3D
+try: # SciPy >= 0.10.0
+ from scipy.constants import Boltzmann as _kB
+except ImportError: # SciPy < 0.10.0
+ from scipy.constants import Bolzmann as _kB
+from scipy.integrate import quad as _quad
+
+
+MANUAL_T_INTEGRAL = True
+MANUAL_GIP_INTEGRAL = True
+MANUAL_GIBBS_INTEGRAL = True
+MANUAL_INTEGRAL_STEPS = 5000
+_T_cache = {}
+
+
+def _integral(func, a, b, args, manual=True, M=MANUAL_INTEGRAL_STEPS):
+ if manual:
+ dx = (b-a)/(M-1)
+ integral = 0
+ for x in _n.linspace(a, b, M):
+ xp = x + dx
+ y = func(x, *args)
+ integral += y*dx
+ else:
+ integral,abserr = _quad(func, a, b, args=args)
+ return integral
+
+
+def _Gip_kernel(theta_prime, theta, Gis, fbkBT, gamma, b):
+ "eqn. 24, kernel"
+ transfer = T(theta=theta, theta_prime=theta_prime, gamma=gamma, b=b)
+ if len(Gis) == 1: # only holding G0
+ return Gis[-1](theta=theta, fbkBT=fbkBT) * transfer
+ else: # holding Gi>0
+ return Gip(
+ theta=theta_prime, Gis=Gis[:-1], fbkBT=fbkBT,
+ gamma=gamma, b=b) * transfer
+
+def Gip(theta, Gis, fbkBT, gamma, b):
+ "eqn. 24 wrapper"
+ integral = _integral(
+ func=_Gip_kernel, a=0, b=2*_n.pi, args=(theta, Gis, fbkBT, gamma, b),
+ manual=MANUAL_GIP_INTEGRAL)
+ return _n.exp(fbkBT*_n.cos(theta)) * abs(integral)
+
+def delta(y, y_p):
+ if (y_p == 0 and y == 0): # double, but we'll get a single too.
+ return 1
+ elif y_p == 0: # single whammy
+ return 1
+ elif y/y_p < 0: # crossing
+ return 1
+ return 0
+
+def delta_h(theta, theta_prime, phi, phi_p, gamma):
+ y = theta_prime - h(theta=theta, phi=phi, gamma=gamma)
+ y_p = theta_prime - h(theta=theta, phi=phi_p, gamma=gamma)
+ return delta(y=y, y_p=y_p)
+
+def T(theta, theta_prime, gamma, b):
+ "eqn. 25"
+ key = (theta, theta_prime, gamma)
+ if key not in _T_cache:
+ assert MANUAL_T_INTEGRAL, MANUAL_T_INTEGRAL
+ #integral,abserr = _quad(
+ # _T_kernel, 0, 2*_n.pi, args=(theta, theta_prime, gamma))
+ M = MANUAL_INTEGRAL_STEPS
+ d_phi = 2*_n.pi/(M-1)
+ integral = 0
+ for phi in _n.linspace(0, 2*_n.pi, M):
+ phi_p = phi + d_phi
+ integral += delta_h(
+ theta=theta, theta_prime=theta_prime, phi=phi, phi_p=phi_p,
+ gamma=gamma)
+ _T_cache[key] = integral
+ return _T_cache[key]
+
+def _h(theta, phi, gamma):
+ "for eqn.s 26 and 31"
+ return _n.cos(theta)*_n.cos(gamma) + _n.sin(theta)*_n.sin(gamma)*_n.cos(phi)
+
+def h(theta, phi, gamma):
+ "eqn. 26"
+ return _n.arccos(_h(theta=theta, phi=phi, gamma=gamma))
+
+def G0(theta, fbkBT):
+ "eqn. 27"
+ return _n.exp(fbkBT*_n.cos(theta)) * abs(_n.sin(theta))
+
+def GN(theta, N, fbkBT, gamma, b=None):
+ "for fig. 4"
+ Gis = [G0] + [Gip]*(N-1)
+ if N == 1:
+ return Gis[-1](theta=theta, fbkBT=fbkBT)
+ return Gis[-1](theta=theta, Gis=Gis[:-1], fbkBT=fbkBT, gamma=gamma, b=b)
+
+def monte_carlo_GN(N=1, fbkBT=1, gamma_degrees=20, b=1, plot=True):
+ "a check for GN()"
+ gamma = gamma_degrees * _n.pi / 180
+ dirs = _n.zeros((N, 3), dtype=float)
+ pos = _n.zeros((N+1, 3), dtype=float)
+ for i in range(N):
+ dice = 1
+ p = 0
+ hits = 0
+ misses = 0
+ scale = 0.001
+ while dice > p:
+ misses += 1
+ if i == 0:
+ new_direction = _n.random.randn(3)
+ new_direction /= _n.linalg.norm(new_direction)
+ else:
+ old_direction = dirs[i-1,:]
+ phi = 2 * _n.pi * _n.random.rand()
+ if old_direction[0] == 0:
+ basis = _n.array([1, 0, 0], dtype=float)
+ else:
+ basis = _n.array([0, 1, 0], dtype=float)
+ perpendicular_1 = _n.cross(old_direction, basis)
+ perpendicular_1 /= _n.linalg.norm(perpendicular_1)
+ perpendicular_2 = _n.cross(old_direction, perpendicular_1)
+ new_direction = (
+ old_direction * _n.cos(gamma) +
+ perpendicular_1 * _n.sin(gamma) * _n.cos(phi) +
+ perpendicular_2 * _n.sin(gamma) * _n.sin(phi))
+ if float(hits)/misses > 0.01:
+ scale /= 2
+ elif float(hits)/misses < 0.001:
+ scale *= 2
+ p = _n.exp(fbkBT*new_direction[0]) * scale
+ dice = _n.random.rand()
+ misses -= 1
+ hits += 1
+ dirs[i,:] = new_direction
+ pos[i+1,:] = pos[i,:] + new_direction
+ print('scale: {}'.format(scale))
+ if plot:
+ figure = _pyplot.figure()
+ axes = figure.add_subplot(1, 1, 1, projection='3d', aspect='equal')
+ axes.plot([x for x,y,z in pos], [y for x,y,z in pos], [z for x,y,z in pos])
+ axes.set_title('monte carlo FRC (gamma={} deg)'.format(gamma_degrees))
+ axes.set_xlabel('x')
+ axes.set_ylabel('y')
+ axes.set_ylabel('z')
+ figure.savefig('mv-frc-g-{}.png'.format(gamma_degrees))
+ return _n.arccos(dirs[-1][0])
+
+def Gibbs(N, fbkBT, gamma, b):
+ "eqn. 28"
+ Gis = [G0] + [Gip]*(N-1)
+ integral = _integral(
+ Gis[-1], 0, 2*_n.pi, args=(Gis[-1], fbkBT, gamma, b),
+ manual=MANUAL_GIBBS_INTEGRAL)
+
+def Rz(N, f, b, kBT, gamma, df):
+ "eqn. 29"
+ G = Gibbs(N=N, fbkBT=f*b/kBT, gamma=gamma)
+ Gp = Gibbs(N=N, fbkBT=(f+df)*b/kBT, gamma=gamma)
+ return (Gp-G)/df
+
+def test_h(gamma_degrees=20):
+ figure = _pyplot.figure()
+ axes = figure.add_subplot(1, 1, 1)
+ axes.hold(True)
+ gamma = gamma_degrees * _n.pi / 180
+ for i in range(100):
+ thetas = _n.linspace(0*_n.pi, 2*_n.pi, 50)
+ phis = _n.linspace(0, 2*_n.pi, 50)
+ data = _n.zeros((len(thetas), len(phis)), dtype=float)
+ for c,phi in enumerate(phis):
+ for r,theta in enumerate(thetas):
+ data[r,c] = h(theta, phi, gamma)
+ X,Y = _n.meshgrid(thetas, phis)
+ pcolor = axes.pcolor(X, Y, data.T)
+ figure.colorbar(pcolor)
+ axes.axis([0, 2*_n.pi, 0, 2*_n.pi])
+ axes.set_title('h(theta, phi, gamma={} deg)'.format(gamma_degrees))
+ axes.set_xlabel('theta (degrees)')
+ axes.set_ylabel('phi (degrees)')
+ figure.savefig('test-h-g-{}.png'.format(gamma_degrees))
+
+def test_T(gamma_degrees=20):
+ figure = _pyplot.figure()
+ axes = figure.add_subplot(1, 1, 1)
+ axes.hold(True)
+ gamma = gamma_degrees * _n.pi / 180
+ for i in range(100):
+ thetas = _n.linspace(0*_n.pi, 2*_n.pi, 50)
+ theta_primes = _n.linspace(0, 2*_n.pi, 50)
+ data = _n.zeros((len(thetas), len(theta_primes)), dtype=float)
+ for c,theta_prime in enumerate(theta_primes):
+ for r,theta in enumerate(thetas):
+ data[r,c] = T(theta, theta_prime, gamma)
+ X,Y = _n.meshgrid(thetas, theta_primes)
+ pcolor = axes.pcolor(X, Y, data.T)
+ figure.colorbar(pcolor)
+ axes.axis([0, 2*_n.pi, 0, 2*_n.pi])
+ axes.set_title("T(theta, theta', gamma={} deg)".format(gamma_degrees))
+ axes.set_xlabel('theta (degrees)')
+ axes.set_ylabel("theta' (degrees)")
+ figure.savefig('test-T-g-{}.png'.format(gamma_degrees))
+
+def figure_4_frame(N=3, gamma_degrees=20):
+ figure = _pyplot.figure()
+ axes = figure.add_subplot(1, 1, 1)
+ axes.hold(True)
+ gamma = gamma_degrees * _n.pi / 180
+ for color,N in [('r', 1), ('b', 2)]: #, ('g', 3)]:
+ fbkBTs = [100] #_n.linspace(0, 100, 5)
+ thetas = _n.linspace(0, 2*_n.pi, 50)
+ data = _n.zeros((len(thetas), len(fbkBTs)), dtype=float)
+ for c,fbkBT in enumerate(fbkBTs):
+ for r,theta in enumerate(thetas):
+ print((r,c))
+ #data[r,c] = G0(theta=theta, fbkBT=fbkBT)
+ data[r,c] = GN(theta=theta, N=N, fbkBT=fbkBT, gamma=gamma)
+ data[:,c] /= data[:,c].sum()
+ degrees = thetas * 180 / _n.pi
+ axes.plot(degrees, data, '{}.-'.format(color))
+ axes.axes.axvline(gamma_degrees/2., color='k')
+ axes.set_xbound([0, 100])
+ axes.set_title('fig. 4 (gamma = {} deg)'.format(gamma_degrees))
+ axes.set_xlabel('theta (degrees)')
+ axes.set_ylabel('probability G_N(theta)')
+ figure.savefig('test-G_N-g-{}.png'.format(gamma_degrees))
+
+###
+#
+#def hz(theta, phi, gamma, b):
+# "eqn. 31"
+# return b * _h(theta=theta, phi=phi, gamma=gamma)
+#
+#def delta_hz(theta, z, z_prime, phi, phi_p, gamma, b):
+# "for eqn. 30"
+# y = z_prime - z - hz(theta=theta, phi=phi, gamma=gamma)
+# y_p = z_prime - hz(theta=theta, phi=phi_p, gamma=gamma)
+# return delta(y=y, y_p=y_p)
+#
+#def Tz(z, theta, z_prime, theta_prime, gamma, b):
+# "eqn. 30"
+# key = (z, theta, z_prime, theta_prime, gamma, b)
+# if key not in _T_cache:
+# assert MANUAL_T_INTEGRAL, MANUAL_T_INTEGRAL
+# #integral,abserr = _quad(
+# # _T_kernel, 0, 2*_n.pi, args=(theta, theta_prime, gamma, b))
+# M = MANUAL_INTEGRAL_STEPS
+# d_phi = 2*_n.pi/(M-1)
+# integral = 0
+# for phi in _n.linspace(0, 2*_n.pi, M):
+# phi_p = phi + d_phi
+# integral += delta_h(
+# theta=theta, theta_prime=theta_prime, phi=phi, phi_p=phi_p,
+# gamma=gamma) * delta_hz(
+# theta=theta, z=z, z_prime=z_prime, phi=phi, phi_p=phi_p,
+# gamma=gamma, b)
+# _T_cache[key] = integral
+# return _T_cache[key]
+#
+#def G0Z(z, theta, gamma, b):
+# return delta(z-b*_n.cos(theta))*_n.sin(theta)
+
+
+
+if __name__ == '__main__':
+ for gamma_degrees in [20, 60]: #, 90]:
+ #test_h(gamma_degrees=gamma_degrees)
+ #test_T(gamma_degrees=gamma_degrees)
+ figure_4_frame(N=100, gamma_degrees=gamma_degrees)
+ #thetas = []
+ #for x in range(500):
+ # print('calc theta {} for {}'.format(x, gamma_degrees))
+ # thetas.append(monte_carlo_GN(N=100, fbkBT=100, gamma_degrees=gamma_degrees, plot=False))
+ #figure = _pyplot.figure()
+ #axes = figure.add_subplot(1, 1, 1)
+ #axes.hist(_n.array(thetas) * 180 / _n.pi)
+ #print('done hist')
+ _pyplot.show()
projects / blog.git / commitdiff
