ls03.gms : Test LS (Least Squares) utility - higher dimensions for GDX

Description

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


Category : GAMS Test library


Main file : ls03.gms

$TITLE Test LS (Least Squares) utility - higher dimensions for GDX (LS03,SEQ=397)

$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 /
  j2 / j2_1 * j2_2 /
  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 /
  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)
  ;

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) = yes;

SCALAR d;
PARAMETER
  b_(p,i4,j4,k4,i5,j5,k5,i6,j6,k6)
  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)
  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)
 * 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);

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;