calibcant/bump_analyze.py

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