Quit[]

<<Combinatorica`

(* Define zeon multiplication *)

Unprotect[Times];(Subscript[ζ, a_])(Subscript[ζ, b_]):=If[Length[a⋂b]>0,0,(Subscript[ζ, a⋃b])];
Protect[Times];
Unprotect[Power];
(((x_/;!FreeQ[x,(Subscript[ζ, a_])]))^(n_Integer)):=Module[{y,f},y=Expand[x];Switch[EvenQ[n],True,If[n==0,Return[1],((Composition[Expand][Distribute[f[x,y]]/.f->Times])^(n/2))],False,If[n==1,Return[x],((Composition[Expand][Distribute[f[y,y]]/.f->Times])^((n-1)/2))x]]];
Protect[Power];Unprotect[Expand];
Expand[x_/;!FreeQ[x,(Subscript[ζ, _])]]:=DeleteCases[Distribute[x,Plus,Times],0.`(Subscript[ζ, _])];Protect[Expand];

GenRandGraph[n_,p_]:=Module[{i,j,A,rw,co},A=Table[0,{i,1,n},{j,1,n}];For[rw=1,rw≤n,rw++,For[co=1,co<rw,co++,If[RandomReal[]<p,A[[rw,co]]=1;A[[co,rw]]=1;,Null];]];Return[A]];

(* Count k-cycles in graph using Combinatorica's HamiltonianCycle algorithm *)
KCycleCount[A_,k_]:=Module[{list, Kset},
list=Table[i,{i,1,Length[A]}];
Kset=Subsets[list,{k}];
Return[Sum[Length[HamiltonianCycle[FromAdjacencyMatrix[A[[Kset[[ell]],Kset[[ell]]]]],All]],{ell,1,Length[Kset]}]]]

n=8;
A=GenRandGraph[n,.5];
GraphPlot[A]
KCycleCount[A,Floor[n/2]]

GraphPlot[A]

(* Count H.C.'s  using approach of E.T. Bax *)
HCycleCount[A_]:=Module[{list, Kset,verts,nulltable,ell},
list=Table[i,{i,1,Length[A]}];
verts=Length[A];
Kset=Subsets[list];
nulltable=Table[0,{i,1,verts},{j,1,verts}];
Return[Sum[(((-1))^((verts-Length[Kset[[ell]]])))MatrixPower[(A ReplacePart[nulltable,CartesianProduct[Kset[[ell]],Kset[[ell]]]->1]),verts],{ell,1,Length[Kset]}]]]

HCycleCount[(A)]//MatrixForm

KCycleCount[A,3]

(* Count all cycles using approach of E.T. Bax *)
AllCycleCount[A_]:=Module[{list, Kset,verts,nulltable,ell},
list=Table[i,{i,1,Length[A]}];
verts=Length[A];
Kset=Subsets[list,{1,verts}];
nulltable=Table[0,{i,1,verts},{j,1,verts}];
Return[Sum[MatrixPower[(A ReplacePart[nulltable,CartesianProduct[Kset[[ell]],Kset[[ell]]]->1]),Length[Kset[[ell]]]].MatrixPower[(IdentityMatrix[verts]-A ReplacePart[nulltable,CartesianProduct[Kset[[ell]],Kset[[ell]]]->1]),verts-Length[Kset[[ell]]]],{ell,1,Length[Kset]}]]]

AllCycleCount[A]//MatrixForm

(* Count k cycles  using approach of E.T. Bax *)
BaxKCycleCount[A_,k_]:=Module[{list, Kset,verts,nulltable,Ksub,ell},
list=Table[i,{i,1,Length[A]}];
verts=Length[A];
Kset=Subsets[list,{1,k}];
nulltable=Table[0,{i,1,verts},{j,1,verts}];
Return[Sum[(((-1))^(k-Length[Kset[[ell]]]))Binomial[n-Length[Kset[[ell]]],k-Length[Kset[[ell]]]]MatrixPower[(A ReplacePart[nulltable,CartesianProduct[Kset[[ell]],Kset[[ell]]]->1]),k],{ell,1,Length[Kset]}]];]

Tr[BaxKCycleCount[A,5]]/5//MatrixForm

BaxKCycleCount[A,5]//MatrixForm

KCycleCount[A,5]

BaxKCycleCount[A,Floor[n/2]]//MatrixForm

(* Randomly-generated Graphs on n vertices *)
(* Cycles of length k satisfying 3≤k≤Floor[n/2]*)
(* are enumerated first by nilpotent adjacency matrix method *)
(* then by Bax HC method applied to all n/2 subgraphs *)
(* then using Tarjan-type approach from Combinatorica package *)
(* Computation times in seconds are averaged over 20 trials and plotted *)

XZ={};
XC={};
XB={};

For[Trials=1,Trials≤20,Trials++,
nstart=6;
nend=12;
Outputtable={{"n","p","Zeon Time","Graph","Bax Time","CombiTarjan Time","cycle size", "#{k-cycles}"}};
EdgeProbability=.2;
Timetable={};
Timetable3=Timetable;
Timetable4=Timetable;

For[n=nstart,n≤nend,n++,
k=Max[{3,Floor[n/2]}];

ClearSystemCache[];
a=GenRandGraph[n,EdgeProbability];
g=FromAdjacencyMatrix[a];

(*Count cycles using zeon approach *)

A=a.DiagonalMatrix[Table[(Subscript[ζ, {i}]),{i,1,n}]];
starttime=AbsoluteTime[];
nilnum=(Simplify[Tr[Expand[MatrixPower[A,k]]]]/k)/.{(Subscript[ζ, _])->1};
endtime=AbsoluteTime[];

elapsed=endtime-starttime;
Timetable=Append[Timetable,elapsed];

ClearSystemCache[];

(* Use Bax method *)
starttime=AbsoluteTime[];
knum=Tr[BaxKCycleCount[a,k]]/k;
endtime=AbsoluteTime[];

(* Check that results agree *)
If[nilnum≠knum,Print[nilnum, " ", knum],];

elapsed3=endtime-starttime;  (* Time to count cycles *)
Timetable3=Append[Timetable3,elapsed3];

(* Combi-Tarjan method *)
ClearSystemCache[];
starttime=AbsoluteTime[];
Cknum=KCycleCount[a,k];
endtime=AbsoluteTime[];

(* Check that results agree *)
If[nilnum≠Cknum,Print[nilnum, " ",knum," ", Cknum],];

elapsed4=endtime-starttime;  (* Time to count cycles *)

Timetable4=Append[Timetable4,elapsed4];

(* Provide occasional runtime feedback *)
Print[n, " of trial ",Trials," done."];

Outputtable=Append[Outputtable, { n, EdgeProbability,elapsed, Framed@GraphPlot[a],elapsed3,elapsed4,k,knum}];];

XB=Append[XB,Timetable3];
XZ=Append[XZ,Timetable];
XC=Append[XC,Timetable4];

(* Provide occasional detailed feedback *)
Print[Grid[Outputtable]];];

(* Given comparison summary Bax v. zeon*)
Print[" --- Complexity Record --- "]
Print["RED Nilpotent adjacency matrix method"];
Print["BLUE Bax method"];

Show[{ListPlot[Mean[XZ], PlotRange->All,PlotStyle->{RGBColor[.6,0,0],AbsolutePointSize[5]},AxesLabel->{"vertices","seconds"},DataRange->{nstart,nend},AxesOrigin->{nstart,0}],ListPlot[Mean[XB],PlotMarkers->{"B"},PlotRange->All, PlotStyle->{RGBColor[0,0,.6], AbsolutePointSize[5]},DataRange->{nstart,nend},AxesOrigin->{nstart,0}] }]

(* Generate random graph on n vertices with e equiprobable edges  *)
EdgeConstrainedGraph[n_,e_]:=Module[{i,j,A,rw,co, Ω},
Ω=Binomial[n,2];
g=CompleteGraph[n];
For[i=1,i≤Ω-e,i++,
g=DeleteEdges[g,Edges[g][[Random[Integer,{1,M[g]}]]]];];
Return[ToAdjacencyMatrix[g]];]

(* Randomly-generated Graphs on n vertices and e ≤ q(Ω/k-1) edges *)
(* Cycles of length 3 ≤ k ≤ Floor[n/2] are enumerated first by zeon method *)
(* then by Bax HC method applied to all k-vertex subgraphs *)
(* then using Tarjan-type approach from Combinatorica package *)
(* Computation times in seconds are averaged over 20 trials and plotted *)

XZ={};
XC={};
XB={};

For[Trials=1,Trials≤20,Trials++,

nstart=10; (* Vertices in smallest graphs *)
nend=16;  (* Vertices in largest graphs *)

q=.95;

Outputtable={{"n","Zeon Time","Graph","Bax Time","CombiTarjan Time","cycle size", "#{k-cycles}"}};



Timetable={};
Timetable3=Timetable;
Timetable4=Timetable;

For[n=nstart,n≤nend,n++,

Ω=Binomial[n,2];


(* Choose cycle length first, then generate graph *)

k=Random[Integer,{3,Max[{3, Floor[n/2]}]}]; (*Cycle length*)

e=Floor[q(Ω/k-1)]; (* Maximum edges *)    

(*

(* Choose number of edges first, then select size of cycles to be counted *)

e=Random[Integer,{3,q(Ω/3-1)}];     (*# edges first*)

k=Random[Integer,{3,q Ω/(e+q)}]; (*Cycle length*)

k=Floor[q Ω/(e+q)]; (* Max Cycle Length *)

*)

ClearSystemCache[];

(*EdgeProbability=q/(k-1);
ps=Append[ps,EdgeProbability];*)

ClearSystemCache[];
a=EdgeConstrainedGraph[n,e];
g=FromAdjacencyMatrix[a];

(*Count cycles using zeon approach *)

A=a.DiagonalMatrix[Table[(Subscript[ζ, {i}]),{i,1,n}]];
starttime=AbsoluteTime[];
nilnum=(Simplify[Tr[Expand[MatrixPower[A,k]]]]/k)/.{(Subscript[ζ, _])->1};
endtime=AbsoluteTime[];

elapsed=endtime-starttime;
Timetable=Append[Timetable,elapsed];

ClearSystemCache[];

(* Use Bax method *)
starttime=AbsoluteTime[];
knum=Tr[BaxKCycleCount[a,k]]/k;
endtime=AbsoluteTime[];

(* Check that results agree *)
If[nilnum≠knum,Print[nilnum, " ", knum],];

elapsed3=endtime-starttime;  (* Time to count cycles *)
Timetable3=Append[Timetable3,elapsed3];

(* Combi-Tarjan method *)
ClearSystemCache[];
starttime=AbsoluteTime[];
Cknum=KCycleCount[a,k];
endtime=AbsoluteTime[];

(* Check that results agree *)
If[nilnum≠Cknum,Print[nilnum, " ",knum," ", Cknum],];

elapsed4=(endtime-starttime);  (* Time to count cycles *)

Timetable4=Append[Timetable4,elapsed4];

(* Provide occasional runtime feedback *)
Print[n, " of trial ",Trials," done."];

Outputtable=Append[Outputtable, { n,elapsed, Framed@GraphPlot[a],elapsed3,elapsed4,k,knum}];];

XB=Append[XB,Timetable3];
XZ=Append[XZ,Timetable];
XC=Append[XC,Timetable4];

(* Provide occasional detailed feedback *)
Print[Grid[Outputtable]];];

(* Given comparison summary *)
Print[" --- Complexity Record --- "]
Print["RED Nilpotent adjacency matrix method"];
Print["BLUE Bax method"];
Print["GREEN Combinatorica method"];
Show[{ListPlot[Mean[XZ], PlotRange->All,PlotStyle->{RGBColor[.6,0,0],AbsolutePointSize[5]},AxesLabel->{"vertices","seconds"},DataRange->{nstart,nend},AxesOrigin->{nstart,0}],ListPlot[Mean[XB],PlotMarkers->{"B"},PlotRange->All, PlotStyle->{RGBColor[0,0,.6], AbsolutePointSize[5]},DataRange->{nstart,nend},AxesOrigin->{nstart,0}],ListPlot[Mean[XC],PlotMarkers->{"C"},PlotRange->All, PlotStyle->{RGBColor[0,.6,0], AbsolutePointSize[5]},DataRange->{nstart,nend},AxesOrigin->{nstart,0}] }]

(* Randomly-generated Graphs on n vertices *)
(* Cycles of random length 3≤ k ≤ n/2 *)(* are enumerated by nilpotent adjacency matrix method *)
(* Constraints on number of edges in graph *)
(* Mean times are compared with theoretical curve (n^4)(((1+q))^n) *)

XZ={};
For[Trials=1,Trials≤100,Trials++,
nstart=6;
nend=26;
q=0.9;
Outputtable={{"n","q","Zeon Time","Graph","cycle size", "#{k-cycles}"}};

s=Table[ (n^4)(((1+q))^n) ,{n,nstart,nend}];
Timetable={};
Timetable3=Timetable;
Timetable4=Timetable;
For[n=nstart,n≤nend,n++,
Ω=Binomial[n,2];


(* Choose cycle length first, then generate graph *)

k=Random[Integer,{3,Max[{3, Floor[n/2]}]}]; (*Cycle length*)

e=Floor[q(Ω/k-1)]; (* Maximum edges *)    

(* Alternate approach available by uncommenting this segment in place of previous *)

(* Choose number of edges first, then select size of cycles to be counted *)

e=Random[Integer,{3,q(Ω/3-1)}];     (*# edges first*)

k=Random[Integer,{3,q Ω/(e+q)}]; (*Cycle length*)

k=Floor[q Ω/(e+q)]; (* Max Cycle Length *)



ClearSystemCache[];

(*EdgeProbability=q/(k-1);
ps=Append[ps,EdgeProbability];*)

ClearSystemCache[];
a=EdgeConstrainedGraph[n,e];
g=FromAdjacencyMatrix[a];

A=a.DiagonalMatrix[Table[(Subscript[ζ, {i}]),{i,1,n}]];

starttime=AbsoluteTime[];
nilnum=(Simplify[Tr[Expand[MatrixPower[A,k]]]]/k)/.{(Subscript[ζ, _])->1};
endtime=AbsoluteTime[];


elapsed=endtime-starttime;

Timetable=Append[Timetable,elapsed];


If[Mod[n,2]==1,Print[n, " done of trial ",Trials],];
Outputtable=Append[Outputtable, { n, q,elapsed, Framed@GraphPlot[a],k,nilnum}];];
XZ=Append[XZ,Timetable];

Print[Grid[Outputtable]];
];
Clear[μ];
mu=Solve[D[Expand[Total[(((μ s-Mean[XZ]))^2)]],μ]==0][[1]];
Print["μ = ",mu];
Show[{ListPlot[Mean[XZ], PlotRange->All,PlotStyle->{RGBColor[.5,0,0],AbsolutePointSize[5]},AxesLabel->{"vertices","seconds"},DataRange->{nstart,nend},AxesOrigin->{nstart,0}],
ListPlot[s μ/.mu,DataRange->{nstart,nend},AxesOrigin->{nstart,0},PlotMarkers->{"*"},PlotStyle->{RGBColor[0,0,.5],AbsolutePointSize[5]}]}]
ListPlot[Mean[XZ], PlotRange->All,PlotStyle->{RGBColor[.5,0,0],AbsolutePointSize[5]},AxesLabel->{"vertices","seconds"},DataRange->{nstart,nend},AxesOrigin->{nstart,0}]

mu

(* Randomly-generated Graphs on n vertices *)
(* Cycles of length 6 are enumerated first by nilpotent adjacency matrix method *)
(* then by Bax HC method applied to all n/2 subgraphs *)
(* then using Tarjan-type approach from Combinatorica package *)
(* Computation times in seconds are averaged over 20 trials and plotted *)

XZ={};
XC={};
XB={};

For[Trials=1,Trials≤20,Trials++,
nstart=6;
nend=15;
Outputtable={{"n","p","Zeon Time","Graph","Bax Time", "CombiTarjan Time","cycle size", "#{k-cycles}"}};
EdgeProbability=.25;
Timetable={};
Timetable3=Timetable;
Timetable4=Timetable;

For[n=nstart,n≤nend,n++,
k=6;

ClearSystemCache[];
a=GenRandGraph[n,EdgeProbability];
g=FromAdjacencyMatrix[a];

(*Count cycles using zeon approach *)

A=a.DiagonalMatrix[Table[(Subscript[ζ, {i}]),{i,1,n}]];
starttime=AbsoluteTime[];
nilnum=(Simplify[Tr[Expand[MatrixPower[A,k]]]]/k)/.{(Subscript[ζ, _])->1};
endtime=AbsoluteTime[];

elapsed=endtime-starttime;
Timetable=Append[Timetable,elapsed];

ClearSystemCache[];

(* Use Bax method *)
starttime=AbsoluteTime[];
knum=Tr[BaxKCycleCount[a,k]]/k;
endtime=AbsoluteTime[];

(* Check that results agree *)
If[nilnum≠knum,Print[nilnum, " ", knum],];

elapsed3=endtime-starttime;  (* Time to count cycles *)
Timetable3=Append[Timetable3,elapsed3];

(* Combi-Tarjan method *)
ClearSystemCache[];
starttime=AbsoluteTime[];
Cknum=KCycleCount[a,k];
endtime=AbsoluteTime[];

(* Check that results agree *)
If[nilnum≠Cknum,Print[nilnum, " ",knum," ", Cknum],];

elapsed4=endtime-starttime;  (* Time to count cycles *)

Timetable4=Append[Timetable4,elapsed4];

(* Provide occasional runtime feedback *)
If[Mod[n,2]==1,Print[n, " of trial ",Trials," done."],];

Outputtable=Append[Outputtable, { n, EdgeProbability,elapsed, Framed@GraphPlot[a],elapsed3,elapsed4, k,knum}];];

XB=Append[XB,Timetable3];
XZ=Append[XZ,Timetable];
XC=Append[XC,Timetable4];

(* Provide occasional detailed feedback *)
If[Mod[Trials,2]==1,Print[Grid[Outputtable]],];];

(* Given comparison summary *)
Print[" --- Complexity Record --- "]
Print["RED Nilpotent adjacency matrix method"];
Print["BLUE Bax method"];
Print["GREEN Combinatorica method"];
Show[{ListPlot[Mean[XZ], PlotRange->All,PlotStyle->{RGBColor[.6,0,0],AbsolutePointSize[5]},AxesLabel->{"vertices","seconds"},DataRange->{nstart,nend},AxesOrigin->{nstart,0}],ListPlot[Mean[XB],PlotMarkers->{"B"},PlotRange->All, PlotStyle->{RGBColor[0,0,.6], AbsolutePointSize[5]},DataRange->{nstart,nend},AxesOrigin->{nstart,0}],ListPlot[Mean[XC],PlotMarkers->{"C"},PlotRange->All, PlotStyle->{RGBColor[0,.6,0], AbsolutePointSize[5]},DataRange->{nstart,nend},AxesOrigin->{nstart,0}] }]

(* Randomly-generated Graphs on n vertices *)
(* Cycles of length 5 are enumerated first by nilpotent adjacency matrix method *)
(* then by Bax HC method applied to all 5 subgraphs *)
(* then using Tarjan-type approach from Combinatorica package *)
(* Computation times in seconds are averaged over 20 trials and plotted *)

XZ={};
XC={};
XB={};

For[Trials=1,Trials≤20,Trials++,
nstart=6;
nend=15;
Outputtable={{"n","p","Zeon Time","Graph","Bax Time","CombiTarjan time","cycle size", "#{k-cycles}"}};
EdgeProbability=.3;
Timetable={};
Timetable3=Timetable;
Timetable4=Timetable;

For[n=nstart,n≤nend,n++,
k=5;

ClearSystemCache[];
a=GenRandGraph[n,EdgeProbability];
g=FromAdjacencyMatrix[a];

(*Count cycles using zeon approach *)

A=a.DiagonalMatrix[Table[(Subscript[ζ, {i}]),{i,1,n}]];
starttime=AbsoluteTime[];
nilnum=(Simplify[Tr[Expand[MatrixPower[A,k]]]]/k)/.{(Subscript[ζ, _])->1};
endtime=AbsoluteTime[];

elapsed=endtime-starttime;
Timetable=Append[Timetable,elapsed];

ClearSystemCache[];

(* Use Bax method *)
starttime=AbsoluteTime[];
knum=Tr[BaxKCycleCount[a,k]]/k;
endtime=AbsoluteTime[];

(* Check that results agree *)
If[nilnum≠knum,Print[nilnum, " ", knum],];

elapsed3=endtime-starttime;  (* Time to count cycles *)
Timetable3=Append[Timetable3,elapsed3];

(* Combi-Tarjan method *)
ClearSystemCache[];
starttime=AbsoluteTime[];
Cknum=KCycleCount[a,k];
endtime=AbsoluteTime[];

(* Check that results agree *)
If[nilnum≠Cknum,Print[nilnum, " ",knum," ", Cknum],];

elapsed4=endtime-starttime;  (* Time to count cycles *)

Timetable4=Append[Timetable4,elapsed4];

(* Provide occasional runtime feedback *)
If[Mod[n,5]==1,Print[n, " of trial ",Trials," done."],];

Outputtable=Append[Outputtable, { n, EdgeProbability,elapsed, Framed@GraphPlot[a],elapsed3,elapsed4,k,knum}];];

XB=Append[XB,Timetable3];
XZ=Append[XZ,Timetable];
XC=Append[XC,Timetable4];

(* Provide occasional detailed feedback *)
If[Mod[Trials,5]==1,Print[Grid[Outputtable]],];];

(* Given comparison summary *)
Print[" --- Complexity Record --- "]
Print["RED Nilpotent adjacency matrix method"];
Print["BLUE Bax method"];
Print["GREEN Combinatorica method"];
Show[{ListPlot[Mean[XZ], PlotRange->All,PlotStyle->{RGBColor[.6,0,0],AbsolutePointSize[5]},AxesLabel->{"vertices","seconds"},DataRange->{nstart,nend},AxesOrigin->{nstart,0}],ListPlot[Mean[XB],PlotMarkers->{"B"},PlotRange->All, PlotStyle->{RGBColor[0,0,.6], AbsolutePointSize[5]},DataRange->{nstart,nend},AxesOrigin->{nstart,0}],ListPlot[Mean[XC],PlotMarkers->{"C"},PlotRange->All, PlotStyle->{RGBColor[0,.5,0], AbsolutePointSize[5]},DataRange->{nstart,nend},AxesOrigin->{nstart,0}] }]