(* :Title: Transfer matrices *)

(* :Context: "GrafPack`TransferMatrices`" *)

(* :Authors:
	D. Andrén, P.H. Lundow (editor), K. Markström
	Department of Mathematics
	Umea University
	S-901 87 Umea
	Sweden
	Bug reports to per.hakan.lundow@math.umu.se
*)
(* :Summary:
	This package contains various routines that returns (compressed)
	transfer matrices for polygraphs.
*)
(* :History: 
	010423 Created.
	020227 Corrected usage message of CompressedTransferMatrix1F.
	030428 Added TransferMatrixIP, CompressedTransferMatrixIP and
		CompressedTransferVectorIP.
*)
(* :Mathematica Version: 3.0 *)
(* :Keywords: *)
(* :Limitations: *)
(* :Discussion: *)

BeginPackage["GrafPack`TransferMatrices`",
	{
		"GrafPack`Combinatorics`",
		"GrafPack`Basics`",
		"GrafPack`GraphPolynomials`"
	}
]

CompressedTransferMatrix1F::usage =
	"CompressedTransferMatrix1F[g, orb] returns the compressed transfer
	matrix for counting 1-factors in G x Pn. The argument orb is the list
	of non-isomorphic 2-colourings returned by Orbits.\n
	Default options: Reduce->True, Verbose->False. \n
	Reduce -> False performs the compression but skips the reduction-part.\n
	Verbose -> True shows the progress during the computations.\n
	Example: \n
	\t This computes the number of 1-factors in C4 x P8 \n
	\t g=Cycle[4];\n
	\t aut=Automorphisms[g];\n
	\t orb=Orbits[aut,2];\n
	\t mat=CompressedTransferMatrix1F[g,orb];\n
	\t MatrixPower[mat,8,1,1]."

CompressedTransferMatrixM::usage =
	"CompressedTransferMatrixM[g, orb] returns the compressed transfer
	matrix for counting matchings in G x Pn. The argument orb is the list
	of non-isomorphic 2-colourings returned by Orbits.  The option
	Verbose->True shows the progress during the computations. See
	CompressedTransferMatrix1F."

CompressedTransferMatrixMP::usage =
	"CompressedTransferMatrixMP[g, orb, x] returns the compressed
	transfer matrix for computing matching polynomials of G x P_n. Argument
	orb is the list of non-isomorphic 2-colourings returned by Orbits. The
	option Verbose->True shows the progress during the computations. See
	CompressedTransferMatrix1F."

CompressedTransferMatrixIP::usage =
	"CompressedTransferMatrixIP[g, orb, x, y] returns the compressed
	transfer matrix for computing the Ising polynomial of G x Pn in
	the variables x and y. Argument orb is the list of non-isomorphic
	2-colourings returned by Orbits. The option Verbose->True shows the
	progress during the computations.
	\nExample: compute polynomial for C4 x Pn
	\n g=Cycle[4];
	\n aut=Automorphisms[g];
	\n orb=Orbits[aut,2];
	\n vec=CompressedTransferVectorIP[g,orb,x,y];
	\n mat=CompressedTransferMatrixIP[g,orb,x,y];
	\n Do[vec=Expand[vec.mat],{n-1}];
	\n pol=Apply[Plus,vec];"

CompressedTransferVectorIP::usage =
	"CompressedTransferVectorIP[g, orb, x, y] returns the compressed
	transfer vector used in the first step for computing the Ising
	polynomial of G x Pn in the varibles x and y. Argument orb is the
	list of non-isomorphic 2-colourings returned by Orbits. The option
	Verbose->True shows the progress during the computations. See
	CompressedTransferMatrixIP."

(* Reduce::usage *)
If[!StringMatchQ[Reduce::usage, "*GrafPack*"],
	Reduce::usage = Reduce::usage <>
	"\n Reduce is also an option for the GrafPack function
	CompressedTransferMatrix1F."
]

TransferMatrix1F::usage =
	"TransferMatrix1F[g, a, b] returns a transfer matrix for counting
	1-factors in a polygraph. Arguments a and b are the binary relations
	of the polygraph at graph g. The option Verbose->True shows the
	progress during the computations.\n
	Example: \n
	\t This computes the number of 1-factors in C4 x C8:\n
	\t g=Cycle[4];\n
	\t a={{1,1},{2,2},{3,3},{4,4}};\n
	\t mat=TransferMatrix1F[g,a,a];\n
	\t Sum[MatrixPower[mat,8,i,i],{i,Length[mat]}]"

TransferMatrixM::usage =
	"TransferMatrixM[g, a, b] returns the transfer matrix for counting
	matchings in a polygraph. The option Verbose->True shows the progress
	during the computations. See TransferMatrix1F."

TransferMatrixMP::usage =
	"TransferMatrixMP[g, a, b, x] returns the transfer matrix for computing
	matching polynomials of a polygraph. Argument x is the variable that
	the polynomials are expressed in. The option Verbose->True shows the
	progress during the computations. See TransferMatrix1F."

TransferMatrixIP::usage =
	"TransferMatrixIP[g, in, out, pos, neg, x, y] returns the general
	transfer matrix for computing the Ising polynomial in variables x
	and y.	Argument in is the set of in-vertices, out is the set of
	out-vertices, pos is the set of vertices fixed to +1, neg is the set
	of vertices fixed to -1, in, pos and neg are 'grey' vertices and do
	not contribute to the magnetization, incident edges contribute to the
	energy though. The union of in, out, pos and neg must be the set of
	vertices of g, pos and neg must be disjoint (obviously), in and out
	must be disjoint from pos and neg, in and out need not be disjoint."

TransferRow1F::usage =
	"TransferRow1F[g, a, b, r, c] returns a row of the transfer matrix
	for counting 1-factors in a polygraph. Arguments a and b are the binary
	relations of the polygraph at graph g. The row returned corresponds to
	the r:th subset of a as in-edges to the graph g. The list c contains the
	subset-numbers of b that should be out-edges.\n
	Example: \n
	\t a={{1,1},{2,2},{3,3},{4,4}};\n
	\t TransferRow1F[Path[4],a,a,5,Table[i,{i,0,15}]]\n
	returns row number 5 (counting from 0) of the transfer matrix for 
	counting 1-factors in P4 x Pn."

TransferRowM::usage=
	"TransferRowM[g, a, b, r, c] returns a row of the transfer matrix for
	counting matchings in a polygraph. See TransferRow1F."

TransferRowMP::usage =
	"TransferRowMP[g, a, b, x, r, c] returns a row of the transfer matrix
	for computing matching polynomials of a polygraph. See TransferRow1F."

Begin["`Private`"]


protected = Unprotect[Reduce]

definitions = {
	CompressedTransferMatrix1F, CompressedTransferMatrixM, 
	CompressedTransferMatrixMP, TransferMatrix1F, TransferMatrixM, 
	TransferMatrixMP, TransferRow1F, TransferRowM, TransferRowMP
}

Scan[Unprotect, definitions]


CompressedTransferMatrix1F[g_Graph, orb_List, opts___?OptionQ] :=
	Module[{a, b, i, orb2, orb3, reduce, verbose, k=0},
    	reduce = Reduce /. {opts} /. Options[CompressedTransferMatrix1F];
		verbose = Verbose /. {opts} /. Options[CompressedTransferMatrix1F];
		a = Table[{i, i}, {i, Order[g]}];
	    Which[

			(*Case 1: do not reduce *)
			!reduce,

			If[verbose, Print[Length[orb]," x ", Length[orb], " matrix"]];
			Map[
				(k++;
				If[verbose, Print["Making row ", k]];
				i = #;
				Map[Apply[Plus, TransferRow1F[g, a, a, i, #]]&, orb])&,
				Map[First, orb]
			],

			(*Case 2: balanced bipartite graph *)
			(b=Bipartition[g])!={} && Length[b[[1]]]==Length[b[[2]]],

			orb2 = balancedClasses[b[[1]], b[[2]], orb];
			If[verbose, Print[Length[orb2]," x ", Length[orb2], " matrix"]];
			Map[
				(k++;
				If[verbose, Print["Making row ", k]];
				i = #;
				Map[Apply[Plus, TransferRow1F[g, a, a, i, #]]&, orb2])&,
				Map[First, orb2]
			],

			(*Case 3: even graph*)
			EvenQ[Order[g]],

			orb2 = evenClasses[orb];
			If[verbose, Print[Length[orb2]," x ", Length[orb2], " matrix"]];
			Map[
				(k++;
				If[verbose, Print["Making row ", k]];
				i = #;
				Map[Apply[Plus, TransferRow1F[g, a, a, i, #]]&, orb2])&,
				Map[First, orb2]
			],

			(*Case default: odd graph *)
			OddQ[Order[g]],

			orb2 = evenClasses[orb];
			orb3 = oddClasses[orb];
			If[verbose, Print[Length[orb2]," x ", Length[orb2], " matrix"]];
			Dot[
				If[verbose, Print["First matrix: ", Length[orb2], " rows"]];
			    Map[
					(k++;
					If[verbose, Print["Making row ", k]];
					i = #;
					Map[Apply[Plus, TransferRow1F[g, a, a, i, #]]&, orb3])&,
					Map[First, orb2]
			    ],
				k=0;
				If[verbose, Print["Second matrix: ", Length[orb3], " rows"]];
			    Map[
					(k++;
					If[verbose, Print["Making row ", k]];
					i = #;
					Map[Apply[Plus, TransferRow1F[g, a, a, i, #]]&, orb2])&,
					Map[First, orb3]
			    ]
			]
	    ]
	]

Options[CompressedTransferMatrix1F] = {Reduce -> True, Verbose -> False}


CompressedTransferMatrixM[g_Graph, orb_List, opts___?OptionQ] :=
	Module[{a, i, verbose, k=0},
		verbose = Verbose /. {opts} /. Options[CompressedTransferMatrixM];
		If[verbose, Print[Length[orb], " x ", Length[orb], " matrix."]];
		a = Table[{i, i}, {i, Order[g]}];
		Map[
			(k++;
			If[verbose, Print["Making row ", k]];
			i = #;
			Map[Apply[Plus, TransferRowM[g, a, a, i, #]]&, orb])&,
			Map[First, orb]
		]
	]

Options[CompressedTransferMatrixM] = {Verbose -> False}


CompressedTransferMatrixMP[g_Graph, orb_List, x_, opts___?OptionQ] :=
	Module[{a, i, verbose, k=0},
		verbose = Verbose /. {opts} /. Options[CompressedTransferMatrixMP];
		If[verbose, Print[Length[orb], " x ", Length[orb], " matrix."]];
		a = Table[{i, i}, {i, Order[g]}];
    	Map[
			(k++;
			If[verbose, Print["Making row ", k]];
			i = #;
			Map[Apply[Plus, TransferRowMP[g, a, a, x, i, #]]&, orb])&,
			Map[First, orb]
	    ]
	]

Options[CompressedTransferMatrixMP] = {Verbose -> False}


CompressedTransferMatrixIP[g_Graph, orb_List, x_, y_, opts___?OptionQ] :=
	Module[{adj, a, b, i, verbose, k=0},
		verbose=Verbose/.{opts}/.Options[CompressedTransferMatrixIP];
		If[verbose,Print[Length[orb]," x ",Length[orb]," matrix."]];
		adj = ToAdjacencyLists[transferGraph[g]];
		b = Range[Order[g]]; (* out-vertices *)
		a = b + Order[g];    (*  in-vertices *)
		Map[(
			k++;
			If[verbose, Print["Making row ", k]];
			i = #;
			Map[
				Apply[Plus,transferRowIP[adj,x,y,a,b,i,#]]&,
				orb
			])&,
			Map[First, orb]
		]
	]

Options[CompressedTransferMatrixIP] = {Verbose -> False}


CompressedTransferVectorIP[g_Graph, orb_List, x_, y_, opts___?OptionQ] :=
	Module[{adj, h, i, u, v=Vertices[g]},
		adj = ToAdjacencyLists[g];
		Map[(
			u = UnrankSet[First[#], v];
			Length[#]*IsingMonomial[adj, u, x, y])&,
			orb
		]
	]


TransferMatrix1F[g_Graph, a:{{_,_}...}, b:{{_,_}...}, opts___?OptionQ] :=
	Module[{rel1, rel2, verbose, k=0},
		verbose = Verbose /. {opts} /. Options[TransferMatrix1F];
		rel1 = goodRelations[a];
		rel2 = goodRelations[b];
		If[verbose, Print[Length[rel1]," x ",Length[rel2]," matrix"]];
		Map[
			(k++;
			If[verbose, Print["Making row ", k]];
			TransferRow1F[g, a, b, #, rel2])&, 
			rel1
		]
	]

Options[TransferMatrix1F] = {Verbose -> False}


TransferMatrixM[g_Graph, a:{{_,_}...}, b:{{_,_}...}, opts___?OptionQ] :=
	Module[{rel1, rel2, verbose, k=0},
		verbose = Verbose /. {opts} /. Options[TransferMatrix1F];
		rel1 = goodRelations[a];
		rel2 = goodRelations[b];
		If[verbose, Print[Length[rel1]," x ",Length[rel2]," matrix"]];
		Map[
			(k++;
			If[verbose, Print["Making row ", k]];
			TransferRowM[g, a, b, #, rel2])&, 
			rel1
		]
	]

Options[TransferMatrixM] = {Verbose -> False}


TransferMatrixMP[g_Graph, a:{{_,_}...}, b:{{_,_}...}, x_, opts___?OptionQ]:=
	Module[{rel1, rel2, verbose, k=0},
		verbose = Verbose /. {opts} /. Options[TransferMatrixMP];
		rel1 = goodRelations[a];
		rel2 = goodRelations[b];
		If[verbose, Print[Length[rel1]," x ",Length[rel2]," matrix"]];
		Map[
			(k++;
			If[verbose, Print["Making row ", k]];
			TransferRowMP[g, a, b, x, #, rel2])&, 
			rel1
		]
	]

Options[TransferMatrixMP] = {Verbose -> False}


TransferMatrixIP[g_Graph, in_List, out_List, pos_List, neg_List, x_, y_, 
	opts___?OptionQ] :=
	Module[{adj, i, j, k, p, q, u, v, verbose},
		verbose = Verbose /. {opts} /. Options[TransferMatrixIP];
		adj = ToAdjacencyLists[g];
		If[verbose,Print[2^Length[in]," x ",2^Length[out]," matrix"]];
		Table[
			If[verbose,Print["Making row ", i]];
			u = UnrankSet[i, in];
			v = Complement[in, u];
			u = Union[u, pos];
			v = Union[v, neg];
			k = Length[v] - Length[u];
			Table[
				p = UnrankSet[j, out];
				q = Complement[out, p];
				p = Union[u, p]; (* positive vertices *)
				q = Union[v, q]; (* negative vertices *)
				If[DisjointQ[p, q],
					IsingMonomial[adj, p, x, y] * y^k,
					0
				],
				{j, 0, 2^Length[out]-1}
			],
			{i, 0, 2^Length[in]-1}
		]
	] /; okArgsQ[g, in, out, pos, neg]

okArgsQ[g_, in_, out_, pos_, neg_] :=
	And[
		DisjointQ[in,  pos, neg],
		DisjointQ[out, pos, neg],
		Union[in, out, pos, neg] == Vertices[g]
	]


TransferRow1F[g_Graph, a:{{_,_}...}, b:{{_,_}...}, r_Integer, c_List] :=
	Module[{x, y},
		x = Map[Last, UnrankSet[r, a]];
		Map[
			(y = Map[First, UnrankSet[#, b]];
			If[Intersection[x, y] == {},
				NumberOfOneFactors[DeleteVertex[g, Union[x, y]]],
			(*Else*)
				0
			])&,
			c
		]
	]


TransferRowM[g_Graph, a:{{_,_}...}, b:{{_,_}...}, r_Integer, c_List] :=
	Module[{x, y},
		x = Map[Last, UnrankSet[r, a]];
		Map[
			(y = Map[First, UnrankSet[#, b]];
			If[Intersection[x,y] == {},
				NumberOfMatchings[DeleteVertex[g, Union[x, y]]],
			(*Else*)
				0
			])&,
			c
		]
	]


TransferRowMP[g_Graph, a:{{_,_}...}, b:{{_,_}...}, z_, r_Integer, c_List] :=
	Module[{s, x, y},
		x = Map[Last, UnrankSet[r, a]];
		Map[
			(y = Map[First, UnrankSet[#, b]];
			s = (-1)^Length[y];
			If[Intersection[x, y] == {},
				s*MatchingPolynomial[DeleteVertex[g, Union[x, y]], z],
			(*Else*)
				0
			])&,
			c
		]
	]


(* Private functions *)

goodRelations[rel_] := (* private function *)
	(* Returns the subsets of rel that are matchings *)
	Flatten[Position[Map[goodRelQ, Subsets[rel]], True]] - 1

goodRelQ[r_] := (* private function *)
	(* Tests whether r is a 'good' relation, ie a matching *)
	Module[{r1=Map[First,r], r2=Map[Last,r]},
		Length[Union[r1]] == Length[r1] && Length[Union[r2]] == Length[r2]
	]

balancedClasses[x_, y_, orb_] :=
	Module[{z},
		Select[
			orb,
			(z=UnrankSet[First[#]];
			Length[Intersection[x, z]] == Length[Intersection[y, z]])&
		]
	]

evenClasses[orb_] := Select[orb, EvenQ[Length[UnrankSet[First[#]]]]&]

oddClasses[orb_] := Select[orb, OddQ[Length[UnrankSet[First[#]]]]&]

transferGraph[g_] :=
	Module[{n=Order[g]},
		AddEdge[AddVertex[g, n], Array[{#, n+#}&, n]]
	]

transferRowIP[adj_, x_, y_, in_List, out_List, row_Integer, col_List] :=
	Module[{a, b, c, k},
		a = UnrankSet[row, in];        (* Positive in *)
		Map[(
			b = UnrankSet[#, out]; (* Positive out *)
			c = Union[a, b];       (* Positive all *)
			k = Length[in] - 2*Length[a];
			IsingMonomial[adj, c, x, y]*y^k)&,
			col
		]
	]

Protect[Evaluate[protected]]

Scan[Protect, definitions]


End[]

EndPackage[]