# Objective: linear
# Constraints: convex quadratic

########################################################################
# Works fine with presolve turned on.  Fails for some reason when 
# presolve turned off.
########################################################################

param n default 12;
param m default  6;	# must be even

param N := if (n>m) then n else m;

set X := {0..n};
set Y := {0..m};

set NODES := X cross Y;		# A lattice of Nodes

set ANCHORS within NODES := { x in X, y in Y : 
			      x == 0 
			      && y >= floor(m/3) && y <= m-floor(m/3) 
			      };

param load {(x,y) in NODES, d in 1..2} default 0;

param gcd {x in -N..N, y in -N..N} :=
    (if x <  0 then gcd[-x,y] else 
    (if x == 0 then y else 
    (if y < x  then gcd[y,x] else 
    (gcd[y mod x, x])
    )
    )
    );

set ARCS := { (xi,yi) in NODES, (xj,yj) in NODES: 
    abs( xj-xi ) <= 3 && abs(yj-yi) <=3 && abs(gcd[ xj-xi, yj-yi ]) == 1 
    && ( xi > xj || (xi == xj && yi > yj) )
    };

param length {(xi,yi,xj,yj) in ARCS} := sqrt( (xj-xi)^2 + (yj-yi)^2 );

var w {NODES, 1..2};

minimize work: 
    -sum {(x,y) in NODES, d in 1..2} load[x,y,d]*w[x,y,d];

subject to element_eqs {(xi,yi,xj,yj) in ARCS}: 
    ( (xj-xi)/length[xi,yi,xj,yj] * (w[xi,yi,1] - w[xj,yj,1]) +
      (yj-yi)/length[xi,yi,xj,yj] * (w[xi,yi,2] - w[xj,yj,2]))^2 
    <= length[xi,yi,xj,yj]^2;

subject to anchor {(x,y) in ANCHORS, d in 1..2}:
    w[x,y,d] = 0;

let load[n,m/2,2] := -1;

#let {(x,y) in NODES, d in 1..2} w[x,y,d] := 0.001;
##let {(x,y) in NODES, d in 1..2} w[x,y,d] := 1;

#option solver loqo;
#option loqo_options "verbose=2";
#option presolve 1;
solve;

param u {NODES, 1..2};
let {(x,y) in NODES, d in 1..2} u[x,y,d] := work * w[x,y,d];

drop work;
drop element_eqs;
drop anchor;

param Vol := 1000;

var sigma {ARCS} >= 0;

maximize zero: 0;

subject to compat {(xi,yi) in NODES, d in 1..2: (xi,yi) not in ANCHORS}:
    sum { (xi,yi,xj,yj) in ARCS }
    ((xj-xi)/length[xi,yi,xj,yj] * (u[xi,yi,1] - u[xj,yj,1]) +
     (yj-yi)/length[xi,yi,xj,yj] * (u[xi,yi,2] - u[xj,yj,2]))
    * sigma[xi,yi,xj,yj] 
    * ( if (d==1) then (xj-xi) else (yj-yi) )/length[xi,yi,xj,yj]^2

    -

    sum { (xk,yk,xi,yi) in ARCS }
    ((xi-xk)/length[xk,yk,xi,yi] * (u[xk,yk,1] - u[xi,yi,1]) +
     (yi-yk)/length[xk,yk,xi,yi] * (u[xk,yk,2] - u[xi,yi,2]))
    * sigma[xk,yk,xi,yi] 
    * ( if (d==1) then (xi-xk) else (yi-yk) )/length[xk,yk,xi,yi]^2

    = 

    load[xi,yi,d];

subject to tot_vol:
    sum { (xi,yi,xj,yj) in ARCS } length[xi,yi,xj,yj]*sigma[xi,yi,xj,yj] = Vol;

option loqo_options "verbose=2 sigfig=6";
solve;

option display_eps 0.0001;

printf: "%d \n", 
    card({(xi,yi,xj,yj) in ARCS: sigma[xi,yi,xj,yj] > 1.0e-8})
    > structure3.out;
printf {(xi,yi,xj,yj) in ARCS: sigma[xi,yi,xj,yj] > 1.0e-8}: 
    "%3d %3d %3d %3d %10.4f \n", xi, yi, xj, yj, sigma[xi,yi,xj,yj]
    > structure3.out;