# 90-45 degree grid flat origami model
# Imagiro 74
# Goran Konjevod, March 2006
# comments appreciated: goran@asu.edu
#
# tested with GLPK

# models the right triangle square grid
set NN := 0..2; # each triangle may have neighbors 0,1 and 2

param n integer default 16; # size of unfolded grid (16=4*4)
# sets: V, E define the input graph/grid
#  V = the set of basic grid polygons
set V := 1..n;
param nside := sqrt(n/4); # number of squares along each side of the grid

param m integer default 4; # size of fold grid (4=2*2)
#  W = the set of basic (unfolded) polygons in the fold projection
#  F = the adjacency relation on W
set W := 1..m;
param mside := sqrt(m/4); # number of squares along each side

# Define neighbors
param NeighborV{(v,i) in V cross NN} default
  if ((i=1) and (v mod 4 = 0 )) then v-3 else
  if ((i=1) and (v mod 4 != 0)) then v+1 else
  if ((i=2) and (v mod 4 = 1 )) then v+3 else
  if ((i=2) and (v mod 4 != 1)) then v-1 else
# if we get to here, then i=0 so no need to check
  if ((v mod 4 = 0) and (v <= n-4*nside)) then v+(4*nside-2) else
  if ((v mod 4 = 1) and (v mod nside != 1)) then v-2 else
  if ((v mod 4 = 2) and (v > 4*nside)) then v-(4*nside-2) else
  if ((v mod 4 = 3) and ((v+1) mod nside != 0)) then v+2
  else 0;

# Do the same thing for W
param NeighborW{(v,i) in W cross NN} default
  if ((i=1) and (v mod 4 = 0 )) then v-3 else
  if ((i=1) and (v mod 4 != 0)) then v+1 else
  if ((i=2) and (v mod 4 = 1 )) then v+3 else
  if ((i=2) and (v mod 4 != 1)) then v-1 else
# if we get to here, then i=0 so no need to check
  if ((v mod 4 = 0) and (v <= m-4*mside)) then v+(4*mside-2) else
  if ((v mod 4 = 1) and (v mod mside != 1)) then v-2 else
  if ((v mod 4 = 2) and (v > 4*mside)) then v-(4*mside-2) else
  if ((v mod 4 = 3) and ((v+1) mod mside != 0)) then v+2
  else 0;

#  E = the adjacency relation on V (list of possible creases)
set E :={(v1,v2) in V cross V: exists {i in NN} v2=NeighborV[v1,i] and v1!=v2};

# similarly for F: (list of adjacent triangle pairs)
set F :={(w1,w2) in W cross W: exists {i in NN} w2=NeighborW[w1,i] and w1!=w2};


#  XF (\subseteq E^2) = the exclusion relation on E
#  (the confusing name stands for excluded folds)
#   At most one of the two creases in each pair of X can be folded.
set XF :={(v1,v2,w1,w2) in E cross E: (v1,v2) in E and (w1,w2) in E and
  (((v1 mod 4 = 1) and (v2 mod 4 = 0) and (w1 mod 4 = 1) and (w2 mod 4 = 2))
   or ((v1 mod 4 = 1) and (v2 mod 4 = 0) and (w1 mod 4 = 3) and (w2 mod 4 = 0))
   or ((v1 mod 4 = 2) and (v2 mod 4 = 3) and (w1 mod 4 = 1) and (w2 mod 4 = 2))
   or ((v1 mod 4 = 2) and (v2 mod 4 = 3) and (w1 mod 4 = 3) and (w2 mod 4 = 0))
  )};
# FF = forced crease pairs (forced to be equal)
set FF :={(v1,v2,w1,w2) in E cross E: (v1,v2) in E and (w1,w2) in E and
  (((v1 mod 4 = 1) and (v2 mod 4 = 0) and (w1 mod 4 = 2) and (w2 mod 4 = 3))
   or ((v1 mod 4 = 1) and (v2 mod 4 = 2) and (w1 mod 4 = 3) and (w2 mod 4 = 0))
  )};


# layers
param Lmax default n/2;
set L := 1..Lmax;

# where which neighbor of v ends up if not folded depends on
# the orientation of v;  Flip thus describes the behavior 
# of the neighborhood ordering under orientation-change
param Flip{i in NN} := if i=0 then 0 else 3-i;


######################################################################
# Variables and constraints
# (This is where the model proper begins)

# creases
var fv{(v1,v2) in V cross V: (v1,v2) in E} integer >=0, <=1;
var fm{(v1,v2) in V cross V: (v1,v2) in E} integer >=0, <=1;

# location assignment
var x{(v,w) in V cross W} integer >=0, <=1;

# flip orientation sigma
var s{V} integer >=0, <=1; 

# layer assignment
var l{V cross L} integer >=0, <=1;

# lambda
var lam{V cross W cross L} integer >=0, <=1;

var ll{V} integer >=0, <=Lmax;

var above{(u,v) in V cross V} integer >=0, <=1;

#
######################################################################
# Constraints

# Edges are undirected:
subject to symmetricfm{(v1,v2) in E}:
  fm[v1,v2] = fm[v2,v1];
subject to symmetricfv{(v1,v2) in E}:
  fv[v1,v2] = fv[v2,v1];
# Each crease has at most one sense:
subject to creasefoldedornot{(v1,v2) in E}:
  fv[v1,v2]+fm[v1,v2] <= 1;
# Some crease pairs are forced equal:
subject to forcedfoldsV{(v1,v2,w1,w2) in FF}:
  fv[v1,v2] = fv[w1,w2];
subject to forcedfoldsM{(v1,v2,w1,w2) in FF}:
  fm[v1,v2] = fm[w1,w2];
# Some crease pairs are exclusive:
subject to excludedfolds{(v1,v2,w1,w2) in XF}:
  fv[v1,v2] +fm[v1,v2] + fv[w1,w2] + fm[w1,w2] <= 1;



# Each v\in V is folded onto some w\in W:
subject to basicgridcellassignment{v in V}:
  sum{w in W} x[v,w] = 1;

# Flip orientation propagates through folds:
subject to sfconsistentA{(u,v) in E}:
  s[u] - s[v] <= fm[u,v] + fv[u,v];
subject to sfconsistentB{(u,v) in E}:
  s[v] - s[u] <= fm[u,v] + fv[u,v];
subject to sfconsistentC{(u,v) in E}:
  s[u] + s[v] >= fm[u,v] + fv[u,v];
subject to sfconsistentD{(u,v) in E}:
  s[u] + s[v] <= 2 - fm[u,v] - fv[u,v];


# Location propagates through folds:
# A0: no fold, no flip
subject to locationconsistentA0{(v,w,i) in V cross W cross NN:
  (exists {vv in V} vv=NeighborV[v,i]) and 
  (exists {ww in W} ww=NeighborW[w,i])}:
    x[NeighborV[v,i],NeighborW[w,i]] 
      >= x[v,w] + 1-s[v]
         - (fm[v,NeighborV[v,i]]+fv[v,NeighborV[v,i]]) - 1;
# A1: no fold, flip
subject to locationconsistentA1{(v,w,i) in V cross W cross NN:
  (exists {vv in V} vv=NeighborV[v,Flip[i]]) and
  (exists {ww in W} ww=NeighborW[w,i])}:
    x[NeighborV[v,Flip[i]],NeighborW[w,i]] 
      >= x[v,w] + s[v]
         - (fm[v,NeighborV[v,Flip[i]]]+fv[v,NeighborV[v,Flip[i]]]) - 1;
# B: fold, doesn't matter if flipped or not (only w is involved)
subject to locationconsistentB{(v,w,i) in V cross W cross NN:
  exists {vv in V} vv=NeighborV[v,i]}:
    x[NeighborV[v,i],w] 
      >= x[v,w] + (fv[v,NeighborV[v,i]] + fm[v,NeighborV[v,i]]) - 1;




# Layers

# First: relate lam to l:
subject to llam1{(i,j,k) in V cross W cross L}:
  lam[i,j,k] <= x[i,j];
subject to llam2{(i,j,k) in V cross W cross L}:
  lam[i,j,k] <= l[i,k];
subject to llam3{(i,j,k) in V cross W cross L}:
  lam[i,j,k] >= x[i,j] + l[i,k] - 1;

# Layers fill up from bottom:
subject to layersfillup{(w,k) in W cross L: k>1}:
  sum{v in V} lam[v,w,k] <= sum{v in V} lam[v,w,k-1];
# Every mapped polygon is at some layer:
subject to layersVcomplete{(v,w) in V cross W}:
  x[v,w] = sum{k in L} lam[v,w,k];
# At most one polygon at layer 1; for higher ones implied by layersfillup
subject to layersunique{w in W}: 
  sum{v in V} lam[v,w,1] <= 1;

# Relate l to ll:
subject to lldefined{v in V}:
  ll[v] = sum{k in L} k*l[v,k];


# above is 0 for (v,v):
subject to fixabove0{v in V}:
  above[v,v] = 0;

# Relate above to ll:
subject to layerorderdefined{(u,v) in V cross V: u!=v}:
  above[u,v] >= (ll[u] - ll[v])/Lmax;

# Above antisymmetry:
subject to abovesymmetry{(u,v) in V cross V : u!=v}:
  above[u,v] + above[v,u] <= 1;

# Layers consistent with folds:
subject to layerfoldconsistency0v{(u,v) in E}:
  above[v,u] >= fv[u,v] - s[u];
subject to layerfoldconsistency1v{(u,v) in E}:
  above[u,v] >= fv[u,v] + s[u] - 1;
subject to layerfoldconsistency0m{(u,v) in E}:
  above[u,v] >= fm[u,v] - s[u];
subject to layerfoldconsistency1m{(u,v) in E}:
  above[v,u] >= fm[u,v] + s[u] - 1;


# Non-crossing constraints
#
# Type X: two unfolded edges coincide: no intersection
# cases: neighbors may flip, depending on s[v1], s[v3]
subject to noXintersection00{(v1,v3,w,i) in 
  V cross V cross W cross NN:
  v1 != v3
  and exists{v2 in V} v2=NeighborV[v1,i] 
  and exists{v4 in V} v4=NeighborV[v3,i]}:
    above[NeighborV[v3,i],NeighborV[v1,i]]
    >= above[v3,v1] + x[v1,w] + x[v3,w]
       + 1 - (fv[v1,NeighborV[v1,i]] + fm[v1,NeighborV[v1,i]])
       + 1 - (fv[v3,NeighborV[v3,i]] + fm[v3,NeighborV[v3,i]])
	   + 1 - s[v1]
	   + 1 - s[v3]
       - 6;
subject to noXintersection11{(v1,v3,w,i) in 
  V cross V cross W cross NN:
  v1 != v3
  and exists{v2 in V} v2=NeighborV[v1,i] 
  and exists{v4 in V} v4=NeighborV[v3,i]}:
    above[NeighborV[v3,i],NeighborV[v1,i]]
    >= above[v3,v1] + x[v1,w] + x[v3,w]
       + 1 - (fv[v1,NeighborV[v1,i]] + fm[v1,NeighborV[v1,i]])
       + 1 - (fv[v3,NeighborV[v3,i]] + fm[v3,NeighborV[v3,i]])
	   + s[v1]
       + s[v3]
       - 6;
subject to noXintersection10{(v1,v3,w,i) in 
  V cross V cross W cross NN:
  v1 != v3
  and exists{v2 in V} v2=NeighborV[v1,Flip[i]] 
  and exists{v4 in V} v4=NeighborV[v3,i]}:
    above[NeighborV[v3,i],NeighborV[v1,Flip[i]]]
    >= above[v3,v1] + x[v1,w] + x[v3,w]
       + 1 - (fv[v1,NeighborV[v1,Flip[i]]] + fm[v1,NeighborV[v1,Flip[i]]])
       + 1 - (fv[v3,NeighborV[v3,i]] + fm[v3,NeighborV[v3,i]])
	   + s[v1]
       + 1 - s[v3]
       - 6;

subject to noXintersection01{(v1,v3,w,i) in 
  V cross V cross W cross NN:
  v1 != v3
  and exists{v2 in V} v2=NeighborV[v1,i] 
  and exists{v4 in V} v4=NeighborV[v3,Flip[i]]}:
    above[NeighborV[v3,Flip[i]],NeighborV[v1,i]]
    >= above[v3,v1] + x[v1,w] + x[v3,w]
       + 1 - (fv[v1,NeighborV[v1,i]] + fm[v1,NeighborV[v1,i]])
       + 1 - (fv[v3,NeighborV[v3,Flip[i]]] + fm[v3,NeighborV[v3,Flip[i]]])
	   + 1 - s[v1]
       + s[v3]
       - 6;


# Type Y: one unfolded edge cannot cross a folded edge
subject to noYintersection00{(v1,v3,w,i) in
  V cross V cross W cross NN:
  v1 != v3
  and exists{v2 in V} v2=NeighborV[v1,i]
  and exists{v4 in V} v4=NeighborV[v3,i]}:
	above[v3,NeighborV[v1,i]]
	>=  above[v3,v1] + x[v1,w] + x[v3,w]
       + fv[v1,NeighborV[v1,i]]
       + 1 - (fv[v3,NeighborV[v3,i]] + fm[v3,NeighborV[v3,i]])
	   + 1 - s[v1]
	   + 1 - s[v3]
       - 6;
subject to noYintersection11{(v1,v3,w,i) in
  V cross V cross W cross NN:
  v1 != v3
  and exists{v2 in V} v2=NeighborV[v1,i]
  and exists{v4 in V} v4=NeighborV[v3,i]}:
	above[v3,NeighborV[v1,i]]
	>=  above[v3,v1] + x[v1,w] + x[v3,w]
       + fm[v1,NeighborV[v1,i]]
       + 1 - (fv[v3,NeighborV[v3,i]] + fm[v3,NeighborV[v3,i]])
	   + s[v1]
	   + s[v3]
       - 6;
subject to noYintersection01{(v1,v3,w,i) in
  V cross V cross W cross NN:
  v1 != v3
  and exists{v2 in V} v2=NeighborV[v1,Flip[i]]
  and exists{v4 in V} v4=NeighborV[v3,i]}:
	above[v3,NeighborV[v1,Flip[i]]]
	>=  above[v3,v1] + x[v1,w] + x[v3,w]
       + fm[v1,NeighborV[v1,Flip[i]]]
       + 1 - (fv[v3,NeighborV[v3,i]] + fm[v3,NeighborV[v3,i]])
	   + s[v1]
	   + 1 - s[v3]
       - 6;
subject to noYintersection10{(v1,v3,w,i) in
  V cross V cross W cross NN:
  v1 != v3
  and exists{v2 in V} v2=NeighborV[v1,i]
  and exists{v4 in V} v4=NeighborV[v3,Flip[i]]}:
	above[v3,NeighborV[v1,i]]
	>=  above[v3,v1] + x[v1,w] + x[v3,w]
       + fv[v1,NeighborV[v1,i]]
       + 1 - (fv[v3,NeighborV[v3,Flip[i]]] + fm[v3,NeighborV[v3,Flip[i]]])
	   + 1 - s[v1]
	   + s[v3]
       - 6;

# Type W: two folded edges not crossing
subject to noWintersection00top{(v1,v3,w,i) in
  V cross V cross W cross NN:
  v1 != v3
  and exists{v2 in V} v2=NeighborV[v1,i] 
  and exists{v4 in V} v4=NeighborV[v3,i]}:
	above[NeighborV[v1,i],NeighborV[v3,i]] 
	>= above[v3,v1] + above[NeighborV[v1,i],v3] + x[v1,w] + x[v3,w]
       + fv[v1,NeighborV[v1,i]]
       + fv[v3, NeighborV[v3,i]] + fm[v3,NeighborV[v3,i]]
       + 1 - s[v1]
	   + 1 - s[v3]
       - 7;
subject to noWintersection00bottom{(v1,v3,w,i) in
  V cross V cross W cross NN:
  v1 != v3
  and exists{v2 in V} v2=NeighborV[v1,i] 
  and exists{v4 in V} v4=NeighborV[v3,i]}:
	above[NeighborV[v3,i],v1] 
	>= above[v3,v1] + above[NeighborV[v1,i], v3] + x[v1,w] + x[v3,w]
       + fv[v1,NeighborV[v1,i]]
       + fv[v3, NeighborV[v3,i]] + fm[v3,NeighborV[v3,i]]
       + 1 - s[v1]
	   + 1 - s[v3]
       - 7;
subject to noWintersection01top{(v1,v3,w,i) in
  V cross V cross W cross NN:
  v1 != v3
  and exists{v2 in V} v2=NeighborV[v1,Flip[i]] 
  and exists{v4 in V} v4=NeighborV[v3,i]}:
	above[NeighborV[v1,Flip[i]],NeighborV[v3,i]] 
	>= above[v1,v3] + above[v3,NeighborV[v1,Flip[i]]] + x[v1,w] + x[v3,w]
       + fv[v1,NeighborV[v1,Flip[i]]]
       + fv[v3, NeighborV[v3,i]] + fm[v3,NeighborV[v3,i]]
       + s[v1]
	   + 1 - s[v3]
       - 7;
subject to noWintersection01bottom{(v1,v3,w,i) in
  V cross V cross W cross NN:
  v1 != v3
  and exists{v2 in V} v2=NeighborV[v1,Flip[i]]
  and exists{v4 in V} v4=NeighborV[v3,i]}:
	above[NeighborV[v3,i],v1]
	>= above[v1,v3] + above[v3,NeighborV[v1,Flip[i]]] + x[v1,w] + x[v3,w]
       + fv[v1,NeighborV[v1,Flip[i]]]
       + fv[v3, NeighborV[v3,i]] + fm[v3,NeighborV[v3,i]]
       + s[v1]
	   + 1 - s[v3]
       - 7;


end;