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 common # common module for the calibcant package
  50 import config # config module for the calibcant package
  51 import data_logger
  52 import z_piezo_utils
  53 import linfit
  54 from splittable_kwargs import splittableKwargsFunction, \
  55     make_splittable_kwargs_function
  56
  57 #@splittableKwargsFunction((bump_plot, 'data'))
  58 # Some of the child functions aren't yet defined, so postpone
  59 # make-splittable until later in the module.
  60 def bump_analyze(data, zpGain=config.zpGain,
  61                  zpSensitivity=config.zpSensitivity,
  62                  Vzp_out2V=config.Vzp_out2V,
  63                  Vphoto_in2V=config.Vphoto_in2V,
  64                  **kwargs) :
  65     """
  66     Return the slope of the bump ;).
  67     Inputs:
  68       data        dictionary of data in DAC/ADC bits
  69       Vzp_out2V   function that converts output DAC bits to Volts
  70       Vphoto_in2V function that converts input ADC bits to Volts
  71       zpGain      zpiezo applied voltage per output Volt
  72       zpSensitivity  nm zpiezo response per applied Volt
  73     Returns:
  74      photoSensitivity (Vphoto/Zcant) in Volts/nm
  75     Checks for strong correlation (r-value) and low randomness chance (p-value)
  76
  77     With the current implementation, the data is regressed in DAC/ADC bits
  78     and THEN converted, so we're assuming that both conversions are LINEAR.
  79     if they aren't, rewrite to convert before the regression.
  80     """
  81     bump_plot_kwargs, = bump_analyze._splitargs(bump_analyze, kwargs)
  82     scale_Vzp_bits2V = Vzp_out2V(1) - Vzp_out2V(0)
  83     scale_Vphoto_bits2V = Vphoto_in2V(1) - Vphoto_in2V(0)
  84     Vphoto2Vzp_out_bit, intercept = \
  85         linfit.linregress(x=data["Z piezo output"],
  86                           y=data["Deflection input"],
  87                           plotVerbose=False)
  88     Vphoto2Vzp_out = Vphoto2Vzp_out_bit *scale_Vphoto_bits2V / scale_Vzp_bits2V
  89
  90     #               1 / (Vzp/Vzp_out  *     Zp/Vzp       * Zcant/Zp )
  91     Vzp_out2Zcant = 1.0/ (zpGain      * zpSensitivity) # * 1
  92     bump_plot(data, **bump_plot_kwargs)
  93     return Vphoto2Vzp_out * Vzp_out2Zcant
  94
  95 @splittableKwargsFunction()
  96 def bump_save(data, log_dir=None) :
  97     "Save the dictionary data, using data_logger.data_log()"
  98     if log_dir != None :
  99         log = data_logger.data_log(log_dir, noclobber_logsubdir=False,
 100                                    log_name="bump")
 101         log.write_dict_of_arrays(data)
 102
 103 def bump_load(datafile) :
 104     "Load the dictionary data, using data_logger.date_load()"
 105     dl = data_logger.data_load()
 106     data = dl.read_dict_of_arrays(datafile)
 107     return data
 108
 109 @splittableKwargsFunction()
 110 def bump_plot(data, plotVerbose=False) :
 111     "Plot the bump (Vphoto vs Vzp) if plotVerbose or PYLAB_VERBOSE == True"
 112     if plotVerbose or config.PYLAB_VERBOSE :
 113         common._import_pylab()
 114         common._pylab.figure(config.BASE_FIGNUM)
 115         common._pylab.plot(data["Z piezo output"], data["Deflection input"],
 116                            '.', label='bump')
 117         common._pylab.title("bump surface")
 118         common._pylab.legend(loc='upper left')
 119         common._flush_plot()
 120
 121 make_splittable_kwargs_function(bump_analyze, (bump_plot, 'data'))
 122
 123 @splittableKwargsFunction((bump_analyze, 'data'))
 124 def bump_load_analyze_tweaked(tweak_file, **kwargs):
 125     "Load the output file of tweak_calib_bump.sh, return an array of slopes"
 126     bump_analyze_kwargs, = \
 127         bump_load_analyze_tweaked._splitargs(bump_load_analyze_tweaked, kwargs)
 128     photoSensitivity = []
 129     for line in file(tweak_file, 'r') :
 130         parsed = line.split()
 131         path = parsed[0].strip()
 132         if path[0] == '#' : # a comment
 133             continue
 134         if config.TEXT_VERBOSE :
 135             print "Reading data from %s with ranges %s" % (path, parsed[1:])
 136         # read the data
 137         full_data = bump_load(path)
 138         if len(parsed) == 1 :
 139             data = full_data # use whole bump
 140         else :
 141             # use the listed sections
 142             zp = []
 143             df = []
 144             for rng in parsed[1:] :
 145                 p = rng.split(':')
 146                 starti = int(p[0])
 147                 stopi = int(p[1])
 148                 zp.extend(full_data['Z piezo output'][starti:stopi])
 149                 df.extend(full_data['Deflection input'][starti:stopi])
 150             data = {'Z piezo output': numpy.array(zp),
 151                     'Deflection input': numpy.array(df)}
 152         pSi = bump_analyze(data, **bump_analyze_kwargs)
 153         photoSensitivity.append(pSi)
 154     return numpy.array(photoSensitivity, dtype=numpy.float)
 155
 156 # commandline interface functions
 157 import scipy.io, sys
 158
 159 def read_data(ifile):
 160     "ifile can be a filename string or open (seekable) file object"
 161     if ifile == None :  ifile = sys.stdin
 162     unlabeled_data=scipy.io.read_array(ifile)
 163     data = {}
 164     data['Z piezo output'] = unlabeled_data[:,0]
 165     data['Deflection input'] = unlabeled_data[:,1]
 166     return data
 167
 168 def remove_further_than(data, zp_crit) :
 169     ndata = {}
 170     ndata['Z piezo output'] = []
 171     ndata['Deflection input'] = []
 172     for zp,df in zip(data['Z piezo output'],data['Deflection input']) :
 173         if zp > zp_crit :
 174             ndata['Z piezo output'].append(zp)
 175             ndata['Deflection input'].append(df)
 176     return ndata
 177
 178 if __name__ == '__main__' :
 179     # command line interface
 180     from optparse import OptionParser
 181
 182     usage_string = ('%prog <input-file>\n'
 183                     '2008, W. Trevor King.\n'
 184                     '\n'
 185                     'There are two operation modes, one to analyze a single bump file,\n'
 186                     'and one to analyze tweak files.\n'
 187                     '\n'
 188                     'Single file mode (the default) :\n'
 189                     'Scales raw DAC/ADC bit data and fits a straight line.\n'
 190                     'Returns photodiode sensitivity Vphotodiode/Zcantilever in V/nm.\n'
 191                     '<input-file> should be whitespace-delimited, 2 column ASCII\n'
 192                     'without a header line.  e.g: "<zp_DAC>\\t<deflection_ADC>\\n"\n'
 193                     '\n'
 194                     'Tweak file mode:\n'
 195                     'Runs the same analysis as in single file mode for each bump in\n'
 196                     'a tweak file.  Each line in the tweak file specifies a single bump.\n'
 197                     'Blank lines and those beginning with a pound sign (#) are ignored.\n'
 198                     'The format of a line is a series of whitespace-separated fields--\n'
 199                     'a base file path followed by optional point index ranges, e.g.:\n'
 200                     '20080919/20080919132500_bump_surface 10:651 1413:2047\n'
 201                     'which only discards all points outside the index ranges [10,651)\n'
 202                     'and [1413,2047) (indexing starts at 0).\n'
 203                     )
 204     parser = OptionParser(usage=usage_string, version='%prog '+common.VERSION)
 205     parser.add_option('-C', '--cut-contact', dest='cut',
 206                       help='bilinear fit to cut out contact region (currently only available in single-file mode)',
 207                       action='store_true', default=False)
 208     parser.add_option('-o', '--output-file', dest='ofilename',
 209                       help='write output to FILE (default stdout)',
 210                       type='string', metavar='FILE')
 211     parser.add_option('-c', '--comma-out', dest='comma_out', action='store_true',
 212                       help='Output comma-seperated values (default %default)',
 213                       default=False)
 214     parser.add_option('-p', '--pylab', dest='pylab', action='store_true',
 215                       help='Produce pylab fit checks during execution',
 216                       default=False)
 217     parser.add_option('-t', '--tweak-mode', dest='tweakmode', action='store_true',
 218                       help='Run in tweak-file mode',
 219                       default=False)
 220     parser.add_option('-v', '--verbose', dest='verbose', action='store_true',
 221                       help='Print lots of debugging information',
 222                       default=False)
 223
 224     options,args = parser.parse_args()
 225     parser.destroy()
 226     assert len(args) >= 1, "Need an input file"
 227
 228     ifilename = args[0]
 229
 230     if options.ofilename != None :
 231         ofile = file(options.ofilename, 'w')
 232     else :
 233         ofile = sys.stdout
 234     config.TEXT_VERBOSE = options.verbose
 235     config.PYLAB_INTERACTIVE = False
 236     config.PYLAB_VERBOSE = options.pylab
 237     config.GNUPLOT_VERBOSE = False
 238
 239     if options.tweakmode == False :
 240         data = read_data(ifilename)
 241         if options.cut :
 242             ddict = {'approach':data} # although it may be any combination of approach and retract
 243             try :
 244                 params = z_piezo_utils.analyzeSurfPosData(ddict, retAllParams=True)
 245                 a,b,c,d = tuple(params) # model : f(x) = x<c ? a + b*x : (a+b*c) + d*(x-c)
 246                 print >> sys.stderr, "fit with", params, ". using zp < %d" % c
 247                 data = remove_further_than(data, c)
 248             except z_piezo_utils.poorFit, s :
 249                 # data not very bilinear, so don't cut anything off.
 250                 print >> sys.stderr, "use everything"
 251
 252         photoSensitivity = bump_analyze(data)
 253
 254         print >> ofile, photoSensitivity
 255     else : # tweak file mode
 256         slopes = bump_load_analyze_tweaked(ifilename)
 257         if options.comma_out :
 258             sep = ','
 259         else :
 260             sep = '\n'
 261         common.write_array(ofile, slopes, sep)
 262
 263     if options.ofilename != None :
 264         ofile.close()