calibcant/bump_analyze.py

   1 # calibcant - tools for thermally calibrating AFM cantilevers
   2 #
   3 # Copyright (C) 2008-2012 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 modify it under
   8 # the terms of the GNU General Public License as published by the Free Software
   9 # Foundation, either version 3 of the License, or (at your option) any later
  10 # version.
  11 #
  12 # calibcant is distributed in the hope that it will be useful, but WITHOUT ANY
  13 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  14 # A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
  15 #
  16 # You should have received a copy of the GNU General Public License along with
  17 # calibcant.  If not, see <http://www.gnu.org/licenses/>.
  18
  19 """Surface bump analysis (measures photodiode sensitivity).
  20
  21 Separate the more general `analyze()` from the other `bump_*()`
  22 functions in calibcant.
  23
  24 The relevant physical quantities are:
  25
  26 * `Vzp_out`  Output z-piezo voltage (what we generate)
  27 * `Vzp`      Applied z-piezo voltage (after external amplification)
  28 * `Zp`       The z-piezo position
  29 * `Zcant`    The cantilever vertical deflection
  30 * `Vphoto`   The photodiode vertical deflection voltage (what we measure)
  31
  32 Which are related by the parameters:
  33
  34 * `zp_gain`           Vzp_out / Vzp
  35 * `zp_sensitivity`    Zp / Vzp
  36 * `photo_sensitivity` Vphoto / Zcant
  37
  38 `photo_sensitivity` is measured by bumping the cantilever against the
  39 surface, where `Zp = Zcant` (see `calibrate.bump_acquire()`).  The
  40 measured slope `Vphoto/Vout` is converted to `photo_sensitivity` with
  41 `analyze()`.
  42
  43 >>> import os
  44 >>> from pprint import pprint
  45 >>> import tempfile
  46 >>> import numpy
  47 >>> from h5config.storage.hdf5 import pprint_HDF5
  48 >>> from pypiezo.config import ChannelConfig, AxisConfig
  49 >>> from .config import BumpConfig
  50
  51 >>> fd,filename = tempfile.mkstemp(suffix='.h5', prefix='calibcant-')
  52 >>> os.close(fd)
  53
  54 >>> config = BumpConfig()
  55 >>> z_channel_config = ChannelConfig()
  56 >>> z_channel_config['name'] = 'z'
  57 >>> z_channel_config['maxdata'] = 200
  58 >>> z_channel_config['conversion-coefficients'] = (0,1)
  59 >>> z_channel_config['conversion-origin'] = 0
  60 >>> z_axis_config = AxisConfig()
  61 >>> z_axis_config['channel'] = z_channel_config
  62 >>> deflection_channel_config = ChannelConfig()
  63 >>> deflection_channel_config['name'] = 'deflection'
  64 >>> deflection_channel_config['maxdata'] = 200
  65 >>> deflection_channel_config['conversion-coefficients'] = (0,1)
  66 >>> deflection_channel_config['conversion-origin'] = 0
  67
  68 >>> raw = {
  69 ...     'z': numpy.arange(100, dtype=numpy.uint16),
  70 ...     'deflection': numpy.arange(100, dtype=numpy.uint16),
  71 ...     }
  72 >>> raw['deflection'][:50] = 50
  73 >>> processed = analyze(
  74 ...     config=config, data=raw, z_axis_config=z_axis_config,
  75 ...     deflection_channel_config=deflection_channel_config)
  76 >>> plot(data=raw)  # TODO: convert to V and m
  77 >>> save(filename=filename, group='/bump/',
  78 ...     config=config, raw=raw, processed=processed)
  79
  80 >>> pprint_HDF5(filename)  # doctest: +ELLIPSIS, +REPORT_UDIFF
  81 /
  82   /bump
  83     /bump/config
  84       /bump/config/bump
  85         <HDF5 dataset "far-steps": shape (), type "<i4">
  86           200
  87         <HDF5 dataset "initial-position": shape (), type "<f8">
  88           -5e-08
  89         <HDF5 dataset "min-slope-ratio": shape (), type "<f8">
  90           10.0
  91         <HDF5 dataset "model": shape (), type "|S9">
  92           quadratic
  93         <HDF5 dataset "push-depth": shape (), type "<f8">
  94           2e-07
  95         <HDF5 dataset "push-speed": shape (), type "<f8">
  96           1e-06
  97         <HDF5 dataset "samples": shape (), type "<i4">
  98           1024
  99         <HDF5 dataset "setpoint": shape (), type "<f8">
 100           2.0
 101     /bump/processed
 102       <HDF5 dataset "data": shape (), type "<f8">
 103         1.00...
 104       <HDF5 dataset "units": shape (), type "|S3">
 105         V/m
 106     /bump/raw
 107       /bump/raw/deflection
 108         <HDF5 dataset "data": shape (100,), type "<u2">
 109           [50 50 ... 50 51 52 ... 97 98 99]
 110         <HDF5 dataset "units": shape (), type "|S4">
 111           bits
 112       /bump/raw/z
 113         <HDF5 dataset "data": shape (100,), type "<u2">
 114           [ 0  1  2  3  ... 97 98 99]
 115         <HDF5 dataset "units": shape (), type "|S4">
 116           bits
 117
 118 >>> data = load(filename=filename, group='/bump/')
 119
 120 >>> pprint(data)  # doctest: +ELLIPSIS, +REPORT_UDIFF
 121 {'config': {'bump': <BumpConfig ...>},
 122  'processed': 1.00...,
 123  'raw': {'deflection': array([50, 50, ..., 52, 53, ..., 98, 99], dtype=uint16),
 124          'z': array([ 0,  1,  2,  ..., 97, 98, 99], dtype=uint16)}}
 125
 126 >>> os.remove(filename)
 127 """
 128
 129 import numpy as _numpy
 130 from scipy.optimize import leastsq as _leastsq
 131
 132 try:
 133     import matplotlib as _matplotlib
 134     import matplotlib.pyplot as _matplotlib_pyplot
 135     import time as _time  # for timestamping lines on plots
 136 except (ImportError, RuntimeError), e:
 137     _matplotlib = None
 138     _matplotlib_import_error = e
 139
 140 from pypiezo.base import convert_bits_to_volts as _convert_bits_to_volts
 141 from pypiezo.base import convert_bits_to_meters as _convert_bits_to_meters
 142 from pypiezo.config import AxisConfig as _AxisConfig
 143 from pypiezo.config import InputChannelConfig as _InputChannelConfig
 144
 145 from . import LOG as _LOG
 146 from . import package_config as _package_config
 147 from .config import Linear as _Linear
 148 from .config import Quadratic as _Quadratic
 149 from .config import BumpConfig as _BumpConfig
 150 from .util import SaveSpec as _SaveSpec
 151 from .util import save as _save
 152 from .util import load as _load
 153
 154
 155 def analyze(config, data, z_axis_config,
 156             deflection_channel_config, plot=False):
 157     """Return the slope of the bump.
 158
 159     Inputs:
 160       data              dictionary of data in DAC/ADC bits
 161       config            `.config._BumpConfig` instance
 162       z_axis_config     z `pypiezo.config.AxisConfig` instance
 163       deflection_channel_config
 164                         deflection `pypiezo.config.InputChannelConfig` instance
 165       plot              boolean overriding matplotlib config setting.
 166     Returns:
 167       photo_sensitivity (Vphoto/Zcant) in Volts/m
 168
 169     Checks for strong correlation (r-value) and low randomness chance
 170     (p-value).
 171     """
 172     z = _convert_bits_to_meters(z_axis_config, data['z'])
 173     deflection = _convert_bits_to_volts(
 174         deflection_channel_config, data['deflection'])
 175     high_voltage_rail = _convert_bits_to_volts(
 176         deflection_channel_config, deflection_channel_config['maxdata'])
 177     if config['model'] == _Linear:
 178         kwargs = {
 179             'param_guesser': limited_linear_param_guess,
 180             'model': limited_linear,
 181             'sensitivity_from_fit_params': limited_linear_sensitivity,
 182             }
 183     else:  # _Quadratic
 184         kwargs = {
 185             'param_guesser': limited_quadratic_param_guess,
 186             'model': limited_quadratic,
 187             'sensitivity_from_fit_params': limited_quadratic_sensitivity,
 188             }
 189     photo_sensitivity = fit(
 190         z, deflection, high_voltage_rail=high_voltage_rail, plot=plot,
 191         **kwargs)
 192     return photo_sensitivity
 193
 194 def limited_linear(x, params, high_voltage_rail):
 195     """
 196     Model the bump as:
 197       flat region (off-surface)
 198       linear region (in-contact)
 199       flat region (high-voltage-rail)
 200     Parameters:
 201       x_contact (x value for the surface-contact kink)
 202       y_contact (y value for the surface-contact kink)
 203       slope (dy/dx at the surface-contact kink)
 204     """
 205     x_contact,y_contact,slope = params
 206     off_surface_mask = x <= x_contact
 207     on_surface_mask = x > x_contact
 208     y = (off_surface_mask * y_contact +
 209          on_surface_mask * (y_contact + slope*(x-x_contact)))
 210     y = _numpy.clip(y, y_contact, high_voltage_rail)
 211     return y
 212
 213 def limited_linear_param_guess(x, y):
 214     """
 215     Guess rough parameters for a limited_linear model.  Assumes the
 216     bump approaches (raising the deflection as it does so) first.
 217     Retracting after the approach is optional.  Approximates the contact
 218     position and an on-surface (high) position by finding first crossings
 219     of thresholds 0.3 and 0.7 of the y value's total range.  Not the
 220     most efficient algorithm, but it seems fairly robust.
 221     """
 222     y_contact = float(y.min())
 223     y_max = float(y.max())
 224     i = 0
 225     y_low  = y_contact + 0.3 * (y_max-y_contact)
 226     y_high = y_contact + 0.7 * (y_max-y_contact)
 227     while y[i] < y_low:
 228         i += 1
 229     i_low = i
 230     while y[i] < y_high:
 231         i += 1
 232     i_high = i
 233     x_contact = float(x[i_low])
 234     x_high = float(x[i_high])
 235     if x_high == x_contact:  # things must be pretty flat
 236         x_contact = (x_contact + x[0]) / 2
 237     slope = (y_high - y_contact) / (x_high - x_contact)
 238     return (x_contact, y_contact, slope)
 239
 240 def limited_linear_sensitivity(params):
 241     """
 242     Return the estimated sensitivity to small deflections according to
 243     limited_linear fit parameters.
 244     """
 245     slope = params[2]
 246     return slope
 247
 248 def limited_quadratic(x, params, high_voltage_rail):
 249     """
 250     Model the bump as:
 251       flat region (off-surface)
 252       quadratic region (in-contact)
 253       flat region (high-voltage-rail)
 254     Parameters:
 255       x_contact (x value for the surface-contact kink)
 256       y_contact (y value for the surface-contact kink)
 257       slope (dy/dx at the surface-contact kink)
 258       quad (d**2 y / dx**2, allow decreasing sensitivity with increased x)
 259     """
 260     x_contact,y_contact,slope,quad = params
 261     off_surface_mask = x <= x_contact
 262     on_surface_mask = x > x_contact
 263     y = (off_surface_mask * y_contact +
 264          on_surface_mask * (
 265             y_contact + slope*(x-x_contact) + quad*(x-x_contact)**2))
 266     y = _numpy.clip(y, y_contact, high_voltage_rail)
 267     return y
 268
 269 def limited_quadratic_param_guess(x, y):
 270     """
 271     Guess rough parameters for a limited_quadratic model.  Assumes the
 272     bump approaches (raising the deflection as it does so) first.
 273     Retracting after the approach is optional.  Approximates the contact
 274     position and an on-surface (high) position by finding first crossings
 275     of thresholds 0.3 and 0.7 of the y value's total range.  Not the
 276     most efficient algorithm, but it seems fairly robust.
 277     """
 278     x_contact,y_contact,linear_slope = limited_linear_param_guess(x,y)
 279     contact_depth = x.max() - x_contact
 280     slope = linear_slope / 2
 281     quad = slope / contact_depth
 282     # The above slope and quad were chosen so
 283     #   x = contact_depth
 284     #   x*slope + x**2 * slope == x * linear_slope
 285     return (x_contact, y_contact, slope, quad)
 286
 287 def limited_quadratic_sensitivity(params):
 288     """
 289     Return the estimated sensitivity to small deflections according to
 290     limited_quadratic fit parameters.
 291     """
 292     slope = params[2]
 293     return slope
 294
 295 def fit(z, deflection, high_voltage_rail,
 296         param_guesser=limited_quadratic_param_guess,
 297         model=limited_quadratic,
 298         sensitivity_from_fit_params=limited_quadratic_sensitivity,
 299         plot=False):
 300     """Fit a aurface bump and return the photodiode sensitivitiy.
 301
 302     Parameters:
 303       z              piezo position in meters
 304       deflection     photodiode deflection in volts
 305       param_guesser  function that guesses initial model parameters
 306       model          parametric model of deflection vs. z
 307       sensitivity_from_fit_params  given fit params, return the sensitivity
 308       plot              boolean overriding matplotlib config setting.
 309     Returns:
 310       photodiode_sensitivity  photodiode volts per piezo meter
 311     """
 312     _LOG.debug('fit bump data with model %s' % model)
 313     def residual(p, deflection, z):
 314         return model(z, p, high_voltage_rail=high_voltage_rail) - deflection
 315     param_guess = param_guesser(z, deflection)
 316     try:
 317         p,cov,info,mesg,ier = _leastsq(
 318             residual, param_guess, args=(deflection, z), full_output=True,
 319             maxfev=int(10e3))
 320     except ValueError:
 321         zd = _numpy.ndarray(list(z.shape) + [2], dtype=z.dtype)
 322         zd[:,0] = z
 323         zd[:,1] = d
 324         _numpy.savetxt('/tmp/z-deflection.dat', zd, delimiter='\t')
 325         raise
 326     _LOG.debug('fitted params: %s' % p)
 327     _LOG.debug('covariance matrix: %s' % cov)
 328     #_LOG.debug('info: %s' % info)
 329     _LOG.debug('message: %s' % mesg)
 330     if ier == 1:
 331         _LOG.debug('solution converged')
 332     else:
 333         _LOG.debug('solution did not converge')
 334     if plot or _package_config['matplotlib']:
 335         yguess = model(z, param_guess, high_voltage_rail=high_voltage_rail)
 336         yfit = model(z, p, high_voltage_rail=high_voltage_rail)
 337         _plot({'z': z, 'deflection': deflection}, yguess=yguess, yfit=yfit)
 338     return sensitivity_from_fit_params(p)
 339
 340 def save(filename=None, group='/', config=None, z_axis_config=None,
 341          deflection_channel_config=None, raw=None, processed=None):
 342     specs = [
 343         _SaveSpec(item=config, relpath='config/bump', config=_BumpConfig),
 344         _SaveSpec(item=z_axis_config, relpath='config/z', config=_AxisConfig),
 345         _SaveSpec(item=deflection_channel_config, relpath='config/deflection',
 346                   config=_InputChannelConfig),
 347         _SaveSpec(item=processed, relpath='processed', units='V/m'),
 348         ]
 349     if raw is not None:
 350         for key in raw.keys():
 351             specs.append(_SaveSpec(
 352                     item=raw[key], relpath='raw/{}'.format(key), array=True,
 353                     units='bits'))
 354     _save(filename=filename, group=group, specs=specs)
 355
 356 def load(filename=None, group='/'):
 357     specs = [
 358         _SaveSpec(key=('config', 'bump'), relpath='config/bump',
 359                   config=_BumpConfig),
 360         _SaveSpec(key=('config', 'z_axis_config'), relpath='config/z',
 361                   config=_AxisConfig),
 362         _SaveSpec(key=('config', 'deflection_channel_config'),
 363                   relpath='config/deflection', config=_InputChannelConfig),
 364         _SaveSpec(key=('raw', 'z'), relpath='raw/z', array=True, units='bits'),
 365         _SaveSpec(key=('raw', 'deflection'), relpath='raw/deflection',
 366                   array=True, units='bits'),
 367         _SaveSpec(key=('processed',), relpath='processed', units='V/m'),
 368         ]
 369     return _load(filename=filename, group=group, specs=specs)
 370
 371 def plot(data, yguess=None, yfit=None):
 372     "Plot the bump (Vphoto vs Vzp)"
 373     if not _matplotlib:
 374         raise _matplotlib_import_error
 375     figure = _matplotlib_pyplot.figure()
 376     if yfit is None:
 377         axes = figure.add_subplot(1, 1, 1)
 378     else:
 379         axes = figure.add_subplot(2, 1, 1)
 380         residual_axes = figure.add_subplot(2, 1, 2)
 381     timestamp = _time.strftime('%H%M%S')
 382
 383     axes.hold(True)
 384     axes.plot(data['z'], data['deflection'], '.', label='bump')
 385     if yguess != None:
 386         axes.plot(data['z'], yguess, 'g-', label='guess')
 387     if yfit != None:
 388         axes.plot(data['z'], yfit, 'r-', label='fit')
 389
 390     axes.set_title('bump surface %s' % timestamp)
 391     #axes.legend(loc='upper left')
 392     axes.set_xlabel('Z piezo (meters)')
 393     axes.set_ylabel('Photodiode (Volts)')
 394     if yfit is not None:
 395         # second subplot for residual
 396         residual_axes.plot(data['z'], data['deflection'] - yfit,
 397                            'r-', label='residual')
 398         #residual_axes.legend(loc='upper right')
 399         residual_axes.set_xlabel('Z piezo (meters)')
 400         residual_axes.set_ylabel('Photodiode (Volts)')
 401     if hasattr(figure, 'show'):
 402         figure.show()
 403 _plot = plot  # alternative name for use inside fit()