From: W. Trevor King Date: Wed, 28 Jan 2009 13:46:10 +0000 (-0500) Subject: Shiny, new, flexible bump fitting framework. X-Git-Tag: 0.4~9 X-Git-Url: http://git.tremily.us/?p=calibcant.git;a=commitdiff_plain;h=8bb4ccc40432d607fe71b272a52fa1216b8f4b01 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 36a6b81..0dc6be2 100755 --- 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 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=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()" @@ -107,18 +188,41 @@ 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.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): @@ -202,9 +306,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') @@ -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('-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) @@ -237,18 +341,10 @@ if __name__ == '__main__' : 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> 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