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 7, 1984 3:55:28 pm PDT

JunoStorage USING [Point, Coords, Frame, nullFrame, ConstrList, Constr,
HorConstr, VerConstr, ParaConstr, CongConstr, PerpConstr, AtConstr, CcwConstr],
JunoMatrix USING [Matrix, GetFrameMatrix, MapCoords],
JunoOldSolver,
Random USING [Next],
JunoGraphics USING [DrawEdge],
Real;

JunoMatrix,
Random,
JunoGraphics

JunoOldSolver

Stor: JunoStorage,
Mat: JunoMatrix,
JunoOldSolver,
Gr: JunoGraphics;

a[i][j] is partial deriv. of ith constraint wrt jth unknown

b[i] is minus the numerical value of the ith constraint

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.

BEGIN

END;

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.

BEGIN

END;

Adds a small random multiple of eps to every non-fixed coordinate that enters in the given constraints.

BEGIN

END;

Displays the motions of the points altered by LinSolve or PerturbPoints.

BEGIN

END;

BEGIN

END;

{n: INT ← 0;
SavePositions[];
[] ← LinSolve[];
WHILE n < 20 AND SomePointMoved[] DO
IF displayingMotions THEN DisplayMotions[];
n ← n + 1;
SavePositions[];
[] ← LinSolve[]
ENDLOOP;
RETURN [n < 20]};

{p: Point ← pointLpad.link;
UNTIL p = pointRpad DO p.oldx ← Real.Fix[p.x]; p.oldy ← Real.Fix[p.y]; p ← p.link ENDLOOP};

{p: Point ← pointLpad.link;
UNTIL p = pointRpad DO
IF p.oldx # Real.Fix[p.x] OR p.oldy # Real.Fix[p.y] THEN RETURN [TRUE];
p ← p.link
ENDLOOP;
RETURN [FALSE]};

{p: Point ← pointLpad.link;
UNTIL p = pointRpad DO
IF p.oldx # Real.Fix[p.x] OR p.oldy # Real.Fix[p.y]
THEN Gr.DrawEdge[p.oldx, p.oldy, p.x, p.y];
p ← p.link
ENDLOOP};

BEGIN

p: Point ← pointLpad.link;
UNTIL p = pointRpad DO p.xheader ← 0; p.yheader ← 0; p ← p.link ENDLOOP;

success ← TRUE;
n ← 0;
m ← 1; -- matrix has n rows and m columns.
-- column 1 is the constant column.

BEGIN -- handle counter-clockwise constraints:
p: CCPtr ← ccLpad.link;
UNTIL p = ccRpad DO
NewtonStep[(p.j.y - p.i.y) * (p.k.x - p.i.x) - (p.j.x - p.i.x) * (p.k.y - p.i.y) + 500.0,
GradCC[p.i, p.j, p.k], geq];
p ← p.link;
ENDLOOP;
END;

BEGIN -- handle string constraints:
p: StringPtr ← stringLpad.link;
UNTIL p = stringRpad DO
NewtonStep [p.b4.x - p.b3.x - p.width, GradVer[p.b4, p.b3], eq];
NewtonStep [p.b4.y - p.b3.y, GradHor[p.b4, p.b3], eq];
p ← p.link;
ENDLOOP;
END;

BEGIN -- handle horizontal constraints:
p: HorPtr ← horLpad.link;
UNTIL p = horRpad DO
NewtonStep[p.i.y - p.j.y, GradHor[p.i, p.j], eq];
p ← p.link
ENDLOOP;
END;

BEGIN -- handle vertical constraints:
p: VerPtr ← verLpad.link;
UNTIL p = verRpad DO
NewtonStep[p.i.x - p.j.x, GradVer[p.i, p.j], eq];
p ← p.link;
ENDLOOP;
END;

BEGIN -- handle parallel constraints:
p: LinPtr ← lineLpad.link;
UNTIL p = lineRpad DO
NewtonStep[(p.j.x - p.i.x)*(p.k.y - p.l.y) + (p.l.x - p.k.x)*(p.j.y - p.i.y),
GradLin[p.i, p.j, p.k, p.l], eq];
p ← p.link;
ENDLOOP;
END;

BEGIN -- handle congruence constraints:
p: CongPtr ← congLpad.link;
UNTIL p = congRpad DO
NewtonStep[0.5 * (Square[p.j.x - p.i.x] + Square[p.j.y - p.i.y]
- Square[p.l.x - p.k.x] - Square[p.l.y - p.k.y]),
  GradCong[p.i, p.j, p.k, p.l], eq];
p ← p.link;
ENDLOOP;
END;

-- Now invert matrix, and move points:

BEGIN
i, j, k, l: INT;
pp: REAL;
-- k is the row in which we are pivoting.
-- l is the column in which we are pivoting.
-- i and j are miscellaneous row and column indices respectively
-- p is the reciprocal of the pivot element; also used as temp for swapping.

FOR k IN [1..n] DO
IF rowhead[k].tag = eq OR rowhead[k].tag = geq AND a[k][1] < -.01
THEN {-- set l so a[k,l] is largest of a[k,1], a[k, 2], a[k, 3], excluding
-- columns in which we have already pivoted.
pp ← 0;
l ← 0;
FOR j IN [2 .. m]
DO IF ABS[a[k][j]] >= pp AND colhead[j].tag IN [x .. y]
  THEN {l ← j; pp ← ABS[a[k][l]]}
ENDLOOP;
-- We will pivot at a[k,l], if it exists and is not too small, and is legal.
-- First check that it is not too small:
IF l = 0 OR ABS[a[k][l]] < .001 THEN
IF ABS[a[k][1]] > .01
  -- then this row is an unsatisfiable constraint so:
  THEN {LOOP}
  ELSE LOOP; -- row is identically zero; constraint redundant.
  -- Next we must check that pivoting at a[k, l] will not ruin any
  -- inequalities that are currently satisfied. If there are any such inequalities,
  -- we change k to the index of one of them (namely, the first that would become
  -- unsatisfied), and pivot there instead. That is the function of the following block:
 
  { -- adjust k to account for inequalities:
  kk: INTEGER ← k; -- save k
  FOR i: INT IN [1..n] DO
  IF i # kk AND rowhead[i].tag = geq AND a[i][1] >= 0.0
  AND a[i][l] < 0.0 AND a[i][1] - a[i][l] * a[k][1] / a[k][l] < 0.0
  THEN k ← i
  ENDLOOP};
 
{temp:TableauLabel ← colhead[l]; colhead[l] ← rowhead[k]; rowhead[k] ← temp};
pp ← 1.0 / a[k][l]; a[k][l] ← 1.0;
-- divide everything in pivot row by the pivot element:
FOR j IN [1..m] DO a[k][j] ← a[k][j] * pp ENDLOOP;

FOR i IN [1..n] DO
IF i # k THEN
  FOR j IN [1..m] DO
  IF j # l THEN -- for each a[i,j] outside the pivot row and column
  a[i][j] ← a[i][j] - a[i][l] * a[k][j]; -- note that a[k,j] was already * pp.
ENDLOOP ENDLOOP;

-- Finally process pivot column:
FOR i IN [1..n] DO IF i # k THEN a[i][l] ← -a[i][l] * pp ENDLOOP;
  }

ENDLOOP;

-- now move the points:

FOR i IN [1..n] DO
IF rowhead[i].tag IN [x .. y]
-- as it will be except for redundant rows
-- (or unsatisfiable rows)
THEN
IF rowhead[i].tag = x
THEN rowhead[i].point.x ← rowhead[i].point.x + a[i][1]
  ELSE rowhead[i].point.y ← rowhead[i].point.y + a[i][1];
ENDLOOP;

END;
END;