# Objective: linear
# Constraints: linear

param m default 26;	# must be even
param n default 39;
#param m default  8;	# must be even
#param n default 12;

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 xload {(x,y) in NODES: (x,y) not in ANCHORS} default 0;
param yload {(x,y) in NODES: (x,y) not in ANCHORS} 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 comp {ARCS} >= 0;
var tens {ARCS} >= 0;

minimize volume: 
    sum {(xi,yi,xj,yj) in ARCS} 
        length[xi,yi,xj,yj]*(comp[xi,yi,xj,yj] + tens[xi,yi,xj,yj]);

subject to Xbalance {(xi,yi) in NODES: (xi,yi) not in ANCHORS}:
    sum { (xi,yi,xj,yj) in ARCS } 
    ((xj-xi)/length[xi,yi,xj,yj])*(comp[xi,yi,xj,yj]-tens[xi,yi,xj,yj])
    +
    sum { (xk,yk,xi,yi) in ARCS } 
    ((xi-xk)/length[xk,yk,xi,yi])*(tens[xk,yk,xi,yi]-comp[xk,yk,xi,yi])
    = 
    xload[xi,yi];
    ;

subject to Ybalance {(xi,yi) in NODES: (xi,yi) not in ANCHORS}:
    sum { (xi,yi,xj,yj) in ARCS } 
    ((yj-yi)/length[xi,yi,xj,yj])*(comp[xi,yi,xj,yj]-tens[xi,yi,xj,yj])
    +
    sum { (xk,yk,xi,yi) in ARCS } 
    ((yi-yk)/length[xk,yk,xi,yi])*(tens[xk,yk,xi,yi]-comp[xk,yk,xi,yi])
    = 
    yload[xi,yi];
    ;

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

option solver loqo;
option loqo_options "verbose=2";
solve;

option display_eps 0.0001;

#printf: "%3d %3d \n", m, n > structure.out;
printf: "%d \n", 
    card({(xi,yi,xj,yj) in ARCS: comp[xi,yi,xj,yj]+tens[xi,yi,xj,yj] > 1.0e-4})
    > structure.out;
printf {(xi,yi,xj,yj) in ARCS: comp[xi,yi,xj,yj] + tens[xi,yi,xj,yj] > 1.0e-4}: 
    "%3d %3d %3d %3d %10.4f \n", 
    xi, yi, xj, yj, tens[xi,yi,xj,yj] - comp[xi,yi,xj,yj]
    > structure.out;