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; Point: TYPE = Stor.Point; Coords: TYPE = Stor.Coords; Item: TYPE = Stor.Item; Frame: TYPE = Stor.Frame; displayingMotions: PUBLIC BOOL _ FALSE; loggingPivots: BOOL _ FALSE; TooManyConstraints: SIGNAL = CODE; TooManyConstrainedPoints: SIGNAL = CODE; 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 PrepareTableau: PROC [items: Stor.ItemList, eps: REAL] RETURNS [success, impossible: BOOLEAN] = BEGIN 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. Constraint: PROC[value: REAL, tag: Tag] = 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; 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; 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; 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; 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; 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]]}; 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; LinSolve: PROC [eps: REAL] RETURNS [moved: BOOLEAN] = BEGIN BEGIN DO r, c: INTEGER _ 0; FOR i: INTEGER IN [1..n] WHILE r = 0 DO champ, sum: REAL _ 0.0; 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; 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; IF c = 0 OR champ < .001*sum THEN LOOP; 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; IF r = 0 THEN GOTO noPivot; BEGIN temp: TableauLabel = colhead[c]; colhead[c] _ rowhead[r]; rowhead[r] _ temp; END; 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]}; {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]}; 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; 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; PerturbPoints: PROC [eps: REAL] = 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; DisplayMotions: PROC = 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; 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. "ė JunoOldSolverImpl.mesa (was OldJunoSolverImpl.mesa) Coded July 1981 by Greg Nelson Last Edited by: Gnelson(?) December 10, 1982 2:25 pm Last Edited by: Gnelson, October 11, 1983 7:00 pm Last Edited by: Stolfi, June 13, 1984 10:06:42 am PDT 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. - - - - IMPORTED TYPES - - - - GLOBAL VARIABLES - - - - SIGNALS - - - - GLOBAL VARIABLES a[i][j] is partial deriv. of ith constraint wrt jth unknown b[i] is minus the numerical value of the ith constraint - - - - 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. Analyzes the constraints and computes tableau for Newton-Raphson step. 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. Procedures to enter constraints in the tableau: (Assuming all calls to Partial are for distinct points). Adds row to matrix for a new constraint whose current value is value. Sets satisfied according to present status of constraint. Also sets unsatisfiable if not trivially satisfied; this flag may be reset by Partial. Adds to the matrix one or two new entries, corresponding to the partial derivatives of the most recent row (constraint) in the tableau, with respect to the x and y coordinates of the given point. Gets column corresponding to the specified coordinate of p. Will assign a new column if not yet assigned. Procedures for specific constraints: frame is irrelevant 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. 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. Other useful procs: Returns positive number telling how much a constraint is violated. Fill in tableau a[][] and constant column b[] by enumerating current constraints: - - - - 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. Does one Newton-Raphson step. 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. Invert matrix: Look for a suitable pivot a[r][c] Ignore dependent variables that are unconstrained or trivially satisfiable 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]]) Is a[i][c] large enough to pivot on it? 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. 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. Exchange row and column labels. Do the pivoting. Adds x+y but returns 0.0 if there was excessive cancellarion 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 substitute xc in other rows: Now move the points, and check if any moved more than eps: - - - - RANDOM TRIALS Adds a small random multiple of eps to every non-fixed coordinate that enters in the given constraints. 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. Perturbs the x and/or y coordinate of the given point. - - - - DEBUGGING Displays the motions of the points altered by LinSolve or PerturbPoints. The displayed points are those in rowhead[] or colhead[] that have p.mark=true. - - - - PUBLIC PROCEDURES Ęļ4JJ4J1J7JņĪb J&&Īk J OWÔJ!!J JEEJJ J JJ (JJ#J)J J JJ6J *J > J;; J77J+J+JĪc<LJ <1;;JĀĀJ~~JzzJjjJJJąąĪnbFFJĐcruJĪrĄBW008I8ãJ78 *?EJĸ 6ĸ ?J, ÚJ .JÃÃ 8Jii ^J ,0ĨJJ J 3 2ĖJ%% ,!8J88 !=J88 &JPPJ ,!8J88"(JPPJ #3J #!#ÜJJJ #3J #!#ÜB,11+æJŋŋJííJ ĩJ 0=J J ĻJčđ+1JĄĄJ #3J #!#ÜH,11+æJŋŋJííJ ļJ 4J "!ÂJ ĻJčđ17JJ &4J U J$$*^J F KG§J$$*^J F KFĻJ /J J22HJ ŋJ33J##JJ39JJQQJ 2JBBJsJJ ]J bJ gJ lĸĸ$RJ2 &K%J #J11J11JGGJGGJGGJ 0 J;;JJ 6J ZJ JJ JJņP ´JŪŪJJ ˇđJ-ŽĘ¸ 6JĄI^JCDJ"# (J LMJ uĨĨJ 1æ((J (JņņĻĻJJ ( J ŨŪJ JMMJ /J