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