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     changed_bump = changed_temperature = changed_vibration = False
 399     for i in range(calibration_config['num-bumps']):
 400         bump_group = '%sbump/%d/' % (group, i)
 401         (raw_bump,bump_config,z_axis_config,
 402          deflection_channel_config,processed_bump) = _bump_load(
 403             filename=filename, group=bump_group)
 404         sensitivity = _bump_analyze(
 405             data=raw_bump, bump_config=bump_config,
 406             z_axis_config=z_axis_config,
 407             deflection_channel_config=deflection_channel_config)
 408         bumps[i] = sensitivity
 409         rel_error = abs(sensitivity - processed_bump)/processed_bump
 410         if rel_error > maximum_relative_error:
 411             _LOG.warn(("new analysis doesn't match for bump %d: %g -> %g "
 412                        "(difference: %g, relative error: %g)")
 413                       % (i, processed_bump, sensitivity,
 414                          sensitivity-processed_bump, rel_error))
 415             changed_bump = True
 416             if not dry_run:
 417                 _bump_save(filename, bump_group, processed_bump=sensitivity)
 418     for i in range(calibration_config['num-temperatures']):
 419         temperature_group = '%stemperature/%d/' % (group, i)
 420         (raw_temperature,temperature_config,processed_temperature
 421          ) = _temperature_load(
 422             filename=filename, group=temperature_group)
 423         temperature = _temperature_analyze(
 424             raw_temperature, temperature_config)
 425         temperatures[i] = temperature
 426         rel_error = abs(temperature - processed_temperature
 427                         )/processed_temperature
 428         if rel_error > maximum_relative_error:
 429             _LOG.warn(("new analysis doesn't match for temperature %d: "
 430                        "%g -> %g (difference: %g, relative error: %g)")
 431                       % (i, processed_temperature, temperature,
 432                          temperature-processed_temperature, rel_error))
 433             changed_temperature = True
 434             if not dry_run:
 435                 _temperature_save(
 436                     filename, temperature_group,
 437                     processed_T=temperature)
 438     for i in range(calibration_config['num-vibrations']):
 439         vibration_group = '%svibration/%d/' % (group, i)
 440         (raw_vibration,vibration_config,deflection_channel_config,
 441          processed_vibration) = _vibration_load(
 442             filename=filename, group=vibration_group)
 443         variance = _vibration_analyze(
 444             deflection=raw_vibration, vibration_config=vibration_config,
 445             deflection_channel_config=deflection_channel_config)
 446         vibrations[i] = variance
 447         rel_error = abs(variance - processed_vibration)/processed_vibration
 448         if rel_error > maximum_relative_error:
 449             _LOG.warn(("new analysis doesn't match for vibration %d: %g -> %g "
 450                        "(difference: %g, relative error: %g)")
 451                       % (i, processed_vibration, variance,
 452                          variance-processed_vibration, rel_error))
 453             changed_vibration = True
 454             if not dry_run:
 455                 _vibration_save(
 456                     filename, vibration_group, processed_vibration=variance)
 457
 458     calib_group = '%scalibration/' % group
 459
 460     if changed_bump and not dry_run:
 461         calib_save(filename, calib_group, bumps=bumps)
 462     if changed_temperature and not dry_run:
 463         calib_save(filename, calib_group, temperatures=temperatures)
 464     if changed_vibration and not dry_run:
 465         calib_save(filename, calib_group, vibrations=vibrations)
 466
 467     new_k,new_k_s = calib_analyze(
 468         bumps=bumps, temperatures=temperatures, vibrations=vibrations)
 469     rel_error = abs(new_k-k)/k
 470     if rel_error > maximum_relative_error:
 471         _LOG.warn(("new analysis doesn't match for k: %g -> %g "
 472                    "(difference: %g, relative error: %g)")
 473                   % (k, new_k, new_k-k, rel_error))
 474         if not dry_run:
 475             calib_save(filename, calib_group, k=new_k)
 476     rel_error = abs(new_k_s-k_s)/k_s
 477     if rel_error > maximum_relative_error:
 478         _LOG.warn(("new analysis doesn't match for k_s: %g -> %g "
 479                    "(difference: %g, relative error: %g)")
 480                   % (k_s, new_k_s, new_k_s-k_s, rel_error))
 481         if not dry_run:
 482             calib_save(filename, calib_group, k_s=new_k_s)
 483     return (new_k, new_k_s)
 484
 485 def calib_plot_all(bumps, bump_details, temperatures, temperature_details,
 486                    vibrations, vibration_details, calibration_config, k, k_s,
 487                    maximum_relative_error=1e-5):
 488     calib_plot(bumps, temperatures, vibrations)
 489     for i,bump in enumerate(bump_details):
 490         sensitivity = _bump_analyze(
 491             data=bump['raw_bump'], bump_config=bump['bump_config'],
 492             z_axis_config=bump['z_axis_config'],
 493             deflection_channel_config=bump['deflection_channel_config'],
 494             plot=True)
 495         rel_error = abs(sensitivity - bump['processed_bump']
 496                         )/bump['processed_bump']
 497         if rel_error > maximum_relative_error:
 498             _LOG.warn(("new analysis doesn't match for bump %d: %g != %g "
 499                        "(difference: %g, relative error: %g)")
 500                       % (i, sensitivity, bump['processed_bump'],
 501                          sensitivity-bump['processed_bump'], rel_error))
 502     # nothing interesting to plot for temperatures...
 503     for i,vibration in enumerate(vibration_details):
 504         variance = _vibration_analyze(
 505             deflection=vibration['raw_vibration'],
 506             vibration_config=vibration['vibration_config'],
 507             deflection_channel_config=vibration['deflection_channel_config'],
 508             plot=True)
 509         rel_error = abs(variance - vibration['processed_vibration']
 510                         )/vibration['processed_vibration']
 511         if rel_error > maximum_relative_error:
 512             _LOG.warn(("new analysis doesn't match for vibration %d: %g != %g "
 513                        "(difference: %g, relative error: %g)")
 514                       % (i, variance, vibration['processed_vibration'],
 515                          variance-vibration['processed_vibration'], rel_error))