(*The idea behind the function "Boundary", is that the user should be able to input the faces of a pure simplicial complex and mathematica should then be able to determine all the information about the boundary maps.*)

union[X_]:=Fold[Union,{},X] (*This function allows us to recover what the set of vertices is for a pure simplicial complex with faces X*)

Ridges[X_,k_]:=Subsets[ union[X],{k}] (* This function then finds the k-simplices of a pure simplicial complex with faces X. *)

Orientation[x_,y_]:=Piecewise[Table[{(-1)^(j) ,Intersection[Delete[Table[x[[i]],{i,1,Length[x]}],j],y]==y},{j,1,Length[x]}]](*This gives the sign of incidence between k-simplex x and (k-1)-simplex y.*)

Incidence[x_,y_,k_]:=Piecewise[{{Orientation[x,y],Length[Intersection[x,y]]==k-1},{0,!{Length[Intersection[x,y]]==k-1}}}](*This function checks whether the k-simplex x and (k-1)-simplex y are incident or not, if they are it returns the sign of the incidence, if they are not it returns 0.*)

Boundary[X_,k_]:=
Piecewise[{
{Piecewise[{{Table[Incidence[Ridges[X,k][[i]],Ridges[X,k-1][[j]],k],{i,1,Length[Ridges[X,k]]},{j,1,Length[Ridges[X,k-1]]}],k<Length[X[[1]]]},{Table[Incidence[X[[i]],Ridges[X,k-1][[j]],k],{i,1,Length[X]},{j,1,Length[Ridges[X,k-1]]}],k==Length[X[[1]]]}}],MatrixQ[X]},{Print["Please input the faces of a pure simplicial complex"],!{MatrixQ[X]}}}](*Given a pure simplicial complex with faces X, this function computes the kth boundary map. It returns an error if the input X is no*)

(*These are just the functions that allow us to compute the tutte polynomial for a matroid with representation w. "tutte" does this when w is a representation over Q. "tutteq" does this when w is a representation over Fq. Then "flow" and "flowq" respectively specialize "tutte" and "tutteq" *)

power[w_]:=Subsets[w] (*This will be what the sum is taken over *)
r[w_]:=Piecewise[{{MatrixRank[w],Length[w]>0},{0,w={}}}] (*This defines the matroid rank*)
term[w_]:=(x-1)^{-r[w]}*(y-1)^{Length[w]-r[w]}(*These are the monomials that the sum is taken over*)
tutte[w_]:=Simplify[Sum[(x-1)^{r[w]}*term[power[w][[i]]],{i,Length[power[w]]}]](*The sum will be taken over the powerset of w.*)

rq[w_,q_]:=Piecewise[{{MatrixRank[w,Modulus->q],Length[w]>0},{0,w={}}}](*This defines the matroid rank over F_q*)
termq[w_,q_]:=(x-1)^{-rq[w,q]}*(y-1)^{Length[w]-rq[w,q]}(*These are the monomials that the sum is taken over*)
tutteq[w_,q_]:=Simplify[Sum[(x-1)^{rq[w,q]}*termq[power[w][[i]],q],{i,Length[power[w]]}]](*The sum will be taken over the powerset of w*)

flow[w_]:=Factor[(-1)^{Length[w]-r[w]}*tutte[w]/.{x->0,y->1-q}](*This specializes tkr to the number of "flows"*)
flowq[w_,p_]:=Factor[(-1)^{Length[w]-rq[w,p]}*tutteq[w,p]/.{x->0,y->1-q}](*This specializes tkrq to the number of "flows"*)

(*The following are just a combination of the above, to allow us to compute the TKR polynomial for a pure simplicial complex with faces X.*)

tkr[X_,k_]:=tutte[Boundary[X,k]](*this computes the TKR polynomial in dimension k with homology coefficients in Z*)
tkrq[X_,k_,q_]:=tutteq[Boundary[X,k],q](*this computes TKR polynomial in dimension k with homology coefficients in Zq*)
simpflow[X_,k_]:=flow[Boundary[X,k]](*this is properly called the Bott polynomial of dimension k.*)
simpflowq[X_,k_,q_]:=flowq[Boundary[X,k],q](*this is the number of q-flows in dimension k*)
TKR[X_]:=Table[tkr[X,i],{i,1,Length[X[[1]]]}](*this gives you the list of each TKR polynomial*)
FLOW[X_]:=Table[simpflow[X,i],{i,1,Length[X[[1]]]}](*this gives you the list of each Bott polynomial*)
FLOWQ[X_]:=Table[Table[simpflowq[X,k,i];Print[i],{i,{0,2,3,5}}],{k,1,Length[X[[1]]]}](*This gives you the number of q-flows for the first 5 3 primes in each dimension.*)

Chaincomplex[X_]:=Table[MatrixForm[Transpose[Boundary[X,i]]],{i,1,Length[X[[1]]]}](*This lists the boundary matrices.*)

Complete[v_,d_]:=Subsets[Table[i,{i,1,v}],{d}](*This is the complete simplicial complex of dimension d on v vertices.*)

Bridgless[X_,k_]:=For[i=0,i<=Length[X],i++,t=True;If[r[Boundary[X,k]]-r[Delete[Boundary[X,k],i]]==1,t=False];If[t==False,Break[]];If[i==Length[X],Return[t]]] (*This checks whether or not the kth boundary matrix of your simplicial complex with face X. It returns a true or false.*)

Rando[v_,d_,f_]:=RandomSample[Complete[v,d],f](*Gives us a random pure complex of dimension d on v vertices with f faces*)
SearchFlows[v_,d_,f_,n_]:=Do[R=Rando[v,d,f];If[TrueQ[Bridgless[R,d]],Print[simpflow[R,d]];Print[R],Print["Not Bridgless"]];R=.,{i,1,n}](*Searches for flows randomly n times*)

(*EXAMPLES*)

Tetrahedron={{1,2,3},{1,2,4},{1,3,4},{2,3,4}}

Chaincomplex[Tetrahedron]
TKR[Tetrahedron]
FLOW[Tetrahedron]

FourBall=Complete[5,5]

Chaincomplex[FourBall]
TKR[FourBall]
FLOW[FourBall]

Kosar={{1,1,1,0,0,0},{0,0,0,1,1,1},{-1,-1,0,-1,0,0},{0,0,-1,0,-1,-1},{1,0,-1,0,-1,0},{0,1,0,1,0,1}}

MatrixForm[Transpose[Kosar]]
tutte[Kosar]
tutteq[Kosar,2]
flow[Kosar]
flowq[Kosar,2]

RP2={{1,2,4},{1,2,6},{1,3,4},{1,3,5},{1,5,6},{2,3,5},{2,3,6},{2,4,5},{3,4,6},{4,5,6}}

MatrixForm[Transpose[Boundary[RP2,3]]]
tkr[RP2,3]
simpflow[RP2,3]
simpflowq[RP2,3,2]

Peterson={{1,2},{1,5},{1,7},{2,3},{2,8},{3,4},{3,9},{4,5},{4,10},{5,6},{6,8},{6,9},{7,9},{7,10},{8,10}}

GraphPlot[{1->2,1->5,1->7,2->3,2->8,3->4,3->9,4->5,4->10,5->6,6->8,6->9,7->9,7->10,8->10}]

MatrixForm[Transpose[Boundary[Peterson,2]]]
tkr[Peterson,2]
simpflow[Peterson,2]

K42=Complete[4,2]
K53=Complete[5,3]
K64=Complete[6,4]

MatrixForm[Transpose[Boundary[K42,2]]]
MatrixForm[Transpose[Boundary[K53,3]]]
MatrixForm[Transpose[Boundary[K64,4]]]

simpflow[K42,2]
simpflow[K53,3]
simpflow[K64,4]

SearchFlows[5,3,8,10]

SearchFlows[6,4,13,10]

SearchFlows[6,2,10,10]