posts/Comparing_velocity_clamp_experiments/crunch.py

   1 #!/usr/bin/env python
   2 # -*- coding: utf-8 -*-
   3 #
   4 # Copyright (C) 2012  W. Trevor King <wking@tremily.us>
   5 #
   6 # This program is free software: you can redistribute it and/or modify
   7 # it under the terms of the GNU General Public License as published by
   8 # the Free Software Foundation, either version 3 of the License, or
   9 # (at your option) any later version.
  10 #
  11 # This program is distributed in the hope that it will be useful,
  12 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14 # GNU General Public License for more details.
  15 #
  16 # You should have received a copy of the GNU General Public License
  17 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
  18
  19 """Calculate all parameters associated with a velocity-clamp pull.
  20
  21 This helps you track down where in the calibration process things may
  22 be going wrong.
  23 """
  24
  25 import collections as _collections
  26 import functools as _functools
  27 import os.path as _os_path
  28
  29 try:  # SciPy >= 0.10.0
  30     from scipy.constants import Boltzmann as _kB
  31 except ImportError:  # SciPy < 0.10.0
  32     from scipy.constants import Bolzmann as _kB
  33
  34 import h5py as _h5py
  35
  36 # for some reason my scatterplot 'pixel' markers were taking up several pixels.
  37 import matplotlib.markers as _markers
  38 def _init(self, marker=None, fillstyle='full'):
  39     self._fillstyle = fillstyle
  40     self.set_marker(marker)
  41     self.set_fillstyle(fillstyle)
  42     if marker == ',':
  43         self._transform = _markers.Affine2D().scale(0.1)
  44 _markers.MarkerStyle.__init__ = _init
  45 # end pixel hack
  46
  47 import matplotlib.pyplot as _pyplot
  48 import matplotlib.widgets as _widgets
  49 import numpy as _numpy
  50 import scipy as _scipy
  51
  52 import calibcant.calibrate as _calibcant_calibrate
  53 import calibcant.config as _calibcant_config
  54 import calibcant.bump_analyze as _calibcant_bump_analyze
  55 import calibcant.vibration_analyze as _calibcant_vibration_analyze
  56 import pypiezo.surface as _pypiezo_surface
  57 import FFT_tools as _FFT_tools
  58
  59 ## your data
  60
  61 piezo_file = '~/tmp/rsrch/calibrate-piezo/data/calibrate-piezo-2012-02-06T16-23-30.dat'
  62
  63 piezo_parameters = {
  64     'pits': [  # (z_bottom_bits, z_top_bits) pairs
  65         (30044, 32736),
  66         (30167, 33033),
  67         ],
  68     'pit_depth': 200e-9,  # meters
  69     'z_gain': 20,  # volts applied at piezo / volts output from DAC
  70     'output_gamma_z': 10./2**15  # volts output from DAC / bit
  71     }
  72
  73 # cantilever calibration
  74
  75 calibcant_file = '~/rsrch/data/calibrate_cantilever/2012-03-26T13-02-19-calibcant-data.h5'
  76 calibcant_bump = 3
  77 calibcant_vibration = 15
  78 deflection_parameters = {
  79     'd_gain': 1,  # volts ADC input / volts photodiode differencer output
  80     'input_gamma_d': 10./2**15,  # volts input at ADC / bit
  81     }
  82
  83 # unfolding pull
  84
  85 pull_file = '~/rsrch/data/unfold-h5/good-2/2012-03-26T13-57-03.h5'
  86 _gamma = 22*_numpy.pi/180
  87 _l_over_b = _numpy.cos(_gamma/2.)/abs(_numpy.log(_numpy.cos(_gamma)))
  88 _b = 0.4e-9/_l_over_b
  89 pull_parameters = {
  90     'drift_slope': 0.06,
  91     'alpha_f_offset': -50,
  92     'smooth_steps': 600,
  93     'start_sawteeth': 50e-9,  # meters from kink to start of sawteeth
  94     'stop_sawteeth': 300e-9,  # meters from kink to end of sawteeth
  95     'dL': 29e-9,  # meters, expected value for I27
  96     'contour_space_model': 'frc',
  97     #'contour_space_model': 'wlc',
  98     #'contour_space_parameters': None,  # use defaults
  99     'contour_space_parameters': {'b': _b, 'gamma': _gamma},
 100     'contour_space_lines': _numpy.arange(75e-9, 350e-9, 29e-9),
 101     }
 102
 103 ## readers for your data
 104
 105 def load_piezo(filename=piezo_file):
 106     data = _numpy.loadtxt(_os_path.expanduser(filename))
 107     x = data[:,0]
 108     z = data[:,1]
 109     deflection = data[:,2]
 110     return {'x': x, 'z': z, 'deflection': deflection}
 111
 112 def load_bump(filename=calibcant_file, index=calibcant_bump):
 113     calibrator,calibcant_data,calibcant_raw_data = (
 114         _calibcant_calibrate.load_all(filename=_os_path.expanduser(filename)))
 115     bump_data = calibcant_raw_data['bump'][index]
 116     z = bump_data['raw']['z']
 117     deflection = bump_data['raw']['deflection']
 118     # HACK, 0 index not necessarily deflection
 119     deflection_high_voltage_rail = (
 120         calibrator.config['afm']['piezo']['inputs'][0]['maxdata'])
 121     return {
 122         'z': z,
 123         'deflection': deflection,
 124         'deflection_high_voltage_rail': deflection_high_voltage_rail,
 125         }
 126
 127 def load_vibration(filename=calibcant_file, index=calibcant_vibration):
 128     calibrator,calibcant_data,calibcant_raw_data = (
 129         _calibcant_calibrate.load_all(filename=_os_path.expanduser(filename)))
 130     vibration_data = calibcant_raw_data['vibration'][index]
 131     temperature = calibcant_data['raw']['temperature'].mean()
 132     deflection = vibration_data['raw']
 133     config = vibration_data['config']['vibration']
 134     psd_kwargs = {
 135         'chunk_size': config['chunk-size'],
 136         'overlap': config['overlap'],
 137         'window': config['window'],
 138         }
 139     fit_kwargs = {
 140         'min_frequency': config['minimum-fit-frequency'],
 141         'max_frequency': config['maximum-fit-frequency'],
 142         'offset': True,
 143         }
 144     #offset=config['model'] == _calibcant_config.OffsetBreitWigner)
 145     return {
 146         'deflection': deflection,
 147         'temperature': temperature,
 148         'frequency': config['frequency'],
 149         'psd_kwargs': psd_kwargs,
 150         'fit_kwargs': fit_kwargs,
 151         }
 152
 153 def load_pull(filename=pull_file, with_temperature=True):
 154     with _h5py.File(_os_path.expanduser(filename)) as f:
 155         alpha_v = _numpy.array(f['unfold']['z'][...].flat, dtype=float)
 156         alpha_f = _numpy.array(
 157             f['unfold']['deflection'][...].flat, dtype=float)
 158         v = float(f['config']['unfold']['velocity'][...])
 159         if with_temperature:
 160             temperature = float(f['environment']['temperature'][...])
 161         else:
 162             temperature = None
 163     return {
 164         'alpha_v': alpha_v,
 165         'alpha_f': alpha_f,
 166         'velocity': v,
 167         'temperature': temperature,
 168         }
 169
 170
 171 ## crunching
 172
 173 class NamedObject (object):
 174     def __init__(self, name=None):
 175         self.name = name
 176
 177     def log(self, message):
 178         print('{}: {}'.format(self.name, message))
 179
 180
 181 class SliderFigure (NamedObject):
 182     def __init__(self, figure, axes=None, parameters=None, data=None,
 183                  callback=None, **kwargs):
 184         super(SliderFigure, self).__init__(**kwargs)
 185         self.figure = figure
 186         if axes is None:
 187             axes = list(figure.axes)
 188         self.axes = axes
 189         self.parameters = parameters
 190         self.data = data
 191         self.callback = callback
 192
 193         self._slider_height = 0.03
 194         self._slider_margin = 0.02
 195         self._total_slider_height = len(parameters)*(
 196             self._slider_height+self._slider_margin) + self._slider_margin
 197
 198         self._shift_original_axes_up()
 199         self._add_sliders()
 200
 201     @staticmethod
 202     def _map_original_axes_y(y, new_bottom):
 203         """map [0,1] -> [new_bottom,1]
 204
 205         >>> SliderFigure._map_original_axes_y(0, 0.5)
 206         0.5
 207         >>> SliderFigure._map_original_axes_y(0.5, 0.4)
 208         0.7
 209         """
 210         return 1 - (1-y)*(1-new_bottom)
 211
 212     def _shift_original_axes_up(self):
 213         for axes in self.figure.axes:
 214             p = axes.get_position()
 215             p.y0 = self._map_original_axes_y(p.y0, self._total_slider_height)
 216             p.y1 = self._map_original_axes_y(p.y1, self._total_slider_height)
 217             axes.set_position(p)
 218
 219     def _add_sliders(self):
 220         self._sliders = {}
 221         n = len(self.parameters)
 222         h = self._slider_height
 223         m = self._slider_margin
 224         left = 0.4
 225         w = 0.4
 226         for i,(k,v) in enumerate(self.parameters.items()):
 227             try:
 228                 x = float(v)
 229             except (TypeError, ValueError):
 230                 continue  # non-float parameter
 231             axes = self.figure.add_axes([left, m+i*(h+m), w, h])
 232             min,max = sorted([v/4., 4*v])  # v might be negative
 233             slider = _widgets.Slider(
 234                 axes, k, min, max, valinit=v, valfmt='%1.5g',
 235                 closedmin=False, closedmax=False)
 236             slider.on_changed(_functools.partial(
 237                     self.on_changed, parameter=k, slider=slider))
 238             self._sliders[k] = slider
 239
 240     def on_changed(self, value, parameter, slider=None):
 241         if slider is None:  # external change (not direct slider movement)
 242             slider = self._sliders[parameter]
 243             slider.set_val(value)
 244             return
 245         self.log('changed {} to {}, recalculating...'.format(
 246                 parameter, value))
 247         self.parameters[parameter] = value
 248         if self.callback:
 249             self.callback(
 250                 figure=self.figure, axes=self.axes, data=self.data,
 251                 parameters=self.parameters)
 252         self.log('finished recalculating change of {} to {}'.format(
 253                 parameter, value))
 254
 255
 256 class Experiment (NamedObject):
 257     def __init__(self, load_data=None, input_parameters=None,
 258                  output_parameters=None, output_parameter_units=None,
 259                  descendants=None, **kwargs):
 260         super(Experiment, self).__init__(**kwargs)
 261         self.load_data = load_data
 262         if input_parameters is None:
 263             input_parameters = tuple()
 264         self.input_parameters = input_parameters
 265         if output_parameters is None:
 266             output_parameters = tuple()
 267         self.output_parameters = output_parameters
 268         if output_parameter_units is None:
 269             output_parameter_units = {}
 270         self.output_parameter_units = output_parameter_units
 271         if descendants is None:
 272             descendants = []
 273         self.descendants = descendants
 274         self._iparameters = self._oparameters = None
 275
 276     def crunch(self, parameters=None):
 277         if self.load_data and (not hasattr(self, 'data') or self.data is None):
 278             self.data = self.load_data()
 279         self._iparameters = dict(
 280             (k,parameters[k]) for k in self.input_parameters)
 281         self._iparameters.update(self.data)
 282         oparameters = self._crunch(**self._iparameters)
 283         return self._crunch_callback(oparameters)
 284
 285     @staticmethod
 286     def _crunch(**kwargs):
 287         return {}
 288
 289     def _crunch_callback(self, oparameters):
 290         if self._oparameters is not None:
 291             for k in self.output_parameters:
 292                 if self._oparameters[k] != oparameters[k]:
 293                     self.log('{} changed to {} during recrunch'.format(
 294                             k, oparameters[k]))
 295                     for d in self.descendants:
 296                         if k in d.input_parameters:
 297                             self.log('{} depends on {}'.format(d.name, k))
 298                             if d.slider_figure:
 299                                 #slider.set_val(value)
 300                                 d.slider_figure.on_changed(
 301                                     value=oparameters[k], parameter=k)
 302         self._oparameters = oparameters
 303         return dict((k,self._oparameters[k]) for k in self.output_parameters)
 304
 305     def plot(self, save=False, slider=False, filename=None, **kwargs):
 306         """Wrapper around plotting functions for consistent handling."""
 307         if save and not filename:
 308             filename = '{}.png'.format(self.name)
 309         assert self._oparameters is not None, 'crunch before plotting!'
 310         kwargs.update(self._iparameters)
 311         kwargs.update(self._oparameters)
 312         kwargs['log'] = self.log
 313         kwargs['slider'] = slider
 314         figure = self._plot(**kwargs)
 315         if filename:
 316             figure.savefig(filename)
 317         if slider:
 318             parameters = dict(
 319                 (k,self._iparameters[k]) for k in self.input_parameters)
 320             parameters['_crunch'] = self._crunch
 321             parameters['_crunch_callback'] = self._crunch_callback
 322             self.slider_figure = SliderFigure(
 323                 name='{}-slider'.format(self.name), figure=figure,
 324                 data=self.data, parameters=parameters, callback=self._replot)
 325         return figure
 326
 327     @staticmethod
 328     def _plot(**kwargs):
 329         raise NotImplementedError()
 330
 331     @staticmethod
 332     def _replot(figure, axes, data=None, parameters=None):
 333         parameters = dict(parameters)
 334         crunch = parameters.pop('_crunch')
 335         crunch_callback = parameters.pop('_crunch_callback')
 336         parameters.update(data)
 337         oparameters = crunch(**parameters)
 338         crunch_callback(oparameters)
 339         return oparameters
 340
 341     # utilities
 342
 343     @staticmethod
 344     def smooth(x, window_len=11):
 345         """
 346         >>> Experiment.smooth([1,2,3,4,5], window_len=2)
 347         array([ 1. ,  1.5,  2.5,  3.5,  4.5])
 348
 349         With the padding, the convolution was over `[1, 1, 2, 3, 4, 5]`.
 350
 351         >>> Experiment.smooth([1,2,3,4,5], window_len=3)
 352         array([ 1.33333333,  2.        ,  3.        ,  4.        ,  4.66666667])
 353
 354         With the padding, the convolution was over `[1, 1, 2, 3, 4, 5, 5]`.
 355         """
 356         assert window_len > 1, window_len
 357         window = _numpy.ones(window_len, 'd')  # moving average
 358         front_pad = int(window_len/2)  # truncate
 359         if int(window_len) % 2 == 0:
 360             back_pad = front_pad - 1
 361         else:
 362             back_pad = front_pad
 363         front_padding = [x[0]] * front_pad
 364         back_padding = [x[-1]] * back_pad
 365         x = _numpy.concatenate((front_padding, x, back_padding))
 366         return _numpy.convolve(window/window.sum(), x, mode='valid')
 367
 368     @staticmethod
 369     def true_set(dictionary):
 370         """
 371         >>> from collections import defaultdict
 372         >>> Experiment.true_set({})
 373         False
 374         >>> Experiment.true_set({'a': True})
 375         True
 376         >>> Experiment.true_set(defaultdict(lambda: True))
 377         True
 378         >>> Experiment.true_set(defaultdict(lambda: 'abc'))
 379         False
 380         """
 381         if True in dictionary.values():
 382             return True
 383         if isinstance(dictionary, _collections.defaultdict):
 384             if dictionary.default_factory() is True:
 385                 return True
 386         return False
 387
 388
 389 class Piezo (Experiment):
 390     def __init__(self, name='piezo', *args, **kwargs):
 391         kwargs['input_parameters'] = [
 392             'pits', 'pit_depth', 'z_gain', 'output_gamma_z']
 393         kwargs['output_parameters'] = [
 394             'alpha_z', 'delta_alpha_z', 'gamma_z', 'sigma_z']
 395         kwargs['output_parameter_units'] = {
 396             'alpha_z': 'bit/m',
 397             'delta_alpha_z': 'bit/m',
 398             'gamma_z': 'V/bit',
 399             'sigma_z': 'V/m',
 400             }
 401         super(Piezo, self).__init__(name=name, *args, **kwargs)
 402
 403     @staticmethod
 404     def _crunch(pits, pit_depth, z_gain, output_gamma_z, **kwargs):
 405         ret = {}
 406         ret['dz'] = _numpy.array([(top-bot) for bot,top in pits], dtype=float)
 407         ret['alpha_z'] = ret['dz'].mean()/pit_depth
 408         ret['delta_alpha_z'] = ret['dz'].std()/pit_depth
 409         ret['gamma_z'] = z_gain * output_gamma_z
 410         ret['sigma_z'] = ret['gamma_z'] * ret['alpha_z']
 411         return ret
 412
 413     @staticmethod
 414     def _plot(x, z, deflection, **kwargs):
 415         figure = _pyplot.figure()
 416         axes = figure.add_subplot(1, 1, 1)
 417         color_dimension = 'deflection'
 418         color_dimension = 'trace'
 419         #color_dimension = 'debug trace'
 420         if color_dimension in ['trace', 'debug trace']:
 421             trace = _numpy.zeros(x.shape, dtype=bool)
 422             direction = None
 423             for i,xi in enumerate(x):
 424                 last_x = x[max(0, i-1)]
 425                 new_direction = xi > last_x
 426                 if direction is None and xi != last_x:  # first movement
 427                     trace[0:i] = direction = new_direction
 428                 else:  # subsequent movement
 429                     if xi != last_x:
 430                         direction = new_direction
 431                     trace[i] = direction
 432         if color_dimension == 'debug trace':
 433             axes.plot(x, 'b.', 2**15*trace, 'r.-')
 434         elif color_dimension == 'trace':
 435             scatter = axes.scatter(
 436                 x, z, c=trace, marker=',', edgecolors='face')
 437         else:  # color_dimension == 'deflection'
 438             scatter = axes.scatter(
 439                 x, z, c=deflection, marker=',',
 440                 edgecolors='face')
 441         if color_dimension in ['trace', 'deflection']:
 442             colorbar = figure.colorbar(scatter)
 443         if color_dimension == 'trace':
 444             colorbar.set_label('trace/retrace')
 445         elif color_dimension == 'deflection':
 446             colorbar.set_label('deflection (bits)')
 447         axes.set_title('piezo calibration')
 448         axes.set_xlabel('x piezo (bits)')
 449         axes.set_ylabel('z piezo (bits)')
 450         return figure
 451
 452
 453 class Bump (Experiment):
 454     def __init__(self, name='bump', *args, **kwargs):
 455         kwargs['input_parameters'] = [
 456             'input_gamma_d', 'd_gain', 'alpha_z']
 457         kwargs['output_parameters'] = [
 458             'alpha_d', 'gamma_d', 'sigma_d']
 459         kwargs['output_parameter_units'] = {
 460             'alpha_d': 'diode_bit/piezo_bit',
 461             'gamma_d': 'V/bit',
 462             'sigma_d': 'V/m',
 463             }
 464         super(Bump, self).__init__(name=name, *args, **kwargs)
 465
 466     @staticmethod
 467     def _crunch(z, deflection, deflection_high_voltage_rail,
 468                 input_gamma_d, d_gain, alpha_z):
 469         parameters = _calibcant_bump_analyze.fit(
 470             z, deflection, high_voltage_rail=deflection_high_voltage_rail,
 471             param_guesser=_calibcant_bump_analyze.limited_linear_param_guess,
 472             model=_calibcant_bump_analyze.limited_linear,
 473             sensitivity_from_fit_params=lambda parameters: parameters)
 474         ret = {}
 475         ret['deflection_fit'] = _calibcant_bump_analyze.limited_linear(
 476             z, parameters, high_voltage_rail=deflection_high_voltage_rail)
 477         ret['alpha_d'] = parameters[2]  # bump diode bits / bump piezo bit
 478         ret['gamma_d'] = input_gamma_d / d_gain
 479         ret['sigma_d'] = ret['gamma_d'] * ret['alpha_d'] * alpha_z
 480         return ret
 481
 482     @staticmethod
 483     def _plot(alpha_d, gamma_d, sigma_d, z, deflection, deflection_fit,
 484               **kwargs):
 485         figure = _pyplot.figure()
 486         axes = figure.add_subplot(1, 1, 1)
 487         axes.plot(z, deflection, 'b.',
 488                   z, deflection_fit, 'r-')
 489         axes.set_title('photodiode calibration')
 490         axes.set_xlabel('z piezo (bits)')
 491         axes.set_ylabel('deflection (bits)')
 492         return figure
 493
 494     @staticmethod
 495     def _replot(figure, axes, data=None, parameters=None):
 496         oparameters = Experiment._replot(
 497             figure=figure, axes=axes, data=data, parameters=parameters)
 498         deflection = data['deflection']
 499         z = data['z']
 500         deflection_fit = oparameters['deflection_fit']
 501         axes[0].lines[0].set_data(z, deflection)
 502         axes[0].lines[1].set_data(z, deflection_fit)
 503         _pyplot.draw()
 504
 505
 506 class Vibration (Experiment):
 507     def __init__(self, name='vibration', *args, **kwargs):
 508         kwargs['input_parameters'] = [
 509             'alpha_z', 'alpha_d']
 510         kwargs['output_parameters'] = [
 511             'vibration_temperature', 'vibration_variance', 'f_0', 'beta_f',
 512             'G_1f', 'Q', 'kappa']
 513         kwargs['output_parameter_units'] = {
 514             'vibration_temperature': 'K',
 515             'vibration_variance': 'bit^2',
 516             'f_0': 'Hz',
 517             'beta_f': 'Hz',
 518             'G_1f': 'bit^2 Hz^3',
 519             'Q': '',
 520             'kappa': 'N/m',
 521             }
 522         super(Vibration, self).__init__(name=name, *args, **kwargs)
 523
 524     @staticmethod
 525     def _crunch(deflection, temperature, frequency, alpha_z, alpha_d,
 526                 psd_kwargs, fit_kwargs):
 527         freq_axis,power = _FFT_tools.unitary_avg_power_spectrum(
 528             (deflection - deflection.mean()), frequency, **psd_kwargs)
 529         f,beta,G,offset = _calibcant_vibration_analyze.fit_psd(
 530             freq_axis, power, **fit_kwargs)
 531         power_fit = _calibcant_vibration_analyze.breit_wigner(
 532             freq_axis, f, beta, G, offset)
 533         variance = _calibcant_vibration_analyze.breit_wigner_area(f, beta, G)
 534         kappa = 1/variance * (alpha_d*alpha_z)**2 * _kB*temperature
 535         return {
 536             'vibration_temperature': temperature,
 537             'freq_axis': freq_axis,
 538             'power': power,
 539             'f_0': f,
 540             'beta_f': beta,
 541             'G_1f': G,
 542             'Q': f/beta,
 543             'offset': offset,
 544             'power_fit': power_fit,
 545             'vibration_variance': variance,
 546             'kappa': kappa,
 547             }
 548
 549     @staticmethod
 550     def _plot(deflection, freq_axis, power, power_fit, frequency, **kwargs):
 551         figure = _pyplot.figure()
 552         axes = figure.add_subplot(2, 1, 1)
 553         max_time = len(deflection)/frequency
 554         time = _numpy.linspace(0, max_time, len(deflection))
 555         axes.plot(time, deflection, '.')
 556         axes.set_xbound([0, max_time])
 557         axes.set_title('cantilever calibration')
 558         axes.set_xlabel('time (s)')
 559         axes.set_ylabel('deflection (bit)')
 560
 561         axes = figure.add_subplot(2, 1, 2)
 562         axes.hold(True)
 563         axes.set_yscale('log')
 564         axes.plot(freq_axis, power, 'b.')
 565         xmin,xmax = axes.get_xbound()
 566         ymin,ymax = axes.get_ybound()
 567         axes.plot(freq_axis, power_fit, 'r-')
 568         axes.axis([xmin, xmax, ymin, ymax])
 569         axes.set_xlabel('frequency (Hz)')
 570         axes.set_ylabel('deflection PSD (bit^2/Hz)')
 571         return figure
 572
 573     @staticmethod
 574     def _replot(figure, axes, data=None, parameters=None):
 575         oparameters = Experiment._replot(
 576             figure=figure, axes=axes, data=data, parameters=parameters)
 577         deflection = data['deflection']
 578         frequency = data['frequency']
 579         freq_axis = oparameters['freq_axis']
 580         power = oparameters['power']
 581         power_fit = oparameters['power_fit']
 582         max_time = len(deflection)/frequency
 583         time = _numpy.linspace(0, max_time, len(deflection))
 584         axes[0].lines[0].set_data(time, deflection)
 585         axes[1].lines[0].set_data(freq_axis, power)
 586         axes[1].lines[1].set_data(freq_axis, power_fit)
 587         _pyplot.draw()
 588
 589
 590 class Pull (Experiment):
 591     def __init__(self, name='pull', *args, **kwargs):
 592         kwargs['input_parameters'] = [
 593             'alpha_z', 'alpha_d', 'vibration_variance',
 594             'vibration_temperature', 'drift_slope', 'alpha_f_offset',
 595             'smooth_steps', 'start_sawteeth', 'stop_sawteeth',
 596             'contour_space_model', 'contour_space_parameters',
 597             'contour_space_lines']
 598         kwargs['output_parameters'] = [
 599             'pull_alpha_d', 'pull_velocity', 'pull_temperature']
 600         kwargs['output_parameter_units'] = {
 601             'pull_alpha_d': 'diode_bit/piezo_bit', 'pull_velocity': 'm/s',
 602             'pull_temperature': 'K',
 603             }
 604         super(Pull, self).__init__(name=name, *args, **kwargs)
 605
 606     @staticmethod
 607     def _crunch(alpha_v, alpha_f, velocity, temperature,
 608                 alpha_z, alpha_d, vibration_variance, vibration_temperature,
 609                 drift_slope, alpha_f_offset, start_sawteeth, stop_sawteeth,
 610                 contour_space_model, contour_space_parameters,
 611                 **kwargs):
 612         alpha_f = alpha_f - (alpha_v-alpha_v.min())*drift_slope
 613         ddict = {
 614             'approach': {  # for _analyze_surface_position_data
 615                 'z': alpha_v[::-1],
 616                 'deflection': alpha_f[::-1],
 617                 }
 618             }
 619         (non_contact_offset,non_contact_slope,kink,alpha_dp
 620          ) = _pypiezo_surface.analyze_surface_position_data(
 621             ddict, return_all_parameters=True)
 622         # The fitted kink offset, but thrown off by protein presense
 623         #alpha_f0 = non_contact_offset + non_contact_slope*(kink-alpha_v.min())
 624         alpha_f0 = alpha_f[-100:].mean()
 625         alpha_v0 = alpha_v.max() - (alpha_f.max()-alpha_f0)/alpha_dp
 626         d_alpha_f = -(alpha_f - alpha_f0 + alpha_f_offset)
 627         d_alpha_v = -(alpha_v - alpha_v0)
 628         kBT = _kB * vibration_temperature
 629         f = d_alpha_f/vibration_variance*(alpha_d**2 * alpha_z)/alpha_dp * kBT
 630         z_piezo = 1./alpha_z * d_alpha_v
 631         z = 1./alpha_z * (d_alpha_v - d_alpha_f/alpha_dp)
 632         z_contour = contour_values(
 633             z, f, temperature=temperature, model=contour_space_model,
 634             parameters=contour_space_parameters)
 635         in_sawteeth = (z > start_sawteeth) * (z < stop_sawteeth)
 636         z_sawteeth = z[in_sawteeth]
 637         f_sawteeth = f[in_sawteeth]
 638         return {
 639             'pull_alpha_d': alpha_d,
 640             'pull_temperature': temperature,
 641             'pull_velocity': velocity,
 642             'z': z,
 643             'f': f,
 644             'z_piezo': z_piezo,
 645             'z_contour': z_contour,
 646             'z_sawteeth': z[in_sawteeth],
 647             'f_sawteeth': f[in_sawteeth],
 648             }
 649
 650     @staticmethod
 651     def _plot(z_piezo, z, z_contour, deflection=None, f=None,
 652               pull_temperature=None, smooth_steps=None,
 653               xlabel='z piezo (m)', ylabel='deflection (N)',
 654               bounds=None, log=None, slider=None,
 655               contour_space_model=None, contour_space_parameters=None,
 656               contour_space_lines=None, **kwargs):
 657         """Generate pull-plots with three possible abscissa.
 658
 659         The contour-space plots are following Puchner et al., 2008.
 660         http://dx.doi.org/10.1529/biophysj.108.129999
 661         """
 662         if deflection is None:
 663             assert f is not None, "must set 'deflection' or 'f'"
 664             deflection = f
 665         else:
 666             assert f is None, "cannot set both 'f' and 'deflection'"
 667         figure = _pyplot.figure()
 668         axes = figure.add_subplot(1, 1, 1)
 669         axes._z_data = {
 670             'z piezo (m)': z_piezo,
 671             'protein extension (m)': z,
 672             'contour-space (m)': z_contour,
 673             }
 674         try:
 675             zp = axes._z_data[xlabel]
 676         except KeyError as e:
 677             raise NotImplementedError(xlabel)  # from e  (for Python 3)
 678         plots = []
 679         plots.extend(
 680             axes.plot(
 681                 zp, deflection,
 682                 color='b', linestyle='none', marker=',', markeredgecolor='b'))
 683         if smooth_steps:
 684             window_len = int(len(z)/smooth_steps)
 685             axes.hold(True)
 686             z_smooth = Experiment.smooth(z, window_len)
 687             d_smooth = Experiment.smooth(deflection, window_len)
 688             try:
 689                 c_smooth = contour_values(
 690                     z_smooth, d_smooth, temperature=pull_temperature,
 691                     model=contour_space_model,
 692                     parameters=contour_space_parameters)
 693             except ValueError:  # avoid crashing figure.py
 694                 c_smooth = None
 695             axes._smooth_z_data = {
 696                 'z piezo (m)': Experiment.smooth(z_piezo, window_len),
 697                 'protein extension (m)': z_smooth,
 698                 'contour-space (m)': c_smooth,
 699                 }
 700             zp_smooth = axes._smooth_z_data[xlabel]
 701             if xlabel in ['z piezo (m)', 'protein extension (m)']:
 702                 plots.extend(axes.plot(
 703                         zp_smooth, d_smooth, color='r', markeredgecolor='r'))
 704             else:  # 'contour-space (m)'
 705                 plots.extend(axes.plot(
 706                         zp_smooth, d_smooth,
 707                         color='r', linestyle='none', marker=',',
 708                         markeredgecolor='r'))
 709         if xlabel == 'contour-space (m)' and contour_space_lines is not None:
 710             for zi in contour_space_lines:
 711                 axes.axvline(zi, color='g')
 712         if bounds:
 713             axes.axis(bounds)
 714         axes.set_title('velocity clamp unfolding')
 715         axes.set_xlabel(xlabel)
 716         axes.set_ylabel(ylabel)
 717
 718         if slider:  # add radio buttons for plot-type selection
 719             figure.subplots_adjust(left=0.5)
 720             left = 0.05
 721             bottom = 0.2
 722             width = 0.3
 723             height = 0.6
 724             rax = figure.add_axes(
 725                 [left, bottom, width, height], axisbg='none')
 726             labels = (
 727                 'z piezo (m)',
 728                 'protein extension (m)',
 729                 'contour-space (m)',
 730                 )
 731             radio = _widgets.RadioButtons(
 732                 rax, labels=labels, active=labels.index(xlabel))
 733             def change_type(label):  # closure; upvalues: log, plots, axes
 734                 log('change plot to {}'.format(label))
 735                 axes.set_xlabel(label)
 736                 for i,plot in enumerate(plots):
 737                     old_z,deflection = plot.get_data()
 738                     if i == 0:  # regular plot
 739                         z = axes._z_data[label]
 740                     elif i == 1:  # smooth plot
 741                         if label in ['z piezo (m)', 'protein extension (m)']:
 742                             plots[-1].set_linestyle('-')
 743                             plots[-1].set_marker(None)
 744                         else:  # 'contour-space (m)'
 745                             plots[-1].set_linestyle('None')
 746                             plots[-1].set_marker(',')
 747                         z = axes._smooth_z_data[label]
 748                     else:
 749                         break
 750                     plot.set_data(z, deflection)
 751                 _pyplot.draw()
 752                 log('changed plot to {}'.format(label))
 753             axes._callback = change_type  # avoid garbage collection
 754             axes._radio = radio           # avoid garbage collection
 755             radio.on_clicked(change_type)
 756         return figure
 757
 758     @staticmethod
 759     def _replot(figure, axes, data=None, parameters=None):
 760         oparameters = Experiment._replot(
 761             figure=figure, axes=axes, data=data, parameters=parameters)
 762         axes[0]._z_data = {
 763             'z piezo (m)': oparameters['z_piezo'],
 764             'protein extension (m)': oparameters['z'],
 765             'contour-space (m)': oparameters['z_contour'],
 766             }
 767         if parameters.get('smooth_steps', None):
 768             window_len = int(len(oparameters['z'])/parameters['smooth_steps'])
 769             z_smooth = Experiment.smooth(oparameters['z'], window_len)
 770             d_smooth = Experiment.smooth(oparameters['f'], window_len)
 771             axes[0]._smooth_z_data = {
 772                 'z piezo (m)': Experiment.smooth(oparameters['z_piezo']),
 773                 'protein extension (m)': z_smooth,
 774                 'contour-space (m)': contour_values(
 775                     z_smooth, d_smooth,
 776                     temperature=oparameters['pull_temperature'],
 777                     model=parameters['contour_space_model'],
 778                     parameters=parameters['contour_space_parameters']),
 779                 }
 780         for i,plot in enumerate(axes[0].lines):
 781             z,deflection = plot.get_data()
 782             if i == 0:  # regular plot
 783                 deflection = oparameters['f']
 784             elif i == 1:  # smooth plot
 785                 deflection = d_smooth
 786             plot.set_data(z, deflection)
 787         axes[0]._callback(axes[0].get_xlabel())
 788
 789
 790 def inverse_wlc(F):
 791     """
 792       f = kB T / p ( 1/4 * (1/(1-z/L)**2 - 1) + z/L )
 793
 794     With the replacements x = z/L and F = fp/(kB T), this gives
 795
 796       F = ( 1/4 * (1/(1-x)**2 - 1) + x )
 797
 798     After some math (see the sawsim manual for details), we get
 799
 800       0 = x**3 - (F+9/4) x**2 + (2F+3/2) x - F
 801
 802     To get x, just solve this cubic equation.
 803
 804     >>> def wlc(x):
 805     ...     return (0.25*(1/(1-x)**2 - 1) + x)
 806     >>> F = wlc(0.8)
 807     >>> '{:g}'.format(F)
 808     '6.8'
 809     >>> '{:g}'.format(inverse_wlc(F))
 810     '0.8'
 811     """
 812     if F < 0:
 813         return 0
 814     roots = _scipy.roots([1, -(F+9./4), (2*F+3./2), -F])
 815     for r in roots:
 816         if r.imag == 0 and r.real > 0 and r.real < 1:
 817             return r.real
 818     raise ValueError((F, roots))
 819
 820 def inverse_frc(F, gamma, c=2, beta=2, g=float('inf')):
 821     """
 822     From Livadaru et al., 2003, equations 22, 48, and 49 [1].
 823
 824       l = b*cos(γ/2)/|log(cos(γ))|                          (eqn. 22)
 825       F_WLC[x] = 3/4 - 1/x + x^2/4                          (eqn. 48)
 826       x = 1 - {(F_WLC^{-1}[Fl/b])^β + (cF)^β}^{-1/β} + F/g  (eqn. 49)
 827
 828     where b is the segment length, γ is the fixed bond angle, x=R_z/L,
 829     and F = fb/kBT.  They use l for persistence length.
 830
 831     They also use f/g for chain extension instead of F/g, but it's
 832     easy enough to make appropriate adjustments to g, and using F
 833     keeps temperature out of the formula.  Combining equations 22 and
 834     49 to minimize constants:
 835
 836       x = 1 - {(F_WLC^{-1}[F cos(γ/2)/|log(cos(γ))|])^β + (cF)^β}^{-1/β} + F/g
 837
 838     The F_WLC form (48) can be used to reconstruct Bustamante's more
 839     familiar version [2] via Livadaru's equation 47:
 840
 841       fl/(kBT) = F_WLC[(1-R_z/L)^{-1}]
 842
 843     I found Livadaru via Puchner et al., 2008, equation 3 [3].  Note
 844     that there are typos in Livadaru's equation 46 which were
 845     propogated in Puchner's equation 3.  It should have read
 846
 847                { fa/(3kBT)               for fb/(kBT) < b/l
 848       R_z/L ~= { 1-(4fl/(kBT))^{-1/2}    for b/l < fb/(kBT) < l/b
 849                { 1-(cfb/(kBT))^{-1}      for l/b < fb/(kBT)
 850
 851     Also, the form's low-force regime is missing a prefactor to allow
 852     it to match up with the mid-force form at the crossover point.
 853     They avoid this issue by using Bustamante's WLC interpolation
 854     formula [2] for both the low- and mid-force portions in equation
 855     49.
 856
 857     Check against values from Livadaru's figure 14.
 858
 859     >>> for gamma in [1, 20, 70]:
 860     ...     for F in [0.01, 0.1, 1, 10]:
 861     ...         x = inverse_frc(F=F, gamma=gamma*_numpy.pi/180)
 862     ...         print(('(1-Rz/L)^(-1) = {:g}    '
 863     ...                '(for gamma={:g} deg, F={:g})').format(
 864     ...             1./(1-x), gamma, F))
 865     (1-Rz/L)^(-1) = 16.1199    (for gamma=1 deg, F=0.01)
 866     (1-Rz/L)^(-1) = 51.2165    (for gamma=1 deg, F=0.1)
 867     (1-Rz/L)^(-1) = 162.053    (for gamma=1 deg, F=1)
 868     (1-Rz/L)^(-1) = 512.833    (for gamma=1 deg, F=10)
 869     (1-Rz/L)^(-1) = 1.11106    (for gamma=20 deg, F=0.01)
 870     (1-Rz/L)^(-1) = 2.26794    (for gamma=20 deg, F=0.1)
 871     (1-Rz/L)^(-1) = 8.05245    (for gamma=20 deg, F=1)
 872     (1-Rz/L)^(-1) = 32.1006    (for gamma=20 deg, F=10)
 873     (1-Rz/L)^(-1) = 1.0053    (for gamma=70 deg, F=0.01)
 874     (1-Rz/L)^(-1) = 1.07101    (for gamma=70 deg, F=0.1)
 875     (1-Rz/L)^(-1) = 2.56046    (for gamma=70 deg, F=1)
 876     (1-Rz/L)^(-1) = 20.6952    (for gamma=70 deg, F=10)
 877
 878     [1]: http://dx.doi.org/10.1021/ma020751g
 879     [2]: http://dx.doi.org/10.1126/science.8079175
 880     [3]: http://dx.doi.org/10.1529/biophysj.108.129999
 881     """
 882     if F < 0:
 883         return 0
 884     l_over_b = _numpy.cos(gamma/2.)/abs(_numpy.log(_numpy.cos(gamma)))
 885     xWLC = inverse_wlc(F*l_over_b)    # R_z/L for a WLC
 886     invFWLC = 1./(1 - xWLC)           # F_WLC^{-1}[Fl/b]  (see eqn.s 47 and 48)
 887     return 1 - (invFWLC**beta + (c*F)**beta)**(-1./beta) + F/g
 888
 889 def contour_values(z, deflection, temperature=300., model='wlc',
 890                    parameters=None):
 891     """Convert (z, f) points to (L, f) space using a polymer model.
 892
 893     Default parameters are from Puchner et al., 2008 [1].  The FJC
 894     parameters are corrected to match the corrected version of
 895     Livadaru et al., 2003 [2].
 896
 897     On Wed, Sep 12, 2012 at 01:37:23AM +0000, Puchner, Georg Elias Michael wrote:
 898     > Hello Trevor,
 899     >
 900     > thanks for pointing this out, I wasn't aware of this typo and
 901     > also think they should have published and erratum.
 902     >
 903     > I dug into my old data and code (I am not doing force
 904     > spectroscopy andy more since some years) and implemented your
 905     > corrected version. I attached a comparis on of my old
 906     > transformation equation (black line with the parameters γ=22°
 907     > and b=0.4 nm as published) and the fixed equation (dashed red
 908     > line with γ=41° and b=0.11 nm)
 909     >
 910     > [snip extra details]
 911
 912     [1]: http://dx.doi.org/10.1529/biophysj.108.129999
 913     [2]: http://dx.doi.org/10.1021/ma020751g
 914     """
 915     if parameters is None:
 916         parameters = {}
 917     if model == 'wlc':
 918         p = parameters.get('p', 0.4e-9)  # meters, persistence length
 919         d_scale = p/(_kB*temperature)
 920     elif model == 'frc':
 921         b = parameters.get('b', 0.11e-9)  # meters, segment length
 922         gamma = parameters.get(
 923             'gamma', 41*_numpy.pi/180)  # radians, fixed bond angle
 924         d_scale = b/(_kB*temperature)
 925     else:
 926         raise NotImplementedError(model)
 927     L_values = _numpy.zeros(z.shape, dtype=float)
 928     for i,(z_i,d_i) in enumerate(zip(z, deflection)):
 929         if model == 'wlc':
 930             x = inverse_wlc(F=d_i*d_scale)
 931         else:
 932             x = inverse_frc(F=d_i*d_scale, gamma=gamma)
 933         if x == 0:
 934             L_values[i] = 0
 935         else:
 936             L_values[i] = z_i/x  # x = z/L
 937     return L_values
 938
 939 def load_experiments(load_piezo=None, piezo_parameters=None,
 940                      load_bump=None, bump_parameters=None,
 941                      load_vibration=None, vibration_parameters=None,
 942                      load_pull=None, pull_parameters=None,
 943                      override_parameters={}):
 944     if load_piezo is None:
 945         load_piezo = globals()['load_piezo']
 946     if piezo_parameters is None:
 947         piezo_parameters = globals()['piezo_parameters']
 948     if load_bump is None:
 949         load_bump = globals()['load_bump']
 950     if bump_parameters is None:
 951         bump_parameters = globals()['deflection_parameters']
 952     if load_vibration is None:
 953         load_vibration = globals()['load_vibration']
 954     if vibration_parameters is None:
 955         vibration_parameters = globals()['deflection_parameters']
 956     if load_pull is None:
 957         load_pull = globals()['load_pull']
 958     if pull_parameters is None:
 959         pull_parameters = globals()['pull_parameters']
 960
 961     parameters = {}
 962     experiments = []
 963     for experiment,external_parameters in [
 964         (Piezo(load_data=load_piezo), piezo_parameters),
 965         (Bump(load_data=load_bump), deflection_parameters),
 966         (Vibration(load_data=load_vibration), deflection_parameters),
 967         (Pull(load_data=load_pull), pull_parameters),
 968         ]:
 969         external_parameters.update(parameters)
 970         parameters.update(experiment.crunch(external_parameters))
 971         parameters.update(override_parameters)
 972         for exp in experiments:
 973             exp.descendants.append(experiment)
 974         experiments.append(experiment)
 975     return (parameters, experiments)
 976
 977 def print_parameters(parameters, experiments):
 978     for experiment in experiments:
 979         for k in experiment.output_parameters:
 980             if k.startswith('delta_'):
 981                 continue
 982             delta_k = 'delta_{}'.format(k)
 983             if delta_k in experiment.output_parameters:
 984                 print('{} = {:g} +/- {:g} {}'.format(
 985                         k, parameters[k], parameters[delta_k],
 986                         experiment.output_parameter_units[k]))
 987             else:
 988                 print('{} = {:g} {}'.format(
 989                         k, parameters[k],
 990                         experiment.output_parameter_units[k]))
 991
 992
 993 if __name__ == '__main__':
 994     import argparse
 995
 996     parser = argparse.ArgumentParser(description='Play a pure tone')
 997     parser.add_argument(
 998         '-f', '--frequency', type=float, help='tone frequency in Hz')
 999
1000     args = parser.parse_args()
1001
1002     sigma_z = 1/9.2e-9  # V/m
1003     gamma_z = piezo_parameters['z_gain'] * piezo_parameters['output_gamma_z']
1004     alpha_z = sigma_z / gamma_z
1005     print(('hacked alpha_z = {:g} to match {:g} m/V piezo calibration'
1006            ).format(alpha_z, 1/sigma_z))
1007     parameters,experiments = load_experiments(
1008         override_parameters={  # decreasing sigma_z spreads the peaks
1009             'sigma_z': sigma_z,
1010             'alpha_z': alpha_z,
1011             'delta_alpha_z': 0,
1012             })
1013     print_parameters(parameters, experiments)
1014     for experiment in experiments:
1015         experiment.plot(slider=True)
1016     _pyplot.show()