def solve_mec_3constr(constraint_fnl, n = 30, initparams = None):
    delta = 1e-3
    if initparams:
        params = N.array(initparams)
    else:
        params = [3.14, 0, 0]
    for i in range(n):
        k, lam1, lam2 = params
        xys, cost, x, y, th = run_elastica(-.5, .5, k, lam1, lam2)
        c1c, c2c, c3c = constraint_fnl(cost, x, y, th)
        print '% constraint_fnl =', c1c, c2c, c3c

        dc1s = []
        dc2s = []
        dc3s = []
        for j in range(len(params)):
            params1 = N.array(params)
            params1[j] += delta
            k, lam1, lam2 = params1
            xys, cost, x, y, th = run_elastica(-.5, .5, k, lam1, lam2)
            c1p, c2p, c3p = constraint_fnl(cost, x, y, th)
            params1 = N.array(params)
            params1[j] -= delta
            k, lam1, lam2 = params1
            xys, cost, x, y, th = run_elastica(-.5, .5, k, lam1, lam2)
            c1m, c2m, c3m = constraint_fnl(cost, x, y, th)
            dc1s.append((c1p - c1m) / (2 * delta))
            dc2s.append((2 * delta) / (c2p - c2m))
            dc3s.append((c3p - c3m) / (2 * delta))

        # Make Jacobian matrix to invert
        jm = N.array([dc1s, dc2s, dc3s])
        #print jm
        ji = la.inverse(jm)

        dp = N.dot(ji, [c1c, c2c, c3c])
        if i < n/2: scale = .25
        else: scale = 1
        params -= scale * dp
        print '%', params
    return params
dataset/ETHPy150Open google/roboto/third_party/spiro/curves/mecsolve.py/solve_mec_3constr/swapped operands of `/`