calibcant/bump_analyze.py

   1 #!/usr/bin/python
   2 #
   3 # calibcant - tools for thermally calibrating AFM cantilevers
   4 #
   5 # Copyright (C) 2008-2010 W. Trevor King <wking@drexel.edu>
   6 #
   7 # This file is part of CalibCant.
   8 #
   9 # CalibCant is free software: you can redistribute it and/or
  10 # modify it under the terms of the GNU Lesser General Public
  11 # License as published by the Free Software Foundation, either
  12 # version 3 of the License, or (at your option) any later version.
  13 #
  14 # CalibCant is distributed in the hope that it will be useful,
  15 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  16 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  17 # GNU Lesser General Public License for more details.
  18 #
  19 # You should have received a copy of the GNU Lesser General Public
  20 # License along with CalibCant.  If not, see
  21 # <http://www.gnu.org/licenses/>.
  22
  23 """
  24 Separate the more general bump_analyze() from the other bump_*()
  25 functions in calibcant.  Also provide a command line interface
  26 for analyzing data acquired through other workflows.
  27
  28 The relevant physical quantities are :
  29  Vzp_out  Output z-piezo voltage (what we generate)
  30  Vzp      Applied z-piezo voltage (after external ZPGAIN)
  31  Zp       The z-piezo position
  32  Zcant    The cantilever vertical deflection
  33  Vphoto   The photodiode vertical deflection voltage (what we measure)
  34
  35 Which are related by the parameters :
  36  zpGain           Vzp_out / Vzp
  37  zpSensitivity    Zp / Vzp
  38  photoSensitivity Vphoto / Zcant
  39
  40 photoSensitivity is measured by bumping the cantilever against the
  41 surface, where Zp = Zcant (see calibrate.bump_aquire()).  The measured
  42 slope Vphoto/Vout is converted to photoSensitivity with bump_analyze().
  43 """
  44
  45 import numpy
  46 import scipy.optimize
  47
  48 import data_logger
  49 from splittable_kwargs import splittableKwargsFunction, \
  50     make_splittable_kwargs_function
  51
  52 from . import common
  53 from . import config
  54
  55
  56 @splittableKwargsFunction()
  57 def Vzp_bits2nm(data_bits, zpGain=config.zpGain,
  58                 zpSensitivity=config.zpSensitivity,
  59                 Vzp_out2V=config.Vzp_out2V):
  60     scale_Vzp_bits2V = Vzp_out2V(1) - Vzp_out2V(0)
  61     data_V = data_bits / scale_Vzp_bits2V
  62     #             bits / (bits/V) = V
  63     data_nm = data_V * zpGain * zpSensitivity
  64     return data_nm
  65
  66 @splittableKwargsFunction()
  67 def Vphoto_bits2V(data_bits, Vphoto_in2V=config.Vphoto_in2V):
  68     scale_Vphoto_bits2V = Vphoto_in2V(1) - Vphoto_in2V(0)
  69     Vphoto_V = data_bits / scale_Vphoto_bits2V
  70     #               bits / (bits/V) = V
  71     return Vphoto_V
  72
  73 @splittableKwargsFunction((Vzp_bits2nm, 'data_bits'),
  74                           (Vphoto_bits2V, 'data_bits'))
  75 def slope_bitspbit2Vpnm(slope_bitspbit, **kwargs):
  76     zp_kwargs,photo_kwargs = slope_bitspbit2Vpnm._splitargs(slope_bitspbit2Vpnm, kwargs)
  77     Vzp_bits = 1.0
  78     Vphoto_bits = slope_bitspbit * Vzp_bits
  79     return Vphoto_bits2V(Vphoto_bits, **photo_kwargs)/Vzp_bits2nm(Vzp_bits, **zp_kwargs)
  80
  81 #@splittableKwargsFunction((bump_fit, 'zpiezo_output_bits',
  82 #                           'deflection_input_bits'),
  83 #                          (slope_bitspbit2Vpnm, 'slope_bitspbit'))
  84 # Some of the child functions aren't yet defined, so postpone
  85 # make-splittable until later in the module.
  86 def bump_analyze(data, **kwargs) :
  87     """
  88     Return the slope of the bump ;).
  89     Inputs:
  90       data        dictionary of data in DAC/ADC bits
  91       Vzp_out2V   function that converts output DAC bits to Volts
  92       Vphoto_in2V function that converts input ADC bits to Volts
  93       zpGain      zpiezo applied voltage per output Volt
  94       zpSensitivity  nm zpiezo response per applied Volt
  95     Returns:
  96      photoSensitivity (Vphoto/Zcant) in Volts/nm
  97     Checks for strong correlation (r-value) and low randomness chance (p-value)
  98
  99     With the current implementation, the data is regressed in DAC/ADC bits
 100     and THEN converted, so we're assuming that both conversions are LINEAR.
 101     If they aren't, rewrite to convert before the regression.
 102     """
 103     bump_fit_kwargs,slope_bitspbit2Vpnm_kwargs = \
 104         bump_analyze._splitargs(bump_analyze, kwargs)
 105     Vphoto2Vzp_out_bit = bump_fit(data['Z piezo output'],
 106                                   data['Deflection input'],
 107                                   **bump_fit_kwargs)
 108     return slope_bitspbit2Vpnm(Vphoto2Vzp_out_bit, **slope_bitspbit2Vpnm_kwargs)
 109
 110 def limited_linear(x, params):
 111     """
 112     Model the bump as:
 113       flat region (off-surface)
 114       linear region (in-contact)
 115       flat region (high-voltage-rail)
 116     Parameters:
 117       x_contact (x value for the surface-contact kink)
 118       y_contact (y value for the surface-contact kink)
 119       slope (dy/dx at the surface-contact kink)
 120     """
 121     high_voltage_rail = 2**16 - 1 # bits
 122     x_contact,y_contact,slope = params
 123     y = slope*(x-x_contact) + y_contact
 124     y = numpy.clip(y, y_contact, high_voltage_rail)
 125     return y
 126
 127 def limited_linear_param_guess(x, y) :
 128     """
 129     Guess rough parameters for a limited_linear model.  Assumes the
 130     bump approaches (raising the deflection as it does so) first.
 131     Retracting after the approach is optional.  Approximates the contact
 132     position and an on-surface (high) position by finding first crossings
 133     of thresholds 0.3 and 0.7 of the y value's total range.  Not the
 134     most efficient algorithm, but it seems fairly robust.
 135     """
 136     y_contact = float(y.min())
 137     y_max = float(y.max())
 138     i = 0
 139     y_low  = y_contact + 0.3 * (y_max-y_contact)
 140     y_high = y_contact + 0.7 * (y_max-y_contact)
 141     while y[i] < y_low :
 142         i += 1
 143     i_low = i
 144     while y[i] < y_high :
 145         i += 1
 146     i_high = i
 147     x_contact = float(x[i_low])
 148     x_high = float(x[i_high])
 149     slope = (y_high - y_contact) / (x_high - x_contact)
 150     return (x_contact, y_contact, slope)
 151
 152 def limited_linear_sensitivity(params):
 153     """
 154     Return the estimated sensitivity to small deflections according to
 155     limited_linear fit parameters.
 156     """
 157     slope = params[2]
 158     return slope
 159
 160 def limited_quadratic(x, params):
 161     """
 162     Model the bump as:
 163       flat region (off-surface)
 164       quadratic region (in-contact)
 165       flat region (high-voltage-rail)
 166     Parameters:
 167       x_contact (x value for the surface-contact kink)
 168       y_contact (y value for the surface-contact kink)
 169       slope (dy/dx at the surface-contact kink)
 170       quad (d**2 y / dx**2, allow decreasing sensitivity with increased x)
 171     """
 172     high_voltage_rail = 2**16 - 1 # bits
 173     x_contact,y_contact,slope,quad = params
 174     y = slope*(x-x_contact) + quad*(x-x_contact)**2+ y_contact
 175     y = numpy.clip(y, y_contact, high_voltage_rail)
 176     return y
 177
 178 def limited_quadratic_param_guess(x, y) :
 179     """
 180     Guess rough parameters for a limited_quadratic model.  Assumes the
 181     bump approaches (raising the deflection as it does so) first.
 182     Retracting after the approach is optional.  Approximates the contact
 183     position and an on-surface (high) position by finding first crossings
 184     of thresholds 0.3 and 0.7 of the y value's total range.  Not the
 185     most efficient algorithm, but it seems fairly robust.
 186     """
 187     x_contact,y_contact,slope = limited_linear_param_guess(x,y)
 188     quad = 0
 189     return (x_contact, y_contact, slope, quad)
 190
 191 def limited_quadratic_sensitivity(params):
 192     """
 193     Return the estimated sensitivity to small deflections according to
 194     limited_quadratic fit parameters.
 195     """
 196     slope = params[2]
 197     return slope
 198
 199 @splittableKwargsFunction()
 200 def bump_fit(zpiezo_output_bits, deflection_input_bits,
 201              param_guesser=limited_quadratic_param_guess,
 202              model=limited_quadratic,
 203              sensitivity_from_fit_params=limited_quadratic_sensitivity,
 204              plotVerbose=False) :
 205     x = zpiezo_output_bits
 206     y = deflection_input_bits
 207     def residual(p, y, x) :
 208         return model(x, p) - y
 209     param_guess = param_guesser(x, y)
 210     p,cov,info,mesg,ier = \
 211         scipy.optimize.leastsq(residual, param_guess, args=(y, x),
 212                                full_output=True, maxfev=int(10e3))
 213     if config.TEXT_VERBOSE :
 214         print "Fitted params:",p
 215         print "Covariance mx:",cov
 216         print "Info:", info
 217         print "mesg:", mesg
 218         if ier == 1 :
 219             print "Solution converged"
 220         else :
 221             print "Solution did not converge"
 222     if plotVerbose or config.PYLAB_VERBOSE :
 223         yguess = model(x, param_guess)
 224         #yguess = None # Don't print the guess, since I'm convinced it's ok ;).
 225         yfit = model(x, p)
 226         bump_plot(data={"Z piezo output":x, "Deflection input":y},
 227                   yguess=yguess, yfit=yfit, plotVerbose=plotVerbose)
 228     return sensitivity_from_fit_params(p)
 229
 230 @splittableKwargsFunction()
 231 def bump_save(data, log_dir=None) :
 232     "Save the dictionary data, using data_logger.data_log()"
 233     if log_dir != None :
 234         log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
 235                                    log_name="bump")
 236         log.write_dict_of_arrays(data)
 237
 238 def bump_load(datafile) :
 239     "Load the dictionary data, using data_logger.date_load()"
 240     dl = data_logger.data_load()
 241     data = dl.read_dict_of_arrays(datafile)
 242     return data
 243
 244 @splittableKwargsFunction()
 245 def bump_plot(data, yguess=None, yfit=None, plotVerbose=False) :
 246     "Plot the bump (Vphoto vs Vzp) if plotVerbose or PYLAB_VERBOSE == True"
 247     if plotVerbose or config.PYLAB_VERBOSE :
 248         common._import_pylab()
 249         common._pylab.figure(config.BASE_FIGNUM)
 250         if yfit != None: # two subplot figure
 251             common._pylab.subplot(211)
 252         common._pylab.hold(False)
 253         common._pylab.plot(data["Z piezo output"], data["Deflection input"],
 254                            '.', label='bump')
 255         common._pylab.hold(True)
 256         if yguess != None:
 257             common._pylab.plot(data["Z piezo output"], yguess,
 258                                'g-', label='guess')
 259         if yfit != None:
 260             common._pylab.plot(data["Z piezo output"], yfit,
 261                                'r-', label='fit')
 262         common._pylab.hold(False)
 263         common._pylab.title("bump surface")
 264         common._pylab.legend(loc='upper left')
 265         common._pylab.xlabel("Z piezo output voltage (bits)")
 266         common._pylab.ylabel("Photodiode input voltage (bits)")
 267         if yfit != None:
 268             # second subplot for residual
 269             common._pylab.subplot(212)
 270             common._pylab.plot(data["Z piezo output"],
 271                                data["Deflection input"] - yfit,
 272                                'r-', label='residual')
 273             common._pylab.legend(loc='upper right')
 274             common._pylab.xlabel("Z piezo output voltage (bits)")
 275             common._pylab.ylabel("Photodiode input voltage (bits)")
 276         common._flush_plot()
 277
 278 make_splittable_kwargs_function(bump_analyze,
 279                                 (bump_fit, 'zpiezo_output_bits',
 280                                  'deflection_input_bits'),
 281                                 (slope_bitspbit2Vpnm, 'slope_bitspbit'))
 282
 283 @splittableKwargsFunction((bump_analyze, 'data'))
 284 def bump_load_analyze_tweaked(tweak_file, **kwargs):
 285     "Load the output file of tweak_calib_bump.sh, return an array of slopes"
 286     bump_analyze_kwargs, = \
 287         bump_load_analyze_tweaked._splitargs(bump_load_analyze_tweaked, kwargs)
 288     photoSensitivity = []
 289     for line in file(tweak_file, 'r') :
 290         parsed = line.split()
 291         path = parsed[0].strip()
 292         if path[0] == '#' : # a comment
 293             continue
 294         if config.TEXT_VERBOSE :
 295             print "Reading data from %s with ranges %s" % (path, parsed[1:])
 296         # read the data
 297         full_data = bump_load(path)
 298         if len(parsed) == 1 :
 299             data = full_data # use whole bump
 300         else :
 301             # use the listed sections
 302             zp = []
 303             df = []
 304             for rng in parsed[1:] :
 305                 p = rng.split(':')
 306                 starti = int(p[0])
 307                 stopi = int(p[1])
 308                 zp.extend(full_data['Z piezo output'][starti:stopi])
 309                 df.extend(full_data['Deflection input'][starti:stopi])
 310             data = {'Z piezo output': numpy.array(zp),
 311                     'Deflection input': numpy.array(df)}
 312         pSi = bump_analyze(data, **bump_analyze_kwargs)
 313         photoSensitivity.append(pSi)
 314     return numpy.array(photoSensitivity, dtype=numpy.float)
 315
 316 # commandline interface functions
 317 import scipy.io, sys
 318
 319 def read_data(ifile):
 320     "ifile can be a filename string or open (seekable) file object"
 321     if ifile == None :  ifile = sys.stdin
 322     unlabeled_data=scipy.io.read_array(ifile)
 323     data = {}
 324     data['Z piezo output'] = unlabeled_data[:,0]
 325     data['Deflection input'] = unlabeled_data[:,1]
 326     return data
 327
 328 def remove_further_than(data, zp_crit) :
 329     ndata = {}
 330     ndata['Z piezo output'] = []
 331     ndata['Deflection input'] = []
 332     for zp,df in zip(data['Z piezo output'],data['Deflection input']) :
 333         if zp > zp_crit :
 334             ndata['Z piezo output'].append(zp)
 335             ndata['Deflection input'].append(df)
 336     return ndata
 337
 338 if __name__ == '__main__' :
 339     # command line interface
 340     from optparse import OptionParser
 341
 342     usage_string = ('%prog <input-file>\n'
 343                     '2008, W. Trevor King.\n'
 344                     '\n'
 345                     'There are two operation modes, one to analyze a single bump file,\n'
 346                     'and one to analyze tweak files.\n'
 347                     '\n'
 348                     'Single file mode (the default) :\n'
 349                     'Scales raw DAC/ADC bit data and fits a bounded quadratic.\n'
 350                     'Returns photodiode sensitivity Vphotodiode/Zcantilever in V/nm, determined by.\n'
 351                     'the slope at the kink between the non-contact region and the contact region.\n'
 352                     '<input-file> should be whitespace-delimited, 2 column ASCII\n'
 353                     'without a header line.  e.g: "<zp_DAC>\\t<deflection_ADC>\\n"\n'
 354                     '\n'
 355                     'Tweak file mode:\n'
 356                     'Runs the same analysis as in single file mode for each bump in\n'
 357                     'a tweak file.  Each line in the tweak file specifies a single bump.\n'
 358                     'Blank lines and those beginning with a pound sign (#) are ignored.\n'
 359                     'The format of a line is a series of whitespace-separated fields--\n'
 360                     'a base file path followed by optional point index ranges, e.g.:\n'
 361                     '20080919/20080919132500_bump_surface 10:651 1413:2047\n'
 362                     'which only discards all points outside the index ranges [10,651)\n'
 363                     'and [1413,2047) (indexing starts at 0).\n'
 364                     )
 365     parser = OptionParser(usage=usage_string, version='%prog '+common.VERSION)
 366     parser.add_option('-o', '--output-file', dest='ofilename',
 367                       help='write output to FILE (default stdout)',
 368                       type='string', metavar='FILE')
 369     parser.add_option('-c', '--comma-out', dest='comma_out', action='store_true',
 370                       help='Output comma-seperated values (default %default)',
 371                       default=False)
 372     parser.add_option('-p', '--pylab', dest='pylab', action='store_true',
 373                       help='Produce pylab fit checks during execution',
 374                       default=False)
 375     parser.add_option('-t', '--tweak-mode', dest='tweakmode', action='store_true',
 376                       help='Run in tweak-file mode',
 377                       default=False)
 378     parser.add_option('-d', '--datalogger-mode', dest='datalogger_mode', action='store_true',
 379                       help='Run input files with datalogger.read_dict_of_arrays().  This is useful, for example, to test a single line from a tweakfile.',
 380                       default=False)
 381     parser.add_option('-q', '--disable-quadratic', dest='quadratic', action='store_false',
 382                       help='Disable quadratic term in fitting (i.e. use bounded linear fits).',
 383                       default=True)
 384     parser.add_option('-v', '--verbose', dest='verbose', action='store_true',
 385                       help='Print lots of debugging information',
 386                       default=False)
 387
 388     options,args = parser.parse_args()
 389     parser.destroy()
 390     assert len(args) >= 1, "Need an input file"
 391
 392     ifilename = args[0]
 393
 394     if options.ofilename != None :
 395         ofile = file(options.ofilename, 'w')
 396     else :
 397         ofile = sys.stdout
 398     config.TEXT_VERBOSE = options.verbose
 399     config.PYLAB_INTERACTIVE = False
 400     config.PYLAB_VERBOSE = options.pylab
 401     config.GNUPLOT_VERBOSE = False
 402     if options.quadratic == True:
 403         param_guesser = limited_quadratic_param_guess
 404         model = limited_quadratic
 405         sensitivity_from_fit_params = limited_quadratic_sensitivity
 406     else:
 407         param_guesser = limited_linear_param_guess
 408         model = limited_linear
 409         sensitivity_from_fit_params = limited_linear_sensitivity
 410
 411     if options.tweakmode == False :
 412         if options.datalogger_mode:
 413             data = bump_load(ifilename)
 414         else:
 415             data = read_data(ifilename)
 416         photoSensitivity = bump_analyze(data,
 417                                         param_guesser=param_guesser,
 418                                         model=model,
 419                                         sensitivity_from_fit_params=sensitivity_from_fit_params)
 420
 421         print >> ofile, photoSensitivity
 422     else : # tweak file mode
 423         slopes = bump_load_analyze_tweaked(ifilename,
 424                                            param_guesser=param_guesser,
 425                                            model=model,
 426                                            sensitivity_from_fit_params=sensitivity_from_fit_params)
 427         if options.comma_out :
 428             sep = ','
 429         else :
 430             sep = '\n'
 431         common.write_array(ofile, slopes, sep)
 432
 433     if options.ofilename != None :
 434         ofile.close()