(* :Title: Graph Polynomials *)

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

(* :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 functions for computing graph polynomials.
*)
(* :History:
	991212 Changed expandexpression.
	001027 Added EstimateNumberOfOneFactors, OneFactorPQ
	010424 Minor bug fixed in CharactersticPolynomial.
	030428 Changed the arguments of IsingMonomial.
*)
(* :Mathematica Version: 3.0 *)
(* :Keywords: *)
(* :Limitations: *)
(* :Discussion: *)

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

CharacteristicPolynomial::usage =
	"CharacteristicPolynomial[mat, x] returns the characteristic polynomial
	of the square matrix mat in variable x. This is Det[x*I-mat] if mat is
	an integer matrix and Apply[Times,x-Eigenvalues[mat]] otherwise. The
	coefficient of x^n (n is the order of the matrix) is always +1.
	CharacteristicPolynomial[g, x] gives the characteristic polynomial of
	the graph g, use option WorkingPrecision if you want a numerical solution.
	Example, for fast computation:
	CharacteristicPolynomial[g, x, WorkingPrecision->$MachinePrecision]"

EstimateNumberOfOneFactors::usage =
	"EstimateNumberOfOneFactors[g, k] returns an estimate of the number of 
	1-factors in graph g by running a Monte Carlo algorithm k times. For
	bipartite graphs EstimatePermanent is used."

IndependencePolynomial::usage =
	"IndependencePolynomial[g, x] returns the independence polynomial of
	graph g in variable x. The coefficient of x^k is the number of
	independent sets on k vertices."

IsingMonomial::usage =
	"IsingMonomial[g, pos] returns the pair {energy-level, magnetization}
	when the vertices in graph g given in list pos are set to +1 and the
	other vertices to -1. IsingMonomial[g, pos, x, y] returns the monomial
	x^energy-level * y^magnetization. Alternatively we may also give g as a
	list of adjacencies."

IsingPolynomial::usage =
	"IsingPolynomial[g, x, y] returns the Ising polynomial of the graph g
	in variables x and y. The function takes several options. The options
	Positive and Negative are the vertices to fix as +1 and -1 respectively
	and are by default set to {}. IndependentSetQ->True (the default) finds
	an independent set to speed up the calculations whereas the setting False
	does not. We may also specify an independent set of our own as in e.g.
	IndependentSetQ->{1,2,3,4}. Verbose->True prints a progress report
	during the calculations but is by default set to False. Example:
	IsingPolynomial[g, x, y, Positive->{1,2}, Negative->{3,4}, Verbose->True,
	IndependentSetQ->{1,3,5,7,9,11}] fixes vertices {1,2} to +1 and vertices
	{3,4} to -1 and takes advantage of that {1,3,5,7,9,11} is an independent
	set in g. Progress is regularly reported."

MatchingPolynomial::usage =
	"MatchingPolynomial[g, x] returns the matching polynomial of graph g
	in variable x."

NumberOfIndependentSets::usage =
	"NumberOfIndependentSets[g] returns the number of independent sets in
	graph g."

NumberOfMatchings::usage =
	"NumberOfMatchings[g] returns the number of matchings in graph g."

NumberOfOneFactors::usage =
	"NumberOfOneFactors[g] returns the number of 1-factors (perfect matchings)
	in graph g."

OneFactorPQ::usage = 
	"OneFactorPQ[g, p] does a probabilistic test of whether the graph g has
	a 1-factor.  The argument p is a small probability.  If the function
	returns True then g has a 1-factor.  Otherwise the probability that g
	has a 1-factor is at most p."


Begin["`Private`"]

protected = Unprotect[CharacteristicPolynomial]

definitions = {
	EstimateNumberOfOneFactors, IndependencePolynomial, IsingMonomial,
	IsingPolynomial, MatchingPolynomial, NumberOfIndependentSets,
	NumberOfMatchings, NumberOfOneFactors, OneFactorPQ
}

Scan[Unprotect, definitions]


Clear[CharacteristicPolynomial]

CharacteristicPolynomial[{}, x_] := 1

CharacteristicPolynomial[mat_, x_] :=
	charpoly1[mat, x] /; MatrixQ[mat, IntegerQ]

CharacteristicPolynomial[mat_, x_] :=
	charpoly2[mat, x] /; MatrixQ[mat, NumberQ]

CharacteristicPolynomial[mat_, x_] :=
	charpoly1[mat, x] /; MatrixQ[mat]

CharacteristicPolynomial[Graph[adj_, _], x_, opts___?OptionQ] :=
	Module[{prec},
		prec = WorkingPrecision/.{opts}/.Options[CharacteristicPolynomial];
		Switch[prec,
			Infinity,
				CharacteristicPolynomial[adj, x],
			_Integer,
				CharacteristicPolynomial[N[adj, prec], x],
			_,
				Return[]
		]
	]

Options[CharacteristicPolynomial] = {WorkingPrecision -> Infinity}

charpoly1[mat_, x_] := (* private function *)
	Module[{p, n=Length[mat]},
		p = Det[x * IdentityMatrix[n] - mat];
		Expand[p * Sign[Coefficient[p, x, n]]]
	]

charpoly2[mat_, x_] := (* private function *)
	Expand[Apply[Times, x - Eigenvalues[mat]]]


EstimateNumberOfOneFactors[g_Graph, k_Integer] := 0 /; OddQ[Order[g]]

EstimateNumberOfOneFactors[g_Graph?BipartiteQ, k_Integer] :=
	EstimatePermanent[ToBiadjacencyMatrix[g], k]

EstimateNumberOfOneFactors[g_Graph, k_Integer] :=
	(* REF: Lovasz & Plummer, Matching theory, ch. 8 *)
	Module[{mat0, mat1, pos=ToEdges[g], sum=0},
		mat0 = Table[0.0, {Order[g]}, {Order[g]}];
		Do[
			mat1 = MapAt[(2.0*Random[Integer]-1.0)&, mat0, pos];
			sum  = sum + Det[mat1 - Transpose[mat1]],
			{k}
		];
		N[sum / k]
	]


IndependencePolynomial[g_Graph, x_] :=
	Module[{p, z},
		p = indpoly[g, z];
		p /. z -> x
	]

indpoly[g_, x_] := (* private function *)
	Module[{maxdeg, p, q, v, deg=DegreeSequence[g]},
		maxdeg = Max[deg];
		Which[
			maxdeg <= 1,
				indpolydeg1[g, x],
			maxdeg <= 2,
				indpolydeg2[g, x, deg],
			True,
				v = Position[deg, maxdeg, {1}, 1][[1, 1]];
				p = indpoly[DeleteVertex[g, v], x];
				q = indpoly[DeleteVertex[g, Neighbours[g, v, 1]], x];
				p + Expand[q * x]
		]
	]

indpolypath[n_,	x_]:= (* private function *)
	Module[{k},
		Sum[Binomial[n + 1 - k, k] * x^k, {k, 0, Ceiling[n / 2]}]
	]

indpolycycle[n_, x_] := (* private function *)
	Module[{k},
		Sum[n * Binomial[n - k, k] / (n - k) * x^k, {k, 0, Floor[n / 2]}]
	]

indpolydeg1[g_, x_] := (* private function *)
	Expand[(1 + 2*x)^Size[g] * (1 + x)^(Order[g] - 2*Size[g])]

indpolydeg2[g_, x_, deg_] := (* private function *)
	Expand[
		Apply[Times,
			Map[
				If[Min[deg[[#]]] < 2,
					indpolypath[Length[#], x],
				(*Else*)
					indpolycycle[Length[#], x]
				]&,
				Component[g, All]
			]
		]
	]


IsingMonomial[g_Graph, pos_List] :=
	IsingMonomial[ToAdjacencyLists[g], pos]

IsingMonomial[g_Graph, pos_List, x_, y_] :=
	Apply[Times, {x, y}^IsingMonomial[ToAdjacencyLists[g], pos]]

IsingMonomial[adj:{___List}, pos_List] :=
	Module[{spin},
		spin = MapAt[1&, Table[-1, {Length[adj]}], Map[List, pos]];
		{
			Apply[Plus, spin * Map[Apply[Plus, spin[[#]]]&, adj]] / 2,
			Apply[Plus, spin]
		}
	]

IsingMonomial[adj:{___List}, pos_List, x_, y_] :=
	Apply[Times, {x, y}^IsingMonomial[adj, pos]]


IsingPolynomial[g_Graph, x_, y_, opts___?OptionQ] :=
	Expand[(y + 1/y)^Order[g]] /; Size[g] == 0

IsingPolynomial[g_Graph, x_, y_, opts___?OptionQ] :=
	Module[{ind, neg, pos},
		pos = Positive /. {opts} /. Options[IsingPolynomial];
		neg = Negative /. {opts} /. Options[IsingPolynomial];
		If[!SubsetQ[pos, Vertices[g]],
			Message[IsingPolynomial::notsubset, pos];
			Return[]
		];
		If[!SubsetQ[neg, Vertices[g]],
			Message[IsingPolynomial::notsubset, neg];
			Return[]
		];
		If[!DisjointQ[pos, neg],
			Message[IsingPolynomial::notdisjoint, pos, neg];
			Return[]
		];
		ind = IndependentSetQ /. {opts} /. Options[IsingPolynomial];
		Block[{verbose},
			verbose = Verbose /. {opts} /. Options[IsingPolynomial];
			Switch[ind,
				True,
					isingpoly1[g, x, y, Union[pos], Union[neg]],
				False,
					isingpoly2[g, x, y, Union[pos], Union[neg]],
				_List,
					If[!IndependentSetQ[g, ind],
						Message[IsingPolynomial::notindset, ind];
						Return[]
					];
					isingpoly3[g, x, y, Union[pos], Union[neg], Union[ind]],
				_,
					Return[]
			]
		]
	]

Options[IsingPolynomial] = {
	Positive->{}, Negative->{}, IndependentSetQ->True, Verbose->False
}

IsingPolynomial::notindset = "`1` is not an independent set"

IsingPolynomial::notsubset = "`1` is not a subset of the vertices."

IsingPolynomial::notdisjoint = "`1` and `2` are not disjoint sets."

isingpoly1[g_, x_, y_, {}, {}] :=
	Module[{ind, poly, v},
		If[verbose, Print["Computing maximum independent set..."]];
		ind = MaximumIndependentSet[g];
		v = First[Complement[Vertices[g], ind]];
		poly = isingpoly3[g, x, y, {v}, {}, ind];
		poly + (poly /. y -> 1/y)
	]

isingpoly1[g_, x_, y_, pos_, neg_] :=
	Module[{a, b, h, ind, p},
	(* 
		Relabels the graph so that the flip vertices comes first,
		then the positive vertices, then the negative vertices and
		finally the independent set which the function finds.
	*)
		p = Join[Complement[Vertices[g], pos, neg], pos, neg];
		h = PermuteGraph[g, p];
		p = InversePermutation[p];
		a = PermuteSet[pos, p];
		b = PermuteSet[neg, p];
		If[verbose, Print["Computing maximum independent set..."]];
		ind = MaximumIndependentSet[
			InduceSubgraph[h, Range[Order[g] - Length[a] - Length[b]]]
		];
		p = Join[Complement[Vertices[g], a, b, ind], a, b, ind];
		h = PermuteGraph[h, p];
		isingpoly5[h, x, y, Length[a], Length[b], Length[ind]]
	]

isingpoly2[g_, x_, y_, {}, {}] :=
	Module[{poly},
		poly = isingpoly2[g, x, y, {1}, {}];
		poly + (poly /. y -> 1/y)
	]

isingpoly2[g_, x_, y_, pos_, neg_] :=
	(*
		Relabels the graph so that the flip vertices comes first,
		then the positive vertices and finally the negative
		vertices. Initializes the variables used in the iterations.
		No independent set is used.
	*)
	Block[{adj, lev, mag, p, spin},
		p = Join[Complement[Vertices[g], pos, neg], pos, neg];
		adj = ToAdjacencyLists[PermuteGraph[g, p]];
		spin = Table[1, {Order[g]}];
		Scan[(spin[[-#]] = -1)&, Range[Length[neg]]];
		mag = Apply[Plus, spin];
		lev = Apply[Plus, spin * Map[Apply[Plus, spin[[#]]]&, adj]] / 2;
		isingpoly4[Order[g] - Length[pos] - Length[neg], x, y]
	]

isingpoly3[g_, x_, y_, {}, {}, ind_] :=
	Module[{poly, v},
		v = First[Complement[Vertices[g], ind]];
		poly = isingpoly3[g, x, y, {v}, {}, ind];
		poly + (poly /. y -> 1/y)
	]

isingpoly3[g_, x_, y_, pos_, neg_, ind_] :=
	(*
		Relabels the graph so that the flip vertices come first,
		then the positive vertices, then the negative vertices
		and finally the independent set.
	*)
	Module[{h, p, set},
		set = Complement[ind, pos, neg];
		p = Join[Complement[Vertices[g], pos, neg, set], pos, neg, set];
		h = PermuteGraph[g, p];
		isingpoly5[h, x, y, Length[pos], Length[neg], Length[set]]
	]

isingpoly4[0, x_, y_] := x^lev * y^mag

isingpoly4[n_, x_, y_] :=
	Module[{code, card, diff, flip, cnt=0, poly=0},
		NextGrayCodeChunk[code, card, flip, n, 12];
		If[verbose, Print["Flipping spins..."]];
		While[True, 
			poly = poly + Apply[Plus,
				Map[
					(
						spin[[#]] = -spin[[#]];
						diff = 2*spin[[#]];
						lev = lev + diff*Apply[Plus, spin[[adj[[#]]]]];
						mag = mag + diff;
						x^lev * y^mag
					)&,
					flip
				]
			];
			If[verbose,
				cnt = cnt + Length[flip];
				Print[N[100 * cnt / 2^n, 3], " %"]
			];
			If[card == 0, Break[]];
			NextGrayCodeChunk[code, card, flip];
		];
		If[verbose, Print["Done."]];
		poly
	]

isingpoly5[g_, x_, y_, a_, b_, c_] :=
	(*
		Initializes the variables used in the iterations
	*)
	Block[{adj, deg, lev, mag, neib, spin, n=Order[g]},
		neib = Map[Neighbours[g, #]&, Range[n + 1 - c, n]];
		spin = Table[1, {n - c}];
		Scan[(spin[[-#]] = -1)&, Range[b]];
		deg = Map[Apply[Plus, spin[[#]]]&, neib];
		neib = Map[Map[First, Position[neib, #]]&, Range[n - c]];
		neib = Map[List, neib, {2}];
		adj = ToAdjacencyLists[InduceSubgraph[g, Range[n - c]]];
		lev = Apply[Plus, spin * Map[Apply[Plus, spin[[#]]]&, adj]] / 2;
		mag = Apply[Plus, spin];
		isingpoly6[n - a - b - c, x, y]
	]

isingpoly6[0, x_, y_] :=
	Module[{poly},
		poly = x^lev * y^mag * Apply[Times, Map[z, deg]];
		poly = (poly /. z[k_] :> x^k * y + 1 / (x^k * y));
		Expand[poly]
	]

isingpoly6[n_, x_, y_] :=
	Module[{code, card, diff, flip, z, cnt=0, poly=0},
		NextGrayCodeChunk[code, card, flip, n, 12];
		If[verbose, Print["Flipping spins..."]];
		While[True,
			poly = poly + Apply[Plus,
				Map[
					(
						spin[[#]] = -spin[[#]];
						diff = 2*spin[[#]];  (* spin difference *)
						lev = lev + diff*Apply[Plus, spin[[adj[[#]]]]];
						mag = mag + diff;
						deg = MapAt[(# + diff)&, deg, neib[[#]]];
						x^lev * y^mag * Apply[Times, Map[z, deg]]
					)&,
					flip
				]
			];
			If[verbose,
				cnt = cnt + Length[flip];
				Print[N[100 * cnt / 2^n, 3], " %"]
			];
			If[card == 0, Break[]];
			NextGrayCodeChunk[code, card, flip];
		];
		poly = (poly /. z[k_] :> x^k * y + 1 / (x^k * y));
		expandexpression[poly, 100]
	]

expandexpression[expr_Plus, k_Integer?Positive] :=
	Module[{i, tmp, len=Length[expr], new=0, old=expr},
		If[verbose, Print["Expanding ", Length[expr], " terms..."]];
		Do[
			If[verbose, Print[N[100 * i / len, 3], " %"]];
			tmp = Take[old, Min[k, Length[old]]];
			old = old - tmp;
			new = new + Expand[tmp],
			{i, 0, len-1, k}
		];
		If[verbose, Print["Done."]];
		new
	]

expandexpression[expr_, k_Integer?Positive] := Expand[expr]


MatchingPolynomial[g_Graph, x_] :=
	Module[{y},
		If[ 2*Size[g] > Binomial[Order[g], 2],
			matchpolydense[g, y] /. y -> x,
		(*Else*)
			matchpolysparse[g, y] /. y -> x
		]
	]

matchpolydense[g_, x_] := (* private function *)
	Module[{k, p, n=Order[g]},
		p = matchpolysparse[GraphComplement[g], x];
		Expand[
			Dot[
				Abs[Reverse[CoefficientList[p, x]][[Range[1, n + 1, 2]]]],
				Table[matchpolycomplete[n - 2*k, x],{k, 0, Floor[n/2]}]
			]
		]
	]

matchpolysparse[g_, x_] := (* private function *)
	Module[{e, maxdeg, deg=DegreeSequence[g]},
		maxdeg = Max[deg];
		Which[
			maxdeg <= 1,
				matchpolydeg1[g, x],
			maxdeg <= 2,
				matchpolydeg2[g, x, deg],
			True,
				e = HighDegreeEdge[g, deg];
				matchpolysparse[DeleteEdge[g, e], x] -
				matchpolysparse[DeleteVertex[g, e], x]
		]
	]

matchpolydeg1[g_, x_] := (* private function *)
	Expand[x^Order[g] * (1 - 1/x^2)^Size[g]]

matchpolydeg2[g_, x_, deg_] :=
	Expand[
		Apply[Times,
			Map[
				If[Min[deg[[#]]] < 2,
					matchpolypath[Length[#], x],
				(*Else*)
					matchpolycycle[Length[#], x]
				]&,
				Component[g, All]
			]
		]
	]

matchpolycomplete[n_, x_] := (* private function *)
	Module[{k},
		Expand[
			Sum[
				(-1)^k * x^(n - 2*k) / ((2*k)!! * (n - 2*k)!),
				{k, 0, Floor[n/2]}
			] * n!
		]
	]

matchpolypath[n_, x_] := (* private function *)
	Module[{k},
		Sum[
			(-1)^k * Binomial[n - k, k] * x^(n - 2*k),
			{k, 0, Floor[n/2]}
		]
	]

matchpolycycle[n_, x_] := (* private function *)
	Module[{k},
		Expand[
			Sum[
				(-1)^k * Binomial[n - k, k] / (n - k) * x^(n - 2*k),
				{k, 0, Floor[n/2]}
			] * n
		]
	]


NumberOfIndependentSets[g_Graph] := IndependencePolynomial[g, 1]


NumberOfMatchings[g_Graph] := Abs[MatchingPolynomial[g, I]]


NumberOfOneFactors[g_Graph] :=
	Module[{bip},
		Which[
			OddQ[Order[g]],
				0,
			MaxDegree[g] <= 2,
				Abs[MatchingPolynomial[g, 0]],
			(bip = Bipartition[g]) != {},
				Permanent[ToAdjacencyMatrix[g][[bip[[1]], bip[[2]]]]],
			True,
				Abs[MatchingPolynomial[g, 0]]
		]
	]


OneFactorPQ[g_Graph, prob_] :=
	(* REF: Lovasz & Plummer, Matching theory, ch. 8 *)
	Module[{mat0, mat1, n=Order[g], pos=ToEdges[g], test=False},
		mat0 = Table[0, {n}, {n}];
		Do[
			mat1 = MapAt[Random[Integer, {-n, n}]&, mat0, pos];
			If[Det[mat1 - Transpose[mat1]] != 0,
				test = True;
				Break[];
			],
			{1 + Ceiling[-Log[2, N[prob]]]}
		];
		test
	]

OneFactorPQ[Graph[{}, ___], ___] := True


Protect[Evaluate[protected]]

Scan[Protect, definitions]

End[]

EndPackage[]