calibcant/analyze.py

   1 # calibcant - tools for thermally calibrating AFM cantilevers
   2 #
   3 # Copyright (C) 2008-2012 W. Trevor King <wking@drexel.edu>
   4 #
   5 # This file is part of calibcant.
   6 #
   7 # calibcant is free software: you can redistribute it and/or modify it under
   8 # the terms of the GNU General Public License as published by the Free Software
   9 # Foundation, either version 3 of the License, or (at your option) any later
  10 # version.
  11 #
  12 # calibcant is distributed in the hope that it will be useful, but WITHOUT ANY
  13 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  14 # A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
  15 #
  16 # You should have received a copy of the GNU General Public License along with
  17 # calibcant.  If not, see <http://www.gnu.org/licenses/>.
  18
  19 """Calculate `k` from arrays of bumps, temperatures, and vibrations.
  20
  21 Separate the more general `analyze()` from the other calibration
  22 functions in calibcant.
  23
  24 The relevent physical quantities are :
  25   Vzp_out  Output z-piezo voltage (what we generate)
  26   Vzp      Applied z-piezo voltage (after external ZPGAIN)
  27   Zp       The z-piezo position
  28   Zcant    The cantilever vertical deflection
  29   Vphoto   The photodiode vertical deflection voltage (what we measure)
  30   Fcant    The force on the cantilever
  31   T        The temperature of the cantilever and surrounding solution
  32            (another thing we measure)
  33   k_b      Boltzmann's constant
  34
  35 Which are related by the parameters:
  36   zp_gain           Vzp_out / Vzp
  37   zp_sensitivity    Zp / Vzp
  38   photo_sensitivity Vphoto / Zcant
  39   k_cant            Fcant / Zcant
  40
  41
  42 >>> import numpy
  43 >>> from .config import CalibrateConfig
  44
  45 >>> config = CalibrateConfig()
  46 >>> bumps = numpy.array((15.9e6, 16.9e6, 16.3e6))
  47 >>> temperatures = numpy.array((295, 295.2, 294.8))
  48 >>> vibrations = numpy.array((2.20e-5, 2.22e-5, 2.21e-5))
  49
  50 >>> k,k_s = analyze(bumps=bumps, temperatures=temperatures,
  51 ...     vibrations=vibrations)
  52 >>> (k, k_s)  # doctest: +ELLIPSIS
  53 (0.0493..., 0.00248...)
  54
  55 Most of the error in this example comes from uncertainty in the
  56 photodiode sensitivity (bumps).
  57
  58 >>> k_s/k  # doctest: +ELLIPSIS
  59 0.0503...
  60 >>> bumps.std()/bumps.mean()  # doctest: +ELLIPSIS
  61 0.0251...
  62 >>> temperatures.std()/temperatures.mean()  # doctest: +ELLIPSIS
  63 0.000553...
  64 >>> vibrations.std()/vibrations.mean()  # doctest: +ELLIPSIS
  65 0.00369...
  66 """
  67
  68 import h5py as _h5py
  69 import numpy as _numpy
  70 try:
  71     from scipy.constants import Boltzmann as _kB  # in J/K
  72 except ImportError:
  73     from scipy.constants import Bolzmann as _kB  # in J/K
  74 # Bolzmann -> Boltzmann patch submitted:
  75 #   http://projects.scipy.org/scipy/ticket/1417
  76 # Fixed in scipy commit 4716d91, Apr 2, 2011, during work after v0.9.0rc5.
  77
  78 try:
  79     import matplotlib as _matplotlib
  80     import matplotlib.pyplot as _matplotlib_pyplot
  81     import time as _time  # for timestamping lines on plots
  82 except (ImportError, RuntimeError), e:
  83     _matplotlib = None
  84     _matplotlib_import_error = e
  85
  86 from h5config.storage.hdf5 import h5_create_group as _h5_create_group
  87 from pypiezo.base import get_axis_name as _get_axis_name
  88
  89 from . import LOG as _LOG
  90 from . import package_config as _package_config
  91
  92 from .bump_analyze import analyze as _bump_analyze
  93 from .bump_analyze import save as _bump_save
  94 from .temperature_analyze import analyze as _temperature_analyze
  95 from .temperature_analyze import save as _temperature_save
  96 from .vibration_analyze import analyze as _vibration_analyze
  97 from .vibration_analyze import save as _vibration_save
  98 from .util import SaveSpec as _SaveSpec
  99 from .util import save as _save
 100
 101
 102 def analyze(bumps, temperatures, vibrations):
 103     """Analyze data from `get_calibration_data()`
 104
 105     Inputs (all are arrays of recorded data):
 106       bumps measured (V_photodiode / nm_tip) proportionality constant
 107       temperatures    measured temperature (K)
 108       vibrations  measured V_photodiode variance in free solution (V**2)
 109     Outputs:
 110       k    cantilever spring constant (in N/m, or equivalently nN/nm)
 111       k_s  standard deviation in our estimate of k
 112
 113     Notes:
 114
 115     We're assuming vib is mostly from thermal cantilever vibrations
 116     (and then only from vibrations in the single vertical degree of
 117     freedom), and not from other noise sources.
 118
 119     If the error is large, check the relative errors
 120     (`x.std()/x.mean()`)of your input arrays.  If one of them is
 121     small, don't bother repeating that measurment too often.  If one
 122     is large, try repeating that measurement more.  Remember that you
 123     need enough samples to have a valid error estimate in the first
 124     place, and that none of this addresses any systematic errors.
 125     """
 126     ps_m = bumps.mean()  # ps for photo-sensitivity
 127     ps_s = bumps.std()
 128     T_m = temperatures.mean()
 129     T_s = temperatures.std()
 130     v2_m = vibrations.mean()  # average voltage variance
 131     v2_s = vibrations.std()
 132
 133     # Vphoto / photo_sensitivity = x
 134     # k = kB T / <x**2> = kB T photo_sensitivity**2 / Vphoto_var
 135     #
 136     # units,  photo_sensitivity =  Vphoto/(Zcant in m),
 137     # so Vphoto/photo_sensitivity = Zcant in m
 138     # so k = J/K * K / m^2 = J / m^2 = N/m
 139     k  = _kB * T_m * ps_m**2 / v2_m
 140
 141     # propogation of errors
 142     # dk/dT = k/T
 143     dk_T = k/T_m * T_s
 144     # dk/dps = 2k/ps
 145     dk_ps = 2*k/ps_m * ps_s
 146     # dk/dv2 = -k/v2
 147     dk_v = -k/v2_m * v2_s
 148
 149     k_s = _numpy.sqrt(dk_T**2 + dk_ps**2 + dk_v**2)
 150
 151     _LOG.info('variable (units)         : '
 152               'mean +/- std. dev. (relative error)')
 153     _LOG.info('cantilever k (N/m)       : %g +/- %g (%g)' % (k, k_s, k_s/k))
 154     _LOG.info('photo sensitivity (V/m)  : %g +/- %g (%g)'
 155               % (ps_m, ps_s, ps_s/ps_m))
 156     _LOG.info('T (K)                    : %g +/- %g (%g)'
 157               % (T_m, T_s, T_s/T_m))
 158     _LOG.info('vibration variance (V^2) : %g +/- %g (%g)'
 159               % (v2_m, v2_s, v2_s/v2_m))
 160
 161     if _package_config['matplotlib']:
 162         plot(bumps, temperatures, vibrations)
 163
 164     return (k, k_s)
 165
 166
 167 def plot(bumps, temperatures, vibrations):
 168     if not _matplotlib:
 169         raise _matplotlib_import_error
 170     figure = _matplotlib_pyplot.figure()
 171
 172     bump_axes = figure.add_subplot(3, 1, 1)
 173     T_axes = figure.add_subplot(3, 1, 2)
 174     vib_axes = figure.add_subplot(3, 1, 3)
 175
 176     timestamp = _time.strftime('%H%M%S')
 177     bump_axes.set_title('cantilever calibration %s' % timestamp)
 178
 179     bump_axes.autoscale(tight=True)
 180     bump_axes.plot(bumps, 'g.-')
 181     bump_axes.set_ylabel('photodiode sensitivity (V/m)')
 182     T_axes.autoscale(tight=True)
 183     T_axes.plot(temperatures, 'r.-')
 184     T_axes.set_ylabel('temperature (K)')
 185     vib_axes.autoscale(tight=True)
 186     vib_axes.plot(vibrations, 'b.-')
 187     vib_axes.set_ylabel('thermal deflection variance (V^2)')
 188
 189     if hasattr(figure, 'show'):
 190         figure.show()
 191     return figure
 192 _plot = plot  # alternative name for use inside analyze_all()
 193
 194 def save_results(filename=None, group='/', bump=None,
 195                  temperature=None, vibration=None, spring_constant=None,
 196                  spring_constant_deviation=None):
 197     specs = [
 198         _SaveSpec(item=bump, relpath='raw/photodiode-sensitivity',
 199                   array=True, units='V/m'),
 200         _SaveSpec(item=temperature, relpath='raw/temperature',
 201                   array=True, units='K'),
 202         _SaveSpec(item=vibration, relpath='raw/vibration',
 203                   array=True, units='V^2/Hz'),
 204         _SaveSpec(item=spring_constant, relpath='processed/spring-constant',
 205                   units='N/m', deviation=spring_constant_deviation),
 206         ]
 207     _save(filename=filename, group=group, specs=specs)
 208
 209 def analyze_all(config, data, raw_data, maximum_relative_error=1e-5,
 210                 filename=None, group=None, plot=False, dry_run=False):
 211     "(Re)analyze (and possibly plot) all data from a `calib()` run."
 212     if not data.get('bump', None):
 213         data['bump'] = _numpy.zeros((config['num-bumps'],), dtype=float)
 214     if not data.get('temperature', None):
 215         data['temperature'] = _numpy.zeros(
 216             (config['num-temperatures'],), dtype=float)
 217     if not data.get('vibrations', None):
 218         data['vibration'] = _numpy.zeros(
 219                 (config['num-vibrations'],), dtype=float)
 220     axis_config = config['afm']['piezo'].select_config(
 221         setting_name='axes',
 222         attribute_value=config['afm']['main-axis'],
 223         get_attribute=_get_axis_name)
 224     input_config = config['afm']['piezo'].select_config(
 225         setting_name='inputs', attribute_value='deflection')
 226     bumps_changed = temperatures_changed = vibrations_changed = False
 227     if not isinstance(group, _h5py.Group) and not dry_run:
 228         f = _h5py.File(filename, mode='a')
 229         group = _h5_create_group(f, group)
 230     else:
 231         f = None
 232     try:
 233         if len(data.get('raw/bump', [])) != len(data['bump']):
 234             bumps_changed = True
 235         for i,bump in enumerate(raw_data['bump']):
 236             data['bump'][i],changed = check_bump(
 237                 index=i, bump=bump, z_axis_config=axis_config,
 238                 deflection_channel_config=input_config, plot=plot,
 239                 maximum_relative_error=maximum_relative_error)
 240             if changed and not dry_run:
 241                 bumps_changed = True
 242                 bump_group = _h5_create_group(group, 'bump/{}'.format(i))
 243                 _bump_save(group=bump_group, processed=data['bump'][i])
 244         if len(data.get('raw/temperature', [])) != len(data['temperature']):
 245             temperatures_changed = True
 246         for i,temperature in enumerate(raw_data['temperature']):
 247             data['temperature'][i],changed = check_temperature(
 248                 index=i, temperature=temperature,
 249                 maximum_relative_error=maximum_relative_error)
 250             if changed and not dry_run:
 251                 temperatures_changed = True
 252                 temperature_group = _h5_create_group(
 253                     group, 'temperature/{}'.format(i))
 254                 _temperature_save(
 255                     group=temerature_group, processed=data['temperature'][i])
 256         if len(data.get('raw/vibration', [])) != len(data['vibration']):
 257             vibrations_changed = True
 258         for i,vibration in enumerate(raw_data['vibration']):
 259             data['vibration'][i],changed = check_vibration(
 260                     index=i, vibration=vibration,
 261                     deflection_channel_config=input_config, plot=plot,
 262                     maximum_relative_error=maximum_relative_error)
 263             if changed and not dry_run:
 264                 vibrations_changed = True
 265                 vibration_group = _h5_create_group(
 266                     group, 'vibration/{}'.format(i))
 267                 _vibration_save(
 268                     group=vibration_group, processed=data['vibration'])
 269         k,k_s,changed = check_calibration(
 270             k=data.get('processed/spring_constant', None),
 271             k_s=data.get('processed/spring_constant_deviation', None),
 272             bumps=data['bump'],
 273             temperatures=data['temperature'], vibrations=data['vibration'],
 274             maximum_relative_error=maximum_relative_error)
 275         if (changed or bumps_changed or temperatures_changed or
 276             vibrations_changed) and not dry_run:
 277             calibration_group = _h5_create_group(group, 'calibration')
 278             if bumps_changed:
 279                 save_results(
 280                     group=calibration_group, bump=data['bump'])
 281             if temperatures_changed:
 282                 save_results(
 283                     group=calibration_group, temperature=data['temperature'])
 284             if vibrations_changed:
 285                 save_results(
 286                     group=calibration_group, vibration=data['vibration'])
 287             if changed:
 288                 save_results(
 289                     group=calibration_group,
 290                     spring_constant=k, spring_constant_deviation=k_s)
 291     finally:
 292         if f:
 293             f.close()
 294     if plot:
 295         _plot(bumps=data['raw']['bump'],
 296              temperatures=data['raw']['temperature'],
 297              vibrations=data['raw']['vibration'])
 298     return (k, k_s)
 299
 300 def check_bump(index, bump, maximum_relative_error, **kwargs):
 301     changed = False
 302     sensitivity = _bump_analyze(
 303         config=bump['config']['bump'], data=bump['raw'], **kwargs)
 304     if bump.get('processed', None) is None:
 305         changed = True
 306         _LOG.warn('new analysis for bump {}: {}'.format(index, sensitivity))
 307     else:
 308         rel_error = abs(sensitivity - bump['processed'])/bump['processed']
 309         if rel_error > maximum_relative_error:
 310             changed = True
 311             _LOG.warn(("new analysis doesn't match for bump {}: {} -> {} "
 312                        "(difference: {}, relative error: {})").format(
 313                     index, bump['processed'], sensitivity,
 314                     sensitivity-bump['processed'], rel_error))
 315     return (sensitivity, changed)
 316
 317 def check_temperature(index, temperature, maximum_relative_error, **kwargs):
 318     changed = False
 319     temp = _temperature_analyze(
 320         config=temperature['config']['temperature'],
 321         temperature=temperature['raw'], **kwargs)
 322     if temperature.get('processed', None) is None:
 323         changed = True
 324         _LOG.warn('new analysis for temperature {}: {}'.format(index, temp))
 325     else:
 326         rel_error = abs(temp - temperature['processed']
 327                         )/temperature['processed']
 328         if rel_error > maximum_relative_error:
 329             changed = True
 330             _LOG.warn(("new analysis doesn't match for temperature "
 331                        "{} -> {} (difference: {}, relative error: {})"
 332                        ).format(
 333                     index, temperature['processed'], temp,
 334                     temp-temperature['processed'], rel_error))
 335     return (temp, changed)
 336
 337 def check_vibration(index, vibration, maximum_relative_error, **kwargs):
 338     changed = False
 339     variance = _vibration_analyze(
 340         config=vibration['config']['vibration'],
 341         deflection=vibration['raw'], **kwargs)
 342     if vibration.get('processed', None) is None:
 343         changed = True
 344         _LOG.warn('new analysis for temperature {}: {}'.format(
 345                 index, variance))
 346     else:
 347         rel_error = abs(variance-vibration['processed'])/vibration['processed']
 348         if rel_error > maximum_relative_error:
 349             _LOG.warn(("new analysis doesn't match for vibration {}: {} != {} "
 350                        "(difference: {}, relative error: {})").format(
 351                     index, variance, vibration['processed'],
 352                     variance-vibration['processed'], rel_error))
 353     return (variance, changed)
 354
 355 def check_calibration(k, k_s, maximum_relative_error, **kwargs):
 356     changed = False
 357     new_k,new_k_s = analyze(**kwargs)
 358     if k is None:
 359         changed = True
 360         _LOG.warn('new analysis for the spring constant: {}'.format(new_k))
 361     else:
 362         rel_error = abs(new_k-k)/k
 363         if rel_error > maximum_relative_error:
 364             changed = True
 365             _LOG.warn(("new analysis doesn't match for the spring constant: "
 366                        "{} != {} (difference: {}, relative error: {})").format(
 367                     new_k, k, new_k-k, rel_error))
 368     if k_s is None:
 369         changed = True
 370         _LOG.warn('new analysis for the spring constant deviation: {}'.format(
 371                 new_k_s))
 372     else:
 373         rel_error = abs(new_k_s-k_s)/k_s
 374         if rel_error > maximum_relative_error:
 375             changed = True
 376             _LOG.warn(
 377                 ("new analysis doesn't match for the spring constant deviation"
 378                  ": {} != {} (difference: {}, relative error: {})").format(
 379                     new_k_s, k_s, new_k_s-k_s, rel_error))
 380     return (new_k, new_k_s, changed)