# Developed by John F. Raffensperger, 2005. Updated 13 March 2006.

# Does your 2-matching model suffer from the heartbreak of fractionality?
# Do you stay awake at night, wishing you could tighten your TSP formulation?
# Do you hurt and ache from the lack of a simple way to learn about comb constraints?
# If so, then your suffering may be over. This AMPL script is just the thing you need!

# The IP formulation for the assignment problem is naturally integer, but, strangely enough,
# the natural IP formulation for 2-matching problem is not.
# 
# The file contains two models, one to solve the 2-matching LP relaxation, and the other to find 2-matching cuts.
# In many ways, this is quite similar to the DFJ algorithm.

param N; set Nodes:= {1..N}; # number of nodes. 
param iteration default 1;
param xcoord {Nodes}; param ycoord {Nodes};
set Arcs := {i in Nodes, j in Nodes: i < j};
param length {i in Nodes, j in Nodes: (i,j) in Arcs}; # := sqrt((xcoord[i] - xcoord[j])^2 + (ycoord[i] - ycoord[j])^2);

set W {(i,j) in Arcs, k in 1..iteration} default {}; # Holds sets of nodes.
set deltaW {(i,j) in Arcs, k in 1..iteration} dimension 2 default {}; # Holds sets of arcs with one end in W[i,j,k].
set F {(i,j) in Arcs, k in 1..iteration} dimension 2 default {}; # Holds sets of arcs.

var x {(i,j) in Arcs} >= 0 <= 1; # binary; # We are studying the LP relaxation.

# The 2-matching problem. ------------------------------------------------------------
minimize Tourlength: sum {(i,j) in Arcs} length[i,j] * x[i,j];

subject to Degree {i in Nodes}: sum {(i,j) in Arcs} x[i,j] + sum {(j,i) in Arcs} x[j,i] = 2;

# You should be able to use any of the following 3 types of comb constraints equivalently.
#subject to Comb {(i,j) in Arcs, t in 1..iteration: card(W[i,j,t])>=1}:
#	sum {(a,b) in Arcs: ((a in W[i,j,t] and b not in W[i,j,t]) 
#		or (a not in W[i,j,t] and b in W[i,j,t])) and (a,b) not in F[i,j,t]} x[a,b]
#		+ card(F[i,j,t]) - sum {(n,o) in F[i,j,t]} x[n,o] >= 1;

#subject to Comb {(i,j) in Arcs, t in 1..iteration: card(W[i,j,t])>=1}:
#	sum {(a,b) in deltaW[i,j,t]: (a,b) not in F[i,j,t]} x[a,b]
#		+ card(F[i,j,t]) - sum {(n,o) in F[i,j,t]} x[n,o] >= 1;

# Nemhauser & Wolsey definition, p. 276.
subject to Comb {(i,j) in Arcs, t in 1..iteration: card(W[i,j,t])>=2}:
	sum {(a,b) in Arcs: a in W[i,j,t] and b in W[i,j,t]} x[a,b]
		+ sum {(n,o) in F[i,j,t]} x[n,o]
		<= card(W[i,j,t]) + floor (card(F[i,j,t])/2);

# The max flow subproblem. ------------------------------------------------------------
set Arcs1 := {i in Nodes, j in Nodes: i<>j};
param bound {Arcs1};
param source symbolic in Nodes; param sink symbolic in Nodes;

var w {(i,j) in Arcs1} >= 0, <= bound[i,j];

maximize Maxflow: sum {(source,j) in Arcs1} w[source,j] - sum {(j,source) in Arcs1} w[j,source];

subject to Conserveflow {i in Nodes: i <> source and i <> sink}: sum {(j,i) in Arcs1} w[j,i] - sum {(i,j) in Arcs1} w[i,j] = 0;