calibcant/bump_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 """Surface bump analysis (measures photodiode sensitivity).
  22
  23 Separate the more general `bump_analyze()` from the other `bump_*()`
  24 functions in calibcant.
  25
  26 The relevant physical quantities are:
  27
  28 * `Vzp_out`  Output z-piezo voltage (what we generate)
  29 * `Vzp`      Applied z-piezo voltage (after external amplification)
  30 * `Zp`       The z-piezo position
  31 * `Zcant`    The cantilever vertical deflection
  32 * `Vphoto`   The photodiode vertical deflection voltage (what we measure)
  33
  34 Which are related by the parameters:
  35
  36 * `zp_gain`           Vzp_out / Vzp
  37 * `zp_sensitivity`    Zp / Vzp
  38 * `photo_sensitivity` Vphoto / Zcant
  39
  40 `photo_sensitivity` is measured by bumping the cantilever against the
  41 surface, where `Zp = Zcant` (see `calibrate.bump_acquire()`).  The
  42 measured slope `Vphoto/Vout` is converted to `photo_sensitivity` with
  43 `bump_analyze()`.
  44
  45 >>> import os
  46 >>> from pprint import pprint
  47 >>> import tempfile
  48 >>> import numpy
  49 >>> from .config import HDF5_BumpConfig
  50 >>> from h5config.hdf5 import pprint_HDF5, HDF5_ChannelConfig, HDF5_AxisConfig
  51
  52 >>> fd,filename = tempfile.mkstemp(suffix='.h5', prefix='calibcant-')
  53 >>> os.close(fd)
  54
  55 >>> bump_config = HDF5_BumpConfig(filename=filename, group='/bump/config/')
  56 >>> bump_config.save()
  57 >>> z_channel_config = HDF5_ChannelConfig(filename=None)
  58 >>> z_channel_config['maxdata'] = 200
  59 >>> z_channel_config['conversion-coefficients'] = (0,1)
  60 >>> z_channel_config['conversion-origin'] = 0
  61 >>> z_axis_config = HDF5_AxisConfig(filename=None)
  62 >>> deflection_channel_config = HDF5_ChannelConfig(filename=None)
  63 >>> deflection_channel_config['maxdata'] = 200
  64 >>> deflection_channel_config['conversion-coefficients'] = (0,1)
  65 >>> deflection_channel_config['conversion-origin'] = 0
  66
  67 >>> raw_bump = {
  68 ...     'z': numpy.arange(100, dtype=numpy.uint16),
  69 ...     'deflection': numpy.arange(100, dtype=numpy.uint16),
  70 ...     }
  71 >>> raw_bump['deflection'][:50] = 50
  72 >>> processed_bump = bump_analyze(
  73 ...     raw_bump, bump_config, z_channel_config, z_axis_config,
  74 ...     deflection_channel_config)
  75 >>> bump_plot(data=raw_bump)  # TODO: convert to V and m
  76 >>> bump_save(filename=filename, group='/bump/', raw_bump=raw_bump,
  77 ...     z_channel_config=z_channel_config, z_axis_config=z_axis_config,
  78 ...     deflection_channel_config=deflection_channel_config,
  79 ...     processed_bump=processed_bump)
  80
  81 >>> pprint_HDF5(filename)  # doctest: +ELLIPSIS, +REPORT_UDIFF
  82 /
  83   /bump
  84     /bump/config
  85       /bump/config/deflection
  86         /bump/config/deflection/channel
  87           <HDF5 dataset "channel": shape (), type "|S1">
  88             0
  89           <HDF5 dataset "conversion-coefficients": shape (), type "|S4">
  90             0, 1
  91           <HDF5 dataset "conversion-origin": shape (), type "|S1">
  92             0
  93           <HDF5 dataset "device": shape (), type "|S12">
  94             /dev/comedi0
  95           <HDF5 dataset "inverse-conversion-coefficients": shape (), type "|S1">
  96 <BLANKLINE>
  97           <HDF5 dataset "inverse-conversion-origin": shape (), type "|S3">
  98             1.0
  99           <HDF5 dataset "maxdata": shape (), type "|S3">
 100             200
 101           <HDF5 dataset "range": shape (), type "|S1">
 102             1
 103           <HDF5 dataset "subdevice": shape (), type "|S2">
 104             -1
 105       <HDF5 dataset "far-steps": shape (), type "|S3">
 106         200
 107       <HDF5 dataset "initial-position": shape (), type "|S6">
 108         -5e-08
 109       <HDF5 dataset "model": shape (), type "|S9">
 110         quadratic
 111       <HDF5 dataset "push-depth": shape (), type "|S5">
 112         2e-07
 113       <HDF5 dataset "push-speed": shape (), type "|S5">
 114         1e-06
 115       <HDF5 dataset "samples": shape (), type "|S4">
 116         1024
 117       <HDF5 dataset "setpoint": shape (), type "|S3">
 118         2.0
 119       /bump/config/z
 120         /bump/config/z/axis
 121           <HDF5 dataset "gain": shape (), type "|S3">
 122             1.0
 123           <HDF5 dataset "maximum": shape (), type "|S4">
 124             None
 125           <HDF5 dataset "minimum": shape (), type "|S4">
 126             None
 127           <HDF5 dataset "sensitivity": shape (), type "|S3">
 128             1.0
 129         /bump/config/z/channel
 130           <HDF5 dataset "channel": shape (), type "|S1">
 131             0
 132           <HDF5 dataset "conversion-coefficients": shape (), type "|S4">
 133             0, 1
 134           <HDF5 dataset "conversion-origin": shape (), type "|S1">
 135             0
 136           <HDF5 dataset "device": shape (), type "|S12">
 137             /dev/comedi0
 138           <HDF5 dataset "inverse-conversion-coefficients": shape (), type "|S1">
 139 <BLANKLINE>
 140           <HDF5 dataset "inverse-conversion-origin": shape (), type "|S3">
 141             1.0
 142           <HDF5 dataset "maxdata": shape (), type "|S3">
 143             200
 144           <HDF5 dataset "range": shape (), type "|S1">
 145             1
 146           <HDF5 dataset "subdevice": shape (), type "|S2">
 147             -1
 148     <HDF5 dataset "processed": shape (), type "<f8">
 149       1.00000002115
 150     /bump/raw
 151       <HDF5 dataset "deflection": shape (100,), type "<u2">
 152         [50 50 ... 50 51 52 ... 97 98 99]
 153       <HDF5 dataset "z": shape (100,), type "<u2">
 154         [ 0  1  2  3  ... 97 98 99]
 155
 156 >>> (raw_bump,bump_config,z_channel_config,z_axis_config,
 157 ...  deflection_channel_config,processed_bump) = bump_load(
 158 ...     filename=filename, group='/bump/')
 159
 160 >>> pprint(raw_bump)  # doctest: +ELLIPSIS
 161 {'deflection': array([50, 50, ... 51, 52, 53, ..., 97, 98, 99], dtype=uint16),
 162  'z': array([ 0,  1,  2,  ..., 97, 98, 99], dtype=uint16)}
 163 >>> processed_bump  # doctest: +ELLIPSIS
 164 1.00...
 165
 166 >>> os.remove(filename)
 167 """
 168
 169 import h5py as _h5py
 170 import numpy as _numpy
 171 from scipy.optimize import leastsq as _leastsq
 172
 173 try:
 174     import matplotlib as _matplotlib
 175     import matplotlib.pyplot as _matplotlib_pyplot
 176     import time as _time  # for timestamping lines on plots
 177 except (ImportError, RuntimeError), e:
 178     _matplotlib = None
 179     _matplotlib_import_error = e
 180
 181 from h5config.hdf5 import h5_create_group as _h5_create_group
 182 from pypiezo.base import convert_bits_to_volts as _convert_bits_to_volts
 183 from pypiezo.base import convert_bits_to_meters as _convert_bits_to_meters
 184 from pypiezo.config import HDF5_ChannelConfig as _HDF5_ChannelConfig
 185 from pypiezo.config import HDF5_AxisConfig as _HDF5_AxisConfig
 186
 187 from . import LOG as _LOG
 188 from . import package_config as _package_config
 189 from .config import Linear as _Linear
 190 from .config import Quadratic as _Quadratic
 191 from .config import HDF5_BumpConfig as _HDF5_BumpConfig
 192
 193
 194 def bump_analyze(data, bump_config, z_channel_config, z_axis_config,
 195                  deflection_channel_config, plot=False):
 196     """Return the slope of the bump.
 197
 198     Inputs:
 199       data              dictionary of data in DAC/ADC bits
 200       bump_config       `.config._BumpConfig` instance
 201       z_channel_config  z `pypiezo.config.ChannelConfig` instance
 202       z_axis_config     z `pypiezo.config.AxisConfig` instance
 203       deflection_channel_config
 204                         deflection `pypiezo.config.ChannelConfig` instance
 205       plot              boolean overriding matplotlib config setting.
 206     Returns:
 207       photo_sensitivity (Vphoto/Zcant) in Volts/m
 208
 209     Checks for strong correlation (r-value) and low randomness chance
 210     (p-value).
 211     """
 212     z = _convert_bits_to_meters(z_channel_config, z_axis_config, data['z'])
 213     deflection = _convert_bits_to_volts(
 214         deflection_channel_config, data['deflection'])
 215     if bump_config['model'] == _Linear:
 216         kwargs = {
 217             'param_guesser': limited_linear_param_guess,
 218             'model': limited_linear,
 219             'sensitivity_from_fit_params': limited_linear_sensitivity,
 220             }
 221     else:  # _Quadratic
 222         kwargs = {
 223             'param_guesser': limited_quadratic_param_guess,
 224             'model': limited_quadratic,
 225             'sensitivity_from_fit_params': limited_quadratic_sensitivity,
 226             }
 227     photo_sensitivity = bump_fit(z, deflection, plot=plot, **kwargs)
 228     return photo_sensitivity
 229
 230 def limited_linear(x, params):
 231     """
 232     Model the bump as:
 233       flat region (off-surface)
 234       linear region (in-contact)
 235       flat region (high-voltage-rail)
 236     Parameters:
 237       x_contact (x value for the surface-contact kink)
 238       y_contact (y value for the surface-contact kink)
 239       slope (dy/dx at the surface-contact kink)
 240     """
 241     high_voltage_rail = 2**16 - 1 # bits
 242     x_contact,y_contact,slope = params
 243     y = slope*(x-x_contact) + y_contact
 244     y = _numpy.clip(y, y_contact, high_voltage_rail)
 245     return y
 246
 247 def limited_linear_param_guess(x, y):
 248     """
 249     Guess rough parameters for a limited_linear model.  Assumes the
 250     bump approaches (raising the deflection as it does so) first.
 251     Retracting after the approach is optional.  Approximates the contact
 252     position and an on-surface (high) position by finding first crossings
 253     of thresholds 0.3 and 0.7 of the y value's total range.  Not the
 254     most efficient algorithm, but it seems fairly robust.
 255     """
 256     y_contact = float(y.min())
 257     y_max = float(y.max())
 258     i = 0
 259     y_low  = y_contact + 0.3 * (y_max-y_contact)
 260     y_high = y_contact + 0.7 * (y_max-y_contact)
 261     while y[i] < y_low:
 262         i += 1
 263     i_low = i
 264     while y[i] < y_high:
 265         i += 1
 266     i_high = i
 267     x_contact = float(x[i_low])
 268     x_high = float(x[i_high])
 269     slope = (y_high - y_contact) / (x_high - x_contact)
 270     return (x_contact, y_contact, slope)
 271
 272 def limited_linear_sensitivity(params):
 273     """
 274     Return the estimated sensitivity to small deflections according to
 275     limited_linear fit parameters.
 276     """
 277     slope = params[2]
 278     return slope
 279
 280 def limited_quadratic(x, params):
 281     """
 282     Model the bump as:
 283       flat region (off-surface)
 284       quadratic region (in-contact)
 285       flat region (high-voltage-rail)
 286     Parameters:
 287       x_contact (x value for the surface-contact kink)
 288       y_contact (y value for the surface-contact kink)
 289       slope (dy/dx at the surface-contact kink)
 290       quad (d**2 y / dx**2, allow decreasing sensitivity with increased x)
 291     """
 292     high_voltage_rail = 2**16 - 1 # bits
 293     x_contact,y_contact,slope,quad = params
 294     y = slope*(x-x_contact) + quad*(x-x_contact)**2+ y_contact
 295     y = _numpy.clip(y, y_contact, high_voltage_rail)
 296     return y
 297
 298 def limited_quadratic_param_guess(x, y):
 299     """
 300     Guess rough parameters for a limited_quadratic model.  Assumes the
 301     bump approaches (raising the deflection as it does so) first.
 302     Retracting after the approach is optional.  Approximates the contact
 303     position and an on-surface (high) position by finding first crossings
 304     of thresholds 0.3 and 0.7 of the y value's total range.  Not the
 305     most efficient algorithm, but it seems fairly robust.
 306     """
 307     x_contact,y_contact,linear_slope = limited_linear_param_guess(x,y)
 308     contact_depth = x.max() - x_contact
 309     slope = linear_slope / 2
 310     quad = slope / contact_depth
 311     # The above slope and quad were chosen so
 312     #   x = contact_depth
 313     #   x*slope + x**2 * slope == x * linear_slope
 314     return (x_contact, y_contact, slope, quad)
 315
 316 def limited_quadratic_sensitivity(params):
 317     """
 318     Return the estimated sensitivity to small deflections according to
 319     limited_quadratic fit parameters.
 320     """
 321     slope = params[2]
 322     return slope
 323
 324 def bump_fit(z, deflection,
 325              param_guesser=limited_quadratic_param_guess,
 326              model=limited_quadratic,
 327              sensitivity_from_fit_params=limited_quadratic_sensitivity,
 328              plot=False):
 329     """Fit a aurface bump and return the photodiode sensitivitiy.
 330
 331     Parameters:
 332       z              piezo position in meters
 333       deflection     photodiode deflection in volts
 334       param_guesser  function that guesses initial model parameters
 335       model          parametric model of deflection vs. z
 336       sensitivity_from_fit_params  given fit params, return the sensitivity
 337       plot              boolean overriding matplotlib config setting.
 338     Returns:
 339       photodiode_sensitivity  photodiode volts per piezo meter
 340     """
 341     _LOG.debug('fit bump data with model %s' % model)
 342     def residual(p, deflection, z):
 343         return model(z, p) - deflection
 344     param_guess = param_guesser(z, deflection)
 345     p,cov,info,mesg,ier = _leastsq(
 346         residual, param_guess, args=(deflection, z), full_output=True,
 347         maxfev=int(10e3))
 348     _LOG.debug('fitted params: %s' % p)
 349     _LOG.debug('covariance matrix: %s' % cov)
 350     #_LOG.debug('info: %s' % info)
 351     _LOG.debug('message: %s' % mesg)
 352     if ier == 1:
 353         _LOG.debug('solution converged')
 354     else:
 355         _LOG.debug('solution did not converge')
 356     if plot or _package_config['matplotlib']:
 357         yguess = model(z, param_guess)
 358         yfit = model(z, p)
 359         bump_plot({'z': z, 'deflection': deflection}, yguess=yguess, yfit=yfit)
 360     return sensitivity_from_fit_params(p)
 361
 362 def bump_save(filename, group='/', raw_bump=None, bump_config=None,
 363               z_channel_config=None, z_axis_config=None,
 364               deflection_channel_config=None, processed_bump=None):
 365     with _h5py.File(filename, 'a') as f:
 366         cwg = _h5_create_group(f, group)
 367         if raw_bump is not None:
 368             for k in ['z', 'deflection']:
 369                 try:
 370                     del cwg['raw/{}'.format(k)]
 371                 except KeyError:
 372                     pass
 373             cwg['raw/z'] = raw_bump['z']
 374             cwg['raw/deflection'] = raw_bump['deflection']
 375         for config,key in [(bump_config, 'config/bump'),
 376                            (z_channel_config, 'config/z/channel'),
 377                            (z_axis_config, 'config/z/axis'),
 378                            (deflection_channel_config,
 379                             'config/deflection/channel')]:
 380             if config is None:
 381                 continue
 382             config_cwg = _h5_create_group(cwg, key)
 383             config.save(group=config_cwg)
 384         if processed_bump is not None:
 385             try:
 386                 del cwg['processed']
 387             except KeyError:
 388                 pass
 389             cwg['processed'] = processed_bump
 390
 391 def bump_load(filename, group='/'):
 392     assert group.endswith('/')
 393     raw_bump = processed_bump = None
 394     configs = []
 395     with _h5py.File(filename, 'a') as f:
 396         try:
 397             raw_bump = {
 398                 'z': f[group+'raw/z'][...],
 399                 'deflection': f[group+'raw/deflection'][...],
 400                 }
 401         except KeyError:
 402             pass
 403         for Config,key in [(_HDF5_BumpConfig, 'config/bump'),
 404                            (_HDF5_ChannelConfig, 'config/z/channel'),
 405                            (_HDF5_AxisConfig, 'config/z/axis'),
 406                            (_HDF5_ChannelConfig, 'config/deflection/channel')]:
 407             config = Config(filename=filename, group=group+key)
 408             configs.append(config)
 409         try:
 410             processed_bump = f[group+'processed'][...]
 411         except KeyError:
 412             pass
 413     ret = [raw_bump]
 414     ret.extend(configs)
 415     ret.append(processed_bump)
 416     for config in configs:
 417         config.load()
 418     return tuple(ret)
 419
 420 def bump_plot(data, yguess=None, yfit=None):
 421     "Plot the bump (Vphoto vs Vzp)"
 422     if not _matplotlib:
 423         raise _matplotlib_import_error
 424     figure = _matplotlib_pyplot.figure()
 425     if yfit is None:
 426         axes = figure.add_subplot(1, 1, 1)
 427     else:
 428         axes = figure.add_subplot(2, 1, 1)
 429         residual_axes = figure.add_subplot(2, 1, 2)
 430     timestamp = _time.strftime('%H%M%S')
 431
 432     axes.hold(True)
 433     axes.plot(data['z'], data['deflection'], '.', label='bump')
 434     if yguess != None:
 435         axes.plot(data['z'], yguess, 'g-', label='guess')
 436     if yfit != None:
 437         axes.plot(data['z'], yfit, 'r-', label='fit')
 438
 439     axes.set_title('bump surface %s' % timestamp)
 440     #axes.legend(loc='upper left')
 441     axes.set_xlabel('Z piezo (meters)')
 442     axes.set_ylabel('Photodiode (Volts)')
 443     if yfit is not None:
 444         # second subplot for residual
 445         residual_axes.plot(data['z'], data['deflection'] - yfit,
 446                            'r-', label='residual')
 447         #residual_axes.legend(loc='upper right')
 448         residual_axes.set_xlabel('Z piezo (meters)')
 449         residual_axes.set_ylabel('Photodiode (Volts)')
 450     figure.show()