From 7c4a36795040ba4ec5b28c18f9cca269a8d15b20 Mon Sep 17 00:00:00 2001 From: Dan White Date: Tue, 28 Aug 2012 15:02:34 -0500 Subject: [PATCH] Add secant optimization algorithm --- python-lib/calibrate.py | 150 +++++++++++++++++++++++++++++--------- python-lib/termplotter.py | 3 +- 2 files changed, 119 insertions(+), 34 deletions(-) diff --git a/python-lib/calibrate.py b/python-lib/calibrate.py index c64d220..9888798 100644 --- a/python-lib/calibrate.py +++ b/python-lib/calibrate.py @@ -1,8 +1,10 @@ +from pprint import pprint from termplotter import TermPlotter -tplot = TermPlotter([-2**11, 2**11]) +tplot = TermPlotter([-2**11, 2**11], width=80) +tpx = TermPlotter([-128, 127], width=80) def bisect(fset, limits=[0.0, 1.0], dgmin=None, up='uk', down='dj'): @@ -52,6 +54,80 @@ def bisect(fset, limits=[0.0, 1.0], dgmin=None, up='uk', down='dj'): +from scipy.optimize import brentq +def secant_opt(func, x0, limits, x1scale=0.5, tol=1e-3, maxiter=50): + """ + Find the function's zero on a bounded interval using the secant method. + + x0 - initial integer guess + limits - tuple of min,max for x + + Evaluates func(x0), guesses x1 by sign of result and x1scale'th of the + range in that direction. The expected function shape is atan()-like. + """ + p0 = x0 + guesses = [] + results = [] + xmin = min(limits) + xmax = max(limits) + + converged = False + + def saturate(x): + return int(round(max(min(x, xmax), xmin))) + + def eval(x): + r = func(x) + guesses.append(x) + results.append(r) + return r + + q0 = eval(p0) + if q0 > 0.0: + p1 = saturate(p0 - x1scale * abs(xmax - p0)) + else: + p1 = saturate(p0 + x1scale * abs(p0 - xmin)) + + q1 = eval(p1) + for niter in range(maxiter): + dp = q1*(p1 - p0)/(q1 - q0) + + #force a +-1 step + if abs(dp) < 1.0: + if dp < 0.0: + dp = -1 + else: + dp = 1 + p = saturate(p1 - dp) + q = eval(p) + + # Convergence criteria + # x to varies by (+1,-1) or (-1,+1) + # or 3 consecutive guesses the same + d0 = p1 - p0 + d1 = p - p1 + if (abs(d0) == 1) and (abs(d1) == 1) and (d0 == -1*d1): + a0 = abs(q0) + a1 = abs(q1) + a = abs(q) + if (a0 <= a1) and (a <= a1): + #print "*** turnaround convergence ***" + #print ''.join(['%5i ' % i for i in guesses]) + #print ''.join(['%9.2e' % i for i in results]) + return p0 + #else: + #print "false alarm, almost there..." + ##return p1 + + p0 = p1 + q0 = q1 + p1 = p + q1 = q + + raise Error("Failed to converge.") + + + def offset2signed(v, bits): return (v - 2**(bits-1)) >> 2 @@ -67,8 +143,9 @@ def mux_a_offset(x): values.append(offset2signed(adc.read(), 16)) sleep(DELAY) mv = mean(values) - print xi, mv - print tplot(mv) + #print xi, mv + #print '%4i'%xi, tplot(mv) + print tpx(xi), tplot(mv) return mv def mux_b_offset(x): @@ -81,8 +158,10 @@ def mux_b_offset(x): values.append(offset2signed(adc.read(), 16)) sleep(DELAY) mv = mean(values) - print xi, mv - print tplot(mv) + #print xi, mv + #print tplot(mv) + #print '%4i'%xi, tplot(mv) + print tpx(xi), tplot(mv) return mv def chain_a_offset(x, n, mux_offset): @@ -100,8 +179,10 @@ def chain_a_offset(x, n, mux_offset): values.append(offset2signed(adc.read(), 16)) sleep(DELAY) mv = mean(values) - print xi, mv - print tplot(mv) + #print xi, mv + #print tplot(mv) + #print '%4i'%xi, tplot(mv) + print tpx(xi), tplot(mv) return mv - mux_offset dac.vina(1.25) @@ -132,37 +213,40 @@ DELAY = 0.05 print print 'Calibrating mux otaA' adc.mux(4) -resA = fopt(mux_a_offset, -128, 127, xtol=0.01, disp=True) -print resA +#resA = fopt(mux_a_offset, -128, 127, xtol=0.01, disp=True) +#def secant_opt(func, x0, limits, x1scale=0.5, tol=1e-3, maxiter=50): +resA = secant_opt(mux_a_offset, 0, [-128, 127]) muxA_offset = mean([offset2signed(adc.read(), 16) for i in range(1000)]) -print '*** mux.otaA.offset =', int(round(resA)) +print '*** mux.otaA.offset =', muxA_offset, ' @', resA + -if 0: +if 1: print print 'Calibrating mux otaB' adc.mux(5) - resB = fopt(mux_b_offset, -128, 127, xtol=0.01, disp=True) - print resB + #resB = fopt(mux_b_offset, -128, 127, xtol=0.01, disp=True) + resB = secant_opt(mux_b_offset, 0, [-128, 127]) muxB_offset = mean([offset2signed(adc.read(), 16) for i in range(1000)]) - print '*** mux.otaB.offset =', int(round(resB)) - - -offset_chain_a = [] -for n in range(48): - print - print 'Calibrating chain %02i A' % n - adc.mux(4) - mux.selA = n - mux.selB = n - mux.otaA.mode = mux.otaA.MUX_BUF - mux.otaB.mode = mux.otaB.MUX_BUF - mux.write() - nos = fopt(lambda x: chain_a_offset(x, n, muxA_offset), - -128, 127, - xtol=0.01, - disp=True) - mean_os = mean([offset2signed(adc.read(), 16) for i in range(1000)]) - print nos, mean_os - print '*** chain.h[%02i].otaA.offset = %+3i' % (n, int(round(nos))) + print '*** mux.otaB.offset =', muxB_offset, ' @', resB + + +if 1: + offset_chain_a = [] + for n in range(48): + print + print 'Calibrating chain %02i A' % n + adc.mux(4) + mux.selA = n + mux.selB = n + mux.otaA.mode = mux.otaA.MUX_BUF + mux.otaB.mode = mux.otaB.MUX_BUF + mux.write() + #nos = fopt(lambda x: chain_a_offset(x, n, muxA_offset), + #-128, 127, + #xtol=0.01, + #disp=True) + nos = secant_opt(lambda x: chain_a_offset(x, n, muxA_offset), 0, (-128, 127)) + mean_os = mean([offset2signed(adc.read(), 16) for i in range(1000)]) + print '*** chain.h[%02i].otaA.offset =' % n, mean_os, ' @', nos diff --git a/python-lib/termplotter.py b/python-lib/termplotter.py index b1a1d28..f1e46f2 100644 --- a/python-lib/termplotter.py +++ b/python-lib/termplotter.py @@ -27,7 +27,8 @@ class TermPlotter(object): fwidth = resolution_len + 2 if mantissa_len == 0: fwidth += 1 - s = '%%%i.0%if' % (fwidth, fraction_len) + s = '%%%i.0%if' % (fwidth, fraction_len if fraction_len>0 else 0) + #print s, fwidth, fraction_len self.leader = s + ' ' else: self.leader = leader -- 2.25.1