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