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             calibration_config.save(group=config_cwg)
 238         if bumps is not None:
 239             try:
 240                 del cwg['raw/photodiode-sensitivity/data']
 241             except KeyError:
 242                 pass
 243             try:
 244                 del cwg['raw/photodiode-sensitivity/units']
 245             except KeyError:
 246                 pass
 247             cwg['raw/photodiode-sensitivity/data'] = bumps
 248             cwg['raw/photodiode-sensitivity/units'] = 'V/m'
 249         if temperatures is not None:
 250             try:
 251                 del cwg['raw/temperature/data']
 252             except KeyError:
 253                 pass
 254             try:
 255                 del cwg['raw/temperature/units']
 256             except KeyError:
 257                 pass
 258             cwg['raw/temperature/data'] = temperatures
 259             cwg['raw/temperature/units'] = 'K'
 260         if vibrations is not None:
 261             try:
 262                 del cwg['raw/thermal-vibration-variance/data']
 263             except KeyError:
 264                 pass
 265             try:
 266                 del cwg['raw/thermal-vibration-variance/units']
 267             except KeyError:
 268                 pass
 269             cwg['raw/thermal-vibration-variance/data'] = vibrations
 270             cwg['raw/thermal-vibration-variance/units'] = 'V^2'
 271         if k is not None:
 272             try:
 273                 del cwg['processed/spring-constant/data']
 274             except KeyError:
 275                 pass
 276             try:
 277                 del cwg['processed/spring-constant/units']
 278             except KeyError:
 279                 pass
 280             cwg['processed/spring-constant/data'] = k
 281             cwg['processed/spring-constant/units'] = 'N/m'
 282         if k_s is not None:
 283             try:
 284                 del cwg['processed/spring-constant/standard-deviation']
 285             except KeyError:
 286                 pass
 287             cwg['processed/spring-constant/standard-deviation'] = k_s
 288
 289 def calib_load(filename, group='/'):
 290     assert group.endswith('/')
 291     bumps = temperatures = vibrations = k = k_s = None
 292     configs = []
 293     with _h5py.File(filename, 'a') as f:
 294         try:
 295             bumps = f[group+'raw/photodiode-sensitivity/data'][...]
 296         except KeyError:
 297             pass
 298         try:
 299             temperatures = f[group+'raw/temperature/data'][...]
 300         except KeyError:
 301             pass
 302         try:
 303             vibrations = f[group+'raw/thermal-vibration-variance/data'][...]
 304         except KeyError:
 305             pass
 306         try:
 307             k = f[group+'processed/spring-constant/data'][...]
 308         except KeyError:
 309             pass
 310         try:
 311             k_s = 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     figure.show()
 339
 340
 341 def calib_load_all(filename, group='/'):
 342     "Load all data from a `calib()` run."
 343     assert group.endswith('/'), group
 344     bumps,temperatures,vibrations,calibration_config,k,k_s = calib_load(
 345         filename, group+'calibration/')
 346     bump_details = []
 347     for i in range(calibration_config['num-bumps']):
 348         (raw_bump,bump_config,z_channel_config,z_axis_config,
 349          deflection_channel_config,processed_bump) = _bump_load(
 350             filename=filename, group='%sbump/%d/' % (group, i))
 351         bump_details.append({
 352                 'raw_bump': raw_bump,
 353                 'bump_config': bump_config,
 354                 'z_channel_config': z_channel_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_channel_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_channel_config=z_channel_config,
 407             z_axis_config=z_axis_config,
 408             deflection_channel_config=deflection_channel_config)
 409         bumps[i] = sensitivity
 410         rel_error = abs(sensitivity - processed_bump)/processed_bump
 411         if rel_error > maximum_relative_error:
 412             _LOG.warn(("new analysis doesn't match for bump %d: %g -> %g "
 413                        "(difference: %g, relative error: %g)")
 414                       % (i, processed_bump, sensitivity,
 415                          sensitivity-processed_bump, rel_error))
 416             changed_bump = True
 417             if not dry_run:
 418                 _bump_save(filename, bump_group, processed_bump=sensitivity)
 419     for i in range(calibration_config['num-temperatures']):
 420         temperature_group = '%stemperature/%d/' % (group, i)
 421         (raw_temperature,temperature_config,processed_temperature
 422          ) = _temperature_load(
 423             filename=filename, group=temperature_group)
 424         temperature = _temperature_analyze(
 425             raw_temperature, temperature_config)
 426         temperatures[i] = temperature
 427         rel_error = abs(temperature - processed_temperature
 428                         )/processed_temperature
 429         if rel_error > maximum_relative_error:
 430             _LOG.warn(("new analysis doesn't match for temperature %d: "
 431                        "%g -> %g (difference: %g, relative error: %g)")
 432                       % (i, processed_temperature, temperature,
 433                          temperature-processed_temperature, rel_error))
 434             changed_temperature = True
 435             if not dry_run:
 436                 _temperature_save(
 437                     filename, temperature_group,
 438                     processed_T=temperature)
 439     for i in range(calibration_config['num-vibrations']):
 440         vibration_group = '%svibration/%d/' % (group, i)
 441         (raw_vibration,vibration_config,deflection_channel_config,
 442          processed_vibration) = _vibration_load(
 443             filename=filename, group=vibration_group)
 444         variance = _vibration_analyze(
 445             deflection=raw_vibration, vibration_config=vibration_config,
 446             deflection_channel_config=deflection_channel_config)
 447         vibrations[i] = variance
 448         rel_error = abs(variance - processed_vibration)/processed_vibration
 449         if rel_error > maximum_relative_error:
 450             _LOG.warn(("new analysis doesn't match for vibration %d: %g -> %g "
 451                        "(difference: %g, relative error: %g)")
 452                       % (i, processed_vibration, variance,
 453                          variance-processed_vibration, rel_error))
 454             changed_vibration = True
 455             if not dry_run:
 456                 _vibration_save(
 457                     filename, vibration_group, processed_vibration=variance)
 458
 459     calib_group = '%scalibration/' % group
 460
 461     if changed_bump and not dry_run:
 462         calib_save(filename, calib_group, bumps=bumps)
 463     if changed_temperature and not dry_run:
 464         calib_save(filename, calib_group, temperatures=temperatures)
 465     if changed_vibration and not dry_run:
 466         calib_save(filename, calib_group, vibrations=vibrations)
 467
 468     new_k,new_k_s = calib_analyze(
 469         bumps=bumps, temperatures=temperatures, vibrations=vibrations)
 470     rel_error = abs(new_k-k)/k
 471     if rel_error > maximum_relative_error:
 472         _LOG.warn(("new analysis doesn't match for k: %g -> %g "
 473                    "(difference: %g, relative error: %g)")
 474                   % (k, new_k, new_k-k, rel_error))
 475         if not dry_run:
 476             calib_save(filename, calib_group, k=new_k)
 477     rel_error = abs(new_k_s-k_s)/k_s
 478     if rel_error > maximum_relative_error:
 479         _LOG.warn(("new analysis doesn't match for k_s: %g -> %g "
 480                    "(difference: %g, relative error: %g)")
 481                   % (k_s, new_k_s, new_k_s-k_s, rel_error))
 482         if not dry_run:
 483             calib_save(filename, calib_group, k_s=new_k_s)
 484     return (new_k, new_k_s)
 485
 486 def calib_plot_all(bumps, bump_details, temperatures, temperature_details,
 487                    vibrations, vibration_details, calibration_config, k, k_s,
 488                    maximum_relative_error=1e-5):
 489     calib_plot(bumps, temperatures, vibrations)
 490     for i,bump in enumerate(bump_details):
 491         sensitivity = _bump_analyze(
 492             data=bump['raw_bump'], bump_config=bump['bump_config'],
 493             z_channel_config=bump['z_channel_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))