Shiny, new, flexible bump fitting framework.

The bump fitting algorithm has been upgraded due to my increased
understanding of the complexity of the photodiode.

Sensitivity (bump slope) decreases as you move away from a zero volt
photodiode output.  Most of the time anyway.  We should be able to
handle data with a gentle camber, and we need to know what voltage
corresponds to the zero-deflection signal.

This prompted the recent change of bump aquisition to begin just _off_
the surface instead of just _on_ the surface.  Now we have to deal with
finding the location of the surface kink, so I split the fitting out
into its own function which least squares fits the data to whatever
model it's passed (currently via Levenburg-Marquardt).

Since the model is now it's own stand-alone entity, I got fancy and
added quadratic fits (to deal with camber) and clipping (to deal with
out-of-range signals).  So there should be no more need of manually
clipping in the tweakfiles or using the --cut-contact option in
bump_analyze (which has been removed).

I also enhanced the output of bump_plot to show the guess, fit, and
residual, in case anyone needs convincing that the fit is working, or is
troubleshooting a new model.

diff --git a/calibcant/bump_analyze.py b/calibcant/bump_analyze.py
index 36a6b81d20edb5b287fdb1e25801ae960dc9b1fc..0dc6be228dd55f15fcdc9e18af18b2ddd425d796 100755 (executable)
--- a/calibcant/bump_analyze.py
+++ b/calibcant/bump_analyze.py
@@ -46,15 +46,16 @@ slope Vphoto/Vout is converted to photoSensitivity with bump_analyze().
 """
 
 import numpy
 """
 
 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 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
 
 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,
 # 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.
     """
     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)
     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
     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
     #               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
 
     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=True) :
+    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()"
 @splittableKwargsFunction()
 def bump_save(data, log_dir=None) :
     "Save the dictionary data, using data_logger.data_log()"


@@ -107,18 +188,41 @@ def bump_load(datafile) :
     return data
 
 @splittableKwargsFunction()
     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)
     "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.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.title("bump surface")
         common._pylab.legend(loc='upper left')
         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()
 
         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):
 
 @splittableKwargsFunction((bump_analyze, 'data'))
 def bump_load_analyze_tweaked(tweak_file, **kwargs):


@@ -202,9 +306,6 @@ if __name__ == '__main__' :
                     'and [1413,2047) (indexing starts at 0).\n'
                     )
     parser = OptionParser(usage=usage_string, version='%prog '+common.VERSION)
                     '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')
     parser.add_option('-o', '--output-file', dest='ofilename',
                       help='write output to FILE (default stdout)',
                       type='string', metavar='FILE')


@@ -217,6 +318,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('-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)
     parser.add_option('-v', '--verbose', dest='verbose', action='store_true',
                       help='Print lots of debugging information',
                       default=False)


@@ -237,18 +341,10 @@ if __name__ == '__main__' :
     config.GNUPLOT_VERBOSE = False
 
     if options.tweakmode == False :
     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
         photoSensitivity = bump_analyze(data)
         
         print >> ofile, photoSensitivity