Get calibcant working with the new load_from_config-based pyafm.
[calibcant.git] / calibcant / bump_analyze.py
old mode 100755 (executable)
new mode 100644 (file)
index 3783936..d5256fc
-#!/usr/bin/python
-#
 # calibcant - tools for thermally calibrating AFM cantilevers
 #
-# Copyright (C) 2007,2008, William Trevor King
+# Copyright (C) 2008-2012 W. Trevor King <wking@drexel.edu>
 #
-# This program is free software; you can redistribute it and/or
-# modify it under the terms of the GNU General Public License as
-# published by the Free Software Foundation; either version 3 of the
-# License, or (at your option) any later version.
+# This file is part of calibcant.
 #
-# This program is distributed in the hope that it will be useful, but
-# WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-# See the GNU General Public License for more details.
+# calibcant is free software: you can redistribute it and/or modify it under
+# the terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
 #
-# You should have received a copy of the GNU General Public License
-# along with this program; if not, write to the Free Software
-# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
-# 02111-1307, USA.
+# calibcant is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+# A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
 #
-# The author may be contacted at <wking@drexel.edu> on the Internet, or
-# write to Trevor King, Drexel University, Physics Dept., 3141 Chestnut St.,
-# Philadelphia PA 19104, USA.
+# You should have received a copy of the GNU General Public License along with
+# calibcant.  If not, see <http://www.gnu.org/licenses/>.
 
-"""
-Separate the more general bump_analyze() from the other bump_*()
-functions in calibcant.  Also provide a command line interface
-for analyzing data acquired through other workflows.
-
-The relevant physical quantities are :
- Vzp_out  Output z-piezo voltage (what we generate)
- Vzp      Applied z-piezo voltage (after external ZPGAIN)
- Zp       The z-piezo position
- Zcant    The cantilever vertical deflection
- Vphoto   The photodiode vertical deflection voltage (what we measure)
-
-Which are related by the parameters :
- zpGain           Vzp_out / Vzp
- zpSensitivity    Zp / Vzp
- photoSensitivity Vphoto / Zcant
-
-photoSensitivity is measured by bumping the cantilever against the
-surface, where Zp = Zcant (see calibrate.bump_aquire()).  The measured
-slope Vphoto/Vout is converted to photoSensitivity with bump_analyze().
+"""Surface bump analysis (measures photodiode sensitivity).
+
+Separate the more general `analyze()` from the other `bump_*()`
+functions in calibcant.
+
+The relevant physical quantities are:
+
+* `Vzp_out`  Output z-piezo voltage (what we generate)
+* `Vzp`      Applied z-piezo voltage (after external amplification)
+* `Zp`       The z-piezo position
+* `Zcant`    The cantilever vertical deflection
+* `Vphoto`   The photodiode vertical deflection voltage (what we measure)
+
+Which are related by the parameters:
+
+* `zp_gain`           Vzp_out / Vzp
+* `zp_sensitivity`    Zp / Vzp
+* `photo_sensitivity` Vphoto / Zcant
+
+`photo_sensitivity` is measured by bumping the cantilever against the
+surface, where `Zp = Zcant` (see `calibrate.bump_acquire()`).  The
+measured slope `Vphoto/Vout` is converted to `photo_sensitivity` with
+`analyze()`.
+
+>>> import os
+>>> from pprint import pprint
+>>> import tempfile
+>>> import numpy
+>>> from h5config.storage.hdf5 import pprint_HDF5
+>>> from pypiezo.config import ChannelConfig, AxisConfig
+>>> from .config import BumpConfig
+
+>>> fd,filename = tempfile.mkstemp(suffix='.h5', prefix='calibcant-')
+>>> os.close(fd)
+
+>>> config = BumpConfig()
+>>> z_channel_config = ChannelConfig()
+>>> z_channel_config['name'] = 'z'
+>>> z_channel_config['maxdata'] = 200
+>>> z_channel_config['conversion-coefficients'] = (0,1)
+>>> z_channel_config['conversion-origin'] = 0
+>>> z_axis_config = AxisConfig()
+>>> z_axis_config['channel'] = z_channel_config
+>>> deflection_channel_config = ChannelConfig()
+>>> deflection_channel_config['name'] = 'deflection'
+>>> deflection_channel_config['maxdata'] = 200
+>>> deflection_channel_config['conversion-coefficients'] = (0,1)
+>>> deflection_channel_config['conversion-origin'] = 0
+
+>>> raw = {
+...     'z': numpy.arange(100, dtype=numpy.uint16),
+...     'deflection': numpy.arange(100, dtype=numpy.uint16),
+...     }
+>>> raw['deflection'][:50] = 50
+>>> processed = analyze(
+...     config=config, data=raw, z_axis_config=z_axis_config,
+...     deflection_channel_config=deflection_channel_config)
+>>> plot(data=raw)  # TODO: convert to V and m
+>>> save(filename=filename, group='/bump/',
+...     config=config, raw=raw, processed=processed)
+
+>>> pprint_HDF5(filename)  # doctest: +ELLIPSIS, +REPORT_UDIFF
+/
+  /bump
+    /bump/config
+      /bump/config/bump
+        <HDF5 dataset "far-steps": shape (), type "<i4">
+          200
+        <HDF5 dataset "initial-position": shape (), type "<f8">
+          -5e-08
+        <HDF5 dataset "min-slope-ratio": shape (), type "<f8">
+          10.0
+        <HDF5 dataset "model": shape (), type "|S9">
+          quadratic
+        <HDF5 dataset "push-depth": shape (), type "<f8">
+          2e-07
+        <HDF5 dataset "push-speed": shape (), type "<f8">
+          1e-06
+        <HDF5 dataset "samples": shape (), type "<i4">
+          1024
+        <HDF5 dataset "setpoint": shape (), type "<f8">
+          2.0
+    /bump/processed
+      <HDF5 dataset "data": shape (), type "<f8">
+        1.00...
+      <HDF5 dataset "units": shape (), type "|S3">
+        V/m
+    /bump/raw
+      /bump/raw/deflection
+        <HDF5 dataset "data": shape (100,), type "<u2">
+          [50 50 ... 50 51 52 ... 97 98 99]
+        <HDF5 dataset "units": shape (), type "|S4">
+          bits
+      /bump/raw/z
+        <HDF5 dataset "data": shape (100,), type "<u2">
+          [ 0  1  2  3  ... 97 98 99]
+        <HDF5 dataset "units": shape (), type "|S4">
+          bits
+
+>>> data = load(filename=filename, group='/bump/')
+
+>>> pprint(data)  # doctest: +ELLIPSIS, +REPORT_UDIFF
+{'config': {'bump': <BumpConfig ...>},
+ 'processed': 1.00...,
+ 'raw': {'deflection': array([50, 50, ..., 52, 53, ..., 98, 99], dtype=uint16),
+         'z': array([ 0,  1,  2,  ..., 97, 98, 99], dtype=uint16)}}
+
+>>> os.remove(filename)
 """
 
-import numpy
-import scipy.optimize
-import common # common module for the calibcant package
-import config # config module for the calibcant package
-import data_logger
-from splittable_kwargs import splittableKwargsFunction, \
-    make_splittable_kwargs_function
-
-
-@splittableKwargsFunction()
-def Vzp_bits2nm(data_bits, zpGain=config.zpGain,
-                zpSensitivity=config.zpSensitivity,
-                Vzp_out2V=config.Vzp_out2V):
-    scale_Vzp_bits2V = Vzp_out2V(1) - Vzp_out2V(0)
-    data_V = data_bits / scale_Vzp_bits2V
-    #             bits / (bits/V) = V
-    data_nm = data_V * zpGain * zpSensitivity
-    return data_nm
-
-@splittableKwargsFunction()
-def Vphoto_bits2V(data_bits, Vphoto_in2V=config.Vphoto_in2V):
-    scale_Vphoto_bits2V = Vphoto_in2V(1) - Vphoto_in2V(0)
-    Vphoto_V = data_bits / scale_Vphoto_bits2V
-    #               bits / (bits/V) = V
-    return Vphoto_V
-
-@splittableKwargsFunction((Vzp_bits2nm, 'data_bits'),
-                          (Vphoto_bits2V, 'data_bits'))
-def slope_bitspbit2Vpnm(slope_bitspbit, **kwargs):
-    zp_kwargs,photo_kwargs = slope_bitspbit2Vpnm._splitargs(slope_bitspbit2Vpnm, kwargs)
-    Vzp_bits = 1.0
-    Vphoto_bits = slope_bitspbit * Vzp_bits
-    return Vphoto_bits2V(Vphoto_bits, **photo_kwargs)/Vzp_bits2nm(Vzp_bits, **zp_kwargs)
-    
-#@splittableKwargsFunction((bump_fit, 'zpiezo_output_bits',
-#                           'deflection_input_bits'),
-#                          (slope_bitspbit2Vpnm, 'slope_bitspbit'))
-# Some of the child functions aren't yet defined, so postpone
-# make-splittable until later in the module.
-def bump_analyze(data, **kwargs) :
-    """
-    Return the slope of the bump ;).
+import numpy as _numpy
+from scipy.optimize import leastsq as _leastsq
+
+try:
+    import matplotlib as _matplotlib
+    import matplotlib.pyplot as _matplotlib_pyplot
+    import time as _time  # for timestamping lines on plots
+except (ImportError, RuntimeError), e:
+    _matplotlib = None
+    _matplotlib_import_error = e
+
+from pypiezo.base import convert_bits_to_volts as _convert_bits_to_volts
+from pypiezo.base import convert_bits_to_meters as _convert_bits_to_meters
+from pypiezo.config import AxisConfig as _AxisConfig
+from pypiezo.config import InputChannelConfig as _InputChannelConfig
+
+from . import LOG as _LOG
+from . import package_config as _package_config
+from .config import Linear as _Linear
+from .config import Quadratic as _Quadratic
+from .config import BumpConfig as _BumpConfig
+from .util import SaveSpec as _SaveSpec
+from .util import save as _save
+from .util import load as _load
+
+
+def analyze(config, data, z_axis_config,
+            deflection_channel_config, plot=False):
+    """Return the slope of the bump.
+
     Inputs:
-      data        dictionary of data in DAC/ADC bits
-      Vzp_out2V   function that converts output DAC bits to Volts
-      Vphoto_in2V function that converts input ADC bits to Volts
-      zpGain      zpiezo applied voltage per output Volt
-      zpSensitivity  nm zpiezo response per applied Volt
+      data              dictionary of data in DAC/ADC bits
+      config            `.config._BumpConfig` instance
+      z_axis_config     z `pypiezo.config.AxisConfig` instance
+      deflection_channel_config
+                        deflection `pypiezo.config.InputChannelConfig` instance
+      plot              boolean overriding matplotlib config setting.
     Returns:
-     photoSensitivity (Vphoto/Zcant) in Volts/nm
-    Checks for strong correlation (r-value) and low randomness chance (p-value)
-    
-    With the current implementation, the data is regressed in DAC/ADC bits
-    and THEN converted, so we're assuming that both conversions are LINEAR.
-    If they aren't, rewrite to convert before the regression.
+      photo_sensitivity (Vphoto/Zcant) in Volts/m
+
+    Checks for strong correlation (r-value) and low randomness chance
+    (p-value).
     """
-    bump_fit_kwargs,slope_bitspbit2Vpnm_kwargs = \
-        bump_analyze._splitargs(bump_analyze, kwargs)
-    Vphoto2Vzp_out_bit = bump_fit(data['Z piezo output'],
-                                  data['Deflection input'],
-                                  **bump_fit_kwargs)
-    return slope_bitspbit2Vpnm(Vphoto2Vzp_out_bit, **slope_bitspbit2Vpnm_kwargs)
-
-def limited_linear(x, params):
+    z = _convert_bits_to_meters(z_axis_config, data['z'])
+    deflection = _convert_bits_to_volts(
+        deflection_channel_config, data['deflection'])
+    high_voltage_rail = _convert_bits_to_volts(
+        deflection_channel_config, deflection_channel_config['maxdata'])
+    if config['model'] == _Linear:
+        kwargs = {
+            'param_guesser': limited_linear_param_guess,
+            'model': limited_linear,
+            'sensitivity_from_fit_params': limited_linear_sensitivity,
+            }
+    else:  # _Quadratic
+        kwargs = {
+            'param_guesser': limited_quadratic_param_guess,
+            'model': limited_quadratic,
+            'sensitivity_from_fit_params': limited_quadratic_sensitivity,
+            }
+    photo_sensitivity = fit(
+        z, deflection, high_voltage_rail=high_voltage_rail, plot=plot,
+        **kwargs)
+    return photo_sensitivity
+
+def limited_linear(x, params, high_voltage_rail):
     """
     Model the bump as:
       flat region (off-surface)
@@ -119,13 +202,15 @@ def limited_linear(x, params):
       y_contact (y value for the surface-contact kink)
       slope (dy/dx at the surface-contact kink)
     """
-    high_voltage_rail = 2**16 - 1 # bits
     x_contact,y_contact,slope = params
-    y = slope*(x-x_contact) + y_contact
-    y = numpy.clip(y, y_contact, high_voltage_rail)
+    off_surface_mask = x <= x_contact
+    on_surface_mask = x > x_contact
+    y = (off_surface_mask * y_contact +
+         on_surface_mask * (y_contact + slope*(x-x_contact)))
+    y = _numpy.clip(y, y_contact, high_voltage_rail)
     return y
 
-def limited_linear_param_guess(x, y) :
+def limited_linear_param_guess(x, y):
     """
     Guess rough parameters for a limited_linear model.  Assumes the
     bump approaches (raising the deflection as it does so) first.
@@ -139,14 +224,16 @@ def limited_linear_param_guess(x, y) :
     i = 0
     y_low  = y_contact + 0.3 * (y_max-y_contact)
     y_high = y_contact + 0.7 * (y_max-y_contact)
-    while y[i] < y_low :
+    while y[i] < y_low:
         i += 1
     i_low = i
-    while y[i] < y_high :
+    while y[i] < y_high:
         i += 1
     i_high = i
     x_contact = float(x[i_low])
     x_high = float(x[i_high])
+    if x_high == x_contact:  # things must be pretty flat
+        x_contact = (x_contact + x[0]) / 2
     slope = (y_high - y_contact) / (x_high - x_contact)
     return (x_contact, y_contact, slope)
 
@@ -158,7 +245,7 @@ def limited_linear_sensitivity(params):
     slope = params[2]
     return slope
 
-def limited_quadratic(x, params):
+def limited_quadratic(x, params, high_voltage_rail):
     """
     Model the bump as:
       flat region (off-surface)
@@ -170,13 +257,16 @@ def limited_quadratic(x, params):
       slope (dy/dx at the surface-contact kink)
       quad (d**2 y / dx**2, allow decreasing sensitivity with increased x)
     """
-    high_voltage_rail = 2**16 - 1 # bits
     x_contact,y_contact,slope,quad = params
-    y = slope*(x-x_contact) + quad*(x-x_contact)**2+ y_contact
-    y = numpy.clip(y, y_contact, high_voltage_rail)
+    off_surface_mask = x <= x_contact
+    on_surface_mask = x > x_contact
+    y = (off_surface_mask * y_contact +
+         on_surface_mask * (
+            y_contact + slope*(x-x_contact) + quad*(x-x_contact)**2))
+    y = _numpy.clip(y, y_contact, high_voltage_rail)
     return y
 
-def limited_quadratic_param_guess(x, y) :
+def limited_quadratic_param_guess(x, y):
     """
     Guess rough parameters for a limited_quadratic model.  Assumes the
     bump approaches (raising the deflection as it does so) first.
@@ -185,8 +275,13 @@ def limited_quadratic_param_guess(x, y) :
     of thresholds 0.3 and 0.7 of the y value's total range.  Not the
     most efficient algorithm, but it seems fairly robust.
     """
-    x_contact,y_contact,slope = limited_linear_param_guess(x,y)
-    quad = 0
+    x_contact,y_contact,linear_slope = limited_linear_param_guess(x,y)
+    contact_depth = x.max() - x_contact
+    slope = linear_slope / 2
+    quad = slope / contact_depth
+    # The above slope and quad were chosen so
+    #   x = contact_depth
+    #   x*slope + x**2 * slope == x * linear_slope
     return (x_contact, y_contact, slope, quad)
 
 def limited_quadratic_sensitivity(params):
@@ -197,239 +292,112 @@ def limited_quadratic_sensitivity(params):
     slope = params[2]
     return slope
 
-@splittableKwargsFunction()
-def bump_fit(zpiezo_output_bits, deflection_input_bits,
-             param_guesser=limited_quadratic_param_guess,
-             model=limited_quadratic,
-             sensitivity_from_fit_params=limited_quadratic_sensitivity,
-             plotVerbose=False) :
-    x = zpiezo_output_bits
-    y = deflection_input_bits
-    def residual(p, y, x) :
-        return model(x, p) - y
-    param_guess = param_guesser(x, y)
-    p,cov,info,mesg,ier = \
-        scipy.optimize.leastsq(residual, param_guess, args=(y, x),
-                               full_output=True, maxfev=int(10e3))
-    if config.TEXT_VERBOSE :
-        print "Fitted params:",p
-        print "Covariance mx:",cov
-        print "Info:", info
-        print "mesg:", mesg
-        if ier == 1 :
-            print "Solution converged"
-        else :
-            print "Solution did not converge"
-    if plotVerbose or config.PYLAB_VERBOSE :
-        yguess = model(x, param_guess)
-        #yguess = None # Don't print the guess, since I'm convinced it's ok ;).
-        yfit = model(x, p)
-        bump_plot(data={"Z piezo output":x, "Deflection input":y},
-                  yguess=yguess, yfit=yfit, plotVerbose=plotVerbose)
+def fit(z, deflection, high_voltage_rail,
+        param_guesser=limited_quadratic_param_guess,
+        model=limited_quadratic,
+        sensitivity_from_fit_params=limited_quadratic_sensitivity,
+        plot=False):
+    """Fit a aurface bump and return the photodiode sensitivitiy.
+
+    Parameters:
+      z              piezo position in meters
+      deflection     photodiode deflection in volts
+      param_guesser  function that guesses initial model parameters
+      model          parametric model of deflection vs. z
+      sensitivity_from_fit_params  given fit params, return the sensitivity
+      plot              boolean overriding matplotlib config setting.
+    Returns:
+      photodiode_sensitivity  photodiode volts per piezo meter
+    """
+    _LOG.debug('fit bump data with model %s' % model)
+    def residual(p, deflection, z):
+        return model(z, p, high_voltage_rail=high_voltage_rail) - deflection
+    param_guess = param_guesser(z, deflection)
+    try:
+        p,cov,info,mesg,ier = _leastsq(
+            residual, param_guess, args=(deflection, z), full_output=True,
+            maxfev=int(10e3))
+    except ValueError:
+        zd = _numpy.ndarray(list(z.shape) + [2], dtype=z.dtype)
+        zd[:,0] = z
+        zd[:,1] = d
+        _numpy.savetxt('/tmp/z-deflection.dat', zd, delimiter='\t')
+        raise
+    _LOG.debug('fitted params: %s' % p)
+    _LOG.debug('covariance matrix: %s' % cov)
+    #_LOG.debug('info: %s' % info)
+    _LOG.debug('message: %s' % mesg)
+    if ier == 1:
+        _LOG.debug('solution converged')
+    else:
+        _LOG.debug('solution did not converge')
+    if plot or _package_config['matplotlib']:
+        yguess = model(z, param_guess, high_voltage_rail=high_voltage_rail)
+        yfit = model(z, p, high_voltage_rail=high_voltage_rail)
+        _plot({'z': z, 'deflection': deflection}, yguess=yguess, yfit=yfit)
     return sensitivity_from_fit_params(p)
 
-@splittableKwargsFunction()
-def bump_save(data, log_dir=None) :
-    "Save the dictionary data, using data_logger.data_log()"
-    if log_dir != None :
-        log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
-                                   log_name="bump")
-        log.write_dict_of_arrays(data)
-
-def bump_load(datafile) :
-    "Load the dictionary data, using data_logger.date_load()"
-    dl = data_logger.data_load()
-    data = dl.read_dict_of_arrays(datafile)
-    return data
-
-@splittableKwargsFunction()
-def bump_plot(data, yguess=None, yfit=None, plotVerbose=False) :
-    "Plot the bump (Vphoto vs Vzp) if plotVerbose or PYLAB_VERBOSE == True"
-    if plotVerbose or config.PYLAB_VERBOSE :
-        common._import_pylab()
-        common._pylab.figure(config.BASE_FIGNUM)
-        if yfit != None: # two subplot figure
-            common._pylab.subplot(211)
-        common._pylab.hold(False)
-        common._pylab.plot(data["Z piezo output"], data["Deflection input"],
-                           '.', label='bump')
-        common._pylab.hold(True)
-        if yguess != None:
-            common._pylab.plot(data["Z piezo output"], yguess,
-                               'g-', label='guess')
-        if yfit != None:
-            common._pylab.plot(data["Z piezo output"], yfit,
-                               'r-', label='fit')
-        common._pylab.hold(False)
-        common._pylab.title("bump surface")
-        common._pylab.legend(loc='upper left')
-        common._pylab.xlabel("Z piezo output voltage (bits)")
-        common._pylab.ylabel("Photodiode input voltage (bits)")
-        if yfit != None:
-            # second subplot for residual
-            common._pylab.subplot(212)
-            common._pylab.plot(data["Z piezo output"],
-                               data["Deflection input"] - yfit,
-                               'r-', label='residual')
-            common._pylab.legend(loc='upper right')
-            common._pylab.xlabel("Z piezo output voltage (bits)")
-            common._pylab.ylabel("Photodiode input voltage (bits)")
-        common._flush_plot()
-
-make_splittable_kwargs_function(bump_analyze,
-                                (bump_fit, 'zpiezo_output_bits',
-                                 'deflection_input_bits'),
-                                (slope_bitspbit2Vpnm, 'slope_bitspbit'))
-
-@splittableKwargsFunction((bump_analyze, 'data'))
-def bump_load_analyze_tweaked(tweak_file, **kwargs):
-    "Load the output file of tweak_calib_bump.sh, return an array of slopes"
-    bump_analyze_kwargs, = \
-        bump_load_analyze_tweaked._splitargs(bump_load_analyze_tweaked, kwargs)
-    photoSensitivity = []
-    for line in file(tweak_file, 'r') :
-        parsed = line.split()
-        path = parsed[0].strip()
-        if path[0] == '#' : # a comment
-            continue
-        if config.TEXT_VERBOSE :
-            print "Reading data from %s with ranges %s" % (path, parsed[1:])
-        # read the data
-        full_data = bump_load(path)
-        if len(parsed) == 1 :
-            data = full_data # use whole bump
-        else :
-            # use the listed sections
-            zp = []
-            df = []
-            for rng in parsed[1:] :
-                p = rng.split(':')
-                starti = int(p[0])
-                stopi = int(p[1])
-                zp.extend(full_data['Z piezo output'][starti:stopi])
-                df.extend(full_data['Deflection input'][starti:stopi])
-            data = {'Z piezo output': numpy.array(zp),
-                    'Deflection input': numpy.array(df)}
-        pSi = bump_analyze(data, **bump_analyze_kwargs)
-        photoSensitivity.append(pSi)
-    return numpy.array(photoSensitivity, dtype=numpy.float)
-
-# commandline interface functions
-import scipy.io, sys
-
-def read_data(ifile):
-    "ifile can be a filename string or open (seekable) file object"
-    if ifile == None :  ifile = sys.stdin
-    unlabeled_data=scipy.io.read_array(ifile)
-    data = {}
-    data['Z piezo output'] = unlabeled_data[:,0]
-    data['Deflection input'] = unlabeled_data[:,1]
-    return data
-
-def remove_further_than(data, zp_crit) :
-    ndata = {}
-    ndata['Z piezo output'] = []
-    ndata['Deflection input'] = []
-    for zp,df in zip(data['Z piezo output'],data['Deflection input']) :
-        if zp > zp_crit :
-            ndata['Z piezo output'].append(zp)
-            ndata['Deflection input'].append(df)
-    return ndata
-
-if __name__ == '__main__' :
-    # command line interface
-    from optparse import OptionParser
-    
-    usage_string = ('%prog <input-file>\n'
-                    '2008, W. Trevor King.\n'
-                    '\n'
-                    'There are two operation modes, one to analyze a single bump file,\n'
-                    'and one to analyze tweak files.\n'
-                    '\n'
-                    'Single file mode (the default) :\n'
-                    'Scales raw DAC/ADC bit data and fits a bounded quadratic.\n'
-                    'Returns photodiode sensitivity Vphotodiode/Zcantilever in V/nm, determined by.\n'
-                    'the slope at the kink between the non-contact region and the contact region.\n'
-                    '<input-file> should be whitespace-delimited, 2 column ASCII\n'
-                    'without a header line.  e.g: "<zp_DAC>\\t<deflection_ADC>\\n"\n'
-                    '\n'
-                    'Tweak file mode:\n'
-                    'Runs the same analysis as in single file mode for each bump in\n'
-                    'a tweak file.  Each line in the tweak file specifies a single bump.\n'
-                    'Blank lines and those beginning with a pound sign (#) are ignored.\n'
-                    'The format of a line is a series of whitespace-separated fields--\n'
-                    'a base file path followed by optional point index ranges, e.g.:\n'
-                    '20080919/20080919132500_bump_surface 10:651 1413:2047\n'
-                    'which only discards all points outside the index ranges [10,651)\n'
-                    'and [1413,2047) (indexing starts at 0).\n'
-                    )
-    parser = OptionParser(usage=usage_string, version='%prog '+common.VERSION)
-    parser.add_option('-o', '--output-file', dest='ofilename',
-                      help='write output to FILE (default stdout)',
-                      type='string', metavar='FILE')
-    parser.add_option('-c', '--comma-out', dest='comma_out', action='store_true',
-                      help='Output comma-seperated values (default %default)',
-                      default=False)
-    parser.add_option('-p', '--pylab', dest='pylab', action='store_true',
-                      help='Produce pylab fit checks during execution',
-                      default=False)
-    parser.add_option('-t', '--tweak-mode', dest='tweakmode', action='store_true',
-                      help='Run in tweak-file mode',
-                      default=False)
-    parser.add_option('-d', '--datalogger-mode', dest='datalogger_mode', action='store_true',
-                      help='Run input files with datalogger.read_dict_of_arrays().  This is useful, for example, to test a single line from a tweakfile.',
-                      default=False)
-    parser.add_option('-q', '--disable-quadratic', dest='quadratic', action='store_false',
-                      help='Disable quadratic term in fitting (i.e. use bounded linear fits).',
-                      default=True)
-    parser.add_option('-v', '--verbose', dest='verbose', action='store_true',
-                      help='Print lots of debugging information',
-                      default=False)
-
-    options,args = parser.parse_args()
-    parser.destroy()
-    assert len(args) >= 1, "Need an input file"
-        
-    ifilename = args[0]
-
-    if options.ofilename != None :
-        ofile = file(options.ofilename, 'w')
-    else :
-        ofile = sys.stdout
-    config.TEXT_VERBOSE = options.verbose
-    config.PYLAB_INTERACTIVE = False
-    config.PYLAB_VERBOSE = options.pylab
-    config.GNUPLOT_VERBOSE = False
-    if options.quadratic == True:
-        param_guesser = limited_quadratic_param_guess
-        model = limited_quadratic
-        sensitivity_from_fit_params = limited_quadratic_sensitivity
+def save(filename=None, group='/', config=None, z_axis_config=None,
+         deflection_channel_config=None, raw=None, processed=None):
+    specs = [
+        _SaveSpec(item=config, relpath='config/bump', config=_BumpConfig),
+        _SaveSpec(item=z_axis_config, relpath='config/z', config=_AxisConfig),
+        _SaveSpec(item=deflection_channel_config, relpath='config/deflection',
+                  config=_InputChannelConfig),
+        _SaveSpec(item=processed, relpath='processed', units='V/m'),
+        ]
+    if raw is not None:
+        for key in raw.keys():
+            specs.append(_SaveSpec(
+                    item=raw[key], relpath='raw/{}'.format(key), array=True,
+                    units='bits'))
+    _save(filename=filename, group=group, specs=specs)
+
+def load(filename=None, group='/'):
+    specs = [
+        _SaveSpec(key=('config', 'bump'), relpath='config/bump',
+                  config=_BumpConfig),
+        _SaveSpec(key=('config', 'z_axis_config'), relpath='config/z',
+                  config=_AxisConfig),
+        _SaveSpec(key=('config', 'deflection_channel_config'),
+                  relpath='config/deflection', config=_InputChannelConfig),
+        _SaveSpec(key=('raw', 'z'), relpath='raw/z', array=True, units='bits'),
+        _SaveSpec(key=('raw', 'deflection'), relpath='raw/deflection',
+                  array=True, units='bits'),
+        _SaveSpec(key=('processed',), relpath='processed', units='V/m'),
+        ]
+    return _load(filename=filename, group=group, specs=specs)
+
+def plot(data, yguess=None, yfit=None):
+    "Plot the bump (Vphoto vs Vzp)"
+    if not _matplotlib:
+        raise _matplotlib_import_error
+    figure = _matplotlib_pyplot.figure()
+    if yfit is None:
+        axes = figure.add_subplot(1, 1, 1)
     else:
-        param_guesser = limited_linear_param_guess
-        model = limited_linear
-        sensitivity_from_fit_params = limited_linear_sensitivity
-    
-    if options.tweakmode == False :
-        if options.datalogger_mode:
-            data = bump_load(ifilename)
-        else:
-            data = read_data(ifilename)
-        photoSensitivity = bump_analyze(data,
-                                        param_guesser=param_guesser,
-                                        model=model,
-                                        sensitivity_from_fit_params=sensitivity_from_fit_params)
-        
-        print >> ofile, photoSensitivity
-    else : # tweak file mode
-        slopes = bump_load_analyze_tweaked(ifilename,
-                                           param_guesser=param_guesser,
-                                           model=model,
-                                           sensitivity_from_fit_params=sensitivity_from_fit_params)
-        if options.comma_out :
-            sep = ','
-        else :
-            sep = '\n'
-        common.write_array(ofile, slopes, sep)
-    
-    if options.ofilename != None :
-        ofile.close()
+        axes = figure.add_subplot(2, 1, 1)
+        residual_axes = figure.add_subplot(2, 1, 2)
+    timestamp = _time.strftime('%H%M%S')
+
+    axes.hold(True)
+    axes.plot(data['z'], data['deflection'], '.', label='bump')
+    if yguess != None:
+        axes.plot(data['z'], yguess, 'g-', label='guess')
+    if yfit != None:
+        axes.plot(data['z'], yfit, 'r-', label='fit')
+
+    axes.set_title('bump surface %s' % timestamp)
+    #axes.legend(loc='upper left')
+    axes.set_xlabel('Z piezo (meters)')
+    axes.set_ylabel('Photodiode (Volts)')
+    if yfit is not None:
+        # second subplot for residual
+        residual_axes.plot(data['z'], data['deflection'] - yfit,
+                           'r-', label='residual')
+        #residual_axes.legend(loc='upper right')
+        residual_axes.set_xlabel('Z piezo (meters)')
+        residual_axes.set_ylabel('Photodiode (Volts)')
+    if hasattr(figure, 'show'):
+        figure.show()
+_plot = plot  # alternative name for use inside fit()