From: W. Trevor King Date: Fri, 14 Sep 2012 22:02:10 +0000 (-0400) Subject: posts:comparing_velocity_clamp_experiments: added post. X-Git-Url: http://git.tremily.us/?a=commitdiff_plain;h=0b1459964fe41b26e71457388845617847005cdb;p=blog.git posts:comparing_velocity_clamp_experiments: added post. --- diff --git a/posts/Comparing_velocity_clamp_experiments.mdwn_itex b/posts/Comparing_velocity_clamp_experiments.mdwn_itex new file mode 100644 index 0000000..4d144cc --- /dev/null +++ b/posts/Comparing_velocity_clamp_experiments.mdwn_itex @@ -0,0 +1,342 @@ +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. + + + +Calibrating the piezo +===================== + + + +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: + + + +\[ + \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 + + + + +\[ +\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. + + + +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. + + + +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. + + + +\[ + \alpha_d = \frac{\text{bump diode bits}}{\text{bump piezo bits}} + = 2.24 \text{ diode bits/piezo bits} \;. +\] + +Related conversion factors are + + + + +\[ +\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} +\] + + + +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} +\] + + + +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]] diff --git a/posts/Comparing_velocity_clamp_experiments/bump.png b/posts/Comparing_velocity_clamp_experiments/bump.png new file mode 100644 index 0000000..5d361eb Binary files /dev/null and b/posts/Comparing_velocity_clamp_experiments/bump.png differ diff --git a/posts/Comparing_velocity_clamp_experiments/crunch.py b/posts/Comparing_velocity_clamp_experiments/crunch.py new file mode 100755 index 0000000..f70829a --- /dev/null +++ b/posts/Comparing_velocity_clamp_experiments/crunch.py @@ -0,0 +1,1013 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# +# Copyright (C) 2012 W. Trevor King +# +# 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 . + +"""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() diff --git a/posts/Comparing_velocity_clamp_experiments/figures.py b/posts/Comparing_velocity_clamp_experiments/figures.py new file mode 100755 index 0000000..29f08ee --- /dev/null +++ b/posts/Comparing_velocity_clamp_experiments/figures.py @@ -0,0 +1,105 @@ +#!/usr/bin/env python +# +# Copyright (C) 2012 W. Trevor King +# +# 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 . + +"""Generate figures for the blog post. +""" + +import numpy as _numpy + +from crunch import * # sloppy import for a one-off hack script. + + +good_piezo_file = '~/tmp/rsrch/calibrate-piezo/data/calibrate-piezo-2012-02-06T16-49-29.dat' + +good_piezo_parameters = { + 'pits': [ # (z_bottom_bits, z_top_bits) pairs + (32036, 32604), + (32070, 32630), + (32175, 32721), + (32264, 32789), + (32292, 32830), + (32413, 32987), + (32565, 33077), + ], + '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 + } + +interference_pull_file = '~/tmp/rsrch/calibrate-piezo/interference/2012-05-23T15-56-50.h5' +interference_pull_bounds = [20000, 35000, 33000, 35000] +interference_alpha_z = 1.56e10 + + +if __name__ == '__main__': + parameters,experiments = load_experiments() + diff = _numpy.array([b-a for a,b in piezo_parameters['pits']], dtype=float) + print('pit_depth = {:g} +/- {:g} bits'.format(diff.mean(), diff.std())) + print('pit_depth = {:g} m'.format(piezo_parameters['pit_depth'])) + print_parameters(parameters, experiments) + for experiment in experiments: + if experiment.name == 'pull': + for filename,kwargs in [ + ('pull-piezo.png', { + 'xlabel': 'z piezo (m)', + 'bounds': [-10e-9, 500e-9, -250e-12, 250e-12], + }), + ('pull-protein.png', { + 'xlabel': 'protein extension (m)', + 'bounds': [-10e-9, 500e-9, -250e-12, 250e-12], + }), + ('pull-contour-space.png', { + 'xlabel': 'contour-space (m)', + 'bounds': [140e-9, 350e-9, 0, 130e-12], + }), + ]: + figure = experiment.plot(**kwargs) + figure.savefig(filename) + else: + experiment.plot(save=True) + load_good_piezo = lambda: load_piezo(filename=good_piezo_file) + piezo = Piezo(load_data=load_good_piezo) + data = piezo.crunch(parameters=good_piezo_parameters) + piezo.plot(filename='piezo-good.png') + print('good piezo:') + for k,v in data.items(): + print(' {} = {:g}'.format(k, v)) + + data = load_pull(filename=interference_pull_file, with_temperature=False) + figure = Pull._plot( + z_piezo=_numpy.array(data['alpha_v'], dtype=float), + z=_numpy.zeros(data['alpha_v'].shape), + z_contour=_numpy.zeros(data['alpha_v'].shape), + deflection=_numpy.array(data['alpha_f'], dtype=float), + pull_temperature=297.15, # only used for contour-space plots + **pull_parameters) + axes = figure.axes[0] + axes.axis(interference_pull_bounds) + axes.set_title('piezo calibration via interference') + axes.set_xlabel('z piezo (bits)') + axes.set_ylabel('deflection (bits)') + figure.savefig('pull-interference.png') + + pull = experiments[-1] + parameters['alpha_z'] = interference_alpha_z + pull_parameters.update(parameters) + pull.crunch(pull_parameters) + pull.plot( + filename='pull-contour-space-interference.png', + xlabel='contour-space (m)', + bounds=[140e-9, 350e-9, 0, 130e-12], + ) diff --git a/posts/Comparing_velocity_clamp_experiments/frc-details.py b/posts/Comparing_velocity_clamp_experiments/frc-details.py new file mode 100755 index 0000000..7131392 --- /dev/null +++ b/posts/Comparing_velocity_clamp_experiments/frc-details.py @@ -0,0 +1,286 @@ +#!/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() diff --git a/posts/Comparing_velocity_clamp_experiments/piezo-good.png b/posts/Comparing_velocity_clamp_experiments/piezo-good.png new file mode 100644 index 0000000..4c4136e Binary files /dev/null and b/posts/Comparing_velocity_clamp_experiments/piezo-good.png differ diff --git a/posts/Comparing_velocity_clamp_experiments/piezo.png b/posts/Comparing_velocity_clamp_experiments/piezo.png new file mode 100644 index 0000000..0610998 Binary files /dev/null and b/posts/Comparing_velocity_clamp_experiments/piezo.png differ diff --git a/posts/Comparing_velocity_clamp_experiments/pull-contour-space-interference.png b/posts/Comparing_velocity_clamp_experiments/pull-contour-space-interference.png new file mode 100644 index 0000000..dc00391 Binary files /dev/null and b/posts/Comparing_velocity_clamp_experiments/pull-contour-space-interference.png differ diff --git a/posts/Comparing_velocity_clamp_experiments/pull-contour-space.png b/posts/Comparing_velocity_clamp_experiments/pull-contour-space.png new file mode 100644 index 0000000..8df385e Binary files /dev/null and b/posts/Comparing_velocity_clamp_experiments/pull-contour-space.png differ diff --git a/posts/Comparing_velocity_clamp_experiments/pull-interference.png b/posts/Comparing_velocity_clamp_experiments/pull-interference.png new file mode 100644 index 0000000..db99331 Binary files /dev/null and b/posts/Comparing_velocity_clamp_experiments/pull-interference.png differ diff --git a/posts/Comparing_velocity_clamp_experiments/pull-piezo.png b/posts/Comparing_velocity_clamp_experiments/pull-piezo.png new file mode 100644 index 0000000..aa0af78 Binary files /dev/null and b/posts/Comparing_velocity_clamp_experiments/pull-piezo.png differ diff --git a/posts/Comparing_velocity_clamp_experiments/pull-protein.png b/posts/Comparing_velocity_clamp_experiments/pull-protein.png new file mode 100644 index 0000000..af0cf52 Binary files /dev/null and b/posts/Comparing_velocity_clamp_experiments/pull-protein.png differ diff --git a/posts/Comparing_velocity_clamp_experiments/vibration.png b/posts/Comparing_velocity_clamp_experiments/vibration.png new file mode 100644 index 0000000..d7cc008 Binary files /dev/null and b/posts/Comparing_velocity_clamp_experiments/vibration.png differ