Set default maxFreq to 25 kHz. Unset pylab.hold for the bump_plot residual.
[calibcant.git] / calibcant / bump_analyze.py
index 11a8f51e719f7f8b96220bec7fca47b4357ddb06..748287967cf8d38d29801e353dd217f162bef608 100755 (executable)
@@ -46,15 +46,16 @@ slope Vphoto/Vout is converted to photoSensitivity with bump_analyze().
 """
 
 import numpy
+import scipy.optimize
 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'))
+
+#@splittableKwargsFunction((bump_fit, 'zpiezo_output_bits',
+#                           'deflection_input_bits'))
 # 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,
@@ -78,20 +79,100 @@ def bump_analyze(data, zpGain=config.zpGain,
     and THEN converted, so we're assuming that both conversions are LINEAR.
     if they aren't, rewrite to convert before the regression.
     """
-    bump_plot_kwargs, = bump_analyze._splitargs(bump_analyze, kwargs)
+    bump_fit_kwargs, = bump_analyze._splitargs(bump_analyze, kwargs)
+    Vphoto2Vzp_out_bit = bump_fit(data['Z piezo output'],
+                                  data['Deflection input'],
+                                  **bump_fit_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
 
+def limited_quadratic(x, params):
+    """
+    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)
+    """
+    high_voltage_rail = 2**16 - 1 # bits
+    x_contact,y_contact,slope,quad = params
+    y = slope*(x-x_contact) + quad*(x-x_contact)**2+ y_contact
+    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.
+    """
+    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])
+    slope = (y_high - y_contact) / (x_high - x_contact)
+    quad = 0
+    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
+
+@splittableKwargsFunction()
+def bump_fit(zpiezo_output_bits, deflection_input_bits,
+             paramGuesser=limited_quadratic_param_guess,
+             model=limited_quadratic,
+             sensitivity_from_fit_params=limited_quadratic_sensitivity,
+             plotVerbose=False) :
+    x = zpiezo_output_bits
+    y = deflection_input_bits
+    def residual(p, y, x) :
+        return model(x, p) - y
+    paramGuess = paramGuesser(x, y)
+    p,cov,info,mesg,ier = \
+        scipy.optimize.leastsq(residual, paramGuess, args=(y, x),
+                               full_output=True, maxfev=int(10e3))
+    if config.TEXT_VERBOSE :
+        print "Fitted params:",p
+        print "Covariance mx:",cov
+        print "Info:", info
+        print "mesg:", mesg
+        if ier == 1 :
+            print "Solution converged"
+        else :
+            print "Solution did not converge"
+    if plotVerbose or config.PYLAB_VERBOSE :
+        yguess = model(x, paramGuess)
+        #yguess = None # Don't print the guess, since I'm convinced it's ok ;).
+        yfit = model(x, p)
+        bump_plot(data={"Z piezo output":x, "Deflection input":y},
+                  yguess=yguess, yfit=yfit, plotVerbose=plotVerbose)
+    return sensitivity_from_fit_params(p)
+
 @splittableKwargsFunction()
 def bump_save(data, log_dir=None) :
     "Save the dictionary data, using data_logger.data_log()"
@@ -107,18 +188,42 @@ def bump_load(datafile) :
     return data
 
 @splittableKwargsFunction()
-def bump_plot(data, plotVerbose=False) :
+def bump_plot(data, yguess=None, yfit=None, 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)
+        if yfit != None: # two subplot figure
+            common._pylab.subplot(211)
+        common._pylab.hold(False)
         common._pylab.plot(data["Z piezo output"], data["Deflection input"],
                            '.', label='bump')
+        common._pylab.hold(True)
+        if yguess != None:
+            common._pylab.plot(data["Z piezo output"], yguess,
+                               'g-', label='guess')
+        if yfit != None:
+            common._pylab.plot(data["Z piezo output"], yfit,
+                               'r-', label='fit')
+        common._pylab.hold(False)
         common._pylab.title("bump surface")
         common._pylab.legend(loc='upper left')
+        common._pylab.xlabel("Z piezo output voltage (bits)")
+        common._pylab.ylabel("Photodiode input voltage (bits)")
+        if yfit != None:
+            # second subplot for residual
+            common._pylab.subplot(212)
+            common._pylab.plot(data["Z piezo output"],
+                               data["Deflection input"] - yfit,
+                               'r-', label='residual')
+            common._pylab.legend(loc='upper right')
+            common._pylab.xlabel("Z piezo output voltage (bits)")
+            common._pylab.ylabel("Photodiode input voltage (bits)")
         common._flush_plot()
 
-make_splittable_kwargs_function(bump_analyze, (bump_plot, 'data'))
+make_splittable_kwargs_function(bump_analyze,
+                                (bump_fit, 'zpiezo_output_bits',
+                                 'deflection_input_bits'))
 
 @splittableKwargsFunction((bump_analyze, 'data'))
 def bump_load_analyze_tweaked(tweak_file, **kwargs):
@@ -128,7 +233,9 @@ def bump_load_analyze_tweaked(tweak_file, **kwargs):
     photoSensitivity = []
     for line in file(tweak_file, 'r') :
         parsed = line.split()
-        path = parsed[0].split('\n')[0]
+        path = parsed[0].strip()
+        if path[0] == '#' : # a comment
+            continue
         if config.TEXT_VERBOSE :
             print "Reading data from %s with ranges %s" % (path, parsed[1:])
         # read the data
@@ -192,6 +299,7 @@ if __name__ == '__main__' :
                     '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'
+                    'Blank lines and those beginning with a pound sign (#) are ignored.\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'
@@ -199,9 +307,6 @@ if __name__ == '__main__' :
                     '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')
@@ -214,6 +319,9 @@ if __name__ == '__main__' :
     parser.add_option('-t', '--tweak-mode', dest='tweakmode', action='store_true',
                       help='Run in tweak-file mode',
                       default=False)
+    parser.add_option('-d', '--datalogger-mode', dest='datalogger_mode', action='store_true',
+                      help='Run input files with datalogger.read_dict_of_arrays().  This is useful, for example, to test a single line from a tweakfile.',
+                      default=False)
     parser.add_option('-v', '--verbose', dest='verbose', action='store_true',
                       help='Print lots of debugging information',
                       default=False)
@@ -232,20 +340,12 @@ if __name__ == '__main__' :
     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<c ? a + b*x : (a+b*c) + d*(x-c)
-                print >> 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"
-                
+        if options.datalogger_mode:
+            data = bump_load(ifilename)
+        else:
+            data = read_data(ifilename)
         photoSensitivity = bump_analyze(data)
         
         print >> ofile, photoSensitivity