Major update: vectorize secant method
authorDan White <dan@whiteaudio.com>
Sat, 1 Sep 2012 03:51:09 +0000 (22:51 -0500)
committerDan White <dan@whiteaudio.com>
Sat, 1 Sep 2012 03:51:09 +0000 (22:51 -0500)
python-lib/calibrate.py
python-lib/mpsse-test.py

index 530d4ff90c256b68bd8e40cc3bbd58f2413c0f08..633f219ab997c42e0866e69df4923b2a99a9ca1d 100644 (file)
@@ -54,7 +54,7 @@ def bisect(fset, limits=[0.0, 1.0], dgmin=None, up='uk', down='dj'):
 
 
 
-def secant_opt(func, x0, limits, x1scale=0.5, tol=1e-3, maxiter=50):
+def secant_opt(func, x0, limits, x1scale=0.1, maxiter=100, k=1.0):
     """
     Find the function's zero on a bounded interval using the secant method.
 
@@ -70,51 +70,90 @@ def secant_opt(func, x0, limits, x1scale=0.5, tol=1e-3, maxiter=50):
     xmin = min(limits)
     xmax = max(limits)
 
-    converged = False
+    converging = ones(x0.shape, dtype=int)
+    evaluations = zeros(x0.shape, dtype=int)
 
     def saturate(x):
-        return int(round(max(min(x, xmax), xmin)))
+        imax = find(x > xmax)
+        imin = find(x < xmin)
+        x.flat[imax] = xmax
+        x.flat[imin] = xmin
+        return int64(around(x))
+        #return int(round(max(min(x, xmax), xmin)))
 
     def eval(x):
+        print 'eval:', x[0:N_CHANNELS]
         r = func(x)
+        print r.shape, x.shape
         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))
+    r = xmax - xmin
+    print type(q0), q0
+    p1 = saturate(p0 - sign(q0)*x1scale*r)
+    #if q0 > 0.0:
+        #p1 = saturate(p0 - x1scale * r)
+    #else:
+        #p1 = saturate(p0 + x1scale * r)
 
     q1 = eval(p1)
     for niter in range(maxiter):
+        iz = find(q1==q0)
+        q1.flat[iz] = q1.flat[iz] + 1e-6
+        #dp = converging*q1*(p1 - p0)/(q1 - q0)
         dp = q1*(p1 - p0)/(q1 - q0)
+        dp = q1*abs(p1 - p0)/abs(q1 - q0)
+        #dp = k*q1*(p1 - p0)/(q1 - q0)
+        #dp = k*sign(q1)*abs(p1-p0)
+
+        ia = find(abs(dp) > (r * 2**(-niter)))
+        dp.flat[ia] = sign(dp.flat[ia]) * (r * 2**(-niter))
 
         #force a +-1 step
-        if abs(dp) < 1.0:
-            if dp < 0.0:
-                dp = -1
-            else:
-                dp = 1
+        ismall = find(abs(dp) <= 1.0)
+        dp.flat[ismall] = sign(dp.flat[ismall])
+
+        #unless the channel is done...
+        dp *= converging
+
+        #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)
-            #ensure we pick the correct value
-            if (a0 <= a1) and (a <= a1):
-                return p0
-            #else:
-                #print "false alarm, almost there..."
-                ##return p1
+        #d0 = p1 - p0
+        #d1 = p - p1
+        #if all(abs(d0) == 1) and all(abs(d1) == 1) and (d0 == -1*d1):
+        for i in range(x0.size):
+            d0 = p1.flat[i] - p0.flat[i]
+            d1 = p.flat[i]- p1.flat[i]
+            if (abs(d0) == 1) and (abs(d1) == 1) and (d0 == -1*d1):
+                a0 = abs(q0.flat[i])
+                a1 = abs(q1.flat[i])
+                a = abs(q.flat[i])
+                #ensure we pick the correct value
+                if (a0 <= a1) and (a <= a1):
+                    converging.flat[i] = 0
+                    evaluations.flat[i] = 3 + niter
+                    #print 'channel done:', i
+                    #print converging
+                    #return p0
+                #else:
+                    #print "false alarm, almost there..."
+                    ##return p1
+        print converging[0:N_CHANNELS]
+        if not (sum(converging[0:N_CHANNELS]) > 0.0): 
+            return (p0, {'nevals': evaluations,
+                         'guesses': guesses,
+                         'results': results})
 
         p0 = p1
         q0 = q1
@@ -131,7 +170,8 @@ def offset2signed(v, bits):
 
 
 def mux_a_offset(x):
-    xi = int(round(x))
+    #xi = int(round(x))
+    xi = int(x[0])
     mux.otaA.offset = xi
     mux.write()
     adc.read()
@@ -141,10 +181,10 @@ def mux_a_offset(x):
         sleep(DELAY)
     mv = mean(values)
     print tpx(xi), tplot(mv)
-    return mv
+    return array((mv,))
 
 def mux_b_offset(x):
-    xi = int(round(x))
+    xi = int(x[0])
     mux.otaB.offset = xi
     mux.write()
     adc.read()
@@ -154,15 +194,16 @@ def mux_b_offset(x):
         sleep(DELAY)
     mv = mean(values)
     print tpx(xi), tplot(mv)
-    return mv
+    return array((mv,))
 
 def chain_a_offset(x, n, mux_offset):
     xi = int(round(x))
     cn = chain.h[n]
     cn.cal = 1
-    cn.otaA.se = 1
+    cn.otaA.se = 0
     cn.otaA.cint = 0
     cn.otaA.offset = xi
+    cn.otaA.gain = 15
     chain.write()
     sleep(0.5)
     adc.read()
@@ -172,7 +213,86 @@ def chain_a_offset(x, n, mux_offset):
         sleep(DELAY)
     mv = mean(values)
     print tpx(xi), tplot(mv)
-    return mv - mux_offset
+    return mv #- mux_offset
+
+
+def chain_outputs(n):
+    adc.triggerMode(adc.MODE_IDLE)
+    adc.mux(4) #mux.otaA
+    adc.fifo(1)
+
+    plotter = TermPlotter([-1000, 1000], width=80)
+
+    #outA = []
+    #outB = []
+    outs = zeros((2, n+1))
+    print outs.shape, outs
+    for i in range(n+1):
+        print 'i:', i
+        mux.selA = i
+        mux.selB = i
+        mux.write()
+
+        for j in range(2):
+            print 'k:', j
+            adc.mux(4+j)
+            sleep(0.1)
+            adc.triggerMode(adc.MODE_MANUAL_MANUAL)
+            #adc.triggerMode(adc.MODE_AUTO_AUTO_SINGLE)
+            adc.read()
+            values = []
+            #print
+            for k in range(N_SAMPLES):
+                v = offset2signed(adc.read(), 16)
+                #print plotter(v)
+                values.append(v)
+                sleep(DELAY)
+            mv = mean(values[-1*min(10, N_SAMPLES):-1])
+
+            #sleep(0.01)
+            outs[j, i] = mv
+            #a = mv
+            #b = offset2signed(adc.read(), 16)
+        #print '%02i:'%i, a, b
+        #outA.append(a)
+        #outB.append(b)
+
+    print outs[:,0:N_CHANNELS]
+    return outs
+    #return array([array(outA), array(outB)])
+
+def chain_offsets(values, mux_offset):
+    #for i,n in enumerate(values):
+        #chain.h[i].otaA.offset = int(n)
+        #chain.h[i].cal = 1
+        #chain.h[i].otaA.se = 1
+        #chain.h[i].otaA.cint = 0
+        #chain.h[i].otaA.zero = 1
+    #chain.write()
+    for i,n in enumerate(values[0]):
+        chain.h[i].otaA.offset = int(n)
+        chain.h[i].cal = 1
+        chain.h[i].otaA.se = 1
+        chain.h[i].otaA.cint = 0
+        chain.h[i].otaA.zero = 0
+
+    for i,n in enumerate(values[1]):
+        chain.h[i].otaB.offset = int(n)
+        chain.h[i].otaB.se = 1
+        chain.h[i].otaB.cint = 0
+        chain.h[i].otaB.zero = 0
+    chain.write()
+    sleep(1.0)
+    outs = chain_outputs(N_CHANNELS-1)
+    for i in range(2):
+        outs[i,:] -= mux_offset[i]
+    #outs = [outA - mux_offset[0], outB - mux_offset[1]]
+    print 'here'
+    return array(outs)
+    a = randint(-10, 10, (48,))
+    #for j in range(N_CHANNELS):
+        #a[j] = outA[j]
+    #return a - mux_offset
 
 dac.vina(1.25)
 dac.vinb(1.25)
@@ -194,39 +314,61 @@ adc.triggerMode(adc.MODE_MANUAL_MANUAL)
 
 
 
+N_CHANNELS = 4
 N_SAMPLES = 10
 DELAY = 0.05
 
-print
-print 'Calibrating mux otaA'
-adc.mux(4)
-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 =', muxA_offset, ' @', resA
-
+if 0:
+    print
+    print 'Calibrating mux otaA'
+    adc.mux(4)
+    resA, tmp = secant_opt(mux_a_offset, zeros((1,), dtype=int), [-128, 127])
+    muxA_offset = mean([offset2signed(adc.read(), 16) for i in range(1000)])
+    print '*** mux.otaA.offset =', muxA_offset, ' @', resA
+else:
+    resA = 53
+    muxA_offset = 0.0
 
 
-if 1:
+if 0:
     print
     print 'Calibrating mux otaB'
     adc.mux(5)
-    resB = secant_opt(mux_b_offset, 0, [-128, 127])
+    #resB = secant_opt(mux_b_offset, 0, [-128, 127])
+    resB, tmp = secant_opt(mux_b_offset, zeros((1,), dtype=int), [-128, 127])
     muxB_offset = mean([offset2signed(adc.read(), 16) for i in range(1000)])
     print '*** mux.otaB.offset =', muxB_offset, ' @', resB
+else:
+    resB = 59
+    muxB_offset = 0.0
 
 
-if 1:
+if 0:
     offset_chain_a = []
     for n in range(48):
         print
         print 'Calibrating chain %02i A' % n
+        adc.triggerMode(adc.MODE_IDLE)
         adc.mux(4)
+        adc.triggerMode(adc.MODE_MANUAL_MANUAL)
         mux.selA = n
         mux.selB = n
         mux.otaA.mode = mux.otaA.MUX_BUF
         mux.otaB.mode = mux.otaB.MUX_BUF
         mux.write()
-        nos = secant_opt(lambda x: chain_a_offset(x, n, muxA_offset), 0, (-128, 127))
+        nos = secant_opt(lambda x: chain_a_offset(x, n, muxA_offset), resA, (-128, 127))
         mean_os = mean([offset2signed(adc.read(), 16) for i in range(1000)])
         print '*** chain.h[%02i].otaA.offset =' % n, mean_os, ' @', nos
+        chain_outputs(n)
+
+if 1:
+    mux.otaA.mode = mux.otaA.MUX_BUF
+    mux.otaB.mode = mux.otaB.MUX_BUF
+    mux.write()
+    chain_os, stats = secant_opt(lambda x: chain_offsets(x, (muxA_offset, muxB_offset)), resA*ones((2,N_CHANNELS), dtype=int), (-128, 127), k=0.5)
+    print '-'*80
+    print 'offsets:', chain_os[:N_CHANNELS]
+    print 'N-evals:', stats['nevals'][:N_CHANNELS]
+    print '-'*80
+    print
 
index 0ca927796a6c579708882977bc22626fb6fd26c3..1cfae60d37071a82259d7b7f7ac6c5e0d780fe5e 100755 (executable)
@@ -281,7 +281,7 @@ for h in chain.h:
     h.cal = 0
     for ota in h.ota:
         ota.se = 1
-        ota.cint = 1
+        ota.cint = 0
         ota.zero = 0
         ota.fast = 1
         ota.gain = 15
@@ -297,8 +297,8 @@ mux.selB = 0
 for ota in mux.ota:
     ota.mode = ota.CAL_CMP
     ota.fast = 1
-    ota.gain = 8
-    ota.offset = 52
+    ota.gain = 15
+    ota.offset = 0
 mux.write()