calibcant/analyze.py

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