$title Test marginals for a scaled MCP problem (MCP11,SEQ=696)

$onText
Start with the simplest MCP model one can imagine:
   f(x) := x-c, x >= 0
   For this model f.lo = c.

Ignoring the degenerate case when c==0, we have:
  A) c > 0: implies
     x = c is the solution, f.l == f.lo, f.m = x.l = c
  B) c < 0: implies
     x = 0 is the solution, f.l == 0, f.m = x.l = 0, x.m = (-c)

For each case A) and B), we consider scaling the variable and the equation,
for a total of 4 cases.

Contributor: Steve Dirkse, May 2016
$offText

$if not set TESTTOL $set TESTTOL 1e-6
scalars tol / %TESTTOL% /;
scalar c;
positive variable x;
equation f,
         fNeg 'f times -1';
f..     x =G= c;
fNeg.. -x =L= -c;
model m    /  f.x /;
model mNeg / -fNeg.x /;

sets
  case / cpos, cneg /
  val  / 'f.F', 'f.m', 'x.L', 'x.m' /
  ;
parameters
  unscaled(case,val) 'results with no scaling'
  scaledX(case,val)  'results with x scaled'
  deltaX(case,val)   'unscaled - scaledX'
  scaledF(case,val)  'results with f scaled'
  deltaF(case,val)   'unscaled - scaledF'
  ;

* case A  c > 0
c = 2;
solve m using mcp;
abort$[abs(f.l-c)  > tol] 'bad f.l==c', f.l, c;
abort$[abs(f.m-c)  > tol] 'bad f.m==c', f.m, c;
abort$[abs(x.l-c)  > tol] 'bad x.l==c', x.l, c;
abort$[abs(x.m-0)  > tol] 'bad x.m==0', x.m;
unscaled('cpos','f.F') = round(f.L - f.lo, 5);
unscaled('cpos','f.m') = round(f.m, 5);
unscaled('cpos','x.L') = round(x.L, 5);
unscaled('cpos','x.m') = round(x.m, 5);

solve mNeg using mcp;
abort$[abs(fNeg.l+f.l)    > tol] 'bad fNeg.l ==-f.l',  fNeg.l,f.l;
abort$[abs(fNeg.up+f.lo)  > tol] 'bad fNeg.up==-f.lo', fNeg.up,f.lo;
abort$[abs(fNeg.lo+f.up)  > tol] 'bad fNeg.lo==-f.up', fNeg.lo,f.up;

* case B  c < 0
c = -4;
solve m using mcp;
abort$[abs(f.l-0)  > tol] 'bad f.l==0', f.l;
abort$[abs(f.m-0)  > tol] 'bad f.m==0', f.m;
abort$[abs(x.l-0)  > tol] 'bad x.l==0', x.l;
abort$[abs(x.m+c)  > tol] 'bad x.m==-c', x.m, c;
unscaled('cneg','f.F') = round(f.L - f.lo, 5);
unscaled('cneg','f.m') = round(f.m, 5);
unscaled('cneg','x.L') = round(x.L, 5);
unscaled('cneg','x.m') = round(x.m, 5);

solve mNeg using mcp;
abort$[abs(fNeg.l+f.l)    > tol] 'bad fNeg.l ==-f.l',  fNeg.l,f.l;
abort$[abs(fNeg.up+f.lo)  > tol] 'bad fNeg.up==-f.lo', fNeg.up,f.lo;
abort$[abs(fNeg.lo+f.up)  > tol] 'bad fNeg.lo==-f.up', fNeg.lo,f.up;

m.scaleopt = 1;
mNeg.scaleopt = 1;

* ========== row scaling ==========
f.scale = 10;
fNeg.scale = 10;
* case A  c > 0
c = 2;
solve m using mcp;
scaledF('cpos','f.F') = round(f.L - f.lo, 5);
scaledF('cpos','f.m') = round(f.m, 5);
scaledF('cpos','x.L') = round(x.L, 5);
scaledF('cpos','x.m') = round(x.m, 5);

solve mNeg using mcp;
abort$[abs(fNeg.l+f.l)    > tol] 'bad fNeg.l ==-f.l',  fNeg.l,f.l;
abort$[abs(fNeg.up+f.lo)  > tol] 'bad fNeg.up==-f.lo', fNeg.up,f.lo;
abort$[abs(fNeg.lo+f.up)  > tol] 'bad fNeg.lo==-f.up', fNeg.lo,f.up;

* case B  c < 0
c = -4;
solve m using mcp;
scaledF('cneg','f.F') = round(f.L - f.lo, 5);
scaledF('cneg','f.m') = round(f.m, 5);
scaledF('cneg','x.L') = round(x.L, 5);
scaledF('cneg','x.m') = round(x.m, 5);

solve mNeg using mcp;
abort$[abs(fNeg.l+f.l)    > tol] 'bad fNeg.l ==-f.l',  fNeg.l,f.l;
abort$[abs(fNeg.up+f.lo)  > tol] 'bad fNeg.up==-f.lo', fNeg.up,f.lo;
abort$[abs(fNeg.lo+f.up)  > tol] 'bad fNeg.lo==-f.up', fNeg.lo,f.up;

f.scale = 1;
fNeg.scale = 1;

* ========== col scaling ==========
x.scale = 100;
* case A  c > 0
c = 2;
solve m using mcp;
scaledX('cpos','f.F') = round(f.L - f.lo, 5);
scaledX('cpos','f.m') = round(f.m, 5);
scaledX('cpos','x.L') = round(x.L, 5);
scaledX('cpos','x.m') = round(x.m, 5);

solve mNeg using mcp;
abort$[abs(fNeg.l+f.l)    > tol] 'bad fNeg.l ==-f.l',  fNeg.l,f.l;
abort$[abs(fNeg.up+f.lo)  > tol] 'bad fNeg.up==-f.lo', fNeg.up,f.lo;
abort$[abs(fNeg.lo+f.up)  > tol] 'bad fNeg.lo==-f.up', fNeg.lo,f.up;

* case B  c < 0
c = -4;
solve m using mcp;
scaledX('cneg','f.F') = round(f.L - f.lo, 5);
scaledX('cneg','f.m') = round(f.m, 5);
scaledX('cneg','x.L') = round(x.L, 5);
scaledX('cneg','x.m') = round(x.m, 5);

solve mNeg using mcp;
abort$[abs(fNeg.l+f.l)    > tol] 'bad fNeg.l ==-f.l',  fNeg.l,f.l;
abort$[abs(fNeg.up+f.lo)  > tol] 'bad fNeg.up==-f.lo', fNeg.up,f.lo;
abort$[abs(fNeg.lo+f.up)  > tol] 'bad fNeg.lo==-f.up', fNeg.lo,f.up;

x.scale = 1;

deltaX(case,val) = unscaled(case,val) - scaledX(case,val);
deltaF(case,val) = unscaled(case,val) - scaledF(case,val);
deltaX(case,val) = round(deltaX(case,val), 5);
deltaF(case,val) = round(deltaF(case,val), 5);

execute_unload 'rrrr';

abort$((card(deltaX) + card(deltaF)) > 0) '** FAIL, check file rrrr.gdx for details **';