# Objective: linear
# Constraints: convex quadratic

# option loqo_options="primal verbose=2 timing=1"

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

var eps >= 1.0e-8;

param n default 12;
param m default 6;	# must be a multiple of 6

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 f {(x,y) in NODES, d in 1..2} default 0;  # external forces 

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 v, default 1000;	# total volume of material to use

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

var x {ARCS} >= 0;
var t {ARCS} >= 0;
var y {ARCS};

minimize sum_of_t: 
    eps + sum {(xi,yi,xj,yj) in ARCS} t[xi,yi,xj,yj];

subject to load_constraints {(xk,yk) in NODES, d in 1..2}:
    sum {(xi,yi,xk,yk) in ARCS} y[xi,yi,xk,yk]*
	(if d==1 then xi-xk else yi-yk)/length[xi,yi,xk,yk]^1.5
    +
    sum {(xk,yk,xj,yj) in ARCS} y[xk,yk,xj,yj]*
	(if d==1 then xj-xk else yj-yk)/length[xk,yk,xj,yj]^1.5
    =
    f[xk,yk,d];

subject to second_order_convex {(xi,yi,xj,yj) in ARCS}:
    sqrt(eps^2 + 4*y[xi,yi,xj,yj]^2 + (t[xi,yi,xj,yj] - x[xi,yi,xj,yj])^2)
    <=
    t[xi,yi,xj,yj] + x[xi,yi,xj,yj];

subject to volume:
    sum  {(xi,yi,xj,yj) in ARCS} length[xi,yi,xj,yj]*x[xi,yi,xj,yj] <= v;

let eps := 1;

let f[n,  m/2,2] := -1;
let f[0,  m/3,2] :=  0.5;
let f[0,2*m/3,2] :=  0.5;

let f[0,  m/3,1] :=  3*n/m;
let f[0,2*m/3,1] := -3*n/m;

let {(xi,yi,xj,yj) in ARCS} t[xi,yi,xj,yj] := 2;
let {(xi,yi,xj,yj) in ARCS} x[xi,yi,xj,yj] := 1;
let {(xi,yi,xj,yj) in ARCS} y[xi,yi,xj,yj] := 1;

solve;

display card({(xi,yi,xj,yj) in ARCS});

display card({(xi,yi,xj,yj) in ARCS: 4*y[xi,yi,xj,yj]^2 + (t[xi,yi,xj,yj] -
x[xi,yi,xj,yj])^2 <= 1.0e-6});

option display_eps 0.0001;

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

let eps := 0;
display sum_of_t;