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     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     axis_config = config['afm']['piezo'].select_config(
 228         setting_name='axes',
 229         attribute_value=config['afm']['main-axis'],
 230         get_attribute=_get_axis_name)
 231     input_config = config['afm']['piezo'].select_config(
 232         setting_name='inputs', attribute_value='deflection')
 233     bumps_changed = temperatures_changed = vibrations_changed = False
 234     calibration_group = None
 235     if not isinstance(group, _h5py.Group) and not dry_run:
 236         f = _h5py.File(filename, mode='a')
 237         group = _h5_create_group(f, group)
 238     else:
 239         f = None
 240     try:
 241         if len(data.get('raw', {}).get('bump', [])) != len(data['bump']):
 242             bumps_changed = True
 243         for i,bump in enumerate(raw_data['bump']):
 244             data['bump'][i],changed = check_bump(
 245                 index=i, bump=bump, config=config, z_axis_config=axis_config,
 246                 deflection_channel_config=input_config, plot=plot,
 247                 maximum_relative_error=maximum_relative_error)
 248             if changed and not dry_run:
 249                 bumps_changed = True
 250                 bump_group = _h5_create_group(group, 'bump/{}'.format(i))
 251                 _bump_save(group=bump_group, processed=data['bump'][i])
 252         if len(data.get('raw', {}).get('temperature', [])
 253                ) != len(data['temperature']):
 254             temperatures_changed = True
 255         for i,temperature in enumerate(raw_data['temperature']):
 256             data['temperature'][i],changed = check_temperature(
 257                 index=i, temperature=temperature, config=config,
 258                 maximum_relative_error=maximum_relative_error)
 259             if changed and not dry_run:
 260                 temperatures_changed = True
 261                 temperature_group = _h5_create_group(
 262                     group, 'temperature/{}'.format(i))
 263                 _temperature_save(
 264                     group=temperature_group, processed=data['temperature'][i])
 265         if len(data.get('raw', {}).get('vibration', [])
 266                ) != len(data['vibration']):
 267             vibrations_changed = True
 268         for i,vibration in enumerate(raw_data['vibration']):
 269             data['vibration'][i],changed = check_vibration(
 270                     index=i, vibration=vibration, config=config,
 271                     deflection_channel_config=input_config, plot=plot,
 272                     maximum_relative_error=maximum_relative_error)
 273             if changed and not dry_run:
 274                 vibrations_changed = True
 275                 vibration_group = _h5_create_group(
 276                     group, 'vibration/{}'.format(i))
 277                 _vibration_save(
 278                     group=vibration_group, processed=data['vibration'])
 279         if (bumps_changed or temperatures_changed or vibrations_changed
 280             ) and not dry_run:
 281             calibration_group = _h5_create_group(group, 'calibration')
 282             if bumps_changed:
 283                 save_results(
 284                     group=calibration_group, bump=data['bump'])
 285             if temperatures_changed:
 286                 save_results(
 287                     group=calibration_group, temperature=data['temperature'])
 288             if vibrations_changed:
 289                 save_results(
 290                     group=calibration_group, vibration=data['vibration'])
 291         if len(raw_data['bump']) != len(data['bump']):
 292             raise ValueError(
 293                 'not enough raw bump data: {} of {}'.format(
 294                     len(raw_data['bump']), len(data['bump'])))
 295         if len(raw_data['temperature']) != len(data['temperature']):
 296             raise ValueError(
 297                 'not enough raw temperature data: {} of {}'.format(
 298                     len(raw_data['temperature']), len(data['temperature'])))
 299         if len(raw_data['vibration']) != len(data['vibration']):
 300             raise ValueError(
 301                 'not enough raw vibration data: {} of {}'.format(
 302                     len(raw_data['vibration']), len(data['vibration'])))
 303         k,k_s,changed = check_calibration(
 304             k=data.get('processed', {}).get('spring_constant', None),
 305             k_s=data.get('processed', {}).get(
 306                 'spring_constant_deviation', None),
 307             bumps=data['bump'],
 308             temperatures=data['temperature'], vibrations=data['vibration'],
 309             maximum_relative_error=maximum_relative_error)
 310         if changed and not dry_run:
 311             if calibration_group is None:
 312                 calibration_group = _h5_create_group(group, 'calibration')
 313             save_results(
 314                 group=calibration_group,
 315                 spring_constant=k, spring_constant_deviation=k_s)
 316     finally:
 317         if f:
 318             f.close()
 319     if plot:
 320         _plot(bumps=data['raw']['bump'],
 321              temperatures=data['raw']['temperature'],
 322              vibrations=data['raw']['vibration'])
 323     return (k, k_s)
 324
 325 def check_bump(index, bump, config=None, maximum_relative_error=0, **kwargs):
 326     changed = False
 327     try:
 328         bump_config = bump['config']['bump']
 329     except KeyError:
 330         bump_config = config['bump']
 331     sensitivity = _bump_analyze(
 332         config=bump_config, data=bump['raw'], **kwargs)
 333     if bump.get('processed', None) is None:
 334         changed = True
 335         _LOG.warn('new analysis for bump {}: {}'.format(index, sensitivity))
 336     else:
 337         rel_error = abs(sensitivity - bump['processed'])/bump['processed']
 338         if rel_error > maximum_relative_error:
 339             changed = True
 340             _LOG.warn(("new analysis doesn't match for bump {}: {} -> {} "
 341                        "(difference: {}, relative error: {})").format(
 342                     index, bump['processed'], sensitivity,
 343                     sensitivity-bump['processed'], rel_error))
 344     return (sensitivity, changed)
 345
 346 def check_temperature(index, temperature, config=None,
 347                       maximum_relative_error=0, **kwargs):
 348     changed = False
 349     try:
 350         temp_config = temperature['config']['temperature']
 351     except KeyError:
 352         temp_config = config['temperature']
 353     temp = _temperature_analyze(
 354         config=temp_config,
 355         temperature=temperature['raw'], **kwargs)
 356     if temperature.get('processed', None) is None:
 357         changed = True
 358         _LOG.warn('new analysis for temperature {}: {}'.format(index, temp))
 359     else:
 360         rel_error = abs(temp - temperature['processed']
 361                         )/temperature['processed']
 362         if rel_error > maximum_relative_error:
 363             changed = True
 364             _LOG.warn(("new analysis doesn't match for temperature "
 365                        "{} -> {} (difference: {}, relative error: {})"
 366                        ).format(
 367                     index, temperature['processed'], temp,
 368                     temp-temperature['processed'], rel_error))
 369     return (temp, changed)
 370
 371 def check_vibration(index, vibration, config=None, maximum_relative_error=0,
 372                     **kwargs):
 373     changed = False
 374     try:
 375         vib_config = vibration['config']['vibration']
 376     except KeyError:
 377         vib_config = config['vibration']
 378     variance = _vibration_analyze(
 379         config=vib_config, deflection=vibration['raw'], **kwargs)
 380     if vibration.get('processed', None) is None:
 381         changed = True
 382         _LOG.warn('new analysis for vibration {}: {}'.format(
 383                 index, variance))
 384     else:
 385         rel_error = abs(variance-vibration['processed'])/vibration['processed']
 386         if rel_error > maximum_relative_error:
 387             _LOG.warn(("new analysis doesn't match for vibration {}: {} != {} "
 388                        "(difference: {}, relative error: {})").format(
 389                     index, variance, vibration['processed'],
 390                     variance-vibration['processed'], rel_error))
 391     return (variance, changed)
 392
 393 def check_calibration(k, k_s, maximum_relative_error, **kwargs):
 394     changed = False
 395     new_k,new_k_s = analyze(**kwargs)
 396     if k is None:
 397         changed = True
 398         _LOG.warn('new analysis for the spring constant: {}'.format(new_k))
 399     else:
 400         rel_error = abs(new_k-k)/k
 401         if rel_error > maximum_relative_error:
 402             changed = True
 403             _LOG.warn(("new analysis doesn't match for the spring constant: "
 404                        "{} != {} (difference: {}, relative error: {})").format(
 405                     new_k, k, new_k-k, rel_error))
 406     if k_s is None:
 407         changed = True
 408         _LOG.warn('new analysis for the spring constant deviation: {}'.format(
 409                 new_k_s))
 410     else:
 411         rel_error = abs(new_k_s-k_s)/k_s
 412         if rel_error > maximum_relative_error:
 413             changed = True
 414             _LOG.warn(
 415                 ("new analysis doesn't match for the spring constant deviation"
 416                  ": {} != {} (difference: {}, relative error: {})").format(
 417                     new_k_s, k_s, new_k_s-k_s, rel_error))
 418     return (new_k, new_k_s, changed)