<> <> DIRECTORY LeastSq, Random, Real; LeastSqImpl: CEDAR PROGRAM IMPORTS Random, Real EXPORTS LeastSq = BEGIN NewPointSet: PUBLIC PROC [ dim: NAT ] RETURNS [ ps: LeastSq.PointSet ] = { ps _ NEW[LeastSq.PointSetRec[dim]]; ps.count _ 0; FOR row: NAT IN [0..dim) DO ps[row] _ NEW[LeastSq.RowRec[dim]]; FOR col: NAT IN [0..dim) DO ps[row][col] _ 0.0; ENDLOOP; ENDLOOP; }; AddPoint: PUBLIC PROC [ ps: LeastSq.PointSet, pt: LeastSq.Row, weight: REAL _ 1.0 ] = { IF pt.colDim # ps.rowDim THEN ERROR; ps.count _ ps.count+1; FOR row: NAT IN [0..ps.rowDim) DO FOR col: NAT IN [0..row] DO ps[row][col] _ ps[row][col]+weight*pt[row]*pt[col]; ENDLOOP; ENDLOOP; }; LeastSqFit: PUBLIC PROC [ ps: LeastSq.PointSet ] RETURNS [ eq: LeastSq.Row ] = { <> <<1 if k=0,>> <0 and k> <> <<>> <> <<>> <> basis: LeastSq.PointSet = NEW[LeastSq.PointSetRec[ps.rowDim]]; FOR row: NAT IN [0..ps.rowDim) DO -- fill ps symmetrically FOR col: NAT IN (row..ps.rowDim) DO ps[row][col] _ ps[col][row]; ENDLOOP; ENDLOOP; FOR row: NAT IN [0..ps.rowDim) DO innerProd: REAL; col: NAT; basis[row] _ NEW[LeastSq.RowRec[ps.rowDim]]; FOR col IN [0..ps.rowDim) DO basis[row][col] _ 0.0 ENDLOOP; FOR prevrow: NAT IN [0..row) DO innerProd _ 0.0; FOR col IN [0..ps.rowDim) DO innerProd _ innerProd+basis[prevrow][col]*ps[row][col]; ENDLOOP; -- col FOR col IN [0..ps.rowDim) DO basis[row][col] _ basis[row][col]-innerProd*basis[prevrow][col]; ENDLOOP; -- col ENDLOOP; -- prevrow basis[row][row] _ 1.0; IF row0.0 THEN { normFact: REAL = 1.0/Real.SqRt[innerProd]; FOR col IN [0..ps.rowDim) DO basis[row][col] _ normFact*basis[row][col]; ENDLOOP; -- col }; END; ENDLOOP; -- row eq _ basis[ps.rowDim-1]; }; NewRow: PROC [ dim: NAT ] RETURNS [ pt: LeastSq.Row ] = { pt _ NEW[LeastSq.RowRec[dim]] }; rs: Random.RandomStream _ Random.Create[]; Test: PROC [ eqIn: LeastSq.Row, n: INT, noise: REAL _ 0 ] RETURNS [ eqOut: LeastSq.Row ] = BEGIN ps: LeastSq.PointSet = NewPointSet[eqIn.colDim]; pt: LeastSq.Row = NEW[LeastSq.RowRec[eqIn.colDim]]; pt[0] _ 1.0; FOR i: INT IN [0..n) DO sum: REAL _ eqIn[0]+noise*rs.ChooseInt[-10000, 10000]; FOR j: NAT IN [1..eqIn.colDim-1) DO pt[j] _ rs.ChooseInt[-10000, 10000]; sum _ sum+eqIn[j]*pt[j]; ENDLOOP; pt[eqIn.colDim-1] _ -sum; AddPoint[ps, pt]; ENDLOOP; eqOut _ LeastSqFit[ps]; END; END.