/* Important note: This program is written in magma, so you need this software on 
your computer. It also makes use of the geng program in
Brendan McKay's graph toolkit and nauty. You can download these two programs from
http://cs.anu.edu.au/people/bdm/nauty/

The following command sets a path variable to point to the authors copy
of the program called geng. You need to change it to point to your copy.
*/

gengpath := "/u/keevash/gtools10/geng";


/* Written by Peter Keevash - January 2002.

1. Purpose of program: find smallest value of a maximal fractional packing of 
monochromatic triangles in a 2-coloured complete graph.

2. Algorithm (description taken from paper): 

Initialisation: Choose a starting value for $n$, the number of vertices, and $d$,
the `search depth'. Check all graphs on $n$ vertices. Find the $d$ smallest values
of $t^*$ that occur. Build a list $L_n$ of all graphs for which $t^*$ is one of the
$d-1$ smallest values. (Our lists will always be `reduced': containing
no isomorphic or complementary pairs.)

Iteration: Check all one vertex extensions of the graphs in $L_n$. By averaging,
any other graph on $n+1$ vertices must have $t^*$ at least $\frac{n+1}{n-1}$ times the
$d$th smallest value for graphs on $n$ vertices. We refer to this lower bound as the 
{\em level} at $n+1$. The smallest values below the level that we find by checking
these extensions are indeed the smallest values for all graphs on $n+1$ vertices.
We find as many of the smallest values as possible up to a maximum of $d$.

If there are $d$ values below the level, then we let $L_{n+1}$ be the graphs with
the $d-1$ smallest values and procede with the iteration. However, there
may not be $d$ values below the level, so then we have to reduce the depth of the
search: we let $L_{n+1}$ be all the graphs we have found below the level, and
calculate the level for $n+2$ to be $\frac{n+2}{n}$ times the level for $n+1$.
Then we procede with the iteration. 

Termination: There may come a point at which no values are found
below the level, at which point nothing further can be said. To get results for
larger $n$ it is then necessary to start again with a larger search depth.

3. Operation: the user is prompted to input 3 starting variables - the starting value
of N, the search depth d, and the value of N to stop at.

4. Remarks on implementation:

The function Fract(G) computes value of maximum fractional packing of triangles
in G and its complement. This is a linear programming problem, for which magma
has a built-in solver. Unfortunately, it will only return (approximate) real solutions
rather than (exact) rational solutions, so it is necessary to use the procedure
Rat, which reconstructs the rational solution under the assumption that is has
denominator dividing 360360, which suffices for most purposes. */



Rat := function(x)

// Reconstructs rational from real approximation assuming denominator divides 360360

	err := false;
	y := 360360*x;
	temp := Round(y);
	if (Abs(y-temp) gt 0.1) then
		err := true;
	end if;
	return temp/360360,err;

end function;

Red := procedure(~glist)

// eliminates isomorphs and complementary pairs from graph list

	length := #glist;
	x := 1;
	while (x lt length) do
		H1 := glist[x];
		H2 := Complement(H1);
		y := x+1;
		while (y le length) do
			H3 := glist[y]; 
			if ( IsIsomorphic(H1,H3 : IgnoreLabels := true) or IsIsomorphic(H2,H3 : IgnoreLabels := true) ) then
				Remove(~glist,y);
				length -:= 1;
			else
				y +:= 1;	
			end if;
		end while;
		x +:= 1;
	end while;

end procedure;

Fract := function(G);

// finds value of maximum fractional packing of triangles in G and its complement

	V := VertexSet(G);
	E := EdgeSet(G);
	N := #V;
	EDG := N*(N-1) div 2;			// number of edges	NB. Slightly wasteful to compute each time!
	TR := N*(N-1)*(N-2) div 6;		// number of triples

	trips := [ x: x in Subsets( {1..N} , 3) ];	// potential triangles	(empty or complete)
	edges := [ x: x in Subsets( {1..N} , 2) ];

	/* We run through all triples. When we find a complete or empty
	triangle we:
	 1. Add its `trips' index to a sequence `triangle'
	 2. For each of its edges we:
	  i) Set the bit of `edge_is_real' corresponding to its index in `edges' to 1
	  ii) Add the number of the triangle to a sequence belonging to the edge.

	Further explanation of 2.ii: `triangle_on_edge' is a sequence of sequences:
	  The first parameter is the index of an edge in `edges'
	  The second parameter is the number of the triangle  */

	ntriang := 0;
	triangle := [];						// list of triangles
	edge_is_real := [ 0: x in [1..EDG]];			// indicates edge in triangle
	triangle_on_edge := [ []: x in [1..EDG]];		// sequence of triangles on edge
	for x in { 1..TR } do
		current_trip := trips[x];
        	H := sub< G | {V!i : i in current_trip} >;
        	if IsComplete(H) or IsEmpty(H) then
			ntriang +:= 1;
			Append(~triangle, x);			// add triangle to list
			current_edges := Subsets(current_trip,2);
			for y in current_edges do
				current_index := Index(edges,y);
				edge_is_real[current_index] := 1;  // mark y as a real edge
				Append( ~(triangle_on_edge[current_index]), ntriang );  // add triangle to y list 
			end for;
		end if;
	end for;

	if (ntriang eq 0) then		// only happens for N<=5
		return 0;
	end if;	

	//  Now we count the edges that belong to triangles

	nreal_edges := 0;
	real_edge := [];
	for x in { 1..EDG } do
		if (edge_is_real[x] eq 1) then
			nreal_edges +:= 1;
			Append(~real_edge, x);
		end if;
	end for;

	//  Now build the edge-triangle incidence matrix 

	R := RealField();

	lhs := Matrix(R, nreal_edges, ntriang, [ 0 : i in [1..nreal_edges*ntriang] ] );
	for x in { 1..nreal_edges } do
		edge_index := real_edge[x];
		for y in triangle_on_edge[edge_index] do	;
			lhs[ x, y ] := 1;
		end for;
	end for;

	//  Solve the linear program

	rhs := Matrix(R, nreal_edges, 1, [ 1 : i in [1..nreal_edges] ] );
	rel := Matrix(R, nreal_edges, 1, [ -1 : i in [1..nreal_edges] ] );	// negative values - leq relation
	obj := Matrix(R, 1, ntriang, [ 1 : i in [1..ntriang] ] );

	vector := MaximalSolution( lhs, rel, rhs, obj )[1];
	value := 0;
	for x in {1..ntriang} do
		value +:= vector[x];
	end for;

	return value;

end function;

exhaust := function(N,k)

/* Checks all graphs on N vertices looking for best k values.
Output is the number of values given (up to max k), the values, and
a list of graphs for up to k-1 values. */

	best := [];
	graphs := [];
	nvalues := 0;

	F := OpenGraphFile("cmd " cat gengpath cat " " cat IntegerToString(N),0,0);

	while true do

		more, G := NextGraph(F);
		if (not more) then
	    		break;		// reached end of F
		end if;
		
		realvalue := Fract(G);
		value, err := Rat(realvalue);		// convert to rational
		if err then
			print "Rounding error: ",realvalue," to ",value;
			print G;
		end if;

		// Find best[index] <= value < best[index+1], with equality flag 

		index := nvalues;
		equality := false;
		for j in {1..nvalues} do
			test := best[j];
			if (value le test) then
			   if (value eq test) then
				index := j;
				equality := true;
				break;
			   else
				index := j-1;
				equality := false;
				break;
			   end if;
			end if;
		end for;

		// update best and graphs
		
		if equality then
			Append(~graphs[index],G);
		else
		   if (index lt k) then
			if (nvalues lt k) then
				nvalues +:= 1;
			end if;
			Insert(~best,index+1,value);
			Insert(~graphs,index+1,[G]);
		   end if;
		end if;
	

	end while;

	returnlist := [];
	nlists := Minimum( {nvalues,k-1} );
	for j in {1..nlists} do
		Red(~graphs[j]);
		returnlist cat:= graphs[j];
	end for;

	return nvalues,best,returnlist;

end function;


fullextend := function(N,list,k,level)

/* Checks one vertex extensions of list (a list of graphs
on N-1 vertices) looking for the best values
<= level, up to a maximum of k.
Output is the number of values given, the values, and
a list of graphs for all but the worst value. */

	best := [];
	graphs := [];
	nvalues := 0;
	
	length := #list;

	for i in { 1..length } do

		G := list[i];

	/* One vertex extensions of G are represented by N-1 bits, indicating adjacency to
	the Nth vertex. We check in binary order. */

		bit := [0 : j in {1..N}];		// Nth bit indicates overflow
		AddVertex(~G,N);

		while (bit[N] eq 0) do
		
			realvalue := Fract(G);
			value, err := Rat(realvalue);		// convert to rational
			if err then
				print "Rounding error: ",realvalue," to ",value;
				print G;
			end if;
			
			if (value le level) then

			  // Find best[index] <= value < best[index+1], with equality flag 

			  index := nvalues;
			  equality := false;
			  for j in {1..nvalues} do
				test := best[j];
				if (value le test) then
				   if (value eq test) then
					index := j;
					equality := true;
					break;
				   else
					index := j-1;
					equality := false;
					break;
				   end if;
				end if;
			  end for;

			  // update best and graphs
			
			  if equality then
				Append(~graphs[index],G);
			  else
			     if (index lt k) then
				if (nvalues lt k) then
					nvalues +:= 1;
				end if;
				Insert(~best,index+1,value);
				Insert(~graphs,index+1,[G]);
		  	     end if;
			  end if;

			end if;

			// move to next graph

			carry := true;
			x := 0;
			while carry do
				x +:= 1;
				bit[x] +:= 1;
				if (bit[x] eq 1) then
					carry := false;
					AddEdge(~G,x,N);
				else
					bit[x] := 0;
					RemoveEdge(~G,x,N);
				end if;
			end while;
		end while;
	end for;

	returnlist := [];
	nlists := Minimum( {nvalues,k-1} );
	for j in {1..nlists} do
		Red(~graphs[j]);
		returnlist cat:= graphs[j];
	end for;

	return nvalues,best,returnlist;

end function;

// start of program

readi d,"Search depth: ";
readi N,"Value to exhaust from: ";
readi Nmax,"Stop at: ";

t := Cputime();

nvals,best,graphs:= exhaust(N,d);
print N," vertices: ",[best[j] : j in [1..nvals]],"  ",Cputime(t)," seconds.";
level := best[nvals];

N +:= 1;
while (N le Nmax) do

	if (nvals eq d) then
		level := best[nvals]*N/(N-2);
	else
		level := level*N/(N-2);
	end if;
	print "level: ",level;		// minimum value if not extension of member of graphs
	t := Cputime();
	nvals,best,graphs := fullextend(N,graphs,d,level);

	if (nvals eq 0) then
		print "No more values!";
		break;
	end if;

	print Cputime(t)," seconds.";
	print N," vertices: ",[best[j] : j in [1..nvals]];	
	
	N +:= 1;

end while;