calibcant/bump_analyze.py

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