option solver cplexamp;
model dfj-simple.mod;
data ..\..\data\20points.txt;

param epsilon := .0001; # tolerance for constraint violation.
param numcut default 0; # number of subtour elimination inequalities appended.
param t {(i,j) in Arcs} >= 0;
param from default 0;
param numInTour default 0;
param inTour {Nodes} default 0;

repeat
{	solve;

	# Find a subtour & add it to the model. This is a weak method - it only finds one subtour cut at a time.
	let {(i,j) in Arcs} t[i,j] := x[i,j];
	let {k in Nodes} inTour[k] := 0;
	let from := first(Nodes);
	let numInTour := 0;

	# Scan for a tour.
	for {k in Nodes}
	{	for {i in Nodes: i <> from}
		{	if (from,i) in Arcs and t[from,i] >= epsilon then
			{	let t[from,i] := 0.0; let inTour [i] := 1; let inTour [from] := 1;
				let from := i; let numInTour := numInTour + 1;
				break;
			}

			if (i, from) in Arcs and t[i, from] >= epsilon then
			{	let t[i,from] := 0.0; let inTour [i] := 1; let inTour [from] := 1;
				let from := i; let numInTour := numInTour + 1;
				break;
	}	}	}

	# If the subtour found is not a full tour, add the cut set. Else quit.
	if numInTour < N then
	{	let numcut := numcut + 1;
		for {k in Nodes: inTour[k] > epsilon} let Cutset [numcut] := Cutset [numcut] union {k};

		let numcut := numcut + 1;
		for {k in Nodes: inTour[k] = 0} let Cutset [numcut] := Cutset [numcut] union {k};
	}
	else break; # We end the loop when the numInTour = N.
};

# If the solution is binary, print the tour and quit.
if forall {(i,j) in Arcs} (x[i,j] <= epsilon or x[i,j] >= 1-epsilon) then
	for {i in Nodes, j in Nodes: (i,j) in Arcs and x[i,j] >= 1 - epsilon} printf "%i, %i\n", i, j;

commands ..\..\makesvg.run;