<<>> <> <> <> DIRECTORY PBasics, QPSetup, QPSolve, RealFns, Rope; QPSetupImpl: CEDAR PROGRAM IMPORTS PBasics, QPSolve, RealFns, Rope EXPORTS QPSetup = BEGIN OPEN QPSolve, QPSetup; SparseMatrix: TYPE ~ REF SparseMatrixRep; SparseMatrixRep: TYPE ~ RECORD [m: NAT, s: SEQUENCE len: NAT OF TermList]; FindFreeVars: PUBLIC PROC [cycleEqns, extraEqns: EqnList, nEqns, nVars: NAT] RETURNS [rowOf: IVector] <> <> ~ { rowOf _ NEW[IVectorRep[nVars]]; FOR col: NAT IN [0..nVars) DO rowOf[col] _ noRow ENDLOOP; FOR eqns: EqnList _ cycleEqns, eqns.rest UNTIL eqns = NIL DO rowOf[eqns.first.first.eData.var] _ 0; -- first var of cycle eqn is back edge ENDLOOP; <<<< Need a QR algorithm here to accomodate extra eqns, ensure independence. >>>> }; FindComponents: PUBLIC PROC [g: Graph] RETURNS [compList: ComponentList _ NIL] <> ~ { edgeStack: EdgeList _ NIL; eqnStack: EqnList _ NIL; Stack: PROC [edge: Edge] RETURNS [EdgeList] ~ {RETURN[edgeStack _ CONS[edge, edgeStack]]}; UnStackComponent: PROC [bottomEdge: EdgeList, bottomEqn: EqnList] ~ { component: Component _ NEW[ComponentRep _ [join: NIL, weight: 1, edges: edgeStack, cycleEqns: eqnStack, extraEqns: NIL]]; edgeStack _ bottomEdge.rest; bottomEdge.rest _ NIL; IF eqnStack = bottomEqn THEN component.cycleEqns _ NIL ELSE { eqnStack _ bottomEqn; FOR eqns: EqnList _ component.cycleEqns, eqns.rest UNTIL eqns = NIL DO IF eqns.rest = bottomEqn THEN {eqns.rest _ NIL; EXIT}; ENDLOOP; }; FOR edges: EdgeList _ component.edges, edges.rest UNTIL edges = NIL DO edge: Edge _ edges.first; edge.eData.join _ component; ENDLOOP; compList _ CONS[component, compList]; }; FormEquation: PROC [e: Edge] RETURNS [eqn: TermList] ~ { v: Vertex _ e.thisEnd; w: Vertex _ e.thatEnd; eqn _ NIL; FOR t: Vertex _ v.tree.thisEnd, v.tree.thisEnd WHILE v.visited > w.visited DO eqn _ CONS[Term[eData: v.tree.eData, coef: v.tree.dir], eqn]; v _ t; ENDLOOP; FOR t: Vertex _ w.tree.thisEnd, w.tree.thisEnd WHILE w.visited > v.visited DO eqn _ CONS[Term[eData: w.tree.eData, coef: -w.tree.dir], eqn]; w _ t; ENDLOOP; eqn _ CONS[Term[eData: e.eData, coef: e.dir], eqn]; -- back edge is first, for FindFreeVars }; BiConn: PROC [g: Graph, treeEdge: Edge] RETURNS [vMin: INT] ~ { from: Vertex _ treeEdge.thisEnd; v: Vertex _ treeEdge.thatEnd; v.tree _ treeEdge; v.visited _ vMin _ g.order _ g.order+1; FOR edges: EdgeList _ v.adj, edges.rest UNTIL edges = NIL DO edge: Edge _ edges.first; w: Vertex _ edge.thatEnd; IF w.visited = 0 THEN { bottomEdge: EdgeList _ Stack[edge]; bottomEqn: EqnList _ eqnStack; wMin: INT _ BiConn[g, edge]; vMin _ MIN[vMin, wMin]; IF wMin >= v.visited THEN UnStackComponent[bottomEdge, bottomEqn]; } ELSE { IF w = from THEN LOOP; vMin _ MIN[vMin, w.visited]; IF v.visited > w.visited THEN { [] _ Stack[edge]; eqnStack _ CONS[FormEquation[edge], eqnStack]; }; }; ENDLOOP; }; FOR vL: VertexList _ g.vList.rest, vL.rest UNTIL vL.first = NIL DO vL.first.visited _ 0 ENDLOOP; g.order _ 0; FOR vL: VertexList _ g.vList.rest, vL.rest UNTIL vL.first = NIL DO IF vL.first.visited = 0 THEN [] _ BiConn[g, NEW[EdgeRep _ [NIL, vL.first, NIL, 1.0]]]; ENDLOOP; }; TangleOf: PROC [x: Component] RETURNS [z: Component] ~ { <> t: Component; FOR z _ x, z.join DO IF z.join = NIL THEN EXIT ENDLOOP; UNTIL x.join = NIL DO t _ x; x _ x.join; t.join _ z ENDLOOP; }; Entangle: PROC [x, y: Component] ~ { <> i, j, z, t: Component; FOR i _ x, i.join DO IF i.join = NIL THEN EXIT ENDLOOP; FOR j _ y, j.join DO IF j.join = NIL THEN EXIT ENDLOOP; IF i # j THEN { IF j.weight > i.weight THEN {j.weight _ j.weight+i.weight; i.join _ z _ j} ELSE {i.weight _ i.weight+j.weight; j.join _ z _ i}; } ELSE z _ i; UNTIL x.join = NIL DO t _ x; x _ x.join; t.join _ z ENDLOOP; UNTIL y.join = NIL DO t _ y; y _ y.join; t.join _ z ENDLOOP; }; CoupleComponents: PUBLIC PROC [compList: ComponentList, extraEqns: EqnList] RETURNS [coupledList: ComponentList] <> ~ { FOR eqns: EqnList _ extraEqns, eqns.rest UNTIL eqns = NIL DO terms: TermList _ eqns.first; comp: Component _ TangleOf[terms.first.eData.join]; FOR terms _ terms.rest, terms.rest UNTIL terms = NIL DO Entangle[comp, terms.first.eData.join]; ENDLOOP; comp.extraEqns _ CONS[eqns.first, comp.extraEqns]; ENDLOOP; FOR comps: ComponentList _ compList, comps.rest UNTIL comps = NIL DO comp: Component _ comps.first; tangle: Component _ TangleOf[comp]; IF tangle = comp THEN coupledList _ CONS[comp, coupledList] ELSE { IF comp.edges # NIL THEN { edges: EdgeList _ comp.edges; UNTIL edges.rest = NIL DO edges _ edges.rest ENDLOOP; edges.rest _ tangle.edges; tangle.edges _ comp.edges; }; IF comp.cycleEqns # NIL THEN { eqnList: EqnList _ comp.cycleEqns; UNTIL eqnList.rest = NIL DO eqnList _ eqnList.rest ENDLOOP; eqnList.rest _ tangle.cycleEqns; tangle.cycleEqns _ comp.cycleEqns; }; IF comp.extraEqns # NIL THEN { eqnList: EqnList _ comp.extraEqns; UNTIL eqnList.rest = NIL DO eqnList _ eqnList.rest ENDLOOP; eqnList.rest _ tangle.extraEqns; tangle.extraEqns _ comp.extraEqns; }; }; ENDLOOP; }; NumberEdges: PROC [edges: EdgeList] RETURNS [nVars: NAT _ 0] <> ~ { UNTIL edges = NIL DO edges.first.eData.var _ nVars _ nVars+1; edges _ edges.rest; ENDLOOP; RETURN[nVars+1]; }; CountEquations: PROC [eqnList: EqnList] RETURNS [nEqns: NAT _ 0] <> ~ {UNTIL eqnList = NIL DO nEqns _ nEqns+1; eqnList _ eqnList.rest ENDLOOP}; SolveComponent: PUBLIC PROC [component: Component] <> ~ { penalty: REAL _ 1.0e6; nVars: NAT _ NumberEdges[component.edges]; nEqns: NAT _ CountEquations[component.cycleEqns]+CountEquations[component.extraEqns]; lobd: RVector _ NEW[RVectorRep[nVars]]; x: RVector _ NEW[RVectorRep[nVars]]; c: RVector _ NEW[RVectorRep[nVars]]; A: Matrix _ NewMatrix[nEqns, nVars]; rowOf: IVector _ FindFreeVars[component.cycleEqns, component.extraEqns, nEqns, nVars]; iVar: IVector _ NEW[IVectorRep[nVars]]; hiMark: INT _ nVars; nFR: NAT _ 0; eqnList: EqnList _ component.cycleEqns; <> lobd.n _ x.n _ c.n _ iVar.n _ nVars; <> FOR j: NAT IN [0..nEqns) DO FOR i: NAT IN [0..nVars) DO A[j][i] _ 0.0 ENDLOOP; FOR termList: TermList _ eqnList.first, termList.rest UNTIL termList = NIL DO i: NAT _ termList.first.eData.var; A[j][i] _ termList.first.coef / termList.first.eData.squeeze; ENDLOOP; eqnList _ eqnList.rest; IF eqnList = NIL THEN eqnList _ component.extraEqns; ENDLOOP; <> FOR edgeList: EdgeList _ component.edges, edgeList.rest UNTIL edgeList = NIL DO eData: EdgeData _ edgeList.first.eData; var: NAT _ eData.var; lobd[var] _ eData.spread*eData.squeeze; c[var] _ -eData.relaxed/eData.squeeze; IF rowOf[var] = noRow THEN {iVar[hiMark _ hiMark-1] _ var; x[var] _ lobd[var]} ELSE {iVar[(nFR _ nFR+1)] _ var; x[var] _ lobd[var]+1.0}; ENDLOOP; <> iVar[0] _ 0; nFR _ nFR+1; x[0] _ 0.0; FOR j: NAT IN [0..nEqns) DO FOR i: NAT IN [1..nVars) DO A[j][0] _ A[j][0] + A[j][i]*x[i] ENDLOOP; x[0] _ x[0] + A[j][0] * A[j][0]; ENDLOOP; x[0] _ RealFns.SqRt[x[0]]; FOR j: NAT IN [0..nEqns) DO A[j][0] _ -A[j][0] / x[0] ENDLOOP; lobd[0] _ 0.0; c[0] _ penalty; <> nFR _ QPSolve[c, A, lobd, x, iVar, nFR]; FOR edgeList: EdgeList _ component.edges, edgeList.rest UNTIL edgeList = NIL DO edgeList.first.eData.val _ x[edgeList.first.eData.var]; ENDLOOP; }; SetValues: PROC [v: Vertex] ~ { v.visited _ 0; FOR edges: EdgeList _ v.adj, edges.rest UNTIL edges = NIL DO edge: Edge _ edges.first; w: Vertex _ edge.thatEnd; IF w.visited # 0 THEN {w.val _ v.val+edge.dir*edge.eData.val/edge.eData.squeeze; SetValues[w]}; ENDLOOP; }; SolveGraph: PUBLIC PROC [g: Graph] <> ~ { biconnList: ComponentList _ FindComponents[g]; coupledList: ComponentList _ CoupleComponents[biconnList, g.eqnList]; FOR cL: ComponentList _ coupledList, cL.rest UNTIL cL = NIL DO edge: Edge _ cL.first.edges.first; IF cL.first.edges.rest = NIL AND cL.first.cycleEqns = NIL AND cL.first.extraEqns = NIL THEN edge.eData.val _ edge.eData.spread*edge.eData.squeeze ELSE SolveComponent[cL.first]; ENDLOOP; FOR vL: VertexList _ g.vList.rest, vL.rest UNTIL vL.first = NIL DO IF vL.first.visited # 0 THEN {vL.first.val _ 0.0; SetValues[vL.first]}; ENDLOOP; }; mem: PBasics.Word _ 12422627672B; RandomBit: PROC RETURNS [on: BOOL] ~ { hiBit: PBasics.Word ~ 20000000000B; bits: PBasics.Word ~ 10000000123B+hiBit; mem _ mem+mem; -- or use PBasics.BITLSHIFT[mem, 1] on _ mem >= hiBit; -- or use PBasics.BITAND[mem, hiBit] # 0 IF on THEN mem _ PBasics.BITXOR[mem, bits]; }; maxLevel: NAT ~ 4; RandomLevel: PROC RETURNS [level: Level _ 1] ~ { WHILE RandomBit[] AND RandomBit[] AND level < maxLevel DO level _ level + 1 ENDLOOP; }; sentinel: VertexList ~ NEW[VertexListRep _ [name: "\377", first: NIL, forward: NIL]]; NewVertexList: PROC [level: Level, name: ROPE, vid: INT _ 0] RETURNS [vL: VertexList] ~ { vL _ NEW[VertexListRep _ [name: name, first: NIL, forward: NEW[SkipsRep[level]]]]; IF vid = 0 THEN {vL.rest _ sentinel; FOR i: Level IN [0 .. level) DO vL.forward[i] _ sentinel ENDLOOP} ELSE vL.first _ NEW[VertexRep _ [name: name, visited: 0, tree: NIL, adj: NIL, val: 0.0, vid: vid]]; }; NewGraph: PUBLIC PROC RETURNS [g: Graph] ~ { g _ NEW[GraphRep _ [order: 0, nV: 0, vList: NewVertexList[maxLevel, NIL, 0], eqnList: NIL]]; }; update: Skips _ NEW[SkipsRep[maxLevel]]; AcquireVertex: PUBLIC PROC [g: Graph, name: ROPE] RETURNS [v: Vertex] ~ { list: VertexList _ g.vList; x: VertexList _ list; FOR i: Level DECREASING IN [0 .. list.forward.level) DO WHILE Rope.Compare[x.forward[i].name, name] = less DO x _ x.forward[i] ENDLOOP; update[i] _ x; ENDLOOP; x _ x.rest; IF NOT Rope.Equal[x.name, name] THEN { level: Level _ RandomLevel[]; x _ NewVertexList[level, name, g.nV _ g.nV+1]; x.rest _ update[0].rest; update[0].rest _ x; FOR i: Level IN [0 .. level) DO x.forward[i] _ update[i].forward[i]; update[i].forward[i] _ x; ENDLOOP; }; RETURN[x.first]; }; AcquireEdge: PUBLIC PROC [g: Graph, thisEnd, thatEnd: Vertex] RETURNS [edge: Edge] ~ { back: Edge; FOR edgeList: EdgeList _ thisEnd.adj, edgeList.rest UNTIL edgeList = NIL DO IF edgeList.first.thatEnd = thatEnd THEN RETURN[edgeList.first]; ENDLOOP; edge _ NEW[EdgeRep _ [thisEnd: thisEnd, thatEnd: thatEnd, eData: NIL, dir: 1.0]]; edge.eData _ NEW[EdgeDataRep]; thisEnd.adj _ CONS[edge, thisEnd.adj]; back _ NEW[EdgeRep _ [thisEnd: thatEnd, thatEnd: thisEnd, eData: edge.eData, dir: -1.0]]; thatEnd.adj _ CONS[back, thatEnd.adj]; }; SpreadEdge: PUBLIC PROC [g: Graph, edge: Edge, spread: REAL] ~ {edge.eData.spread _ spread}; SqueezeEdge: PUBLIC PROC [g: Graph, edge: Edge, squeeze: REAL] ~ {edge.eData.squeeze _ RealFns.SqRt[ABS[squeeze]]}; RelaxEdge: PUBLIC PROC [g: Graph, edge: Edge, relaxed: REAL] ~ {edge.eData.relaxed _ relaxed}; ConstrainEdges: PUBLIC PROC [g: Graph, termList: TermList] ~ {g.eqnList _ CONS[termList, g.eqnList]}; ListGraph: PUBLIC PROC [g: Graph, full: BOOL _ FALSE] RETURNS [list: LIST OF REF _ NIL] ~ { FOR vL: VertexList _ g.vList.rest, vL.rest UNTIL vL.first = NIL DO v: Vertex _ vL.first; adj: LIST OF REF _ NIL; FOR vAdj: EdgeList _ v.adj, vAdj.rest UNTIL vAdj = NIL DO IF vAdj.first.dir > 0 THEN adj _ CONS[LIST[vAdj.first.thatEnd.name, vAdj.first.eData], adj] ELSE { w: Vertex _ vAdj.first.thatEnd; FOR wAdj: EdgeList _ w.adj, wAdj.rest UNTIL wAdj = NIL DO IF wAdj.first.thatEnd = v THEN EXIT; REPEAT FINISHED => adj _ CONS[LIST[Rope.Concat["?", w.name], vAdj.first.eData], adj]; ENDLOOP; }; ENDLOOP; list _ CONS[LIST[vL.name, adj], list]; ENDLOOP; }; ESee: PUBLIC PROC [edge: Edge, full: BOOL _ FALSE] RETURNS [list: LIST OF REF _ NIL] ~ { rope: ROPE; IF edge = NIL THEN RETURN[LIST["<>"]]; rope _ Rope.Cat[edge.thisEnd.name, IF edge.dir > 0 THEN ">" ELSE "<", edge.thatEnd.name]; IF full THEN list _ LIST[rope, edge.eData] ELSE list _ LIST[rope]; }; ELSee: PUBLIC PROC [edges: EdgeList, full: BOOL _ FALSE] RETURNS [list: LIST OF REF _ NIL] ~ { UNTIL edges = NIL DO list _ CONS[ESee[edges.first, full], list]; edges _ edges.rest ENDLOOP; }; ListComp: PUBLIC PROC [comp: Component, full: BOOL _ FALSE] RETURNS [list: LIST OF REF] ~ { list _ LIST[ELSee[comp.edges, full], comp.cycleEqns, comp.extraEqns]; }; END.. <> <> <<~ {>> <> <> <> <<~ {>> <> <> <> <> <> <> <> <> <> <> <> <> <> <<};>> <> <> <> <> <> <<};>> <<>> <> <<~ {>> <> <> <> <> <> <> <> <> <> <> <> <<};>> <> <> <> <> <> <> <> <<~ {>> <> <> <> <> <> <> <> <> <> <<};>> <<>>