"""
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,
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()"
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):
'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('-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)
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
projects / calibcant.git / commitdiff