<< JunoOldSolverImpl.mesa (was OldJunoSolverImpl.mesa)>> << >> <> <> <> <> << This is the solver for an interactive graphics program that allows a user to specify the positions of points by declaring geometrical contraints that the positions are to satisfy. The program solves constraints (which can be viewed as simultaneous polynomial equations) by an n-dimensional version of Newton's method for finding the root of a differentiable function.>> << TO FIX: >> << handling of inequalities: check code.>> DIRECTORY JunoStorage USING [Point, Coords, Frame, nullFrame, ItemList, Item, ConstrKind], JunoMatrix USING [Matrix, GetFrameMatrix, MapCoords], JunoOldSolver, Random USING [Next], JunoGraphics USING [DrawEdge], Real; JunoOldSolverImpl: PROGRAM IMPORTS JunoMatrix, Random, JunoGraphics EXPORTS JunoOldSolver = BEGIN OPEN Stor: JunoStorage, Mat: JunoMatrix, JunoOldSolver, Gr: JunoGraphics; << - - - - IMPORTED TYPES >> Point: TYPE = Stor.Point; Coords: TYPE = Stor.Coords; Item: TYPE = Stor.Item; Frame: TYPE = Stor.Frame; << - - - - GLOBAL VARIABLES >> displayingMotions: PUBLIC BOOL _ FALSE; loggingPivots: BOOL _ FALSE; << - - - - SIGNALS >> TooManyConstraints: SIGNAL = CODE; TooManyConstrainedPoints: SIGNAL = CODE; << - - - - GLOBAL VARIABLES >> rowmax: INTEGER = 100; colmax: INTEGER = 100; Tag: TYPE = {x, y, eq, geq}; TableauLabel: TYPE = RECORD [point: Point, tag: Tag]; Column: TYPE = ARRAY [1..rowmax] OF REAL; Array: TYPE = ARRAY [1..rowmax] OF ARRAY [1..colmax] OF REAL; a: REF Array _ NEW[Array]; <> b: REF Column _ NEW[Column]; <> rowhead: ARRAY[1..rowmax] OF TableauLabel; colhead: ARRAY[1..colmax] OF TableauLabel; n, m: INTEGER; -- the tableau a has n relevant rows and m relevant columns. mFrame: REF Mat.Matrix _ NEW [Mat.Matrix]; -- a work matrix << - - - - CONSTRAINT AND DERIVATIVE EVALUATION >> << Solving geometric constraints by the Newton-Raphson method>> << The geometric constraints given to the solver can be translated into a set of equalities and inequalities of the form fi(v1,v2,...,vm) = 0 or fi(v1,v2,...,vm) geq 0, i=1..n. To solve these equations, Juno uses the Newton-Raphson method, extended to handle inequalities. Given an approximate solution X=(X1, X2, ..., Xm), the next approximation is computed by replacing each fi by its Taylor expansion around X, truncated to the first-order terms:>> << fi(X1+x1, X2+x2, ..., Xm+xm) = fi(X1, X2, ..., Xm) + dfi1*x1 + dfi2*x2 + ... + dfim*xm + (higher order terms)>> << where dfij is the patial derivative of fi with respect to its jth argument vj, computed at the point X=(X1, X2, ..., Xm).>> << The next approximation (X1',X2',...,Xm') is then found by solving the system of equations and inequations>> << dfi1*x1 + dfi2*x2 + ... + dfim*xm + fi(X1, X2, ..., Xm) = 0 or dfi1*x1 + dfi2*x2 + ... + dfim*xm + fi(X1, X2, ..., Xm) geq 0>> << and letting Xj' = Xj+xj.>> << The procedure PrepareTableau below constructs the matrix a[i][j] = dfij(X1, X2, ..., Xm) and b[i] = fi(X1, X2, ..., Xm), that is input to the linear equation/inequation solver.>> PrepareTableau: PROC [items: Stor.ItemList, eps: REAL] RETURNS [success, impossible: BOOLEAN] = <> << Returns impossible=TRUE if it concludes it is impossible to satisfy the current constraints given the current positions of the fixed points.>> << Returns success=TRUE if all constraints are satisfied with the points as they are now.>> BEGIN << Procedures to enter constraints in the tableau:>> constant: BOOL; -- TRUE if last constraint depends only on fixed points. maxc: REAL; -- max point coord. (fixed or not) * its partial der. in last constraint maxp: REAL; -- max partial deriv (fixed or not) for last constraint. <<(Assuming all calls to Partial are for distinct points).>> Constraint: PROC[value: REAL, tag: Tag] = <> << Sets satisfied according to present status of constraint. Also sets unsatisfiable if not trivially satisfied; this flag may be reset by Partial.>> BEGIN n _ n + 1; -- add new row IF n > rowmax THEN SIGNAL TooManyConstraints[]; rowhead[n] _ [point: NIL, tag: tag]; FOR j: INT IN [1..m] DO a[n][j] _ 0.0 ENDLOOP; b[n] _ value; constant _ TRUE; -- for now maxp _ maxc _ 0.0 END; Partial: PROC[p: Point, xCoef, yCoef: REAL] = <> BEGIN GetCol: PROC[which: Tag[x..y]] RETURNS [col: INTEGER] = <> BEGIN col _ IF which=x THEN p.xCol ELSE p.yCol; IF col NOT IN [1..m] OR colhead[col].point # p THEN {-- add new column: m _ m + 1; IF m > colmax THEN SIGNAL TooManyConstrainedPoints[]; FOR i: INT IN [1..n] DO a[i][m] _ 0.0 ENDLOOP; IF which=x THEN {p.xCol _ m} ELSE {p.yCol _ m}; p.old _ p.coords; p.mark _ FALSE; -- used to plot motions colhead[m] _ [point: p, tag: which]; col _ m}; END; maxp _ MAX [maxp, ABS[xCoef], ABS[yCoef]]; maxc _ MAX [maxc, ABS [xCoef * p.coords.x]]; maxc _ MAX [maxc, ABS [yCoef * p.coords.y]]; IF NOT p.fixed THEN {constant _ FALSE; IF ABS[xCoef] # 0.0 THEN {col: INTEGER _ GetCol[x]; a[n][col] _ a[n][col] + xCoef}; IF ABS[yCoef] # 0.0 THEN {col: INTEGER _ GetCol[y]; a[n][col] _ a[n][col] + yCoef}}; END; << Procedures for specific constraints:>> EnterVer: PROC[frame: Frame, i, j: Point] = BEGIN IF frame.ver#NIL THEN -- same as (i, j) para (org, ver) {EnterPara[Stor.nullFrame, i, j, frame.org, frame.ver]} ELSE IF frame.hor#NIL THEN -- same as (i, j) perp (org, hor) {EnterPerp[Stor.nullFrame, i, j, frame.org, frame.ver]} ELSE -- simple verticality constraint {Constraint[i.coords.x - j.coords.x, eq]; Partial[i, 1, 0]; Partial[j, -1, 0]} END; EnterHor: PROC[frame: Frame, i, j: Point] = BEGIN IF frame.hor#NIL THEN -- same as (i, j) para (org, hor) {EnterPara[Stor.nullFrame, i, j, frame.org, frame.hor]} ELSE -- simple horizontality constraint {Constraint[i.coords.y - j.coords.y, eq]; Partial[i, 0, 1]; Partial[j, 0, -1]} END; EnterPara: PROC[frame: Frame, i, j, k, l: Point] = BEGIN ijx: REAL = j.coords.x - i.coords.x; -- displacement vector from i to j. ijy: REAL = j.coords.y - i.coords.y; klx: REAL = l.coords.x - k.coords.x; -- displacement vector from k to l. kly: REAL = l.coords.y - k.coords.y; << frame is irrelevant>> Constraint[-ijx*kly + klx*ijy, eq]; Partial[i, kly, -klx]; Partial[j, -kly, klx]; Partial[k, -ijy, ijx]; Partial[l, ijy, -ijx] END; EnterCong: PROC[frame: Frame, i, j, k, l: Point] = BEGIN ijx: REAL = j.coords.x - i.coords.x; -- displacement vector from i to j. ijy: REAL = j.coords.y - i.coords.y; klx: REAL = l.coords.x - k.coords.x; -- displacement vector from k to l. kly: REAL = l.coords.y - k.coords.y; IF frame.hor#NIL AND frame.ver#NIL THEN -- relativized congruence {ux: REAL = frame.hor.coords.x - frame.org.coords.x; -- frame vectors. uy: REAL = frame.hor.coords.y - frame.org.coords.y; vx: REAL = frame.ver.coords.x - frame.org.coords.x; vy: REAL = frame.ver.coords.y - frame.org.coords.y; << The following computations essentially multiply the coordinates of the constrained points by the inverse of the frame matrix, except that the division by the frame's determinant is left out.>> << This still preserves the zeros of the constraint, and is MUCH easier to compute and differentiate. It also avoids division by zero if the frame becomes singular. However, in that case it will consider the constraint as being satisfied.>> iju: REAL = vy*ijx - vx*ijy; -- coords of (j-i) rel frame. ijv: REAL = -uy*ijx + ux*ijy; klu: REAL = vy*klx - vx*kly; -- coords of (l-k) rel frame klv: REAL = -uy*klx + ux*kly; value: REAL = 0.5 * (iju*iju + ijv*ijv- klu*klu - klv*klv); dijx: REAL = vy*iju - uy*ijv; -- der. of value w.r.t. (j.x-i.x) dijy: REAL = -vx*iju + ux*ijv; dklx: REAL = -vy*klu + uy*klv; dkly: REAL = vx*klu - ux*klv; dvx: REAL = -iju*ijy + klu*kly; -- der. w.r.t. vx dvy: REAL = iju*ijx - klu*klx; dux: REAL = ijv*ijy - klv*kly; -- der. w.r.t. ux duy: REAL = -ijv*ijx + klv*klx; Constraint[value, eq]; Partial[i, -dijx, -dijy]; Partial[j, dijx, dijy]; Partial[k, -dklx, -dkly]; Partial[l, dklx, dkly]; Partial[frame.hor, dux, duy]; Partial[frame.ver, dvx, dvy]; Partial[frame.org, -dux-dvx, -duy-dvy]} -- Phew! ELSE -- frame is orthogonal; absolute congruence {Constraint[0.5 * (ijx*ijx + ijy*ijy - klx*klx - kly*kly), eq]; Partial[i, -ijx, -ijy]; Partial[j, ijx, ijy]; Partial[k, klx, kly]; Partial[l, -klx, -kly]} END; EnterPerp: PROC[frame: Frame, i, j, k, l: Point] = BEGIN ijx: REAL = j.coords.x - i.coords.x; -- displacement vector from i to j. ijy: REAL = j.coords.y - i.coords.y; klx: REAL = l.coords.x - k.coords.x; -- displacement vector from k to l. kly: REAL = l.coords.y - k.coords.y; IF frame.hor#NIL AND frame.ver#NIL THEN -- relativized perpendicularity {ux: REAL = frame.hor.coords.x - frame.org.coords.x; -- frame vectors. uy: REAL = frame.hor.coords.y - frame.org.coords.y; vx: REAL = frame.ver.coords.x - frame.org.coords.x; vy: REAL = frame.ver.coords.y - frame.org.coords.y; << The following computations essentially multiply the coordinates of the constrained points by the inverse of the frame matrix, except that the division by the frame's determinant is left out.>> << This still preserves the zeros of the constraint, and is MUCH easier to compute and differentiate. It also avoids division by zero if the frame becomes singular. However, in that case it will consider the constraint as being satisfied.>> iju: REAL = vy*ijx - vx*ijy; -- coords of (j-i) rel frame. ijv: REAL = -uy*ijx + ux*ijy; klu: REAL = vy*klx - vx*kly; -- coords of (l-k) rel frame klv: REAL = -uy*klx + ux*kly; value: REAL = iju*klu+ijv*klv; -- constraint value dijx: REAL = vy*klu - uy*klv; -- der. of value w.r.t. (j.x-i.x) dijy: REAL = -vx*klu + ux*klv; dklx: REAL = iju*vy - ijv*uy; -- der. of value w.r.t. (l.x-k.x) dkly: REAL = -iju*vx + ijv*ux; dvx: REAL = -ijy*klu - iju*kly; -- der. w.r.t. vx dvy: REAL = ijx*klu + iju*klx; dux: REAL = -ijy*klv - ijv*kly; -- der. w.r.t. ux duy: REAL = ijx*klv + ijv*klx; Constraint[value, eq]; Partial[i, -dijx, -dijy]; Partial[j, dijx, dijy]; Partial[k, -dklx, -dkly]; Partial[l, dklx, dkly]; Partial[frame.hor, dux, duy]; Partial[frame.ver, dvx, dvy]; Partial[frame.org, -dux-dvx, -duy-dvy]} -- Phew! ELSE -- frame is orthogonal; absolute perpendicularity {Constraint[ijx*klx + ijy*kly, eq]; Partial[i, -klx, -kly]; Partial[j, klx, kly]; Partial[k, -ijx, -ijy]; Partial[l, ijx, ijy]} END; EnterAt: PROC[frame: Frame, pp: Point, c: Coords] = BEGIN cT: Coords; IF frame # Stor.nullFrame THEN {mFrame _ Mat.GetFrameMatrix [frame, mFrame]; cT _ Mat.MapCoords[c, mFrame]} ELSE {cT _ c}; Constraint[cT.x - pp.coords.x, eq]; Partial[pp, -1, 0]; IF frame # Stor.nullFrame THEN -- enter derivatives with respect to frame {Partial[frame.org, 1, 0]; IF frame.hor # NIL THEN {Partial[frame.hor, c.x, 0]; Partial[frame.org, -c.x, 0]; IF frame.ver # NIL THEN {Partial[frame.ver, c.y, 0]; Partial[frame.org, -c.y, 0]} ELSE {Partial[frame.hor, 0, -c.y]; Partial[frame.org, 0, c.y]}}}; Constraint[cT.y - pp.coords.y, eq]; Partial[pp, 0, -1]; IF frame # Stor.nullFrame THEN -- enter derivatives with respect to frame {Partial[frame.org, 0, 1]; IF frame.hor # NIL THEN {Partial[frame.hor, 0, c.x]; Partial[frame.org, 0, -c.x]; IF frame.ver # NIL THEN {Partial[frame.ver, 0, c.y]; Partial[frame.org, 0, -c.y]} ELSE {Partial[frame.hor, c.y, 0]; Partial[frame.org, -c.y, 0]}}} END; EnterCcw: PROC[frame: Frame, i, j, k: Point] = BEGIN ix: REAL = i.coords.x; iy: REAL = i.coords.y; jx: REAL = j.coords.x; jy: REAL = j.coords.y; kx: REAL = k.coords.x; ky: REAL = k.coords.y; detijk: REAL = (jy-iy)*(kx-ix) - (jx-ix)*(ky-iy); IF frame.hor#NIL AND frame.ver#NIL THEN -- relativized right-handedness {ox: REAL = frame.org.coords.x; oy: REAL =frame.org.coords.y; hx: REAL = frame.hor.coords.x; hy: REAL = frame.hor.coords.y; vx: REAL = frame.ver.coords.x; vy: REAL = frame.ver.coords.y; detframe: REAL = (hy-oy)*(vx-ox)-(hx-ox)*(vy-oy); Constraint[detframe*detijk, geq]; Partial[i, detframe*(ky-jy), detframe*(jx-kx)]; Partial[j, detframe*(iy-ky), detframe*(kx-ix)]; Partial[k, detframe*(jy-iy), detframe*(ix-jx)]; Partial[frame.org, detijk*(vy-hy), detijk*(hx-vx)]; Partial[frame.hor, detijk*(oy-ky), detijk*(vx-ox)]; Partial[frame.ver, detijk*(hy-oy), detijk*(ox-hx)]} ELSE -- frame is right-handed, absolute right-handedness {Constraint[detijk, geq]; Partial[i, ky-jy, jx-kx]; Partial[j, iy-ky, kx-ix]; Partial[k, jy-iy, ix-jx]} END; << Other useful procs:>> Error: PROC[row: INTEGER] RETURNS [error: REAL] = <> BEGIN RETURN[SELECT rowhead[row].tag FROM eq => ABS[b[row]], geq => MAX [0.0, -b[row]], ENDCASE => ERROR] END; P1: PROC[args: LIST OF REF ANY] RETURNS [p: Point] = INLINE {RETURN[NARROW[args.first]]}; P2: PROC[args: LIST OF REF ANY] RETURNS [p: Point] = INLINE {RETURN[NARROW[args.rest.first]]}; P3: PROC[args: LIST OF REF ANY] RETURNS [p: Point] = INLINE {RETURN[NARROW[args.rest.rest.first]]}; P4: PROC[args: LIST OF REF ANY] RETURNS [p: Point] = INLINE {RETURN[NARROW[args.rest.rest.rest.first]]}; << Fill in tableau a[][] and constant column b[] by enumerating current constraints:>> success _ TRUE; impossible _ FALSE; n _ 0; m _ 0; FOR item: Item _ items.first, item.link UNTIL item = NIL OR impossible DO IF item.kind IN Stor.ConstrKind THEN BEGIN args: LIST OF REF ANY = item.args; SELECT item.kind FROM hor => EnterHor[item.frame, P1[args], P2[args]]; ver => EnterVer[item.frame, P1[args], P2[args]]; para => EnterPara[item.frame, P1[args], P2[args], P3[args], P4[args]]; cong => EnterCong[item.frame, P1[args], P2[args], P3[args], P4[args]]; perp => EnterPerp[item.frame, P1[args], P2[args], P3[args], P4[args]]; at => EnterAt [item.frame, P1[args], [NARROW[args.rest.first, REF REAL]^, NARROW[args.rest.rest.first, REF REAL]^]]; ccw => EnterCcw[item.frame, P1[args], P2[args], P3[args]]; ENDCASE => ERROR; IF ABS[b[n]] < MAX[1e-28, 1e-6*maxc] THEN b[n] _ 0.0; IF Error[n] > eps * maxp THEN {success _ FALSE; IF constant THEN impossible _ TRUE}; IF constant THEN {n _ n-1} END; ENDLOOP; END; << - - - - LINEAR SYSTEM SOLVER >> << The linear system solver:>> << Before and after each pivoting step, the tableau (represented by the arrray a[][], the column vector b[], and the header vectors rowhad[] and colhead[]) determines a set of linear equalities and inequalities among n+m scalar variables: the independent ones x1, x2, ..., xm, each corresponding to one of the m columns of a, and the dependent ones, y1, y2, ..., yn, associated with the rows of a. The LinSolve procedure tries to determine the values of those variables so as to satisfy those equalities and inequalities. >> << Each variable can be either a point coordinate or a "slack variable". The former is the x- or y-coordinate of some undetermined point p, and is denoted by the corresponding element of rowhead or colhead being of the form [p, x] or [p, y]. A slack variable represents the numerical value of some constraint; for example, the constraint hor (p, q) has numerical value p.y-q.y. The TableauLabel for a slack variable is either [NIL, eq] or [NIL, geq], depending on whether its value must be zero or non-negative for the constraint to be satisfied. For now, geq labels are generated only by the ccw constraint.>> << At the start of the Newton-Raphson step, all point coordinates are independent variables, and all slack variables are dependent, but this will in general not be true during and after the algorithm. At all times, the tableau can be read as follows: for each row i, we have one equation of the form yi = a[i][1].x1 + a[i][2].x2 + ... + a[i][m].xm + b[i]. In addition, for each variable z (dependent or independent) with label [NIL, eq] we have the equation z=0, and for each one with label [NIL, geq] we have the inequality z geq 0.>> << The goal of the algorithm is to transform the original tableau into an equivalent one with the property that every row constraint is trivially satisfiable. That is, for every dependent slack variable yi we must have b[i]=0 if rowhead[i].tag=eq, and b[i] geq 0 if rowhead[i].tag=geq. Such a tableau admits the trivial solution xj=0 for j=1..m, and yi=b[i] for i=1..n.>> << The algorihm achieves this by a sequence of pivoting steps. A pivoting step consists in rewriting the system so that some independent variable xj becomes dependent, and some yi becomes independent (the element a[i][j] is called the pivot for that step). This changes the matrix a and the vector b, so it may change the satisfiablity of other row constraints. The trick is to choose the pivot so that none of the the rows that are already trivially satisfiable ceases to be so, and so that their number is bound to increase after a sufficient number of pivoting steps.>> LinSolve: PROC [eps: REAL] RETURNS [moved: BOOLEAN] = <> << Returns moved = TRUE if the coordinates of one or more points have changed by more than eps. >> << Th points that moved so are those in rowhead[] with p.mark = true. >> BEGIN << Invert matrix:>> BEGIN DO r, c: INTEGER _ 0; << Look for a suitable pivot a[r][c] >> FOR i: INTEGER IN [1..n] WHILE r = 0 DO champ, sum: REAL _ 0.0; << Ignore dependent variables that are unconstrained or trivially satisfiable >> IF rowhead[i].tag = eq AND b[i] = 0.0 OR rowhead[i].tag = geq AND b[i] >= 0.0 OR rowhead[i].tag IN [x..y] THEN LOOP; << Find independent variable xc that can be changed so as to relieve the constraint on yi, and that requires the minimum such change (i.e. maximizes ABS[a[i][c]]) >> FOR j: INTEGER IN [1 .. m] DO tent: REAL = ABS[a[i][j]]; IF tent >= champ AND (colhead[j].tag IN [x .. y] OR colhead[j].tag = geq AND a[i][j] > 0.0) THEN {c _ j; champ _ tent}; sum _ sum+tent ENDLOOP; << Is a[i][c] large enough to pivot on it?>> IF c = 0 OR champ < .001*sum THEN LOOP; << We must check that pivoting at a[i][c] will not ruin any row constraints that are currently satisfied. If there are any such rows, we find the first that would become unsatisfied if xc were to increase, and pivot there instead of on row i.>> << In any case, a row with an equality constraint should be pivoted upon: such slack variables never become dependent again, so pivoting on them is definite progress. >> r _ i; FOR k: INTEGER IN [1..n] WHILE rowhead[r].tag # eq DO IF k # i AND (rowhead[k].tag = eq AND a[k][c] # 0.0 OR rowhead[k].tag = geq AND b[k] >= 0.0 AND a[k][c] < 0.0 AND b[k] - a[k][c] * b[r] / a[r][c] < 0.0) THEN {r _ k} ENDLOOP ENDLOOP; << Have we got a pivot? If not, then any unsatisfied constraints that remain are usatisfiable. We just ignore them (and hope that non-linearities or random perturbations will improve things in a future Newton-Raphson step. >> IF r = 0 THEN GOTO noPivot; << Exchange row and column labels.>> BEGIN temp: TableauLabel = colhead[c]; colhead[c] _ rowhead[r]; rowhead[r] _ temp; END; << Do the pivoting.>> BEGIN Add: PROC [x, y: REAL] RETURNS [REAL] = INLINE <> {RETURN [IF ABS [x+y] < 1.0e-5*(ABS[x]+ABS[y]) THEN 0.0 ELSE x+y]}; << rewrite the equation>> << yr = ar1.x1 + ar2.x2... + arc.xc + ... + arm.xm + br as>> << xc = - (ar1/arc)x1 - (ar2/arc)x2 - ... + (1/arc)yr - ... - (arm/arc)xm - br/arc>> {coef: REAL = - 1.0/a[r][c]; a[r][c] _ -1.0; FOR j: INTEGER IN [1..m] DO a[r][j] _ coef * a[r][j] ENDLOOP; b[r] _ coef * b[r]}; << substitute xc in other rows:>> FOR i: INTEGER IN [1..n] DO IF i # r THEN {coef: REAL = a[i][c]; a[i][c] _ 0.0; FOR j: INTEGER IN [1..m] DO a[i][j] _ Add[a[i][j], coef*a[r][j]] ENDLOOP; b[i] _ Add[b[i], coef * b[r]]} ENDLOOP; END; REPEAT noPivot => NULL ENDLOOP; END; << Now move the points, and check if any moved more than eps:>> moved _ FALSE; BEGIN FOR i: INTEGER IN [1..n] DO IF rowhead[i].tag IN [x..y] THEN {p: Point = rowhead[i].point; big: REAL = MAX [4*(ABS[p.coords.x]+ABS[p.coords.y]), 100*eps]; b[i] _ MIN[big, MAX[-big, b[i]]]; -- don't move too much at once IF rowhead[i].tag = x THEN {p.coords.x _ p.coords.x + b[i]} ELSE {p.coords.y _ p.coords.y + b[i]}; IF ABS[b[i]] > eps THEN {moved _ TRUE; p.mark _ TRUE}}; ENDLOOP; END; END; << - - - - RANDOM TRIALS >> PerturbPoints: PROC [eps: REAL] = <> << To be called after PrepareTableau, in lieu of LinSolve or after it.>> << The points to be perturbed are those in rowhead[] and colhead[]. They will have their marks set, for the benfit of DisplayMotions.>> BEGIN Perturb: PROC[p: Point, which: Tag[x..y]] = <> BEGIN pert: REAL _ (10.0*eps*Random.Next[])/LAST[INT]; p.mark _ TRUE; IF which=x THEN {p.coords.x _ p.coords.x+pert} ELSE {p.coords.y _ p.coords.y+pert} END; FOR i: INTEGER IN [1..n] DO IF rowhead[i].tag=x OR rowhead[i].tag=y THEN {Perturb[rowhead[i].point, rowhead[i].tag]}; ENDLOOP; FOR j: INTEGER IN [1..m] DO IF colhead[j].tag=x OR colhead[j].tag=y THEN {Perturb[colhead[j].point, colhead[j].tag]}; ENDLOOP; END; << - - - - DEBUGGING >> DisplayMotions: PROC = <> << The displayed points are those in rowhead[] or colhead[] that have p.mark=true. >> BEGIN DrawMotion: PROC[p: Point] = INLINE {IF p.mark THEN {Gr.DrawEdge[p.old, p.coords]; p.mark _ FALSE}}; FOR i: INTEGER IN [1..n] DO IF rowhead[i].tag=x OR rowhead[i].tag=y THEN {DrawMotion[rowhead[i].point]}; ENDLOOP; FOR j: INTEGER IN [1..m] DO IF colhead[j].tag=x OR colhead[j].tag=y THEN {DrawMotion[colhead[j].point]}; ENDLOOP; END; << - - - - PUBLIC PROCEDURES >> Solve: PUBLIC PROC [items: Stor.ItemList, eps: REAL] RETURNS [outcome: Outcome] = BEGIN success, impossible: BOOL _ FALSE; moved: BOOL; THROUGH [1..15] DO [success, impossible] _ PrepareTableau[items, eps]; IF success THEN RETURN [true]; IF impossible THEN RETURN [false]; moved _ LinSolve[eps]; IF NOT moved THEN PerturbPoints [eps]; IF displayingMotions THEN DisplayMotions[] ENDLOOP; RETURN [uncertain] END; END.