Forum > Beginners
Need help/explanation - please
(1/1)
noob:
i have a unit that contains the following procedure, that crashes on exit of procedure with
EAccessViolation: Access violation
--- Code: ---{****************************************************************
* LEAST SQUARES POLYNOMIAL FITTING PROCEDURE *
* ------------------------------------------------------------- *
* This program least squares fits a polynomial to input data. *
* Forsythe orthogonal polynomials are used in the fitting. *
* The number of data points is n. *
* The data is input to the subroutine in x[i], y[i] pairs. *
* The coefficients are returned in c[i], *
* the smoothed data is returned in v[i], *
* the order of the fit is specified by m. *
* The standard deviation of the fit is returned in d. *
* There are two options available by use of the parameter e: *
* 1. if e = 0, the fit is to order m, *
* 2. if e > 0, the order of fit increases towards m, but will *
* stop if the relative standard deviation does not decrease *
* by more than e between successive fits. *
* The order of the fit then obtained is l. *
****************************************************************}
PROCEDURE LS_POLY(m: INTEGER; e1: EXTENDED; x, y: ARRAY OF EXTENDED; VAR outcoeff: ARRAY OF EXTENDED; VAR orderFit: INTEGER; VAR stddev: EXTENDED);
LABEL
10, 15, 20, 30, 50, fin;
VAR
i, l2, n1, l, n: INTEGER;
a1, a2, b1, b2, c1, d1, f1, f2, v1, v2, w, dd, vv: EXTENDED;
v, a, b, c, d, c2, e, f: ARRAY OF EXTENDED;
BEGIN
n := High(x);
l := 0;
n1 := m + 1;
v1 := 1e7;
{Initialize the arrays}
SetLength(v, n); SetLength(c, n); SetLength(d, n); SetLength(c2, n);
SetLength(e, n); SetLength(a, n1); SetLength(b, n1); SetLength(f, n1);
FOR i := 0 TO n1 DO
BEGIN
a[i] := 0;
b[i] := 0;
f[i] := 0;
END;
FOR i := 0 TO n DO
BEGIN
v[i] := 0;
d[i] := 0;
END;
d1 := SQRT(n);
w := d1;
FOR i := 0 TO n DO
e[i] := 1 / w;
f1 := d1;
a1 := 0;
FOR i := 0 TO n DO
a1 := a1 + x[i] * e[i] * e[i];
c1 := 0;
FOR i := 0 TO n DO
c1 := c1 + y[i] * e[i];
b[1] := 1 / f1;
f[1] := b[1] * c1;
FOR i := 0 TO n DO
v[i] := v[i] + e[i] * c1;
m := 1;
10: {Save latest results}
FOR i := 0 TO l DO
c2[i] := c[i];
l2 := l;
v2 := v1;
f2 := f1;
a2 := a1;
f1 := 0;
FOR i := 0 TO n DO
BEGIN
b1 := e[i];
e[i] := (x[i] - a2) * e[i] - f2 * d[i];
d[i] := b1;
f1 := f1 + e[i] * e[i];
END;
f1 := SQRT(f1);
FOR i := 0 TO n DO
e[i] := e[i] / f1;
a1 := 0;
FOR i := 0 TO n DO
a1 := a1 + x[i] * e[i] * e[i];
c1 := 0;
FOR i := 0 TO n DO
c1 := c1 + e[i] * y[i];
m := m + 1;
i := 0;
15:
l := m - i;
b2 := b[l];
d1 := 0;
IF l > 1 THEN
d1 := b[l - 1];
d1 := d1 - a2 * b[l] - f2 * a[l];
b[l] := d1 / f1;
a[l] := b2;
i := i + 1;
IF i <> m THEN
GOTO 15;
FOR i := 0 TO n DO
v[i] := v[i] + e[i] * c1;
FOR i := 0 TO n1 DO
BEGIN
f[i] := f[i] + b[i] * c1;
c[i] := f[i];
END;
vv := 0;
FOR i := 0 TO n DO
vv := vv + (v[i] - y[i]) * (v[i] - y[i]);
{Note the division is by the number of degrees of freedom}
vv := SQRT(vv / (n - l - 1));
l := m;
IF e1 = 0 THEN
GOTO 20;
{Test for minimal improvement}
IF ABS(v1 - vv) / vv < e1 THEN
GOTO 50;
{If error is larger, quit}
IF e1 * vv > e1 * v1 THEN
GOTO 50;
v1 := vv;
20:
IF m = n1 THEN
GOTO 30;
GOTO 10;
30: {Shift the c[i] down, so c(0) is the constant term}
FOR i := 0 TO l DO
c[i - 1] := c[i];
c[l] := 0;
{l is the order of the polynomial fitted}
l := l - 1;
dd := vv;
GOTO fin;
50: {Aborted sequence, recover last values}
l := l2;
vv := v2;
FOR i := 0 TO l DO
c[i] := c2[i];
GOTO 30;
fin:
FOR i := 0 TO l DO
outcoeff[i] := c[i];
orderFit := l;
stddev := dd;
SetLength(v, 0); SetLength(c, 0); SetLength(d, 0); SetLength(c2, 0);
SetLength(e, 0); SetLength(a, 0); SetLength(b, 0); SetLength(f, 0);
exit;
END;
--- End code ---
i call the procedure with:
--- Code: ---VAR
i,j: INTEGER;
m: INTEGER;
e1, stddev:EXTENDED;
x, y: ARRAY OF EXTENDED;
outcoeff: ARRAY OF EXTENDED;
orderfitobtained: INTEGER;
begin
i:= HIGH(orgdat);
SetLength(x, i);
for j := 0 TO i do
x[j] := j+1;
m := 3;
SetLength(outcoeff, m+1);
e1 := 0; orderfitobtained := 0; stddev := 0.0;
LS_POLY(m, e1, x, orgdat{y}, outcoeff,orderfitobtained,stddev);
--- End code ---
any help possible?
i use the newest lazarus ()64bit, win8.1
marcov:
Did you compile and run with range and overflow checks on?
howardpc:
There may be other faults, but this is clearly wrong:
--- Code: ---...
SetLength(x, i);
for j := 0 TO i do
x[j] := j+1;
...
--- End code ---
It should of course be
--- Code: --- ...
for j:=0 to i-1 do
...
--- End code ---
Edit: looking briefly at LS_Poly(), you make the same mistake with the upper bound of every dynamic array you access in your FOR loops, since the upper index you attempt to access is beyond the range of the array.
Also the line towards the end of LS_Poly
--- Code: --- for i:=0 to ... do
begin
...
c[i]:=f[i];
end;
--- End code ---
assumes that c[] and f[] have identical dimensions. Yet you initialise them with differing lengths.
noob:
i am kind of ashamed now.
this was the fix!
Thank You!!!http://forum.lazarus.freepascal.org/Smileys/ExcellentSmileys1/biggrin.gif
Navigation
[0] Message Index