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         _changed_bump = False
 410         bump_group = '%sbump/%d/' % (group, i)
 411         (raw_bump,bump_config,z_axis_config,
 412          deflection_channel_config,processed_bump) = _bump_load(
 413             filename=filename, group=bump_group)
 414         sensitivity = _bump_analyze(
 415             data=raw_bump, bump_config=bump_config,
 416             z_axis_config=z_axis_config,
 417             deflection_channel_config=deflection_channel_config)
 418         bumps[i] = sensitivity
 419         if processed_bump is None:
 420             _changed_bump = True
 421             _LOG.warn('new analysis for bump %d: %g' % (i, sensitivity))
 422         else:
 423             rel_error = abs(sensitivity - processed_bump)/processed_bump
 424             if rel_error > maximum_relative_error:
 425                 _changed_bump = True
 426                 _LOG.warn(("new analysis doesn't match for bump %d: %g -> %g "
 427                            "(difference: %g, relative error: %g)")
 428                           % (i, processed_bump, sensitivity,
 429                              sensitivity-processed_bump, rel_error))
 430         if _changed_bump and not dry_run:
 431             changed_bump = True
 432             _bump_save(filename, bump_group, processed_bump=sensitivity)
 433     for i in range(calibration_config['num-temperatures']):
 434         _changed_temperature = False
 435         temperature_group = '%stemperature/%d/' % (group, i)
 436         (raw_temperature,temperature_config,processed_temperature
 437          ) = _temperature_load(
 438             filename=filename, group=temperature_group)
 439         temperature = _temperature_analyze(
 440             raw_temperature, temperature_config)
 441         temperatures[i] = temperature
 442         if processed_temperature is None:
 443             _changed_temperature = True
 444             _LOG.warn('new analysis for temperature %d: %g' % (i, temperature))
 445         else:
 446             rel_error = abs(temperature - processed_temperature
 447                             )/processed_temperature
 448             if rel_error > maximum_relative_error:
 449                 _changed_temperature = True
 450                 _LOG.warn(("new analysis doesn't match for temperature %d: "
 451                            "%g -> %g (difference: %g, relative error: %g)")
 452                           % (i, processed_temperature, temperature,
 453                              temperature-processed_temperature, rel_error))
 454         if _changed_temperature and not dry_run:
 455             changed_temperature = True
 456             _temperature_save(
 457                 filename, temperature_group,
 458                 processed_T=temperature)
 459     for i in range(calibration_config['num-vibrations']):
 460         _changed_vibration = False
 461         vibration_group = '%svibration/%d/' % (group, i)
 462         (raw_vibration,vibration_config,deflection_channel_config,
 463          processed_vibration) = _vibration_load(
 464             filename=filename, group=vibration_group)
 465         variance = _vibration_analyze(
 466             deflection=raw_vibration, vibration_config=vibration_config,
 467             deflection_channel_config=deflection_channel_config)
 468         vibrations[i] = variance
 469         if processed_vibration is None:
 470             _changed_vibration = True
 471             _LOG.warn('new analysis for vibration %d: %g' % (i, variance))
 472         else:
 473             rel_error = abs(variance - processed_vibration)/processed_vibration
 474             if rel_error > maximum_relative_error:
 475                 _changed_vibration = True
 476                 _LOG.warn(("new analysis doesn't match for vibration %d: "
 477                            "%g -> %g (difference: %g, relative error: %g)")
 478                           % (i, processed_vibration, variance,
 479                              variance-processed_vibration, rel_error))
 480         if _changed_vibration and not dry_run:
 481             changed_vibration = True
 482             _vibration_save(
 483                 filename, vibration_group, processed_vibration=variance)
 484
 485     calib_group = '%scalibration/' % group
 486
 487     if changed_bump and not dry_run:
 488         calib_save(filename, calib_group, bumps=bumps)
 489     if changed_temperature and not dry_run:
 490         calib_save(filename, calib_group, temperatures=temperatures)
 491     if changed_vibration and not dry_run:
 492         calib_save(filename, calib_group, vibrations=vibrations)
 493
 494     new_k,new_k_s = calib_analyze(
 495         bumps=bumps, temperatures=temperatures, vibrations=vibrations)
 496     new_calib_k = False
 497     if k is None:
 498         new_calib_k = True
 499         _LOG.warn('new analysis for k: %g' % new_k)
 500     else:
 501         rel_error = abs(new_k-k)/k
 502         if rel_error > maximum_relative_error:
 503             new_calib_k = True
 504             _LOG.warn(("new analysis doesn't match for k: %g -> %g "
 505                        "(difference: %g, relative error: %g)")
 506                       % (k, new_k, new_k-k, rel_error))
 507     if new_calib_k and not dry_run:
 508         calib_save(filename, calib_group, k=new_k)
 509     new_calib_k_s = False
 510     if k_s is None:
 511         new_calib_k_s = True
 512         _LOG.warn('new analysis for k_s: %g' % new_k_s)
 513     else:
 514         rel_error = abs(new_k_s-k_s)/k_s
 515         if rel_error > maximum_relative_error:
 516             new_calib_k_s = True
 517             _LOG.warn(("new analysis doesn't match for k_s: %g -> %g "
 518                        "(difference: %g, relative error: %g)")
 519                       % (k_s, new_k_s, new_k_s-k_s, rel_error))
 520     if new_calib_k_s and not dry_run:
 521         calib_save(filename, calib_group, k_s=new_k_s)
 522     return (new_k, new_k_s)
 523
 524 def calib_plot_all(bumps, bump_details, temperatures, temperature_details,
 525                    vibrations, vibration_details, calibration_config, k, k_s,
 526                    maximum_relative_error=1e-5):
 527     calib_plot(bumps, temperatures, vibrations)
 528     for i,bump in enumerate(bump_details):
 529         sensitivity = _bump_analyze(
 530             data=bump['raw_bump'], bump_config=bump['bump_config'],
 531             z_axis_config=bump['z_axis_config'],
 532             deflection_channel_config=bump['deflection_channel_config'],
 533             plot=True)
 534         rel_error = abs(sensitivity - bump['processed_bump']
 535                         )/bump['processed_bump']
 536         if rel_error > maximum_relative_error:
 537             _LOG.warn(("new analysis doesn't match for bump %d: %g != %g "
 538                        "(difference: %g, relative error: %g)")
 539                       % (i, sensitivity, bump['processed_bump'],
 540                          sensitivity-bump['processed_bump'], rel_error))
 541     # nothing interesting to plot for temperatures...
 542     for i,vibration in enumerate(vibration_details):
 543         variance = _vibration_analyze(
 544             deflection=vibration['raw_vibration'],
 545             vibration_config=vibration['vibration_config'],
 546             deflection_channel_config=vibration['deflection_channel_config'],
 547             plot=True)
 548         rel_error = abs(variance - vibration['processed_vibration']
 549                         )/vibration['processed_vibration']
 550         if rel_error > maximum_relative_error:
 551             _LOG.warn(("new analysis doesn't match for vibration %d: %g != %g "
 552                        "(difference: %g, relative error: %g)")
 553                       % (i, variance, vibration['processed_vibration'],
 554                          variance-vibration['processed_vibration'], rel_error))