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
  99
 100 def analyze(bumps, temperatures, vibrations):
 101     """Analyze data from `get_calibration_data()`
 102
 103     Inputs (all are arrays of recorded data):
 104       bumps measured (V_photodiode / nm_tip) proportionality constant
 105       temperatures    measured temperature (K)
 106       vibrations  measured V_photodiode variance in free solution (V**2)
 107     Outputs:
 108       k    cantilever spring constant (in N/m, or equivalently nN/nm)
 109       k_s  standard deviation in our estimate of k
 110
 111     Notes:
 112
 113     We're assuming vib is mostly from thermal cantilever vibrations
 114     (and then only from vibrations in the single vertical degree of
 115     freedom), and not from other noise sources.
 116
 117     If the error is large, check the relative errors
 118     (`x.std()/x.mean()`)of your input arrays.  If one of them is
 119     small, don't bother repeating that measurment too often.  If one
 120     is large, try repeating that measurement more.  Remember that you
 121     need enough samples to have a valid error estimate in the first
 122     place, and that none of this addresses any systematic errors.
 123     """
 124     ps_m = bumps.mean()  # ps for photo-sensitivity
 125     ps_s = bumps.std()
 126     T_m = temperatures.mean()
 127     T_s = temperatures.std()
 128     v2_m = vibrations.mean()  # average voltage variance
 129     v2_s = vibrations.std()
 130
 131     # Vphoto / photo_sensitivity = x
 132     # k = kB T / <x**2> = kB T photo_sensitivity**2 / Vphoto_var
 133     #
 134     # units,  photo_sensitivity =  Vphoto/(Zcant in m),
 135     # so Vphoto/photo_sensitivity = Zcant in m
 136     # so k = J/K * K / m^2 = J / m^2 = N/m
 137     k  = _kB * T_m * ps_m**2 / v2_m
 138
 139     # propogation of errors
 140     # dk/dT = k/T
 141     dk_T = k/T_m * T_s
 142     # dk/dps = 2k/ps
 143     dk_ps = 2*k/ps_m * ps_s
 144     # dk/dv2 = -k/v2
 145     dk_v = -k/v2_m * v2_s
 146
 147     k_s = _numpy.sqrt(dk_T**2 + dk_ps**2 + dk_v**2)
 148
 149     _LOG.info('variable (units)         : '
 150               'mean +/- std. dev. (relative error)')
 151     _LOG.info('cantilever k (N/m)       : %g +/- %g (%g)' % (k, k_s, k_s/k))
 152     _LOG.info('photo sensitivity (V/m)  : %g +/- %g (%g)'
 153               % (ps_m, ps_s, ps_s/ps_m))
 154     _LOG.info('T (K)                    : %g +/- %g (%g)'
 155               % (T_m, T_s, T_s/T_m))
 156     _LOG.info('vibration variance (V^2) : %g +/- %g (%g)'
 157               % (v2_m, v2_s, v2_s/v2_m))
 158
 159     if _package_config['matplotlib']:
 160         plot(bumps, temperatures, vibrations)
 161
 162     return (k, k_s)
 163
 164
 165 def plot(bumps, temperatures, vibrations):
 166     if not _matplotlib:
 167         raise _matplotlib_import_error
 168     figure = _matplotlib_pyplot.figure()
 169
 170     bump_axes = figure.add_subplot(3, 1, 1)
 171     T_axes = figure.add_subplot(3, 1, 2)
 172     vib_axes = figure.add_subplot(3, 1, 3)
 173
 174     timestamp = _time.strftime('%H%M%S')
 175     bump_axes.set_title('cantilever calibration %s' % timestamp)
 176
 177     bump_axes.autoscale(tight=True)
 178     bump_axes.plot(bumps, 'g.-')
 179     bump_axes.set_ylabel('photodiode sensitivity (V/m)')
 180     T_axes.autoscale(tight=True)
 181     T_axes.plot(temperatures, 'r.-')
 182     T_axes.set_ylabel('temperature (K)')
 183     vib_axes.autoscale(tight=True)
 184     vib_axes.plot(vibrations, 'b.-')
 185     vib_axes.set_ylabel('thermal deflection variance (V^2)')
 186
 187     if hasattr(figure, 'show'):
 188         figure.show()
 189     return figure
 190 _plot = plot  # alternative name for use inside analyze_all()
 191
 192
 193 def analyze_all(config, data, raw_data, maximum_relative_error=1e-5,
 194                 filename=None, group=None, plot=False, dry_run=False):
 195     "(Re)analyze (and possibly plot) all data from a `calib()` run."
 196     if not data.get('bump', None):
 197         data['bump'] = _numpy.zeros((config['num-bumps'],), dtype=float)
 198     if not data.get('temperature', None):
 199         data['temperature'] = _numpy.zeros(
 200             (config['num-temperatures'],), dtype=float)
 201     if not data.get('vibrations', None):
 202         data['vibration'] = _numpy.zeros(
 203                 (config['num-vibrations'],), dtype=float)
 204     axis_config = config['afm']['piezo'].select_config(
 205         setting_name='axes',
 206         attribute_value=config['afm']['main-axis'],
 207         get_attribute=_get_axis_name)
 208     input_config = config['afm']['piezo'].select_config(
 209         setting_name='inputs', attribute_value='deflection')
 210     bumps_changed = temperatures_changed = vibrations_changed = False
 211     if not isinstance(group, _h5py.Group) and not dry_run:
 212         f = _h5py.File(filename, mode)
 213         group = _h5_create_group(f, group)
 214     else:
 215         f = None
 216     try:
 217         for i,bump in enumerate(raw_data['bump']):
 218             data['bump'][i],changed = check_bump(
 219                 index=i, bump=bump, z_axis_config=axis_config,
 220                 deflection_channel_config=input_config, plot=plot,
 221                 maximum_relative_error=maximum_relative_error)
 222             if changed and not dry_run:
 223                 bumps_changed = True
 224                 bump_group = _h5_create_group(group, 'bump/{}'.format(i))
 225                 _bump_save(group=bump_group, processed=data['bump'][i])
 226         for i,temperature in enumerate(raw_data['temperature']):
 227             data['temperature'][i],changed = check_temperature(
 228                 index=i, temperature=temperature,
 229                 maximum_relative_error=maximum_relative_error)
 230             if changed and not dry_run:
 231                 temperatures_changed = True
 232                 temperature_group = _h5_create_group(
 233                     group, 'temperature/{}'.format(i))
 234                 _temperature_save(
 235                     group=temerature_group, processed=data['temperature'][i])
 236         for i,vibration in enumerate(raw_data['vibration']):
 237             data['vibration'][i],changed = check_vibration(
 238                     index=i, vibration=vibration,
 239                     deflection_channel_config=input_config, plot=plot,
 240                     maximum_relative_error=maximum_relative_error)
 241             if changed and not dry_run:
 242                 vibrations_changed = True
 243                 vibration_group = _h5_create_group(
 244                     group, 'vibration/{}'.format(i))
 245                 _vibration_save(
 246                     group=vibration_group, processed=data['vibration'])
 247         k,k_s,changed = check_calibration(
 248             k=data['processed']['spring_constant'],
 249             k_s=data['processed']['spring_constant_deviation'],
 250             bumps=data['bump'],
 251             temperatures=data['temperature'], vibrations=data['vibration'],
 252             maximum_relative_error=maximum_relative_error)
 253         if (changed or bumps_changed or temperatures_changed or
 254             vibrations_changed) and not dry_run:
 255             calibration_group = _h5_create_group(group, 'calibration')
 256             if bumps_changed:
 257                 calib_save(group=calibration_group, bump=data['bump'])
 258             if temperatures_changed:
 259                 calib_save(
 260                     group=calibration_group, temperature=data['temperature'])
 261             if vibrations_changed:
 262                 calib_save(
 263                     group=calibration_group, vibration=data['vibration'])
 264             if changed:
 265                 calib_save(group=calibration_group, k=k, k_s=k_s)
 266     finally:
 267         if f:
 268             f.close()
 269     if plot:
 270         _plot(bumps=data['raw']['bump'],
 271              temperatures=data['raw']['temperature'],
 272              vibrations=data['raw']['vibration'])
 273     return (k, k_s)
 274
 275 def check_bump(index, bump, maximum_relative_error, **kwargs):
 276     changed = False
 277     sensitivity = _bump_analyze(
 278         config=bump['config']['bump'], data=bump['raw'], **kwargs)
 279     if bump.get('processed', None) is None:
 280         changed = True
 281         _LOG.warn('new analysis for bump {}: {}'.format(index, sensitivity))
 282     else:
 283         rel_error = abs(sensitivity - bump['processed'])/bump['processed']
 284         if rel_error > maximum_relative_error:
 285             changed = True
 286             _LOG.warn(("new analysis doesn't match for bump {}: {} -> {} "
 287                        "(difference: {}, relative error: {})").format(
 288                     index, bump['processed'], sensitivity,
 289                     sensitivity-bump['processed'], rel_error))
 290     return (sensitivity, changed)
 291
 292 def check_temperature(index, temperature, maximum_relative_error, **kwargs):
 293     changed = False
 294     temp = _temperature_analyze(
 295         config=temperature['config']['temperature'],
 296         temperature=temperature['raw'], **kwargs)
 297     if temperature.get('processed', None) is None:
 298         changed = True
 299         _LOG.warn('new analysis for temperature {}: {}'.format(index, temp))
 300     else:
 301         rel_error = abs(temp - temperature['processed']
 302                         )/temperature['processed']
 303         if rel_error > maximum_relative_error:
 304             changed = True
 305             _LOG.warn(("new analysis doesn't match for temperature "
 306                        "{} -> {} (difference: {}, relative error: {})"
 307                        ).format(
 308                     index, temperature['processed'], temp,
 309                     temp-temperature['processed'], rel_error))
 310     return (temp, changed)
 311
 312 def check_vibration(index, vibration, maximum_relative_error, **kwargs):
 313     changed = False
 314     variance = _vibration_analyze(
 315         config=vibration['config']['vibration'],
 316         deflection=vibration['raw'], **kwargs)
 317     if vibration.get('processed', None) is None:
 318         changed = True
 319         _LOG.warn('new analysis for temperature {}: {}'.format(
 320                 index, variance))
 321     else:
 322         rel_error = abs(variance-vibration['processed'])/vibration['processed']
 323         if rel_error > maximum_relative_error:
 324             _LOG.warn(("new analysis doesn't match for vibration {}: {} != {} "
 325                        "(difference: {}, relative error: {})").format(
 326                     index, variance, vibration['processed'],
 327                     variance-vibration['processed'], rel_error))
 328     return (variance, changed)
 329
 330 def check_calibration(k, k_s, maximum_relative_error, **kwargs):
 331     changed = False
 332     new_k,new_k_s = analyze(**kwargs)
 333     if k is None:
 334         changed = True
 335         _LOG.warn('new analysis for the spring constant: {}'.format(new_k))
 336     else:
 337         rel_error = abs(new_k-k)/k
 338         if rel_error > maximum_relative_error:
 339             _LOG.warn(("new analysis doesn't match for the spring constant: "
 340                        "{} != {} (difference: {}, relative error: {})").format(
 341                     new_k, k, new_k-k, rel_error))
 342     if k_s is None:
 343         changed = True
 344         _LOG.warn('new analysis for the spring constant deviation: {}'.format(
 345                 new_k_s))
 346     else:
 347         rel_error = abs(new_k-k)/k
 348         if rel_error > maximum_relative_error:
 349             _LOG.warn(
 350                 ("new analysis doesn't match for the spring constant deviation"
 351                  ": {} != {} (difference: {}, relative error: {})").format(
 352                     new_k_s, k_s, new_k_s-k_s, rel_error))
 353     return (new_k, new_k_s, changed)