bigmat.py 21 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662
  1. # Solver based on direct Newton solving of 4 parameters for each curve
  2. # segment
  3. import sys
  4. from math import *
  5. from Numeric import *
  6. import LinearAlgebra as la
  7. import poly3
  8. import band
  9. class Seg:
  10. def __init__(self, chord, th):
  11. self.ks = [0., 0., 0., 0.]
  12. self.chord = chord
  13. self.th = th
  14. def compute_ends(self, ks):
  15. chord, ch_th = poly3.integ_chord(ks)
  16. l = chord / self.chord
  17. thl = ch_th - (-.5 * ks[0] + .125 * ks[1] - 1./48 * ks[2] + 1./384 * ks[3])
  18. thr = (.5 * ks[0] + .125 * ks[1] + 1./48 * ks[2] + 1./384 * ks[3]) - ch_th
  19. k0l = l * (ks[0] - .5 * ks[1] + .125 * ks[2] - 1./48 * ks[3])
  20. k0r = l * (ks[0] + .5 * ks[1] + .125 * ks[2] + 1./48 * ks[3])
  21. l2 = l * l
  22. k1l = l2 * (ks[1] - .5 * ks[2] + .125 * ks[3])
  23. k1r = l2 * (ks[1] + .5 * ks[2] + .125 * ks[3])
  24. l3 = l2 * l
  25. k2l = l3 * (ks[2] - .5 * ks[3])
  26. k2r = l3 * (ks[2] + .5 * ks[3])
  27. return (thl, k0l, k1l, k2l), (thr, k0r, k1r, k2r), l
  28. def set_ends_from_ks(self):
  29. self.endl, self.endr, self.l = self.compute_ends(self.ks)
  30. def fast_pderivs(self):
  31. l = self.l
  32. l2 = l * l
  33. l3 = l2 * l
  34. return [((.5, l, 0, 0), (.5, l, 0, 0)),
  35. ((-1./12, -l/2, l2, 0), (1./12, l/2, l2, 0)),
  36. ((1./48, l/8, -l2/2, l3), (1./48, l/8, l2/2, l3)),
  37. ((-1./480, -l/48, l2/8, -l3/2), (1./480, l/48, l2/8, l3/2))]
  38. def compute_pderivs(self):
  39. rd = 2e6
  40. delta = 1./rd
  41. base_ks = self.ks
  42. base_endl, base_endr, dummy = self.compute_ends(base_ks)
  43. result = []
  44. for i in range(4):
  45. try_ks = base_ks[:]
  46. try_ks[i] += delta
  47. try_endl, try_endr, dummy = self.compute_ends(try_ks)
  48. deriv_l = (rd * (try_endl[0] - base_endl[0]),
  49. rd * (try_endl[1] - base_endl[1]),
  50. rd * (try_endl[2] - base_endl[2]),
  51. rd * (try_endl[3] - base_endl[3]))
  52. deriv_r = (rd * (try_endr[0] - base_endr[0]),
  53. rd * (try_endr[1] - base_endr[1]),
  54. rd * (try_endr[2] - base_endr[2]),
  55. rd * (try_endr[3] - base_endr[3]))
  56. result.append((deriv_l, deriv_r))
  57. return result
  58. class Node:
  59. def __init__(self, x, y, ty, th):
  60. self.x = x
  61. self.y = y
  62. self.ty = ty
  63. self.th = th
  64. def continuity(self):
  65. if self.ty == 'o':
  66. return 4
  67. elif self.ty in ('c', '[', ']'):
  68. return 2
  69. else:
  70. return 0
  71. def mod_2pi(th):
  72. u = th / (2 * pi)
  73. return 2 * pi * (u - floor(u + 0.5))
  74. def setup_path(path):
  75. segs = []
  76. nodes = []
  77. nsegs = len(path)
  78. if path[0][2] == '{':
  79. nsegs -= 1
  80. for i in range(nsegs):
  81. i1 = (i + 1) % len(path)
  82. x0, y0, t0 = path[i]
  83. x1, y1, t1 = path[i1]
  84. s = Seg(hypot(y1 - y0, x1 - x0), atan2(y1 - y0, x1 - x0))
  85. segs.append(s)
  86. for i in range(len(path)):
  87. x0, y0, t0 = path[i]
  88. if t0 in ('{', '}', 'v'):
  89. th = 0
  90. else:
  91. s0 = segs[(i + len(path) - 1) % len(path)]
  92. s1 = segs[i]
  93. th = mod_2pi(s1.th - s0.th)
  94. n = Node(x0, y0, t0, th)
  95. nodes.append(n)
  96. return segs, nodes
  97. def count_vec(nodes):
  98. jincs = []
  99. n = 0
  100. for i in range(len(nodes)):
  101. i1 = (i + 1) % len(nodes)
  102. t0 = nodes[i].ty
  103. t1 = nodes[i1].ty
  104. if t0 in ('{', '}', 'v', '[') and t1 in ('{', '}', 'v', ']'):
  105. jinc = 0
  106. elif t0 in ('{', '}', 'v', '[') and t1 == 'c':
  107. jinc = 1
  108. elif t0 == 'c' and t1 in ('{', '}', 'v', ']'):
  109. jinc = 1
  110. elif t0 == 'c' and t1 == 'c':
  111. jinc = 2
  112. else:
  113. jinc = 4
  114. jincs.append(jinc)
  115. n += jinc
  116. return n, jincs
  117. thscale, k0scale, k1scale, k2scale = 1, 1, 1, 1
  118. def inversedot_woodbury(m, v):
  119. a = zeros((n, 11), Float)
  120. for i in range(n):
  121. for j in range(max(-7, -i), min(4, n - i)):
  122. a[i, j + 7] = m[i, i + j]
  123. print a
  124. al, indx, d = band.bandec(a, 7, 3)
  125. VtZ = identity(4, Float)
  126. Z = zeros((n, 4), Float)
  127. for i in range(4):
  128. u = zeros(n, Float)
  129. for j in range(4):
  130. u[j] = m[j, n - 4 + i]
  131. band.banbks(a, 7, 3, al, indx, u)
  132. for k in range(n):
  133. Z[k, i] = u[k]
  134. #Z[:,i] = u
  135. for j in range(4):
  136. VtZ[j, i] += u[n - 4 + j]
  137. print Z
  138. print VtZ
  139. H = la.inverse(VtZ)
  140. print H
  141. band.banbks(a, 7, 3, al, indx, v)
  142. return(v - dot(Z, dot(H, v[n - 4:])))
  143. def inversedot(m, v):
  144. return dot(la.inverse(m), v)
  145. n, nn = m.shape
  146. if 1:
  147. for i in range(n):
  148. sys.stdout.write('% ')
  149. for j in range(n):
  150. if m[i, j] > 0: sys.stdout.write('+ ')
  151. elif m[i, j] < 0: sys.stdout.write('- ')
  152. else: sys.stdout.write(' ')
  153. sys.stdout.write('\n')
  154. cyclic = False
  155. for i in range(4):
  156. for j in range(n - 4, n):
  157. if m[i, j] != 0:
  158. cyclic = True
  159. print '% cyclic:', cyclic
  160. if not cyclic:
  161. a = zeros((n, 11), Float)
  162. for i in range(n):
  163. for j in range(max(-5, -i), min(6, n - i)):
  164. a[i, j + 5] = m[i, i + j]
  165. for i in range(n):
  166. sys.stdout.write('% ')
  167. for j in range(11):
  168. if a[i, j] > 0: sys.stdout.write('+ ')
  169. elif a[i, j] < 0: sys.stdout.write('- ')
  170. else: sys.stdout.write(' ')
  171. sys.stdout.write('\n')
  172. al, indx, d = band.bandec(a, 5, 5)
  173. print a
  174. band.banbks(a, 5, 5, al, indx, v)
  175. return v
  176. else:
  177. #return inversedot_woodbury(m, v)
  178. bign = 3 * n
  179. a = zeros((bign, 11), Float)
  180. u = zeros(bign, Float)
  181. for i in range(bign):
  182. u[i] = v[i % n]
  183. for j in range(-7, 4):
  184. a[i, j + 7] = m[i % n, (i + j + 7 * n) % n]
  185. #print a
  186. if 1:
  187. for i in range(bign):
  188. sys.stdout.write('% ')
  189. for j in range(11):
  190. if a[i, j] > 0: sys.stdout.write('+ ')
  191. elif a[i, j] < 0: sys.stdout.write('- ')
  192. else: sys.stdout.write(' ')
  193. sys.stdout.write('\n')
  194. #print u
  195. al, indx, d = band.bandec(a, 5, 5)
  196. band.banbks(a, 5, 5, al, indx, u)
  197. #print u
  198. return u[n + 2: 2 * n + 2]
  199. def iter(segs, nodes):
  200. n, jincs = count_vec(nodes)
  201. print '%', jincs
  202. v = zeros(n, Float)
  203. m = zeros((n, n), Float)
  204. for i in range(len(segs)):
  205. segs[i].set_ends_from_ks()
  206. j = 0
  207. j0 = 0
  208. for i in range(len(segs)):
  209. i1 = (i + 1) % len(nodes)
  210. t0 = nodes[i].ty
  211. t1 = nodes[i1].ty
  212. seg = segs[i]
  213. derivs = seg.compute_pderivs()
  214. print '%derivs:', derivs
  215. jinc = jincs[i] # the number of params on this seg
  216. print '%', t0, t1, jinc, j0
  217. # The constraints are laid out as follows:
  218. # constraints that cross the node on the left
  219. # constraints on the left side
  220. # constraints on the right side
  221. # constraints that cross the node on the right
  222. jj = j0 # the index into the constraint row we're writing
  223. jthl, jk0l, jk1l, jk2l = -1, -1, -1, -1
  224. jthr, jk0r, jk1r, jk2r = -1, -1, -1, -1
  225. # constraints crossing left
  226. if t0 == 'o':
  227. jthl = jj + 0
  228. jk0l = jj + 1
  229. jk1l = jj + 2
  230. jk2l = jj + 3
  231. jj += 4
  232. elif t0 in ('c', '[', ']'):
  233. jthl = jj + 0
  234. jk0l = jj + 1
  235. jj += 2
  236. # constraints on left
  237. if t0 in ('[', 'v', '{') and jinc == 4:
  238. jk1l = jj
  239. jj += 1
  240. if t0 in ('[', 'v', '{', 'c') and jinc == 4:
  241. jk2l = jj
  242. jj += 1
  243. # constraints on right
  244. if t1 in (']', 'v', '}') and jinc == 4:
  245. jk1r = jj
  246. jj += 1
  247. if t1 in (']', 'v', '}', 'c') and jinc == 4:
  248. jk2r = jj
  249. jj += 1
  250. # constraints crossing right
  251. jj %= n
  252. j1 = jj
  253. if t1 == 'o':
  254. jthr = jj + 0
  255. jk0r = jj + 1
  256. jk1r = jj + 2
  257. jk2r = jj + 3
  258. jj += 4
  259. elif t1 in ('c', '[', ']'):
  260. jthr = jj + 0
  261. jk0r = jj + 1
  262. jj += 2
  263. print '%', jthl, jk0l, jk1l, jk2l, jthr, jk0r, jk1r, jk2r
  264. if jthl >= 0:
  265. v[jthl] += thscale * (nodes[i].th - seg.endl[0])
  266. if jinc == 1:
  267. m[jthl][j] += derivs[0][0][0]
  268. elif jinc == 2:
  269. m[jthl][j + 1] += derivs[0][0][0]
  270. m[jthl][j] += derivs[1][0][0]
  271. elif jinc == 4:
  272. m[jthl][j + 2] += derivs[0][0][0]
  273. m[jthl][j + 3] += derivs[1][0][0]
  274. m[jthl][j + 0] += derivs[2][0][0]
  275. m[jthl][j + 1] += derivs[3][0][0]
  276. if jk0l >= 0:
  277. v[jk0l] += k0scale * seg.endl[1]
  278. if jinc == 1:
  279. m[jk0l][j] -= derivs[0][0][1]
  280. elif jinc == 2:
  281. m[jk0l][j + 1] -= derivs[0][0][1]
  282. m[jk0l][j] -= derivs[1][0][1]
  283. elif jinc == 4:
  284. m[jk0l][j + 2] -= derivs[0][0][1]
  285. m[jk0l][j + 3] -= derivs[1][0][1]
  286. m[jk0l][j + 0] -= derivs[2][0][1]
  287. m[jk0l][j + 1] -= derivs[3][0][1]
  288. if jk1l >= 0:
  289. v[jk1l] += k1scale * seg.endl[2]
  290. m[jk1l][j + 2] -= derivs[0][0][2]
  291. m[jk1l][j + 3] -= derivs[1][0][2]
  292. m[jk1l][j + 0] -= derivs[2][0][2]
  293. m[jk1l][j + 1] -= derivs[3][0][2]
  294. if jk2l >= 0:
  295. v[jk2l] += k2scale * seg.endl[3]
  296. m[jk2l][j + 2] -= derivs[0][0][3]
  297. m[jk2l][j + 3] -= derivs[1][0][3]
  298. m[jk2l][j + 0] -= derivs[2][0][3]
  299. m[jk2l][j + 1] -= derivs[3][0][3]
  300. if jthr >= 0:
  301. v[jthr] -= thscale * seg.endr[0]
  302. if jinc == 1:
  303. m[jthr][j] += derivs[0][1][0]
  304. elif jinc == 2:
  305. m[jthr][j + 1] += derivs[0][1][0]
  306. m[jthr][j + 0] += derivs[1][1][0]
  307. elif jinc == 4:
  308. m[jthr][j + 2] += derivs[0][1][0]
  309. m[jthr][j + 3] += derivs[1][1][0]
  310. m[jthr][j + 0] += derivs[2][1][0]
  311. m[jthr][j + 1] += derivs[3][1][0]
  312. if jk0r >= 0:
  313. v[jk0r] -= k0scale * seg.endr[1]
  314. if jinc == 1:
  315. m[jk0r][j] += derivs[0][1][1]
  316. elif jinc == 2:
  317. m[jk0r][j + 1] += derivs[0][1][1]
  318. m[jk0r][j + 0] += derivs[1][1][1]
  319. elif jinc == 4:
  320. m[jk0r][j + 2] += derivs[0][1][1]
  321. m[jk0r][j + 3] += derivs[1][1][1]
  322. m[jk0r][j + 0] += derivs[2][1][1]
  323. m[jk0r][j + 1] += derivs[3][1][1]
  324. if jk1r >= 0:
  325. v[jk1r] -= k1scale * seg.endr[2]
  326. m[jk1r][j + 2] += derivs[0][1][2]
  327. m[jk1r][j + 3] += derivs[1][1][2]
  328. m[jk1r][j + 0] += derivs[2][1][2]
  329. m[jk1r][j + 1] += derivs[3][1][2]
  330. if jk2r >= 0:
  331. v[jk2r] -= k2scale * seg.endr[3]
  332. m[jk2r][j + 2] += derivs[0][1][3]
  333. m[jk2r][j + 3] += derivs[1][1][3]
  334. m[jk2r][j + 0] += derivs[2][1][3]
  335. m[jk2r][j + 1] += derivs[3][1][3]
  336. j += jinc
  337. j0 = j1
  338. #print m
  339. dk = inversedot(m, v)
  340. dkmul = 1
  341. j = 0
  342. for i in range(len(segs)):
  343. jinc = jincs[i]
  344. if jinc == 1:
  345. segs[i].ks[0] += dkmul * dk[j]
  346. elif jinc == 2:
  347. segs[i].ks[0] += dkmul * dk[j + 1]
  348. segs[i].ks[1] += dkmul * dk[j + 0]
  349. elif jinc == 4:
  350. segs[i].ks[0] += dkmul * dk[j + 2]
  351. segs[i].ks[1] += dkmul * dk[j + 3]
  352. segs[i].ks[2] += dkmul * dk[j + 0]
  353. segs[i].ks[3] += dkmul * dk[j + 1]
  354. j += jinc
  355. norm = 0.
  356. for i in range(len(dk)):
  357. norm += dk[i] * dk[i]
  358. return sqrt(norm)
  359. def plot_path(segs, nodes, tol = 1.0, show_cpts = False):
  360. if show_cpts:
  361. cpts = []
  362. j = 0
  363. cmd = 'moveto'
  364. for i in range(len(segs)):
  365. i1 = (i + 1) % len(nodes)
  366. n0 = nodes[i]
  367. n1 = nodes[i1]
  368. x0, y0, t0 = n0.x, n0.y, n0.ty
  369. x1, y1, t1 = n1.x, n1.y, n1.ty
  370. ks = segs[i].ks
  371. abs_ks = abs(ks[0]) + abs(ks[1] / 2) + abs(ks[2] / 8) + abs(ks[3] / 48)
  372. n_subdiv = int(ceil(.001 + tol * abs_ks))
  373. n_subhalf = (n_subdiv + 1) / 2
  374. if n_subdiv > 1:
  375. n_subdiv = n_subhalf * 2
  376. ksp = (ks[0] * .5, ks[1] * .25, ks[2] * .125, ks[3] * .0625)
  377. pside = poly3.int_3spiro_poly(ksp, n_subhalf)
  378. ksm = (ks[0] * -.5, ks[1] * .25, ks[2] * -.125, ks[3] * .0625)
  379. mside = poly3.int_3spiro_poly(ksm, n_subhalf)
  380. mside.reverse()
  381. for j in range(len(mside)):
  382. mside[j] = (-mside[j][0], -mside[j][1])
  383. if n_subdiv > 1:
  384. pts = mside + pside[1:]
  385. else:
  386. pts = mside[:1] + pside[1:]
  387. chord_th = atan2(y1 - y0, x1 - x0)
  388. chord_len = hypot(y1 - y0, x1 - x0)
  389. rot = chord_th - atan2(pts[-1][1] - pts[0][1], pts[-1][0] - pts[0][0])
  390. scale = chord_len / hypot(pts[-1][1] - pts[0][1], pts[-1][0] - pts[0][0])
  391. u, v = scale * cos(rot), scale * sin(rot)
  392. xt = x0 - u * pts[0][0] + v * pts[0][1]
  393. yt = y0 - u * pts[0][1] - v * pts[0][0]
  394. s = -.5
  395. for x, y in pts:
  396. xp, yp = xt + u * x - v * y, yt + u * y + v * x
  397. thp = (((ks[3]/24 * s + ks[2]/6) * s + ks[1] / 2) * s + ks[0]) * s + rot
  398. up, vp = scale / (1.5 * n_subdiv) * cos(thp), scale / (1.5 * n_subdiv) * sin(thp)
  399. if s == -.5:
  400. if cmd == 'moveto':
  401. print xp, yp, 'moveto'
  402. cmd = 'curveto'
  403. else:
  404. if show_cpts:
  405. cpts.append((xlast + ulast, ylast + vlast))
  406. cpts.append((xp - up, yp - vp))
  407. print xlast + ulast, ylast + vlast, xp - up, yp - vp, xp, yp, 'curveto'
  408. xlast, ylast, ulast, vlast = xp, yp, up, vp
  409. s += 1. / n_subdiv
  410. if t1 == 'v':
  411. j += 2
  412. else:
  413. j += 1
  414. print 'stroke'
  415. if show_cpts:
  416. for x, y in cpts:
  417. print 'gsave 0 0 1 setrgbcolor', x, y, 'translate circle fill grestore'
  418. def plot_ks(segs, nodes, xo, yo, xscale, yscale):
  419. j = 0
  420. cmd = 'moveto'
  421. x = xo
  422. for i in range(len(segs)):
  423. i1 = (i + 1) % len(nodes)
  424. n0 = nodes[i]
  425. n1 = nodes[i1]
  426. x0, y0, t0 = n0.x, n0.y, n0.ty
  427. x1, y1, t1 = n1.x, n1.y, n1.ty
  428. ks = segs[i].ks
  429. chord, ch_th = poly3.integ_chord(ks)
  430. l = chord/segs[i].chord
  431. k0 = l * poly3.eval_cubic(ks[0], ks[1], .5 * ks[2], 1./6 * ks[3], -.5)
  432. print x, yo + yscale * k0, cmd
  433. cmd = 'lineto'
  434. k3 = l * poly3.eval_cubic(ks[0], ks[1], .5 * ks[2], 1./6 * ks[3], .5)
  435. k1 = k0 + l/3 * (ks[1] - 0.5 * ks[2] + .125 * ks[3])
  436. k2 = k3 - l/3 * (ks[1] + 0.5 * ks[2] + .125 * ks[3])
  437. print x + xscale / l / 3., yo + yscale * k1
  438. print x + 2 * xscale / l / 3., yo + yscale * k2
  439. print x + xscale / l, yo + yscale * k3, 'curveto'
  440. x += xscale / l
  441. if t1 == 'v':
  442. j += 2
  443. else:
  444. j += 1
  445. print 'stroke'
  446. print xo, yo, 'moveto', x, yo, 'lineto stroke'
  447. def plot_nodes(nodes, segs):
  448. for i in range(len(nodes)):
  449. n = nodes[i]
  450. print 'gsave', n.x, n.y, 'translate'
  451. if n.ty in ('c', '[', ']'):
  452. th = segs[i].th - segs[i].endl[0]
  453. if n.ty == ']': th += pi
  454. print 180 * th / pi, 'rotate'
  455. if n.ty == 'o':
  456. print 'circle fill'
  457. elif n.ty == 'c':
  458. print '3 4 poly fill'
  459. elif n.ty in ('v', '{', '}'):
  460. print '45 rotate 3 4 poly fill'
  461. elif n.ty in ('[', ']'):
  462. print '0 -3 moveto 0 0 3 90 270 arc fill'
  463. else:
  464. print '5 3 poly fill'
  465. print 'grestore'
  466. def prologue():
  467. print '/ss 2 def'
  468. print '/circle { ss 0 moveto currentpoint exch ss sub exch ss 0 360 arc } bind def'
  469. print '/poly {'
  470. print ' dup false exch {'
  471. print ' 0 3 index 2 index { lineto } { moveto } ifelse pop true'
  472. print ' 360.0 2 index div rotate'
  473. print ' } repeat pop pop pop'
  474. print '} bind def'
  475. def run_path(path, show_iter = False, n = 10, xo = 36, yo = 550, xscale = .25, yscale = 2000, pl_nodes = True):
  476. segs, nodes = setup_path(path)
  477. print '.5 setlinewidth'
  478. for i in range(n):
  479. if i == n - 1:
  480. print '0 0 0 setrgbcolor 1 setlinewidth'
  481. elif i == 0:
  482. print '1 0 0 setrgbcolor'
  483. elif i == 1:
  484. print '0 0.5 0 setrgbcolor'
  485. elif i == 2:
  486. print '0.3 0.3 0.3 setrgbcolor'
  487. norm = iter(segs, nodes)
  488. print '% norm =', norm
  489. if show_iter or i == n - 1:
  490. #print '1 0 0 setrgbcolor'
  491. #plot_path(segs, nodes, tol = 5)
  492. #print '0 0 0 setrgbcolor'
  493. plot_path(segs, nodes, tol = 1.0)
  494. plot_ks(segs, nodes, xo, yo, xscale, yscale)
  495. if pl_nodes:
  496. plot_nodes(nodes, segs)
  497. if __name__ == '__main__':
  498. if 1:
  499. path = [(100, 350, 'o'), (225, 350, 'o'), (350, 450, 'o'),
  500. (450, 400, 'o'), (315, 205, 'o'), (300, 200, 'o'),
  501. (285, 205, 'o')]
  502. if 1:
  503. path = [(100, 350, 'o'), (175, 375, '['), (250, 375, ']'), (325, 425, '['),
  504. (325, 450, ']'),
  505. (400, 475, 'o'), (350, 200, 'c')]
  506. if 0:
  507. ecc, ty, ty1 = 0.56199, 'c', 'c'
  508. ecc, ty, ty1 = 0.49076, 'o', 'o',
  509. ecc, ty, ty1 = 0.42637, 'o', 'c'
  510. path = [(300 - 200 * ecc, 300, ty), (300, 100, ty1),
  511. (300 + 200 * ecc, 300, ty), (300, 500, ty1)]
  512. # difficult path #3
  513. if 0:
  514. path = [(100, 300, '{'), (225, 350, 'o'), (350, 500, 'o'),
  515. (450, 500, 'o'), (450, 450, 'o'), (300, 200, '}')]
  516. if 0:
  517. path = [(100, 100, '{'), (200, 180, 'v'), (250, 215, 'o'),
  518. (300, 200, 'o'), (400, 100, '}')]
  519. if 0:
  520. praw = [
  521. (134, 90, 'o'),
  522. (192, 68, 'o'),
  523. (246, 66, 'o'),
  524. (307, 107, 'o'),
  525. (314, 154, '['),
  526. (317, 323, ']'),
  527. (347, 389, 'o'),
  528. (373, 379, 'v'),
  529. (386, 391, 'v'),
  530. (365, 409, 'o'),
  531. (335, 419, 'o'),
  532. (273, 390, 'v'),
  533. (251, 405, 'o'),
  534. (203, 423, 'o'),
  535. (102, 387, 'o'),
  536. (94, 321, 'o'),
  537. (143, 276, 'o'),
  538. (230, 251, 'o'),
  539. (260, 250, 'v'),
  540. (260, 220, '['),
  541. (258, 157, ']'),
  542. (243, 110, 'o'),
  543. (159, 99, 'o'),
  544. (141, 121, 'o'),
  545. (156, 158, 'o'),
  546. (121, 184, 'o'),
  547. (106, 117, 'o')]
  548. if 0:
  549. praw = [
  550. (275, 56, 'o'),
  551. (291, 75, 'v'),
  552. (312, 61, 'o'),
  553. (344, 50, 'o'),
  554. (373, 72, 'o'),
  555. (356, 91, 'o'),
  556. (334, 81, 'o'),
  557. (297, 85, 'v'),
  558. (306, 116, 'o'),
  559. (287, 180, 'o'),
  560. (236, 213, 'o'),
  561. (182, 212, 'o'),
  562. (157, 201, 'v'),
  563. (149, 209, 'o'),
  564. (143, 230, 'o'),
  565. (162, 246, 'c'),
  566. (202, 252, 'o'),
  567. (299, 260, 'o'),
  568. (331, 282, 'o'),
  569. (341, 341, 'o'),
  570. (308, 390, 'o'),
  571. (258, 417, 'o'),
  572. (185, 422, 'o'),
  573. (106, 377, 'o'),
  574. (110, 325, 'o'),
  575. (133, 296, 'o'),
  576. (147, 283, 'v'),
  577. (117, 238, 'o'),
  578. (133, 205, 'o'),
  579. (147, 191, 'v'),
  580. (126, 159, 'o'),
  581. (128, 94, 'o'),
  582. (167, 50, 'o'),
  583. (237, 39, 'o')]
  584. path = []
  585. for x, y, ty in praw:
  586. #if ty == 'o': ty = 'c'
  587. path.append((x, 550 - y, ty))
  588. if 0:
  589. path = [(100, 300, 'o'), (300, 100, 'o'), (300, 500, 'o')]
  590. if 0:
  591. # Woodford data set
  592. ty = 'o'
  593. praw = [(0, 0, '{'), (1, 1.9, ty), (2, 2.7, ty), (3, 2.6, ty),
  594. (4, 1.6, ty), (5, .89, ty), (6, 1.2, '}')]
  595. path = []
  596. for x, y, t in praw:
  597. path.append((100 + 80 * x, 100 + 80 * y, t))
  598. if 0:
  599. ycen = 32
  600. yrise = 0
  601. path = []
  602. ty = '{'
  603. for i in range(10):
  604. path.append((50 + i * 30, 250 + (10 - i) * yrise, ty))
  605. ty = 'o'
  606. path.append((350, 250 + ycen, ty))
  607. for i in range(1, 10):
  608. path.append((350 + i * 30, 250 + i * yrise, ty))
  609. path.append((650, 250 + 10 * yrise, '}'))
  610. prologue()
  611. run_path(path, show_iter = True, n=5)