(* ::Package:: *)

(* :Title: Zedhat *)

(* :Author: Nikita A. Sopenko *)

(* :Summary:
   <summary>
 *)

(* :Copyright: \[Copyright] July 2019 by California Institute of Technology *)

(* :Package Version: 1.0 *)

(* :Mathematica Version: 12.0.0.0 *)

(* :History:
   1.0 for Programming in Mathematica, 1st ed.
*)

(* :Keywords: Chern-Simons, Z-hat, ... *)

(* :Warnings:
   <description of global effects, incompatibilities>
*)

(* :Limitations:
	SU(2) invariants only
   <special cases not handled, known problems>
*)

(* :Discussion:
   <description of algorithm, information for experts>
*)

(* :Requirements:
   <requirements>
*)

(* :Examples:
   <sample input that demonstrates the features of this package>
*)


BeginPackage["Zedhat`"]


	(* Chern-Simons S and T *)

	Snm::usage = "Snm[n,m] gives matrix elements for S-matrix; Snm[n,m,k] is a specialization of k."
	Tn::usage = "Tn[n] gives diagonal elements for T-matrix; Tn[n,k] is a specialization of k."
	ff::usage = "ff[n] gives the n-th power of framing multiplier; ff[n,k] is a specialization of k."
	SM::usage = "SM[k]: S matrix k x k."
	TM::usage = "TM[k]: T matrix k x k."
	TMp::usage = "TMp[k,p]: T^p matrix k x k."
	
	
	(* False theta-funcitons *)
	
	Psi::usage = "Psi[q,p,a,nmax] is a false theta \[CapitalPsi]_{p,a}(q) up to order n"
	Psir::usage = "Psir[p,a,k] is a false theta \[CapitalPsi]_{p,a}(q) for q = Exp[2 \[Pi] i / k]"
	
	
	(* Graph data *)
	
	MforGraph::usage = "MforGraph[graph,signs]: M for graph with choice signs."
	PlotGraph::usage = "PlotGraph[graph,sign] plots a graph with a choive of sings."
	ZLMZL::usage = "ZLMZL[M] gives a list of representatives of Z^L/M Z^L for a matrix M."
	blables::usage = "blables[graph,sign] gives a set of b lables for a graph for a choice of signs"
	
		
	(* WRT invariant *)
	
	ZWRT::usage = "ZWRT[k,graph] computes a numeric value of WRT for a graph at level k."
	CRound::usage = "CRound[N,\[Epsilon]]: sets Re(N) and Im(N) to zero if less than \[Epsilon]."
	
	
	(* Laplace transformation *)
	
	LaplTr::usage = "LaplTr[q,M,zs,b,expression] (in expression) replaces z^n by q^(n M^(-1) n / 4) if n = 2 M l + b and by 0 otherwise."
	
	
	(* Zed-hat *)
	
	Zhat::usage = "Zhat[q,graph,b,kmax,qmax,side,sign] invariant for a graph with label b, choice of signs for edges, and side of the unit disk."
	
	
(* error messages for the exported objects *)



Begin["`Private`"]


(* definition of auxiliary functions and local (static) variables *)

	Eij[i_,j_,n_]:=Table[KroneckerDelta[i,ii]KroneckerDelta[j,jj],{ii,n},{jj,n}];
	\[Sigma]M[graph_,sign_]:=Total[Sign[Eigenvalues[MforGraph[graph,sign]]]]
	bplus[graph_,sign_]:=(\[Sigma]M[graph,sign]+Length[graph[[1]]])/2
	degV[graph_]:=With[{L=graph[[1]]//Length},Total/@Table[If[i==j,2-2Count[graph[[2]],i\[UndirectedEdge]i],-(Count[graph[[2]],i\[UndirectedEdge]j]+Count[graph[[2]],j\[UndirectedEdge]i])],{i,L},{j,L}]];
	seVanderM[dg_,arg_,kmax_]:=1/2 (Normal[Series[(x-x^-1)^dg,{x,0,kmax}]]+Normal[Series[(x-x^-1)^dg,{x,\[Infinity],kmax}]])/.{x->arg}//Expand;
	LaplTr[q_,M_,zs_,b_,expr_]:=With[{MM=Inverse[M]},(((With[{l=Exponent[#,zs]},If[And@@(IntegerQ/@Flatten[1/2 (MM.(l-b))]),q^(l.MM.l/4) #/Times@@Flatten[zs^l],0]])&/@(List@@(Expand[expr+xxx-1])))/.{xxx->1}//Total)]



(* definition of the exported functions *)
	
	(* Chern-Simons S and T *)

	Snm[n_,m_]:=-(I/Sqrt[2k])(q^(n m/2)-q^(-n m/2));
	Snm[n_,m_,kk_]:=-(I/Sqrt[2 kk])(q^(n m/2)-q^(-n m/2))/.{q->E^((2 \[Pi] I)/kk)};
	
	Tn[n_]:= q^((n^2-1)/4);
	Tn[n_,kk_] := q^((n^2-1)/4)/.{q->E^((2 \[Pi] I)/kk)};
	
	ff[fr_]=(E^(-((\[Pi] I)/4)) q^(1/4))^fr;
	ff[fr_,kk_]=(E^(-((\[Pi] I)/4)) q^(1/4))^fr/.{q->E^((2 \[Pi] I)/kk)};

	SM[kn_]:=Table[Snm[n,m]/.{k->kn}/.{q->E^((2 \[Pi] I)/kn)},{n,1,kn-1},{m,1,kn-1}]//N
	TM[kn_]:=Table[KroneckerDelta[n,m]Tn[n]/.{q->E^((2 \[Pi] I)/kn)},{n,1,kn-1},{m,1,kn-1}]//N
	TMp[kn_,p_]:=Table[KroneckerDelta[n,m]Tn[n]^p/.{q->E^((2 \[Pi] I)/kn)},{n,1,kn-1},{m,1,kn-1}]//N
	


	(* False theta-functions *)
	
	Psi[q_,p_,a_,nmax_]:=Sum[If[Mod[n,2p]==a,1,If[Mod[n,2p]==2p-a,-1,0]]q^(n^2/(4 p)),{n,0,nmax}]
	Psir[p_,a_,kk_]:=Limit[Sum[If[Mod[n,2p]==a,1,If[Mod[n,2p]==2p-a,-1,0]] E^((2\[Pi] I)/kk n^2/(4 p)) t^n/(1-t^(2p kk)),{n,0,2p kk-1}],t -> 1]


	(* Graph data *)
	
	MforGraph[graph_,sign_]:=With[{L=graph[[1]]//Length},DiagonalMatrix[graph[[1]]]+Sum[sign[[i]]With[{ed=graph[[2,i]]},(Eij[ed[[1]],ed[[2]],L]+Eij[ed[[2]],ed[[1]],L])],{i,Length[graph[[2]]]}]];
	MforGraph[graph_]:=MforGraph[graph,ConstantArray[1,Length[graph[[2]]]]]
	PlotGraph[graph_,sign_]:=Graph[Table[Labeled[i,graph[[1,i]]],{i,Length[graph[[1]]]}],Table[Labeled[graph[[2,e]],If[sign[[e]]==1,"+","-"]],{e,Length[graph[[2]]]}]]
	PlotGraph[graph_]:=PlotGraph[graph,ConstantArray[1,Length[graph[[2]]]]]
	ZLMZL[M_]:=With[{\[CapitalKappa]=Transpose[M],\[CapitalKappa]i=Inverse[Transpose[M]]},(#[[Round[Length[#]/2]]])&/@Gather[Tuples[Range[-#,#]&/@(Total[Abs[\[CapitalKappa]]])],(And@@(IntegerQ/@(\[CapitalKappa]i.(#1-#2))))&]]
	blables[graph_,sign_]:=(#+Mod[degV[graph],2])&/@(2 ZLMZL[MforGraph[graph,sign]])
	blables[graph_]:=blables[graph,ConstantArray[1,Length[graph[[2]]]]]
	


	(* WRT invariant *)
	
	ZWRT[kk_,graph_]:=With[{deg=Total/@Table[If[i==j,2-2Count[graph[[2]],i\[UndirectedEdge]i],-(Count[graph[[2]],i\[UndirectedEdge]j]+Count[graph[[2]],j\[UndirectedEdge]i])],{i,Length[graph[[1]]]},{j,Length[graph[[1]]]}],s=SM[kk],t=TM[kk]}, Sum[(Product[t[[n[[i]],n[[i]]]]^graph[[1,i]] s[[1,n[[i]]]]^deg[[i]],{i,Length[graph[[1]]]}]Product[s[[n[[e[[1]]]],n[[e[[2]]]]]],{e,graph[[2]]}]),{n,Tuples[Range[1,kk-1],Length[graph[[1]]]]}]]	
	CRound[cnum_,\[Epsilon]_]:=If[Abs[Re[cnum]]<\[Epsilon],0,Re[cnum]]+I If[Abs[Im[cnum]]<\[Epsilon],0,Im[cnum]]	
	


	(* Laplace transformation *)
	
	LaplTr[q_,M_,zs_,b_,expr_]:=With[{MM=Inverse[M]},(((With[{l=Exponent[#,zs]},If[And@@(IntegerQ/@Flatten[1/2 (MM.(l-b))]),q^(l.MM.l/4) #/Times@@Flatten[zs^l],0]])&/@(List@@(Expand[expr+xxx-1])))/.{xxx->1}//Total)]
	


	(* Zed-hat *)
	
	Zhat[q_,graph_,b_,kmax_,qmax_,side_,sign_]:=Module[{L,M,sign0,zs,deg,ser,ser2},
		L=graph[[1]]//Length;
		\[CapitalEpsilon]=graph[[2]]//Length;
		M=MforGraph[graph,sign];
		sign0=ConstantArray[1,\[CapitalEpsilon]];
		zs=(z[#])&/@Range[L];
		deg=degV[graph];
		ser=((-1)^(bplus[graph,sign0]) (Times@@sign) E^((\[Pi] I)/4 (\[Sigma]M[graph,sign]-\[Sigma]M[graph,sign0])) q^((3 \[Sigma]M[graph,sign0]-Total[graph[[1]]])/4))  LaplTr[q,-M,zs,b,Product[seVanderM[deg[[i]],z[i],kmax],{i,L}]];
		ser2=Simplify[(Series[ser,{q,0,qmax}]//Normal)/.{q->q^side}, Assumptions->q>0]
	]
	Zhat[q_,graph_,b_,kmax_,qmax_,side_]:=Zhat[q,graph,b,kmax,qmax,side,ConstantArray[1,Length[graph[[2]]]]]
	Zhat[q_,graph_,b_,kmax_,qmax_]:=Zhat[q,graph,b,kmax,qmax,1,ConstantArray[1,Length[graph[[2]]]]]
	


End[ ]
	(* Protect[ ... ] *)    (* protect exported symbols *)


EndPackage[ ] 
