X-Git-Url: http://git.tremily.us/?p=calibcant.git;a=blobdiff_plain;f=calibcant%2Fbump_analyze.py;h=d5256fc84bd1ea660c05c81b6bfc73d188d35de9;hp=11a8f51e719f7f8b96220bec7fca47b4357ddb06;hb=560f9f94abbdf396b9f624e03bdb2a1f6d840bd1;hpb=71236f900c6b4276a7a93cd010f403213f8e2f18 diff --git a/calibcant/bump_analyze.py b/calibcant/bump_analyze.py old mode 100755 new mode 100644 index 11a8f51..d5256fc --- a/calibcant/bump_analyze.py +++ b/calibcant/bump_analyze.py @@ -1,261 +1,403 @@ -#!/usr/bin/python -# # calibcant - tools for thermally calibrating AFM cantilevers # -# Copyright (C) 2007,2008, William Trevor King +# Copyright (C) 2008-2012 W. Trevor King # -# 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 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 . -""" -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 + + 200 + + -5e-08 + + 10.0 + + quadratic + + 2e-07 + + 1e-06 + + 1024 + + 2.0 + /bump/processed + + 1.00... + + V/m + /bump/raw + /bump/raw/deflection + + [50 50 ... 50 51 52 ... 97 98 99] + + bits + /bump/raw/z + + [ 0 1 2 3 ... 97 98 99] + + bits + +>>> data = load(filename=filename, group='/bump/') + +>>> pprint(data) # doctest: +ELLIPSIS, +REPORT_UDIFF +{'config': {'bump': }, + '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 common # common module for the calibcant package -import config # config module for the calibcant package -import data_logger -import z_piezo_utils -import linfit -from splittable_kwargs import splittableKwargsFunction, \ - make_splittable_kwargs_function - -#@splittableKwargsFunction((bump_plot, 'data')) -# Some of the child functions aren't yet defined, so postpone -# make-splittable until later in the module. -def bump_analyze(data, zpGain=config.zpGain, - zpSensitivity=config.zpSensitivity, - Vzp_out2V=config.Vzp_out2V, - Vphoto_in2V=config.Vphoto_in2V, - **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). + """ + 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): """ - bump_plot_kwargs, = bump_analyze._splitargs(bump_analyze, kwargs) - scale_Vzp_bits2V = Vzp_out2V(1) - Vzp_out2V(0) - scale_Vphoto_bits2V = Vphoto_in2V(1) - Vphoto_in2V(0) - Vphoto2Vzp_out_bit, intercept = \ - linfit.linregress(x=data["Z piezo output"], - y=data["Deflection input"], - plotVerbose=False) - Vphoto2Vzp_out = Vphoto2Vzp_out_bit *scale_Vphoto_bits2V / scale_Vzp_bits2V - - # 1 / (Vzp/Vzp_out * Zp/Vzp * Zcant/Zp ) - Vzp_out2Zcant = 1.0/ (zpGain * zpSensitivity) # * 1 - bump_plot(data, **bump_plot_kwargs) - return Vphoto2Vzp_out * Vzp_out2Zcant - -@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, 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) - common._pylab.plot(data["Z piezo output"], data["Deflection input"], - '.', label='bump') - common._pylab.title("bump surface") - common._pylab.legend(loc='upper left') - common._flush_plot() - -make_splittable_kwargs_function(bump_analyze, (bump_plot, 'data')) - -@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].split('\n')[0] - 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 \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 straight line.\n' - 'Returns photodiode sensitivity Vphotodiode/Zcantilever in V/nm.\n' - ' should be whitespace-delimited, 2 column ASCII\n' - 'without a header line. e.g: "\\t\\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' - '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('-C', '--cut-contact', dest='cut', - help='bilinear fit to cut out contact region (currently only available in single-file mode)', - action='store_true', default=False) - 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('-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.tweakmode == False : - data = read_data(ifilename) - if options.cut : - ddict = {'approach':data} # although it may be any combination of approach and retract - try : - params = z_piezo_utils.analyzeSurfPosData(ddict, retAllParams=True) - a,b,c,d = tuple(params) # model : f(x) = x> sys.stderr, "fit with", params, ". using zp < %d" % c - data = remove_further_than(data, c) - except z_piezo_utils.poorFit, s : - # data not very bilinear, so don't cut anything off. - print >> sys.stderr, "use everything" - - photoSensitivity = bump_analyze(data) - - print >> ofile, photoSensitivity - else : # tweak file mode - slopes = bump_load_analyze_tweaked(ifilename) - if options.comma_out : - sep = ',' - else : - sep = '\n' - common.write_array(ofile, slopes, sep) - - if options.ofilename != None : - ofile.close() + Model the bump as: + flat region (off-surface) + linear region (in-contact) + flat region (high-voltage-rail) + Parameters: + x_contact (x value for the surface-contact kink) + y_contact (y value for the surface-contact kink) + slope (dy/dx at the surface-contact kink) + """ + x_contact,y_contact,slope = params + 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): + """ + Guess rough parameters for a limited_linear model. Assumes the + bump approaches (raising the deflection as it does so) first. + Retracting after the approach is optional. Approximates the contact + position and an on-surface (high) position by finding first crossings + of thresholds 0.3 and 0.7 of the y value's total range. Not the + most efficient algorithm, but it seems fairly robust. + """ + y_contact = float(y.min()) + y_max = float(y.max()) + 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: + i += 1 + i_low = i + 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) + +def limited_linear_sensitivity(params): + """ + Return the estimated sensitivity to small deflections according to + limited_linear fit parameters. + """ + slope = params[2] + return slope + +def limited_quadratic(x, params, high_voltage_rail): + """ + Model the bump as: + flat region (off-surface) + quadratic region (in-contact) + flat region (high-voltage-rail) + Parameters: + x_contact (x value for the surface-contact kink) + y_contact (y value for the surface-contact kink) + slope (dy/dx at the surface-contact kink) + quad (d**2 y / dx**2, allow decreasing sensitivity with increased x) + """ + x_contact,y_contact,slope,quad = params + 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): + """ + Guess rough parameters for a limited_quadratic model. Assumes the + bump approaches (raising the deflection as it does so) first. + Retracting after the approach is optional. Approximates the contact + position and an on-surface (high) position by finding first crossings + 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,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): + """ + Return the estimated sensitivity to small deflections according to + limited_quadratic fit parameters. + """ + slope = params[2] + return slope + +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) + +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: + 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()