Add secant optimization algorithm

diff --git a/python-lib/calibrate.py b/python-lib/calibrate.py
index c64d2200516bb1cc49e29737a956154e4151e30c..9888798d49968a306a19b660d5658281129186c6 100644 (file)
--- 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 b1a1d28e7983333357d737c97db666a1727fdf05..f1e46f2a3e821fde4dbde7579190684bc4f08213 100644 (file)
--- 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