From: Dan White Date: Sat, 1 Sep 2012 03:51:09 +0000 (-0500) Subject: Major update: vectorize secant method X-Git-Tag: calibrations~13 X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=711f7b165b8bbc5e274cab6699522d3ddf0e134c;p=430.git Major update: vectorize secant method --- diff --git a/python-lib/calibrate.py b/python-lib/calibrate.py index 530d4ff..633f219 100644 --- a/python-lib/calibrate.py +++ b/python-lib/calibrate.py @@ -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 diff --git a/python-lib/mpsse-test.py b/python-lib/mpsse-test.py index 0ca9277..1cfae60 100755 --- a/python-lib/mpsse-test.py +++ b/python-lib/mpsse-test.py @@ -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()