(* :Title: Boltzmann *)

(* :Context: "Boltzmann`" *)

(* :Author: P.H. Lundow
	Bug reports to phl@kth.se
*)
(* :Summary:
	Functions for converting between the list of K-values and 
	the Boltzmann distribution of energies. Functions that return
	the traditional quantities of statistical physics.
*)
(* :History:
	031107 Created.
	031119 Much changed, added Susceptibility, deleted UC and UCM.
	040302 New normalisation of InternalEnergy and SpecificHeat.
	040512 Added Expectation, Moment, Entropy.
	040602 Expectation overloaded.
	040603 Expectation and Moment.
	041015 Clean-up. Improved listIntegrate. Deleted InternalEnergy, 
	       SpecificHeat, Magnetisation.
	041216 Added some comments.
	050404 listIntegrate0, BoltzmannDistribution0.
	071010 BoltzmannDistribution for a partition function.
	071015 NormalList01, NormalFunction01.
	100505 Added Cumulant.
	110902 Adapted to version 8: Renamed Cumulant, Expectation, Moment, 
	       Fractile.
*)
(* :Mathematica Version: 4.0 *)
(* :Keywords: *)
(* :Limitations: *)
(* :Discussion:	
*)

BeginPackage["Boltzmann`"]

BoltzmannDistribution::usage =
	"BoltzmannDistribution[K, listK, m] returns the probability
	for each energy in listK at coupling K for a graph with m
	edges. The argument listK is on the form {{e1,K1},..,{en,Kn}}
	where Ki is the coupling at relative energy ei. NOTE:
	e1<e2<...<en. The returned list is on the form
	{{e1,p1},...,{en,pn}} where pi is the probability for energy
	ei. Note that we must have equal distance between energies
	so that e2-e1=e3-e2=... If this is not the case, use 
	BoltzmannDistribution1 instead.\n
	BoltzmannDistribution[K, z, x, m] returns the probabilities
	for a partition function z in the formal variable x."

BoltzmannDistribution0::usage =
	"Does not assume equal distance between energies. Uses
	left-point rule for integration and adds an extra point.
	This point can also be provided as an optional argument."

BoltzmannDistribution1::usage =
	"Older and safer version of BoltzmannDistribution. Does not
	assume equal distance between energies. Uses trapezoidal rule
	for integration instead of Simpson's rule. Less accurate."

CouplingList::usage =
	"CouplingList[K, dist, m] returns the list of K-values at each
	energy for a graph with m edges given the coupling K and the
	distribution dist. This is the inverse of
	BoltzmannDistribution."

CumulantD::usage =
	"CumulantD[dist, k] returns the k:th cumulant in terms of central
	moments. CumulantD[dist, 1] is the mean."

EntropyD::usage =
	"EntropyD[dist] returns the entropy of distribution dist."

ExpectationD::usage =
	"ExpectationD[dist] returns the expected value of distribution
	dist. ExpectationD[dist, func] returns the expected value of
	the function func which takes the energies from dist as
	argument.  As a special case of function func
	ExpectationD[dist, k] returns the expected k:th power (moment)
	of the energy. See CentralMomentD."

FractileD::usage =
	"FractileD[dist, frac] returns the first x-value where the sum
	of probabilities is at least frac."

CentralMomentD::usage =
	"CentralMomentD[dist, k] returns the k:th central moment of
	distribution dist."

Normal01::usage =
	"Normal01[dist] returns an interpolating function
	corresponding to the distribution dist normalised to mean 0
	and variance 1."

NormalFunction01::usage =
	"Normal01[dist] returns an interpolating function
	corresponding to the distribution dist normalised to mean 0
	and variance 1. Same as Normal01."

NormalList01::usage =
	"NormalList01[dist] returns the distribution dist normalised to 
	mean 0 and variance 1."

	
Begin["`Private`"]

definitions = {BoltzmannDistribution, BoltzmannDistribution0,
	BoltzmannDistribution1, CouplingList, CumulantD, EntropyD, ExpectationD,
	FractileD, CentralMomentD, Normal01, NormalFunction01, NormalList01
}

Scan[Unprotect, definitions]


listDerivative[list_]:= (* private function *)
	Module[{x, y, z},
		{x, y} = Transpose[list];
		z = Drop[x, 1] - Drop[x, -1];
		y = Drop[y, 1] - Drop[y, -1];
		x = Drop[x, 1] + Drop[x, -1];
		Transpose[{x/2, y/z}]
	]


listIntegrate0[list_, False]:= (* Private function. Left-point rule. *)
	Module[{x, y, h},
		{x, y} = Transpose[list];
		h = Drop[x, 1] - Drop[x, -1];
		Transpose[{x, FoldList[Plus, 0, h*Drop[y,-1]]}]
	]

listIntegrate0[list_, True]:= (* Private function. Guess right limit point.*)
	Module[{x, y, d, h},
		{x, y} = Transpose[list];
		d = Drop[x, 1] - Drop[x, -1];
		h = Min[d];
		AppendTo[d, h];
		AppendTo[x, Last[x]+h];
		Transpose[{x, FoldList[Plus, 0, d*y]}]
	]


listIntegrate0[list_, r_]:= (* Private function. Give right limit point.*)
	Module[{x, y, d},
		{x, y} = Transpose[list];
		AppendTo[x, r];
		d = Drop[x, 1] - Drop[x, -1];
		Transpose[{x, FoldList[Plus, 0, d*y]}]
	]


listAverage[list_]:= (* private function *)
	Module[{x, y, z},
		{x, y} = Transpose[list];
		x = Drop[x, 1] + Drop[x, -1];
		y = Drop[y, 1] + Drop[y, -1];
		Transpose[{x/2, y/2}]
	]


(* OLD VERSION: trapezoidal rule *)
listIntegrate1[list_]:= (* private function *)
	Module[{x, y, xd, ys},
		{x, y} = Transpose[list];
		xd = Drop[x, 1] - Drop[x, -1];
		ys = Drop[y, 1] + Drop[y, -1];
		Transpose[{x, FoldList[Plus, 0, xd*ys/2]}]
	]


(* NEW VERSION: modified Simpson's rule *)

listIntegrate[list_]:= (* private function *)
	{{list[[1,1]],0}} /; Length[list] == 1

listIntegrate[list_]:= (* private function *)
	Module[{x, y, h, tmp},
		{x, y} = Transpose[list];
		h = Min[Drop[x, 1] - Drop[x, -1]];
		tmp = Flatten[
			Map[{
				(5*#[[1]] + 8*#[[2]] - #[[3]])/12,
				(-#[[1]] + 8*#[[2]] + 5*#[[3]])/12}&,
				Partition[y, 3, 2]
			]
		];
		If[EvenQ[Length[list]], 
			AppendTo[tmp, (y[[-2]]+y[[-1]])/2]
		];
		Transpose[{x, h*FoldList[Plus, 0, tmp]}]
	]


BoltzmannDistribution[K_, listK_List, m_]:=
	Module[{sum, list},
		list = Map[{#[[1]], m*(K - #[[2]])}&, listK];
		list = listIntegrate[list];
		list = Map[{#[[1]], Exp[#[[2]]]}&, list];
		sum  = Apply[Plus, Map[Last, list]];
		Map[{#[[1]], #[[2]]/sum}&, list]
	]

BoltzmannDistribution[K_, z_, x_Symbol, m_]:=
	Module[{zz=z/.x->Exp[K]},
		Map[
			{#/m, Coefficient[z, x, #]*Exp[#*K]/zz}&,
			Sort[Exponent[z, x, List]]
		]
	]


BoltzmannDistribution0[K_, listK_List, m_, r_:False]:=
	Module[{sum, list},
		list = Map[{#[[1]], m*(K - #[[2]])}&, listK];
		list = listIntegrate0[list, r];
		list = Map[{#[[1]], Exp[#[[2]]]}&, list];
		sum  = Apply[Plus, Map[Last, list]];
		Map[{#[[1]], #[[2]]/sum}&, list]
	]


BoltzmannDistribution1[K_, listK_List, m_]:=
	Module[{sum, list},
		list = Map[{#[[1]], m*(K - #[[2]])}&, listK];
		list = listIntegrate1[list];
		list = Map[{#[[1]], Exp[#[[2]]]}&, list];
		sum  = Apply[Plus, Map[Last, list]];
		Map[{#[[1]], #[[2]]/sum}&, list]
	]


CouplingList[K_, dist_List, m_]:=
	Module[{list1, list2},
		list1 = listDerivative[dist];
		list2 = listAverage[dist];
		MapThread[
			{#1[[1]], K - #1[[2]]/#2[[2]]/m}&, 
			{list1, list2}
		]
	]


CumulantD[dist_, 1] := ExpectationD[dist]

CumulantD[dist_, k_] :=
	Module[{c, m, n, z},
		c = D[Log[z[x]], {x, k}];
		c = c/.Derivative[n_][z][x] -> m[n]*z[x];
		c = c/.m[1] -> 0;
		Do[m[n] = CentralMomentD[dist, n], {n, 1, k}];
		c
	]



EntropyD[dist:{__List}] :=
	-Apply[Plus, Map[(#[[2]]*Log[#[[2]]])&, dist]]


ExpectationD[dist:{__List}] :=
	Apply[Plus, Apply[Times, dist, {1}]]

ExpectationD[dist:{__List}, k_Integer]:=
	Apply[Plus, Map[(#[[2]] * #[[1]]^k)&, dist]]

ExpectationD[dist:{__List}, func_]:=
	Apply[Plus, Map[(#[[2]] * func[#[[1]]])&, dist]]


FractileD[dist:{__List}, frac_] :=
	Module[{sum, pos},
		sum = Rest[FoldList[Plus, 0, Map[Last, dist]]];
		pos = Position[sum, x_/;x>=frac, Infinity, 1][[1, 1]];
		dist[[pos, 1]]
	]


CentralMomentD[dist:{__List}, 0] := 1

CentralMomentD[dist:{__List}, 1] := 0

CentralMomentD[dist:{__List}, k_Integer] :=
	Dot[
		(Map[First, dist] - ExpectationD[dist])^k,
		Map[Last, dist]
	]


Normal01[dist:{__List}] := NormalFunction01[dist]


NormalFunction01[dist:{__List}] := Interpolation[NormalList01[dist]]


NormalList01[dist:{__List}] :=
	Module[{a, d, n, s, tmp},
		tmp = Map[First, dist];
		d = Min[Drop[tmp, 1] - Drop[tmp, -1]];
		n = Apply[Plus, Map[Last, dist]];
		tmp = Map[{#[[1]], #[[2]]/n}&, dist];
		a = ExpectationD[tmp];
		s = Sqrt[CentralMomentD[tmp,2]];
		Map[{(#[[1]]-a)/s, s*#[[2]]/d}&, tmp]
	]


Scan[Protect, definitions]

End[]

EndPackage[]