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 `calib_analyze()` from the other `calib_*()`
  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 os
  43 >>> import tempfile
  44 >>> import numpy
  45 >>> from h5config.storage.hdf5 import HDF5_Storage, pprint_HDF5
  46 >>> from .config import CalibrationConfig
  47
  48 >>> fd,filename = tempfile.mkstemp(suffix='.h5', prefix='calibcant-')
  49 >>> os.close(fd)
  50
  51 >>> calibration_config = CalibrationConfig(storage=HDF5_Storage(
  52 ...         filename=filename, group='/calib/config/'))
  53 >>> bumps = numpy.array((15.9e6, 16.9e6, 16.3e6))
  54 >>> temperatures = numpy.array((295, 295.2, 294.8))
  55 >>> vibrations = numpy.array((2.20e-5, 2.22e-5, 2.21e-5))
  56
  57 >>> k,k_s = calib_analyze(bumps=bumps, temperatures=temperatures,
  58 ...     vibrations=vibrations)
  59 >>> (k, k_s)  # doctest: +ELLIPSIS
  60 (0.0493..., 0.00248...)
  61
  62 Most of the error in this example comes from uncertainty in the
  63 photodiode sensitivity (bumps).
  64
  65 >>> k_s/k  # doctest: +ELLIPSIS
  66 0.0503...
  67 >>> bumps.std()/bumps.mean()  # doctest: +ELLIPSIS
  68 0.0251...
  69 >>> temperatures.std()/temperatures.mean()  # doctest: +ELLIPSIS
  70 0.000553...
  71 >>> vibrations.std()/vibrations.mean()  # doctest: +ELLIPSIS
  72 0.00369...
  73
  74 >>> calib_save(filename=filename, group='/calib/',
  75 ...     bumps=bumps, temperatures=temperatures, vibrations=vibrations,
  76 ...     calibration_config=calibration_config, k=k, k_s=k_s)
  77 >>> pprint_HDF5(filename)  # doctest: +ELLIPSIS, +REPORT_UDIFF
  78 /
  79   /calib
  80     /calib/config
  81       <HDF5 dataset "bump": shape (), type "|S1">
  82 <BLANKLINE>
  83       <HDF5 dataset "num-bumps": shape (), type "<i4">
  84         10
  85       <HDF5 dataset "num-temperatures": shape (), type "<i4">
  86         10
  87       <HDF5 dataset "num-vibrations": shape (), type "<i4">
  88         20
  89       <HDF5 dataset "temperature": shape (), type "|S1">
  90 <BLANKLINE>
  91       <HDF5 dataset "temperature-sleep": shape (), type "<i4">
  92         1
  93       <HDF5 dataset "vibration": shape (), type "|S1">
  94 <BLANKLINE>
  95       <HDF5 dataset "vibration-spacing": shape (), type "<f8">
  96         5e-05
  97     /calib/processed
  98       /calib/processed/spring-constant
  99         <HDF5 dataset "data": shape (), type "<f8">
 100           0.0493...
 101         <HDF5 dataset "standard-deviation": shape (), type "<f8">
 102           0.00248...
 103         <HDF5 dataset "units": shape (), type "|S3">
 104           N/m
 105     /calib/raw
 106       /calib/raw/photodiode-sensitivity
 107         <HDF5 dataset "data": shape (3,), type "<f8">
 108           [ 15900000.  16900000.  16300000.]
 109         <HDF5 dataset "units": shape (), type "|S3">
 110           V/m
 111       /calib/raw/temperature
 112         <HDF5 dataset "data": shape (3,), type "<f8">
 113           [ 295.   295.2  294.8]
 114         <HDF5 dataset "units": shape (), type "|S1">
 115           K
 116       /calib/raw/thermal-vibration-variance
 117         <HDF5 dataset "data": shape (3,), type "<f8">
 118           [  2.20...e-05   2.220...e-05   2.210...e-05]
 119         <HDF5 dataset "units": shape (), type "|S3">
 120           V^2
 121
 122 >>> bumps,temperatures,vibrations,calibration_config,k,k_s = calib_load(
 123 ...     filename=filename, group='/calib/')
 124 >>> (k, k_s)  # doctest: +ELLIPSIS
 125 (0.0493..., 0.00248...)
 126
 127 >>> os.remove(filename)
 128 """
 129
 130 import h5py as _h5py
 131 import numpy as _numpy
 132 try:
 133     from scipy.constants import Boltzmann as _kB  # in J/K
 134 except ImportError:
 135     from scipy.constants import Bolzmann as _kB  # in J/K
 136 # Bolzmann -> Boltzmann patch submitted:
 137 #   http://projects.scipy.org/scipy/ticket/1417
 138 # Fixed in scipy commit 4716d91, Apr 2, 2011, during work after v0.9.0rc5.
 139
 140 try:
 141     import matplotlib as _matplotlib
 142     import matplotlib.pyplot as _matplotlib_pyplot
 143     import time as _time  # for timestamping lines on plots
 144 except (ImportError, RuntimeError), e:
 145     _matplotlib = None
 146     _matplotlib_import_error = e
 147
 148 from h5config.storage.hdf5 import HDF5_Storage as _HDF5_Storage
 149 from h5config.storage.hdf5 import h5_create_group as _h5_create_group
 150
 151 from . import LOG as _LOG
 152 from . import package_config as _package_config
 153 from .bump_analyze import bump_analyze as _bump_analyze
 154 from .bump_analyze import bump_load as _bump_load
 155 from .bump_analyze import bump_save as _bump_save
 156 from .config import CalibrationConfig as _CalibrationConfig
 157 from .T_analyze import T_analyze as _temperature_analyze
 158 from .T_analyze import T_load as _temperature_load
 159 from .T_analyze import T_save as _temperature_save
 160 from .vib_analyze import vib_analyze as _vibration_analyze
 161 from .vib_analyze import vib_load as _vibration_load
 162 from .vib_analyze import vib_save as _vibration_save
 163
 164
 165 def calib_analyze(bumps, temperatures, vibrations):
 166     """Analyze data from `get_calibration_data()`
 167
 168     Inputs (all are arrays of recorded data):
 169       bumps measured (V_photodiode / nm_tip) proportionality constant
 170       temperatures    measured temperature (K)
 171       vibrations  measured V_photodiode variance in free solution (V**2)
 172     Outputs:
 173       k    cantilever spring constant (in N/m, or equivalently nN/nm)
 174       k_s  standard deviation in our estimate of k
 175
 176     Notes:
 177
 178     We're assuming vib is mostly from thermal cantilever vibrations
 179     (and then only from vibrations in the single vertical degree of
 180     freedom), and not from other noise sources.
 181
 182     If the error is large, check the relative errors
 183     (`x.std()/x.mean()`)of your input arrays.  If one of them is
 184     small, don't bother repeating that measurment too often.  If one
 185     is large, try repeating that measurement more.  Remember that you
 186     need enough samples to have a valid error estimate in the first
 187     place, and that none of this addresses any systematic errors.
 188     """
 189     ps_m = bumps.mean()  # ps for photo-sensitivity
 190     ps_s = bumps.std()
 191     T_m = temperatures.mean()
 192     T_s = temperatures.std()
 193     v2_m = vibrations.mean()  # average voltage variance
 194     v2_s = vibrations.std()
 195
 196     # Vphoto / photo_sensitivity = x
 197     # k = kB T / <x**2> = kB T photo_sensitivity**2 / Vphoto_var
 198     #
 199     # units,  photo_sensitivity =  Vphoto/(Zcant in m),
 200     # so Vphoto/photo_sensitivity = Zcant in m
 201     # so k = J/K * K / m^2 = J / m^2 = N/m
 202     k  = _kB * T_m * ps_m**2 / v2_m
 203
 204     # propogation of errors
 205     # dk/dT = k/T
 206     dk_T = k/T_m * T_s
 207     # dk/dps = 2k/ps
 208     dk_ps = 2*k/ps_m * ps_s
 209     # dk/dv2 = -k/v2
 210     dk_v = -k/v2_m * v2_s
 211
 212     k_s = _numpy.sqrt(dk_T**2 + dk_ps**2 + dk_v**2)
 213
 214     _LOG.info('variable (units)         : '
 215               'mean +/- std. dev. (relative error)')
 216     _LOG.info('cantilever k (N/m)       : %g +/- %g (%g)' % (k, k_s, k_s/k))
 217     _LOG.info('photo sensitivity (V/m)  : %g +/- %g (%g)'
 218               % (ps_m, ps_s, ps_s/ps_m))
 219     _LOG.info('T (K)                    : %g +/- %g (%g)'
 220               % (T_m, T_s, T_s/T_m))
 221     _LOG.info('vibration variance (V^2) : %g +/- %g (%g)'
 222               % (v2_m, v2_s, v2_s/v2_m))
 223
 224     if _package_config['matplotlib']:
 225         calib_plot(bumps, temperatures, vibrations)
 226
 227     return (k, k_s)
 228
 229 def calib_save(filename, group='/', bumps=None, temperatures=None,
 230                vibrations=None, calibration_config=None, k=None, k_s=None):
 231     with _h5py.File(filename, 'a') as f:
 232         cwg = _h5_create_group(f, group)
 233         if calibration_config is not None:
 234             config_cwg = _h5_create_group(cwg, 'config')
 235             storage = _HDF5_Storage()
 236             storage.save(config=calibration_config, group=config_cwg)
 237         if bumps is not None:
 238             try:
 239                 del cwg['raw/photodiode-sensitivity/data']
 240             except KeyError:
 241                 pass
 242             try:
 243                 del cwg['raw/photodiode-sensitivity/units']
 244             except KeyError:
 245                 pass
 246             cwg['raw/photodiode-sensitivity/data'] = bumps
 247             cwg['raw/photodiode-sensitivity/units'] = 'V/m'
 248         if temperatures is not None:
 249             try:
 250                 del cwg['raw/temperature/data']
 251             except KeyError:
 252                 pass
 253             try:
 254                 del cwg['raw/temperature/units']
 255             except KeyError:
 256                 pass
 257             cwg['raw/temperature/data'] = temperatures
 258             cwg['raw/temperature/units'] = 'K'
 259         if vibrations is not None:
 260             try:
 261                 del cwg['raw/thermal-vibration-variance/data']
 262             except KeyError:
 263                 pass
 264             try:
 265                 del cwg['raw/thermal-vibration-variance/units']
 266             except KeyError:
 267                 pass
 268             cwg['raw/thermal-vibration-variance/data'] = vibrations
 269             cwg['raw/thermal-vibration-variance/units'] = 'V^2'
 270         if k is not None:
 271             try:
 272                 del cwg['processed/spring-constant/data']
 273             except KeyError:
 274                 pass
 275             try:
 276                 del cwg['processed/spring-constant/units']
 277             except KeyError:
 278                 pass
 279             cwg['processed/spring-constant/data'] = k
 280             cwg['processed/spring-constant/units'] = 'N/m'
 281         if k_s is not None:
 282             try:
 283                 del cwg['processed/spring-constant/standard-deviation']
 284             except KeyError:
 285                 pass
 286             cwg['processed/spring-constant/standard-deviation'] = k_s
 287
 288 def calib_load(filename, group='/'):
 289     assert group.endswith('/')
 290     bumps = temperatures = vibrations = k = k_s = None
 291     configs = []
 292     with _h5py.File(filename, 'a') as f:
 293         try:
 294             bumps = f[group+'raw/photodiode-sensitivity/data'][...]
 295         except KeyError:
 296             pass
 297         try:
 298             temperatures = f[group+'raw/temperature/data'][...]
 299         except KeyError:
 300             pass
 301         try:
 302             vibrations = f[group+'raw/thermal-vibration-variance/data'][...]
 303         except KeyError:
 304             pass
 305         try:
 306             k = float(f[group+'processed/spring-constant/data'][...])
 307         except KeyError:
 308             pass
 309         try:
 310             k_s = float(
 311                 f[group+'processed/spring-constant/standard-deviation'][...])
 312         except KeyError:
 313             pass
 314     calibration_config = _CalibrationConfig(storage=_HDF5_Storage(
 315             filename=filename, group=group+'config/'))
 316     calibration_config.load()
 317     return (bumps, temperatures, vibrations, calibration_config, k, k_s)
 318
 319 def calib_plot(bumps, temperatures, vibrations):
 320     if not _matplotlib:
 321         raise _matplotlib_import_error
 322     figure = _matplotlib_pyplot.figure()
 323
 324     bump_axes = figure.add_subplot(3, 1, 1)
 325     T_axes = figure.add_subplot(3, 1, 2)
 326     vib_axes = figure.add_subplot(3, 1, 3)
 327
 328     timestamp = _time.strftime('%H%M%S')
 329     bump_axes.set_title('cantilever calibration %s' % timestamp)
 330
 331     bump_axes.plot(bumps, 'g.-')
 332     bump_axes.set_ylabel('photodiode sensitivity (V/m)')
 333     T_axes.plot(temperatures, 'r.-')
 334     T_axes.set_ylabel('temperature (K)')
 335     vib_axes.plot(vibrations, 'b.-')
 336     vib_axes.set_ylabel('thermal deflection variance (V^2)')
 337
 338     if hasattr(figure, 'show'):
 339         figure.show()
 340
 341
 342 def calib_load_all(filename, group='/'):
 343     "Load all data from a `calib()` run."
 344     assert group.endswith('/'), group
 345     bumps,temperatures,vibrations,calibration_config,k,k_s = calib_load(
 346         filename, group+'calibration/')
 347     bump_details = []
 348     for i in range(calibration_config['num-bumps']):
 349         (raw_bump,bump_config,z_axis_config,deflection_channel_config,
 350          processed_bump) = _bump_load(
 351             filename=filename, group='%sbump/%d/' % (group, i))
 352         bump_details.append({
 353                 'raw_bump': raw_bump,
 354                 'bump_config': bump_config,
 355                 'z_axis_config': z_axis_config,
 356                 'deflection_channel_config': deflection_channel_config,
 357                 'processed_bump': processed_bump,
 358                 })
 359     temperature_details = []
 360     for i in range(calibration_config['num-temperatures']):
 361         (raw_temperature,temperature_config,processed_temperature
 362          ) = _temperature_load(
 363             filename=filename, group='%stemperature/%d/' % (group, i))
 364         temperature_details.append({
 365                 'raw_temperature': raw_temperature,
 366                 'temperature_config': temperature_config,
 367                 'processed_temperature': processed_temperature,
 368                 })
 369     vibration_details = []
 370     for i in range(calibration_config['num-vibrations']):
 371         (raw_vibration,vibration_config,deflection_channel_config,
 372          processed_vibration) = _vibration_load(
 373             filename=filename, group='%svibration/%d/' % (group, i))
 374         vibration_details.append({
 375                 'raw_vibration': raw_vibration,
 376                 'vibration_config': vibration_config,
 377                 'deflection_channel_config': deflection_channel_config,
 378                 'processed_vibration': processed_vibration,
 379                 })
 380     return {
 381         'bumps': bumps,
 382         'bump_details': bump_details,
 383         'temperatures': temperatures,
 384         'temperature_details': temperature_details,
 385         'vibrations': vibrations,
 386         'vibration_details': vibration_details,
 387         'calibration_config': calibration_config,
 388         'k': k,
 389         'k_s': k_s,
 390         }
 391
 392 def calib_analyze_all(filename, group='/', maximum_relative_error=1e-5,
 393                       dry_run=False):
 394     "(Re)analyze all data from a `calib()` run."
 395     assert group.endswith('/'), group
 396     bumps,temperatures,vibrations,calibration_config,k,k_s = calib_load(
 397         filename, group+'calibration/')
 398     if bumps is None:
 399         bumps = _numpy.zeros(
 400             (calibration_config['num-bumps'],), dtype=float)
 401     if temperatures is None:
 402         temperatures = _numpy.zeros(
 403             (calibration_config['num-temperatures'],), dtype=float)
 404     if vibrations is None:
 405         vibrations = _numpy.zeros(
 406             (calibration_config['num-vibrations'],), dtype=float)
 407     changed_bump = changed_temperature = changed_vibration = False
 408     for i in range(calibration_config['num-bumps']):
 409         bump_group = '%sbump/%d/' % (group, i)
 410         (raw_bump,bump_config,z_axis_config,
 411          deflection_channel_config,processed_bump) = _bump_load(
 412             filename=filename, group=bump_group)
 413         sensitivity = _bump_analyze(
 414             data=raw_bump, bump_config=bump_config,
 415             z_axis_config=z_axis_config,
 416             deflection_channel_config=deflection_channel_config)
 417         bumps[i] = sensitivity
 418         rel_error = abs(sensitivity - processed_bump)/processed_bump
 419         if rel_error > maximum_relative_error:
 420             _LOG.warn(("new analysis doesn't match for bump %d: %g -> %g "
 421                        "(difference: %g, relative error: %g)")
 422                       % (i, processed_bump, sensitivity,
 423                          sensitivity-processed_bump, rel_error))
 424             changed_bump = True
 425             if not dry_run:
 426                 _bump_save(filename, bump_group, processed_bump=sensitivity)
 427     for i in range(calibration_config['num-temperatures']):
 428         temperature_group = '%stemperature/%d/' % (group, i)
 429         (raw_temperature,temperature_config,processed_temperature
 430          ) = _temperature_load(
 431             filename=filename, group=temperature_group)
 432         temperature = _temperature_analyze(
 433             raw_temperature, temperature_config)
 434         temperatures[i] = temperature
 435         rel_error = abs(temperature - processed_temperature
 436                         )/processed_temperature
 437         if rel_error > maximum_relative_error:
 438             _LOG.warn(("new analysis doesn't match for temperature %d: "
 439                        "%g -> %g (difference: %g, relative error: %g)")
 440                       % (i, processed_temperature, temperature,
 441                          temperature-processed_temperature, rel_error))
 442             changed_temperature = True
 443             if not dry_run:
 444                 _temperature_save(
 445                     filename, temperature_group,
 446                     processed_T=temperature)
 447     for i in range(calibration_config['num-vibrations']):
 448         vibration_group = '%svibration/%d/' % (group, i)
 449         (raw_vibration,vibration_config,deflection_channel_config,
 450          processed_vibration) = _vibration_load(
 451             filename=filename, group=vibration_group)
 452         variance = _vibration_analyze(
 453             deflection=raw_vibration, vibration_config=vibration_config,
 454             deflection_channel_config=deflection_channel_config)
 455         vibrations[i] = variance
 456         rel_error = abs(variance - processed_vibration)/processed_vibration
 457         if rel_error > maximum_relative_error:
 458             _LOG.warn(("new analysis doesn't match for vibration %d: %g -> %g "
 459                        "(difference: %g, relative error: %g)")
 460                       % (i, processed_vibration, variance,
 461                          variance-processed_vibration, rel_error))
 462             changed_vibration = True
 463             if not dry_run:
 464                 _vibration_save(
 465                     filename, vibration_group, processed_vibration=variance)
 466
 467     calib_group = '%scalibration/' % group
 468
 469     if changed_bump and not dry_run:
 470         calib_save(filename, calib_group, bumps=bumps)
 471     if changed_temperature and not dry_run:
 472         calib_save(filename, calib_group, temperatures=temperatures)
 473     if changed_vibration and not dry_run:
 474         calib_save(filename, calib_group, vibrations=vibrations)
 475
 476     new_k,new_k_s = calib_analyze(
 477         bumps=bumps, temperatures=temperatures, vibrations=vibrations)
 478     rel_error = abs(new_k-k)/k
 479     if rel_error > maximum_relative_error:
 480         _LOG.warn(("new analysis doesn't match for k: %g -> %g "
 481                    "(difference: %g, relative error: %g)")
 482                   % (k, new_k, new_k-k, rel_error))
 483         if not dry_run:
 484             calib_save(filename, calib_group, k=new_k)
 485     rel_error = abs(new_k_s-k_s)/k_s
 486     if rel_error > maximum_relative_error:
 487         _LOG.warn(("new analysis doesn't match for k_s: %g -> %g "
 488                    "(difference: %g, relative error: %g)")
 489                   % (k_s, new_k_s, new_k_s-k_s, rel_error))
 490         if not dry_run:
 491             calib_save(filename, calib_group, k_s=new_k_s)
 492     return (new_k, new_k_s)
 493
 494 def calib_plot_all(bumps, bump_details, temperatures, temperature_details,
 495                    vibrations, vibration_details, calibration_config, k, k_s,
 496                    maximum_relative_error=1e-5):
 497     calib_plot(bumps, temperatures, vibrations)
 498     for i,bump in enumerate(bump_details):
 499         sensitivity = _bump_analyze(
 500             data=bump['raw_bump'], bump_config=bump['bump_config'],
 501             z_axis_config=bump['z_axis_config'],
 502             deflection_channel_config=bump['deflection_channel_config'],
 503             plot=True)
 504         rel_error = abs(sensitivity - bump['processed_bump']
 505                         )/bump['processed_bump']
 506         if rel_error > maximum_relative_error:
 507             _LOG.warn(("new analysis doesn't match for bump %d: %g != %g "
 508                        "(difference: %g, relative error: %g)")
 509                       % (i, sensitivity, bump['processed_bump'],
 510                          sensitivity-bump['processed_bump'], rel_error))
 511     # nothing interesting to plot for temperatures...
 512     for i,vibration in enumerate(vibration_details):
 513         variance = _vibration_analyze(
 514             deflection=vibration['raw_vibration'],
 515             vibration_config=vibration['vibration_config'],
 516             deflection_channel_config=vibration['deflection_channel_config'],
 517             plot=True)
 518         rel_error = abs(variance - vibration['processed_vibration']
 519                         )/vibration['processed_vibration']
 520         if rel_error > maximum_relative_error:
 521             _LOG.warn(("new analysis doesn't match for vibration %d: %g != %g "
 522                        "(difference: %g, relative error: %g)")
 523                       % (i, variance, vibration['processed_vibration'],
 524                          variance-vibration['processed_vibration'], rel_error))