/* Find a minimum spanning 1-tree with Held & Karp subgradient optimization.
   
   For your viewing pleasure, this script writes an SVG file every time it constructs a 1-tree.
   Set the processing priority of AMPL to be "below normal"
   (using Ctrl-Alt-Delete, Task List, Processes, right-click on AMPL). Then start the script.
   After AMPL has started, open the SVG file with your browser. (You'll need Adobe's free SVG plug-in.)
   Then you can refresh the screen frequently, to see the 1tree as AMPL runs.
   If your browser reports an error with the SVG file, try refreshing the page.

   In Internet Explorer, when you print or view the SVG file in 'Print Preview', the SVG graph may be distorted.
   Solution: When your SVG is distorted in the Print Preview window,
   close it with the 'X' button in the topright corner NOT the 'Close' button.
   View the Print Preview again and the image will display as intended.
   
   Great for 50 or 100 city problems. Slow for 200 nodes.
*/

# Script for the Held & Karp subgradient optimisation algorithm.
param N; # number of nodes.
set Nodes := {1..N} ordered;

set InTree within Nodes; 
set NotInTree within Nodes; 
set Arcs := {i in Nodes, j in Nodes: ord(i) < ord(j)};

param xcoord {Nodes};
param ycoord {Nodes};

data ..\data\50points.txt;

# Use this definition for Euclidean models, where the data gives point coordinates.
param length {(i,j) in Arcs} := sqrt((xcoord[i] - xcoord[j])^2 + (ycoord[i] - ycoord[j])^2);

# Use this definition for asymmetric models, where the data gives length explicitly.
# param length {(i,j) in Arcs};
param bitesize default 1; # Bitesize will gradually get smaller, to ensure that our stepsize goes to zero.
param degree {Nodes} integer default 0;
param iteration integer;
param lambda {Nodes} default 0;
param maxiterations integer;

param slacksq;
param smallest_arc;
param stepsize default 1; # The stepsize depends on the bitesize and the slacks. It will go up and down, but tends to 0.
param Tourlength default 0;
param upperbound default N*(sum {(i,j) in Arcs} length[i,j])/card(Arcs); # Could be any feasible tour length.
param weight {Arcs} default Infinity;

param u symbolic in Nodes;
param v symbolic in Nodes;
param w symbolic in Nodes;

# x[i,j]=1 if the arc (i,j) is in the tree, else 0.
var x {(i,j) in Arcs} binary; # Actually, this could be a parameter, since AMPL never calls a solver.
let iteration := 0;
let maxiterations := 400;

repeat while stepsize > 0.001 and iteration <= maxiterations
{	let iteration := iteration + 1;
	printf "Iteration %i.\n", iteration;

	# Update the stepsize based on the slacks.	
	let slacksq := 0; for {i in Nodes} let slacksq := slacksq + (degree[i] - 2)^2;
	let stepsize := bitesize*(upperbound - Tourlength )/slacksq; display stepsize;
	let bitesize := bitesize*0.96;

	# Now update the lambdas.
	for {i in Nodes} let lambda [i] := lambda [i] + stepsize*(degree[i] - 2);
	
	# Reset the degrees to zero, to prepare for the next loop, and decrease the stepsize.
	for {i in Nodes} let degree [i] := 0;
	
	# Calculate the weight on each arc as a function of its length and the current lambdas.
	for {(i,j) in Arcs} let weight [i,j] := length[i,j] + lambda [i] + lambda [j];

	# Build a minimum spanning tree.
	# Leave off first, since we start with it in the set InTree.
	# Leave off last, since we will use it to make a 1-tree at the end.
	# First find the MST for nodes 1 thru n-2.
	for {(i,j) in Arcs}	let x[i,j] := 0;
	let u := first(Nodes);
	let v := last(Nodes);
	let InTree := {v};
	let NotInTree := Nodes diff {v, u};
	
	for {1..card(Nodes) - 2}
	{	let smallest_arc := min ( min {i in InTree, j in NotInTree: (i,j) in Arcs} weight[i,j],
			min {i in InTree, j in NotInTree: (j,i) in Arcs} weight[j,i] );
	
		if smallest_arc = Infinity then
		{	printf "Disconnected graph, iteration %i. No spanning tree.\n", iteration; 
			break all;
		};
	
		for find {i in InTree, j in NotInTree: (i,j) in Arcs or (j,i) in Arcs} 
		{	if (i,j) in Arcs and weight[i,j] = smallest_arc then
			{	let InTree := InTree union {j}; 
				let NotInTree := NotInTree diff {j}; 
				let x[i,j] := 1; 
				break find;
			}
			if (j,i) in Arcs and weight[j,i] = smallest_arc then
			{	let InTree := InTree union {j};
				let NotInTree := NotInTree diff {j};
				let x[j,i] := 1;
				break find;
			} 
		}
	};

	# Add arc 1 to the last node u (still not in NotInTree and not in InTree) to complete the MST.
	let smallest_arc := min (min {i in InTree: (i,u) in Arcs} weight[i,u], min {i in InTree: (u,i) in Arcs} weight[u,i]);
	if smallest_arc = Infinity then
	{	printf "Could not add last node. No spanning tree.\n"; 
		break all;
	};

	for find {i in InTree: (i,u) in Arcs or (u,i) in Arcs} 
	{	if (i,u) in Arcs and weight[i,u] = smallest_arc then
		{	let InTree := InTree union {u}; 
			let x[i,u] := 1;
			# printf "Added %i, %i to 1-tree.\n", i, u;
			let w := i;
			break find;
		}
		if (u,i) in Arcs and weight[u,i] = smallest_arc then
		{	let InTree := InTree union {u};
			let x[u,i] := 1;
			# printf "Added %i, %i to 1-tree.\n", u, i;
			let w := i;
			break find;
		}
	}
	
	# Add arc 2 to the last node u (now in InTree) to make a 1-tree.
	let smallest_arc := min ( min {i in InTree: i <> w and (i,u) in Arcs} weight[i,u],
	                  min {i in InTree: i <> w and (u,i) in Arcs} weight[u,i] );
	if smallest_arc = Infinity then
	{	printf "Could not add last node. No 1-tree.\n"; 
		break all;
	};
	
	for find {i in InTree: i <> w and ((i,u) in Arcs or (u,i) in Arcs)} 
	{	if (i,u) in Arcs and weight[i,u] = smallest_arc then
		{	let InTree := InTree union {u}; 
			let x[i,u] := 1;
			# printf "Added %i, %i to 1-tree.\n", i, u;
			break find;
		}
		if (u,i) in Arcs and weight[u,i] = smallest_arc then
		{	let InTree := InTree union {u};
			let x[u,i] := 1;
			# printf "Added %i, %i to 1-tree.\n", u, i;
			break find;
		}
	}
	
	let Tourlength := sum {(i,j) in Arcs} length[i,j]*x[i,j];
	printf "Tree length is %f\n", Tourlength;
	
	# Calculate the degree of each node. (Student version may be too small for this.)
	for {(i,j) in Arcs: x[i,j]=1}
	{	let degree [i] := degree [i] + 1;
		let degree [j] := degree [j] + 1;
	}
	
	# If all nodes have degree 2, then we have a tour, so stop.
	if iteration >= maxiterations or min {i in Nodes} degree[i] = max {i in Nodes} degree[i] then break;

	commands ..\makesvg.run;
}

for {(i,j) in Arcs: x[i,j]=1} printf " %i %i\n", i, j;

for {i in 1..card(Nodes)} printf "Node %i has degree %i.\n", i, degree[i];

printf "Tree length is %f.\n", Tourlength;
printf "Sum of lambdas is %f.\n", sum {i in Nodes} lambda [i];
reset;