# Solver based on direct Newton solving of 4 parameters for each curve
# segment

import sys
from math import *

from Numeric import *
import LinearAlgebra as la

import poly3
import band

class Seg:
    def __init__(self, chord, th):
        self.ks = [0., 0., 0., 0.]
        self.chord = chord
        self.th = th
    def compute_ends(self, ks):
        chord, ch_th = poly3.integ_chord(ks)
        l = chord / self.chord
        thl = ch_th - (-.5 * ks[0] + .125 * ks[1] - 1./48 * ks[2] + 1./384 * ks[3])
        thr = (.5 * ks[0] + .125 * ks[1] + 1./48 * ks[2] + 1./384 * ks[3]) - ch_th
        k0l = l * (ks[0] - .5 * ks[1] + .125 * ks[2] - 1./48 * ks[3])
        k0r = l * (ks[0] + .5 * ks[1] + .125 * ks[2] + 1./48 * ks[3])
        l2 = l * l
        k1l = l2 * (ks[1] - .5 * ks[2] + .125 * ks[3])
        k1r = l2 * (ks[1] + .5 * ks[2] + .125 * ks[3])
        l3 = l2 * l
        k2l = l3 * (ks[2] - .5 * ks[3])
        k2r = l3 * (ks[2] + .5 * ks[3])
        return (thl, k0l, k1l, k2l), (thr, k0r, k1r, k2r), l
    def set_ends_from_ks(self):
        self.endl, self.endr, self.l = self.compute_ends(self.ks)
    def fast_pderivs(self):
        l = self.l
        l2 = l * l
        l3 = l2 * l
        return [((.5, l, 0, 0), (.5, l, 0, 0)),
                ((-1./12, -l/2, l2, 0), (1./12, l/2, l2, 0)),
                ((1./48, l/8, -l2/2, l3), (1./48, l/8, l2/2, l3)),
                ((-1./480, -l/48, l2/8, -l3/2), (1./480, l/48, l2/8, l3/2))]
    def compute_pderivs(self):
        rd = 2e6
        delta = 1./rd
        base_ks = self.ks
        base_endl, base_endr, dummy = self.compute_ends(base_ks)
        result = []
        for i in range(4):
            try_ks = base_ks[:]
            try_ks[i] += delta
            try_endl, try_endr, dummy = self.compute_ends(try_ks)
            deriv_l = (rd * (try_endl[0] - base_endl[0]),
                       rd * (try_endl[1] - base_endl[1]),
                       rd * (try_endl[2] - base_endl[2]),
                       rd * (try_endl[3] - base_endl[3]))
            deriv_r = (rd * (try_endr[0] - base_endr[0]),
                       rd * (try_endr[1] - base_endr[1]),
                       rd * (try_endr[2] - base_endr[2]),
                       rd * (try_endr[3] - base_endr[3]))
            result.append((deriv_l, deriv_r))
        return result

class Node:
    def __init__(self, x, y, ty, th):
        self.x = x
        self.y = y
        self.ty = ty
        self.th = th
    def continuity(self):
        if self.ty == 'o':
            return 4
        elif self.ty in ('c', '[', ']'):
            return 2
        else:
            return 0

def mod_2pi(th):
    u = th / (2 * pi)
    return 2 * pi * (u - floor(u + 0.5))

def setup_path(path):
    segs = []
    nodes = []
    nsegs = len(path)
    if path[0][2] == '{':
        nsegs -= 1
    for i in range(nsegs):
        i1 = (i + 1) % len(path)
        x0, y0, t0 = path[i]
        x1, y1, t1 = path[i1]
        s = Seg(hypot(y1 - y0, x1 - x0), atan2(y1 - y0, x1 - x0))
        segs.append(s)
    for i in range(len(path)):
        x0, y0, t0 = path[i]

        if t0 in ('{', '}', 'v'):
            th = 0
        else:
            s0 = segs[(i + len(path) - 1) % len(path)]
            s1 = segs[i]
            th = mod_2pi(s1.th - s0.th)

        n = Node(x0, y0, t0, th)
        nodes.append(n)
    return segs, nodes

def count_vec(nodes):
    jincs = []
    n = 0
    for i in range(len(nodes)):
        i1 = (i + 1) % len(nodes)
        t0 = nodes[i].ty
        t1 = nodes[i1].ty
        if t0 in ('{', '}', 'v', '[') and t1 in ('{', '}', 'v', ']'):
            jinc = 0
        elif t0 in ('{', '}', 'v', '[') and t1 == 'c':
            jinc = 1
        elif t0 == 'c' and t1 in ('{', '}', 'v', ']'):
            jinc = 1
        elif t0 == 'c' and t1 == 'c':
            jinc = 2
        else:
            jinc = 4
        jincs.append(jinc)
        n += jinc
    return n, jincs

thscale, k0scale, k1scale, k2scale = 1, 1, 1, 1

def inversedot_woodbury(m, v):
    a = zeros((n, 11), Float)
    for i in range(n):
        for j in range(max(-7, -i), min(4, n - i)):
            a[i, j + 7] = m[i, i + j]
    print a
    al, indx, d = band.bandec(a, 7, 3)
    VtZ = identity(4, Float)
    Z = zeros((n, 4), Float)
    for i in range(4):
        u = zeros(n, Float)
        for j in range(4):
            u[j] = m[j, n - 4 + i]
        band.banbks(a, 7, 3, al, indx, u)
        for k in range(n):
            Z[k, i] = u[k]
        #Z[:,i] = u
        for j in range(4):
            VtZ[j, i] += u[n - 4 + j]
    print Z
    print VtZ
    H = la.inverse(VtZ)
    print H
    band.banbks(a, 7, 3, al, indx, v)
    return(v - dot(Z, dot(H, v[n - 4:])))

def inversedot(m, v):
    return dot(la.inverse(m), v)
    n, nn = m.shape
    if 1:
        for i in range(n):
            sys.stdout.write('% ')
            for j in range(n):
                if m[i, j] > 0: sys.stdout.write('+ ')
                elif m[i, j] < 0: sys.stdout.write('- ')
                else: sys.stdout.write('  ')
            sys.stdout.write('\n')

    cyclic = False
    for i in range(4):
        for j in range(n - 4, n):
            if m[i, j] != 0:
                cyclic = True
    print '% cyclic:', cyclic
    if not cyclic:
        a = zeros((n, 11), Float)
        for i in range(n):
            for j in range(max(-5, -i), min(6, n - i)):
                a[i, j + 5] = m[i, i + j]
        for i in range(n):
            sys.stdout.write('% ')
            for j in range(11):
                if a[i, j] > 0: sys.stdout.write('+ ')
                elif a[i, j] < 0: sys.stdout.write('- ')
                else: sys.stdout.write('  ')
            sys.stdout.write('\n')
        al, indx, d = band.bandec(a, 5, 5)
        print a
        band.banbks(a, 5, 5, al, indx, v)
        return v
    else:
        #return inversedot_woodbury(m, v)
        bign = 3 * n
        a = zeros((bign, 11), Float)
        u = zeros(bign, Float)
        for i in range(bign):
            u[i] = v[i % n]
            for j in range(-7, 4):
                a[i, j + 7] = m[i % n, (i + j + 7 * n) % n]
        #print a
        if 1:
            for i in range(bign):
                sys.stdout.write('% ')
                for j in range(11):
                    if a[i, j] > 0: sys.stdout.write('+ ')
                    elif a[i, j] < 0: sys.stdout.write('- ')
                    else: sys.stdout.write('  ')
                sys.stdout.write('\n')
        #print u
        al, indx, d = band.bandec(a, 5, 5)
        band.banbks(a, 5, 5, al, indx, u)
        #print u
        return u[n + 2: 2 * n + 2]

def iter(segs, nodes):
    n, jincs = count_vec(nodes)
    if n == 0: return
    print '%', jincs
    v = zeros(n, Float)
    m = zeros((n, n), Float)
    for i in range(len(segs)):
        segs[i].set_ends_from_ks()
    j = 0
    j0 = 0
    for i in range(len(segs)):
        i1 = (i + 1) % len(nodes)
        t0 = nodes[i].ty
        t1 = nodes[i1].ty
        seg = segs[i]

        derivs = seg.compute_pderivs()
        print '%derivs:', derivs

        jinc = jincs[i]   # the number of params on this seg
        print '%', t0, t1, jinc, j0

        # The constraints are laid out as follows:
        #   constraints that cross the node on the left
        #   constraints on the left side
        #   constraints on the right side
        #   constraints that cross the node on the right

        jj = j0    # the index into the constraint row we're writing
        jthl, jk0l, jk1l, jk2l = -1, -1, -1, -1
        jthr, jk0r, jk1r, jk2r = -1, -1, -1, -1

        # constraints crossing left

        if t0 == 'o':
            jthl = jj + 0
            jk0l = jj + 1
            jk1l = jj + 2
            jk2l = jj + 3
            jj += 4
        elif t0 in ('c', '[', ']'):
            jthl = jj + 0
            jk0l = jj + 1
            jj += 2

        # constraints on left

        if t0 in ('[', 'v', '{') and jinc == 4:
            jk1l = jj
            jj += 1
        if t0 in ('[', 'v', '{', 'c') and jinc == 4:
            jk2l = jj
            jj += 1

        # constraints on right

        if t1 in (']', 'v', '}') and jinc == 4:
            jk1r = jj
            jj += 1
        if t1 in (']', 'v', '}', 'c') and jinc == 4:
            jk2r = jj
            jj += 1

        # constraints crossing right

        jj %= n
        j1 = jj

        if t1 == 'o':
            jthr = jj + 0
            jk0r = jj + 1
            jk1r = jj + 2
            jk2r = jj + 3
            jj += 4
        elif t1 in ('c', '[', ']'):
            jthr = jj + 0
            jk0r = jj + 1
            jj += 2

        print '%', jthl, jk0l, jk1l, jk2l, jthr, jk0r, jk1r, jk2r

        if jthl >= 0:
            v[jthl] += thscale * (nodes[i].th - seg.endl[0])
            if jinc == 1:
                m[jthl][j] += derivs[0][0][0]
            elif jinc == 2:
                m[jthl][j + 1] += derivs[0][0][0]
                m[jthl][j] += derivs[1][0][0]
            elif jinc == 4:
                m[jthl][j + 2] += derivs[0][0][0]
                m[jthl][j + 3] += derivs[1][0][0]
                m[jthl][j + 0] += derivs[2][0][0]
                m[jthl][j + 1] += derivs[3][0][0]
        if jk0l >= 0:
            v[jk0l] += k0scale * seg.endl[1]
            if jinc == 1:
                m[jk0l][j] -= derivs[0][0][1]
            elif jinc == 2:
                m[jk0l][j + 1] -= derivs[0][0][1]
                m[jk0l][j] -= derivs[1][0][1]
            elif jinc == 4:
                m[jk0l][j + 2] -= derivs[0][0][1]
                m[jk0l][j + 3] -= derivs[1][0][1]
                m[jk0l][j + 0] -= derivs[2][0][1]
                m[jk0l][j + 1] -= derivs[3][0][1]
        if jk1l >= 0:
            v[jk1l] += k1scale * seg.endl[2]
            m[jk1l][j + 2] -= derivs[0][0][2]
            m[jk1l][j + 3] -= derivs[1][0][2]
            m[jk1l][j + 0] -= derivs[2][0][2]
            m[jk1l][j + 1] -= derivs[3][0][2]
        if jk2l >= 0:
            v[jk2l] += k2scale * seg.endl[3]
            m[jk2l][j + 2] -= derivs[0][0][3]
            m[jk2l][j + 3] -= derivs[1][0][3]
            m[jk2l][j + 0] -= derivs[2][0][3]
            m[jk2l][j + 1] -= derivs[3][0][3]

        if jthr >= 0:
            v[jthr] -= thscale * seg.endr[0]
            if jinc == 1:
                m[jthr][j] += derivs[0][1][0]
            elif jinc == 2:
                m[jthr][j + 1] += derivs[0][1][0]
                m[jthr][j + 0] += derivs[1][1][0]
            elif jinc == 4:
                m[jthr][j + 2] += derivs[0][1][0]
                m[jthr][j + 3] += derivs[1][1][0]
                m[jthr][j + 0] += derivs[2][1][0]
                m[jthr][j + 1] += derivs[3][1][0]
        if jk0r >= 0:
            v[jk0r] -= k0scale * seg.endr[1]
            if jinc == 1:
                m[jk0r][j] += derivs[0][1][1]
            elif jinc == 2:
                m[jk0r][j + 1] += derivs[0][1][1]
                m[jk0r][j + 0] += derivs[1][1][1]
            elif jinc == 4:
                m[jk0r][j + 2] += derivs[0][1][1]
                m[jk0r][j + 3] += derivs[1][1][1]
                m[jk0r][j + 0] += derivs[2][1][1]
                m[jk0r][j + 1] += derivs[3][1][1]
        if jk1r >= 0:
            v[jk1r] -= k1scale * seg.endr[2]
            m[jk1r][j + 2] += derivs[0][1][2]
            m[jk1r][j + 3] += derivs[1][1][2]
            m[jk1r][j + 0] += derivs[2][1][2]
            m[jk1r][j + 1] += derivs[3][1][2]
        if jk2r >= 0:
            v[jk2r] -= k2scale * seg.endr[3]
            m[jk2r][j + 2] += derivs[0][1][3]
            m[jk2r][j + 3] += derivs[1][1][3]
            m[jk2r][j + 0] += derivs[2][1][3]
            m[jk2r][j + 1] += derivs[3][1][3]

        j += jinc
        j0 = j1
    #print m
    dk = inversedot(m, v)
    dkmul = 1
    j = 0
    for i in range(len(segs)):
        jinc = jincs[i]
        if jinc == 1:
            segs[i].ks[0] += dkmul * dk[j]
        elif jinc == 2:
            segs[i].ks[0] += dkmul * dk[j + 1]
            segs[i].ks[1] += dkmul * dk[j + 0]
        elif jinc == 4:
            segs[i].ks[0] += dkmul * dk[j + 2]
            segs[i].ks[1] += dkmul * dk[j + 3]
            segs[i].ks[2] += dkmul * dk[j + 0]
            segs[i].ks[3] += dkmul * dk[j + 1]
        j += jinc

    norm = 0.
    for i in range(len(dk)):
        norm += dk[i] * dk[i]
    return sqrt(norm)


def plot_path(segs, nodes, tol = 1.0, show_cpts = False):
    if show_cpts:
        cpts = []
    j = 0
    cmd = 'moveto'
    for i in range(len(segs)):
        i1 = (i + 1) % len(nodes)
        n0 = nodes[i]
        n1 = nodes[i1]
        x0, y0, t0 = n0.x, n0.y, n0.ty
        x1, y1, t1 = n1.x, n1.y, n1.ty
        ks = segs[i].ks
        abs_ks = abs(ks[0]) + abs(ks[1] / 2) + abs(ks[2] / 8) + abs(ks[3] / 48)
        n_subdiv = int(ceil(.001 + tol * abs_ks))
        n_subhalf = (n_subdiv + 1) / 2
        if n_subdiv > 1:
            n_subdiv = n_subhalf * 2
        ksp = (ks[0] * .5, ks[1] * .25, ks[2] * .125, ks[3] * .0625)
        pside = poly3.int_3spiro_poly(ksp, n_subhalf)
        ksm = (ks[0] * -.5, ks[1] * .25, ks[2] * -.125, ks[3] * .0625)
        mside = poly3.int_3spiro_poly(ksm, n_subhalf)
        mside.reverse()
        for j in range(len(mside)):
            mside[j] = (-mside[j][0], -mside[j][1])
        if n_subdiv > 1:
            pts = mside + pside[1:]
        else:
            pts = mside[:1] + pside[1:]
        chord_th = atan2(y1 - y0, x1 - x0)
        chord_len = hypot(y1 - y0, x1 - x0)
        rot = chord_th - atan2(pts[-1][1] - pts[0][1], pts[-1][0] - pts[0][0])
        scale = chord_len / hypot(pts[-1][1] - pts[0][1], pts[-1][0] - pts[0][0])
        u, v = scale * cos(rot), scale * sin(rot)
        xt = x0 - u * pts[0][0] + v * pts[0][1]
        yt = y0 - u * pts[0][1] - v * pts[0][0]
        s = -.5
        for x, y in pts:
            xp, yp = xt + u * x - v * y, yt + u * y + v * x
            thp = (((ks[3]/24 * s + ks[2]/6) * s + ks[1] / 2) * s + ks[0]) * s + rot
            up, vp = scale / (1.5 * n_subdiv) * cos(thp), scale / (1.5 * n_subdiv) * sin(thp)
            if s == -.5:
                if cmd == 'moveto':
                    print xp, yp, 'moveto'
                    cmd = 'curveto'
            else:
                if show_cpts:
                    cpts.append((xlast + ulast, ylast + vlast))
                    cpts.append((xp - up, yp - vp))
                print xlast + ulast, ylast + vlast, xp - up, yp - vp, xp, yp, 'curveto'
            xlast, ylast, ulast, vlast = xp, yp, up, vp
            s += 1. / n_subdiv
        if t1 == 'v':
            j += 2
        else:
            j += 1
    print 'stroke'
    if show_cpts:
        for x, y in cpts:
            print 'gsave 0 0 1 setrgbcolor', x, y, 'translate circle fill grestore'

def plot_ks(segs, nodes, xo, yo, xscale, yscale):
    j = 0
    cmd = 'moveto'
    x = xo
    for i in range(len(segs)):
        i1 = (i + 1) % len(nodes)
        n0 = nodes[i]
        n1 = nodes[i1]
        x0, y0, t0 = n0.x, n0.y, n0.ty
        x1, y1, t1 = n1.x, n1.y, n1.ty
        ks = segs[i].ks
        chord, ch_th = poly3.integ_chord(ks)
        l = chord/segs[i].chord
        k0 = l * poly3.eval_cubic(ks[0], ks[1], .5 * ks[2], 1./6 * ks[3], -.5)
        print x, yo + yscale * k0, cmd
        cmd = 'lineto'
        k3 = l * poly3.eval_cubic(ks[0], ks[1], .5 * ks[2], 1./6 * ks[3], .5)
        k1 = k0 + l/3 * (ks[1] - 0.5 * ks[2] + .125 * ks[3])
        k2 = k3 - l/3 * (ks[1] + 0.5 * ks[2] + .125 * ks[3])
        print x + xscale / l / 3., yo + yscale * k1
        print x + 2 * xscale / l / 3., yo + yscale * k2
        print x + xscale / l, yo + yscale * k3, 'curveto'
        x += xscale / l
        if t1 == 'v':
            j += 2
        else:
            j += 1
    print 'stroke'
    print 'gsave 0.25 setlinewidth'
    print xo, yo, 'moveto', x, yo, 'lineto stroke'
    print 'grestore'

def plot_nodes(nodes, segs):
    for i in range(len(nodes)):
        n = nodes[i]
        print 'gsave', n.x, n.y, 'translate'
        if n.ty in ('c', '[', ']'):
            th = segs[i].th - segs[i].endl[0]
            if n.ty == ']': th += pi
            print 180 * th / pi, 'rotate'
        if n.ty == 'o':
            print 'circle fill'
        elif n.ty == 'c':
            #print '3 4 poly fill'
            print '0 -3 moveto 0 3 lineto stroke'
        elif n.ty in ('v', '{', '}'):
            print '45 rotate 3 4 poly fill'
        elif n.ty in ('[', ']'):
            print '0 -3 moveto 0 0 3 90 270 arc fill'
        else:
            print '5 3 poly fill'
        print 'grestore'

def prologue():
    print '%!PS-Adobe-3.0 EPSF'
    print '/ss 2 def'
    print '/circle { ss 0 moveto currentpoint exch ss sub exch ss 0 360 arc } bind def'
    print '/poly {'
    print '    dup false exch {'
    print '        0 3 index 2 index { lineto } { moveto } ifelse pop true'
    print '        360.0 2 index div rotate'
    print '    } repeat pop pop pop'
    print '} bind def'

def run_path(path, show_iter = False, n = 10, xo = 36, yo = 550, xscale = .25, yscale = 2000, pl_nodes = True, pl_path = True):
    segs, nodes = setup_path(path)
    print '.5 setlinewidth'
    for i in range(n):
        if i == n - 1 and pl_path:
            print '0 0 0 setrgbcolor 1 setlinewidth'
        elif i == 0:
            print '1 0 0 setrgbcolor'
        elif i == 1:
            print '0 0.5 0 setrgbcolor'
        elif i == 2:
            print '0.3 0.3 0.3 setrgbcolor'
        norm = iter(segs, nodes)
        print '% norm =', norm
        if show_iter or i == n - 1:
            #print '1 0 0 setrgbcolor'
            #plot_path(segs, nodes, tol = 5)
            #print '0 0 0 setrgbcolor'
	    for i in range(len(nodes)):
		node = nodes[i]
		print '% node', node.x, node.y
		if i < len(segs):
		    seg = segs[i]
		    print '% seg', seg.ks
            if pl_path:
                plot_path(segs, nodes, tol = 1.0)
            if xo != None:
                plot_ks(segs, nodes, xo, yo, xscale, yscale)
    if pl_nodes:
        plot_nodes(nodes, segs)

if __name__ == '__main__':
    if 1:
        path = [(100, 350, 'o'), (225, 350, 'o'), (350, 450, 'o'),
                (450, 400, 'o'), (315, 205, 'o'), (300, 200, 'o'),
                (285, 205, 'o')]

    if 1:
        path = [(100, 350, 'o'), (175, 375, '['), (250, 375, ']'), (325, 425, '['),
                (325, 450, ']'),
                (400, 475, 'o'), (350, 200, 'c')]

    if 0:
        ecc, ty, ty1 = 0.56199, 'c', 'c'
        ecc, ty, ty1 = 0.49076, 'o', 'o',
        ecc, ty, ty1 = 0.42637, 'o', 'c'
        path = [(300 - 200 * ecc, 300, ty), (300, 100, ty1),
                (300 + 200 * ecc, 300, ty), (300, 500, ty1)]

    # difficult path #3
    if 0:
        path = [(100, 300, '{'), (225, 350, 'o'), (350, 500, 'o'),
                (450, 500, 'o'), (450, 450, 'o'), (300, 200, '}')]

    if 0:
        path = [(100, 100, '{'), (200, 180, 'v'), (250, 215, 'o'),
                (300, 200, 'o'), (400, 100, '}')]

    if 0:
        praw = [
        (134, 90, 'o'),
        (192, 68, 'o'),
        (246, 66, 'o'),
        (307, 107, 'o'),
        (314, 154, '['),
        (317, 323, ']'),
        (347, 389, 'o'),
        (373, 379, 'v'),
        (386, 391, 'v'),
        (365, 409, 'o'),
        (335, 419, 'o'),
        (273, 390, 'v'),
        (251, 405, 'o'),
        (203, 423, 'o'),
        (102, 387, 'o'),
        (94, 321, 'o'),
        (143, 276, 'o'),
        (230, 251, 'o'),
        (260, 250, 'v'),
        (260, 220, '['),
        (258, 157, ']'),
        (243, 110, 'o'),
        (159, 99, 'o'),
        (141, 121, 'o'),
        (156, 158, 'o'),
        (121, 184, 'o'),
        (106, 117, 'o')]
        if 0:
            praw = [
            (275, 56, 'o'),
            (291, 75, 'v'),
            (312, 61, 'o'),
            (344, 50, 'o'),
            (373, 72, 'o'),
            (356, 91, 'o'),
            (334, 81, 'o'),
            (297, 85, 'v'),
            (306, 116, 'o'),
            (287, 180, 'o'),
            (236, 213, 'o'),
            (182, 212, 'o'),
            (157, 201, 'v'),
            (149, 209, 'o'),
            (143, 230, 'o'),
            (162, 246, 'c'),
            (202, 252, 'o'),
            (299, 260, 'o'),
            (331, 282, 'o'),
            (341, 341, 'o'),
            (308, 390, 'o'),
            (258, 417, 'o'),
            (185, 422, 'o'),
            (106, 377, 'o'),
            (110, 325, 'o'),
            (133, 296, 'o'),
            (147, 283, 'v'),
            (117, 238, 'o'),
            (133, 205, 'o'),
            (147, 191, 'v'),
            (126, 159, 'o'),
            (128, 94, 'o'),
            (167, 50, 'o'),
            (237, 39, 'o')]
            
        path = []
        for x, y, ty in praw:
            #if ty == 'o': ty = 'c'
            path.append((x, 550 - y, ty))

    if 0:
        path = [(100, 300, 'o'), (300, 100, 'o'), (300, 500, 'o')]

    if 0:
        # Woodford data set
        ty = 'o'
        praw = [(0, 0, '{'), (1, 1.9, ty), (2, 2.7, ty), (3, 2.6, ty),
                (4, 1.6, ty), (5, .89, ty), (6, 1.2, '}')]
        path = []
        for x, y, t in praw:
            path.append((100 + 80 * x, 100 + 80 * y, t))

    if 0:
        ycen = 32
        yrise = 0
        path = []
        ty = '{'
        for i in range(10):
            path.append((50 + i * 30, 250 + (10 - i) * yrise, ty))
            ty = 'o'
        path.append((350, 250 + ycen, ty))
        for i in range(1, 10):
            path.append((350 + i * 30, 250 + i * yrise, ty))
        path.append((650, 250 + 10 * yrise, '}'))

    prologue()

    run_path(path, show_iter = True, n=5)