ls04.gms : Test LS (Least Squares) utility - higher dimensions for GDX
Run the LS solver on made-up data to see that it works when the
GDX output file contains params with many dimensions.
Contributor: Steve Dirkse, July 2008.
Small Model of Type: GAMS
$TITLE Test LS (Least Squares) utility - higher dimensions for GDX (LS04,SEQ=398)
$ontext
Run the LS solver on made-up data to see that it works when the
GDX output file contains params with many dimensions.
Contributor: Steve Dirkse, July 2008.
$offtext
SET
p 'unknown coefficients' / p0, p1 /
i1 / i1_1 * i1_2 /
j1 / j1_1 * j1_2 /
k1 / k1_1 * k1_2 /
i2 / i2_1 * i2_2 /
j2 / j2 /
k2 / k2_1 * k2_2 /
i3 / i3_1 * i3_2 /
j3 / j3_1 * j3_2 /
k3 / k3_1 * k3_2 /
i4 / i4 /
j4 / j4 /
k4 / k4 /
i5 / i5 /
j5 / j5 /
k5 / k5 /
i6 / i6 /
j6 / j6 /
k6 / k6 /
m 'over the edge' / m /
entire(i1,j1,k1,i2,j2,k2,i3,j3,k3,i4,j4,k4,i5,j5,k5,i6,j6,k6)
pp(p,i4,j4,k4,i5,j5,k5,i6,j6,k6,m)
;
entire(i1,j1,k1,i2,j2,k2,i3,j3,k3,i4,j4,k4,i5,j5,k5,i6,j6,k6) = yes;
pp(p,i4,j4,k4,i5,j5,k5,i6,j6,k6,m) = yes;
SCALAR d;
PARAMETER
b_(p,i4,j4,k4,i5,j5,k5,i6,j6,k6,m)
x (i1,j1,k1,i2,j2,k2,i3,j3,k3,i4,j4,k4,i5,j5,k5,i6,j6,k6) 'observed independent variable'
y (i1,j1,k1,i2,j2,k2,i3,j3,k3,i4,j4,k4,i5,j5,k5,i6,j6,k6) 'observed dependent variable'
fitted(i1,j1,k1,i2,j2,k2,i3,j3,k3,i4,j4,k4,i5,j5,k5,i6,j6,k6) 'fitted values for dependent variable'
;
b_(pp) = 1;
x(i1,j1,k1,i2,j2,k2,i3,j3,k3,i4,j4,k4,i5,j5,k5,i6,j6,k6) =
256 * (ord(i1)-1) + 128 * (ord(j1)-1) + 64 * (ord(k1)-1)
+ 32 * (ord(i2)-1) + 16 * (ord(j2)-1) + 8 * (ord(k2)-1)
+ 4 * (ord(i3)-1) + 2 * (ord(j3)-1) + 1 * (ord(k3)-1);
y(entire) = 1 * x(entire) + 1;
variables
b(p,i4,j4,k4,i5,j5,k5,i6,j6,k6,m)
sse 'sum of squared errors'
;
equation
fit(i1,j1,k1,i2,j2,k2,i3,j3,k3,i4,j4,k4,i5,j5,k5,i6,j6,k6) 'equation to fit'
sumsq
;
sumsq.. sse =n= 0;
fit(i1,j1,k1,i2,j2,k2,i3,j3,k3,i4,j4,k4,i5,j5,k5,i6,j6,k6)..
y(i1,j1,k1,i2,j2,k2,i3,j3,k3,i4,j4,k4,i5,j5,k5,i6,j6,k6) =e=
b('p1',i4,j4,k4,i5,j5,k5,i6,j6,k6,'m')
* x(i1,j1,k1,i2,j2,k2,i3,j3,k3,i4,j4,k4,i5,j5,k5,i6,j6,k6)
+ b('p0',i4,j4,k4,i5,j5,k5,i6,j6,k6,'m');
option lp = ls;
model someReallyBigModelNameToTestWithJustToTryIt /fit,sumsq/;
solve someReallyBigModelNameToTestWithJustToTryIt using lp minimizing sse;
option decimals=8;
display b.l;
abort$(someReallyBigModelNameToTestWithJustToTryIt.solvestat <> %solvestat.NormalCompletion%) 'solver return not normal';
d = smax{pp, abs(b.l(pp)-b_(pp))};
abort$[d > 1e-8] 'bad solution b.l', b.l, b_, d;
execute_load 'ls', fitted;
abort$[execerror > 0] 'Could not load statistics from GDX';
d = smax{entire, abs(y(entire) - fitted(entire))};
abort$[d > 1e-8] 'bad fitted', y, fitted, d;