calibcant/analyze.py

   1 # calibcant - tools for thermally calibrating AFM cantilevers
   2 #
   3 # Copyright (C) 2008-2013 W. Trevor King <wking@tremily.us>
   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     if ps_m == 0:
 134         raise ValueError('invalid bumps: {}'.format(bumps))
 135     if T_m == 0:
 136         raise ValueError('invalid temperatures: {}'.format(temperatures))
 137     if v2_m == 0:
 138         raise ValueError('invalid vibrations: {}'.format(vibrations))
 139
 140     # Vphoto / photo_sensitivity = x
 141     # k = kB T / <x**2> = kB T photo_sensitivity**2 / Vphoto_var
 142     #
 143     # units,  photo_sensitivity =  Vphoto/(Zcant in m),
 144     # so Vphoto/photo_sensitivity = Zcant in m
 145     # so k = J/K * K / m^2 = J / m^2 = N/m
 146     k  = _kB * T_m * ps_m**2 / v2_m
 147
 148     # propogation of errors
 149     # dk/dT = k/T
 150     dk_T = k/T_m * T_s
 151     # dk/dps = 2k/ps
 152     dk_ps = 2*k/ps_m * ps_s
 153     # dk/dv2 = -k/v2
 154     dk_v = -k/v2_m * v2_s
 155
 156     k_s = _numpy.sqrt(dk_T**2 + dk_ps**2 + dk_v**2)
 157
 158     _LOG.info('variable (units)         : '
 159               'mean +/- std. dev. (relative error)')
 160     _LOG.info('cantilever k (N/m)       : %g +/- %g (%g)' % (k, k_s, k_s/k))
 161     _LOG.info('photo sensitivity (V/m)  : %g +/- %g (%g)'
 162               % (ps_m, ps_s, ps_s/ps_m))
 163     _LOG.info('T (K)                    : %g +/- %g (%g)'
 164               % (T_m, T_s, T_s/T_m))
 165     _LOG.info('vibration variance (V^2) : %g +/- %g (%g)'
 166               % (v2_m, v2_s, v2_s/v2_m))
 167
 168     if _package_config['matplotlib']:
 169         plot(bumps, temperatures, vibrations)
 170
 171     return (k, k_s)
 172
 173
 174 def plot(bumps, temperatures, vibrations):
 175     if not _matplotlib:
 176         raise _matplotlib_import_error
 177     figure = _matplotlib_pyplot.figure()
 178
 179     bump_axes = figure.add_subplot(3, 1, 1)
 180     T_axes = figure.add_subplot(3, 1, 2)
 181     vib_axes = figure.add_subplot(3, 1, 3)
 182
 183     timestamp = _time.strftime('%H%M%S')
 184     bump_axes.set_title('cantilever calibration %s' % timestamp)
 185
 186     bump_axes.autoscale(tight=True)
 187     bump_axes.plot(bumps, 'g.-')
 188     bump_axes.set_ylabel('photodiode sensitivity (V/m)')
 189     T_axes.autoscale(tight=True)
 190     T_axes.plot(temperatures, 'r.-')
 191     T_axes.set_ylabel('temperature (K)')
 192     vib_axes.autoscale(tight=True)
 193     vib_axes.plot(vibrations, 'b.-')
 194     vib_axes.set_ylabel('thermal deflection variance (V^2)')
 195
 196     if hasattr(figure, 'show'):
 197         figure.show()
 198     return figure
 199 _plot = plot  # alternative name for use inside analyze_all()
 200
 201 def save_results(filename=None, group='/', bump=None,
 202                  temperature=None, vibration=None, spring_constant=None,
 203                  spring_constant_deviation=None):
 204     specs = [
 205         _SaveSpec(item=bump, relpath='raw/photodiode-sensitivity',
 206                   array=True, units='V/m'),
 207         _SaveSpec(item=temperature, relpath='raw/temperature',
 208                   array=True, units='K'),
 209         _SaveSpec(item=vibration, relpath='raw/vibration',
 210                   array=True, units='V^2/Hz'),
 211         _SaveSpec(item=spring_constant, relpath='processed/spring-constant',
 212                   units='N/m', deviation=spring_constant_deviation),
 213         ]
 214     _save(filename=filename, group=group, specs=specs)
 215
 216 def analyze_all(config, data, raw_data, maximum_relative_error=1e-5,
 217                 filename=None, group=None, plot=False, dry_run=False):
 218     "(Re)analyze (and possibly plot) all data from a `calib()` run."
 219     if not data.get('bump', None):
 220         data['bump'] = _numpy.zeros((config['num-bumps'],), dtype=float)
 221     if not data.get('temperature', None):
 222         data['temperature'] = _numpy.zeros(
 223             (config['num-temperatures'],), dtype=float)
 224     if not data.get('vibrations', None):
 225         data['vibration'] = _numpy.zeros(
 226                 (config['num-vibrations'],), dtype=float)
 227     if 'raw' not in data:
 228         data['raw'] = {}
 229     if 'bump' not in data['raw']:
 230         data['raw']['bump'] = _numpy.zeros((config['num-bumps'],), dtype=float)
 231     if 'temperature' not in data['raw']:
 232         data['raw']['temperature'] = _numpy.zeros(
 233         (config['num-temperatures'],), dtype=float)
 234     if 'vibration' not in data['raw']:
 235         data['raw']['vibration'] = _numpy.zeros(
 236         (config['num-vibrations'],), dtype=float)
 237     axis_config = config['afm']['piezo'].select_config(
 238         setting_name='axes',
 239         attribute_value=config['afm']['main-axis'],
 240         get_attribute=_get_axis_name)
 241     input_config = config['afm']['piezo'].select_config(
 242         setting_name='inputs', attribute_value='deflection')
 243     calibration_group = None
 244     if not isinstance(group, _h5py.Group) and not dry_run:
 245         f = _h5py.File(filename, mode='a')
 246         group = _h5_create_group(f, group)
 247     else:
 248         f = None
 249     try:
 250         bumps_changed = len(data['raw']['bump']) != len(data['bump'])
 251         for i,bump in enumerate(raw_data.get('bump', [])):  # compare values
 252             data['bump'][i],changed = check_bump(
 253                 index=i, bump=bump, config=config, z_axis_config=axis_config,
 254                 deflection_channel_config=input_config, plot=plot,
 255                 maximum_relative_error=maximum_relative_error)
 256             if changed and not dry_run:
 257                 bumps_changed = True
 258                 bump_group = _h5_create_group(group, 'bump/{}'.format(i))
 259                 _bump_save(group=bump_group, processed=data['bump'][i])
 260         temperatures_changed = len(data['raw']['temperature']) != len(
 261             data['temperature'])
 262         for i,temperature in enumerate(raw_data.get('temperature', [])):
 263             data['temperature'][i],changed = check_temperature(
 264                 index=i, temperature=temperature, config=config,
 265                 maximum_relative_error=maximum_relative_error)
 266             if changed and not dry_run:
 267                 temperatures_changed = True
 268                 temperature_group = _h5_create_group(
 269                     group, 'temperature/{}'.format(i))
 270                 _temperature_save(
 271                     group=temperature_group, processed=data['temperature'][i])
 272         vibrations_changed = len(data['raw']['vibration']) != len(
 273             data['vibration'])
 274         for i,vibration in enumerate(raw_data.get('vibration', [])):
 275             data['vibration'][i],changed = check_vibration(
 276                     index=i, vibration=vibration, config=config,
 277                     deflection_channel_config=input_config, plot=plot,
 278                     maximum_relative_error=maximum_relative_error)
 279             if changed and not dry_run:
 280                 vibrations_changed = True
 281                 vibration_group = _h5_create_group(
 282                     group, 'vibration/{}'.format(i))
 283                 _vibration_save(
 284                     group=vibration_group, processed=data['vibration'][i])
 285         if (bumps_changed or temperatures_changed or vibrations_changed
 286             ) and not dry_run:
 287             calibration_group = _h5_create_group(group, 'calibration')
 288             if bumps_changed:
 289                 save_results(
 290                     group=calibration_group, bump=data['bump'])
 291             if temperatures_changed:
 292                 save_results(
 293                     group=calibration_group, temperature=data['temperature'])
 294             if vibrations_changed:
 295                 save_results(
 296                     group=calibration_group, vibration=data['vibration'])
 297         if len(raw_data.get('bump', [])) != len(data['bump']):
 298             raise ValueError(
 299                 'not enough raw bump data: {} of {}'.format(
 300                     len(raw_data.get('bump', [])), len(data['bump'])))
 301         if len(raw_data.get('temperature', [])) != len(data['temperature']):
 302             raise ValueError(
 303                 'not enough raw temperature data: {} of {}'.format(
 304                     len(raw_data.get('temperature', [])),
 305                     len(data['temperature'])))
 306         if len(raw_data['vibration']) != len(data['vibration']):
 307             raise ValueError(
 308                 'not enough raw vibration data: {} of {}'.format(
 309                     len(raw_data.get('vibration', [])),
 310                     len(data['vibration'])))
 311         k,k_s,changed = check_calibration(
 312             k=data.get('processed', {}).get('spring_constant', None),
 313             k_s=data.get('processed', {}).get(
 314                 'spring_constant_deviation', None),
 315             bumps=data['bump'],
 316             temperatures=data['temperature'], vibrations=data['vibration'],
 317             maximum_relative_error=maximum_relative_error)
 318         if changed and not dry_run:
 319             if calibration_group is None:
 320                 calibration_group = _h5_create_group(group, 'calibration')
 321             save_results(
 322                 group=calibration_group,
 323                 spring_constant=k, spring_constant_deviation=k_s)
 324     finally:
 325         if f:
 326             f.close()
 327     if plot:
 328         _plot(bumps=data['bump'],
 329               temperatures=data['temperature'],
 330               vibrations=data['vibration'])
 331     return (k, k_s)
 332
 333 def check_bump(index, bump, config=None, maximum_relative_error=0, **kwargs):
 334     changed = False
 335     try:
 336         bump_config = bump['config']['bump']
 337     except KeyError:
 338         bump_config = config['bump']
 339     sensitivity = _bump_analyze(
 340         config=bump_config, data=bump['raw'], **kwargs)
 341     if bump.get('processed', None) is None:
 342         changed = True
 343         _LOG.warn('new analysis for bump {}: {}'.format(index, sensitivity))
 344     else:
 345         rel_error = abs(sensitivity - bump['processed'])/bump['processed']
 346         if rel_error > maximum_relative_error:
 347             changed = True
 348             _LOG.warn(("new analysis doesn't match for bump {}: {} -> {} "
 349                        "(difference: {}, relative error: {})").format(
 350                     index, bump['processed'], sensitivity,
 351                     sensitivity-bump['processed'], rel_error))
 352     return (sensitivity, changed)
 353
 354 def check_temperature(index, temperature, config=None,
 355                       maximum_relative_error=0, **kwargs):
 356     changed = False
 357     try:
 358         temp_config = temperature['config']['temperature']
 359     except KeyError:
 360         temp_config = config['temperature']
 361     temp = _temperature_analyze(
 362         config=temp_config,
 363         temperature=temperature['raw'], **kwargs)
 364     if temperature.get('processed', None) is None:
 365         changed = True
 366         _LOG.warn('new analysis for temperature {}: {}'.format(index, temp))
 367     else:
 368         rel_error = abs(temp - temperature['processed']
 369                         )/temperature['processed']
 370         if rel_error > maximum_relative_error:
 371             changed = True
 372             _LOG.warn(("new analysis doesn't match for temperature "
 373                        "{} -> {} (difference: {}, relative error: {})"
 374                        ).format(
 375                     index, temperature['processed'], temp,
 376                     temp-temperature['processed'], rel_error))
 377     return (temp, changed)
 378
 379 def check_vibration(index, vibration, config=None, maximum_relative_error=0,
 380                     **kwargs):
 381     changed = False
 382     try:
 383         vib_config = vibration['config']['vibration']
 384     except KeyError:
 385         vib_config = config['vibration']
 386     variance = _vibration_analyze(
 387         config=vib_config, deflection=vibration['raw'], **kwargs)
 388     if vibration.get('processed', None) is None:
 389         changed = True
 390         _LOG.warn('new analysis for vibration {}: {}'.format(
 391                 index, variance))
 392     else:
 393         rel_error = abs(variance-vibration['processed'])/vibration['processed']
 394         if rel_error > maximum_relative_error:
 395             changed = True
 396             _LOG.warn(("new analysis doesn't match for vibration {}: {} != {} "
 397                        "(difference: {}, relative error: {})").format(
 398                     index, variance, vibration['processed'],
 399                     variance-vibration['processed'], rel_error))
 400     return (variance, changed)
 401
 402 def check_calibration(k, k_s, maximum_relative_error, **kwargs):
 403     changed = False
 404     new_k,new_k_s = analyze(**kwargs)
 405     if k is None:
 406         changed = True
 407         _LOG.warn('new analysis for the spring constant: {}'.format(new_k))
 408     else:
 409         rel_error = abs(new_k-k)/k
 410         if rel_error > maximum_relative_error:
 411             changed = True
 412             _LOG.warn(("new analysis doesn't match for the spring constant: "
 413                        "{} != {} (difference: {}, relative error: {})").format(
 414                     new_k, k, new_k-k, rel_error))
 415     if k_s is None:
 416         changed = True
 417         _LOG.warn('new analysis for the spring constant deviation: {}'.format(
 418                 new_k_s))
 419     else:
 420         rel_error = abs(new_k_s-k_s)/k_s
 421         if rel_error > maximum_relative_error:
 422             changed = True
 423             _LOG.warn(
 424                 ("new analysis doesn't match for the spring constant deviation"
 425                  ": {} != {} (difference: {}, relative error: {})").format(
 426                     new_k_s, k_s, new_k_s-k_s, rel_error))
 427     return (new_k, new_k_s, changed)