pivots[A_]:=Flatten[FirstPosition[1]/@Transpose[A]]

pMatrix[n_]:=Table[Which[i==2&&j==1||i==1&&j==2,-1/2,i==j&&i>2,1,True,0],{i, n},{j,n}]

gramMatrix[W_]:=Simplify[Transpose[W].pMatrix[Length[W]].W]

linearRelationsMatrix[G_]:=RowReduce[G][[1;;MatrixRank[G]]]//Transpose

linearRelations[L_]:=Select[Table[Equal[(Subscript[b, i]),L[[i]].Table[(Subscript[b, i]),{i,pivots[L]}]],{i,Length[L]}],Not[BooleanQ[#]]&]

linearRelationsFromW=linearRelations@*linearRelationsMatrix@*gramMatrix;

inverseLinearRelations[G_]:=Table[If[i==#,1,0],{i,Length[G[[1]]]}]&/@Flatten[FirstPosition[#,1]&/@(RowReduce[G]//Select[Not[AllTrue[#,#==0&]]&])]

quadraticFormMatrix[G_,L_]:=Inverse[Table[G[[i, j]],{i,pivots[L]},{j,pivots[L]}]]

quadraticForm[Q_,L_]:=Module[{v=Table[(Subscript[b, i]),{i,pivots[L]}]},Simplify[v.Q.v==0]]

quadform[W_]:=Module[{G=gramMatrix[W]},Module[{L=linearRelationsMatrix[G]},quadraticForm[quadraticFormMatrix[G,L],L]]]

quadformFromG[G_]:=Module[{L=linearRelationsMatrix[G]},quadraticForm[quadraticFormMatrix[G,L],L]]

graphFromG[G_]:=PlanarGraph[AdjacencyGraph[If[#==-1,1,0]&/@#&/@G,VertexLabels->"Name"]]

nextCandidate[s_,t_,adj_]:=Block[{length,pos},length=Length[adj];
pos=Mod[Position[adj,s][[1,1]]+1,length,1];
{t,adj[[pos]]}];

FindFace[g_?PlanarGraphQ]:=Block[{emb},emb=GraphEmbedding[g,"PlanarEmbedding"];
FindFace[g,emb]];

FindFace[g_?PlanarGraphQ,emb_]:=Block[{m,orderings,pAdj,rightF,s,t,initial,face},m=AdjacencyMatrix[g];
Table[pAdj[v]=SortBy[Pick[VertexList[g],m[[v]],1],ArcTan@@(emb[[v]]-emb[[#]])&],{v,VertexList[g]}];
rightF[_]:=False;
Reap[Table[If[!rightF[e],s=e[[1]];
t=e[[2]];
initial=s;
face={s};
While[t=!=initial,rightF[UndirectedEdge[s,t]]=True;
{s,t}=nextCandidate[s,t,pAdj[t]];
face=Join[face,{s}];];
Sow[face];],{e,EdgeList[g]}]][[2,1]]]

findGenerators[faces_,vertices_,G_]:=Table[Block[{face,otherVertex,verticesOnFace,allVertices,Gp,P,Gnew,L,Q,Qform,sol1,sol2,sols,coeff,coeffTable,small,Li,P2},
face=faces[[i]];
verticesOnFace=Sort[RandomSample[face,3]];
otherVertex=RandomChoice[Select[vertices,Not[MemberQ[face,#]]&]];
allVertices=Append[verticesOnFace,otherVertex];
P=Transpose[IdentityMatrix[Length[vertices]][[#]]&/@DeleteDuplicates[Join[allVertices,vertices]]];
Gnew=Transpose[P].G.P;
L=linearRelationsMatrix[Gnew];
Gp=Table[G[[i,j]],{i,allVertices},{j,allVertices}];
Q=Inverse[Gp];
Qform=(#.Q.#&[Table[(Subscript[b, j]),{j,allVertices}]])==0;
sols=Expand[FullSimplify[Total[#[[1]][[2]]&/@Solve[Qform,(Subscript[b, otherVertex])]]]];
coeff=Table[sols/.Table[(Subscript[b, i])->If[i==curr,1,0],{i,allVertices}],{curr,allVertices}];
small=IdentityMatrix[Length[allVertices]];
small[[-1]]=Table[If[j==4,-1,If[MemberQ[verticesOnFace,allVertices[[j]]],coeff[[j]],0]],{j,Length[allVertices]}];
Li=Table[Table[If[j==k,1,0],{j,Length[vertices]}],{k,allVertices}];
P2=Table[If[j==k,1,0],{j,Length[vertices]},{k,DeleteDuplicates[Join[allVertices,vertices]]}];
P2.L.small.Li
],{i,Length[faces]}]

findGeneratorsN[faces_,vertices_,G_]:=Table[Block[{face,otherVertex,verticesOnFace,allVertices,Gp,P,Gnew,L,Q,Qform,sol1,sol2,sols,coeff,coeffTable,small,Li,P2},
face=faces[[i]];
verticesOnFace=Sort[RandomSample[face,3]];
otherVertex=RandomChoice[Select[vertices,Not[MemberQ[face,#]]&]];
allVertices=Append[verticesOnFace,otherVertex];
P=Transpose[IdentityMatrix[Length[vertices]][[#]]&/@DeleteDuplicates[Join[allVertices,vertices]]];
Gnew=Transpose[P].G.P;
L=linearRelationsMatrix[Gnew];
Gp=Table[G[[i,j]],{i,allVertices},{j,allVertices}];
Q=Inverse[Gp];
Qform=(#.Q.#&[Table[(Subscript[b, j]),{j,allVertices}]])==0;
sols=Expand[FullSimplify[Total[#[[1]][[2]]&/@NSolve[Qform,(Subscript[b, otherVertex])]]]];
coeff=Table[sols/.Table[(Subscript[b, i])->If[i==curr,1,0],{i,allVertices}],{curr,allVertices}];
small=IdentityMatrix[Length[allVertices]];
small[[-1]]=Table[If[j==4,-1,If[MemberQ[verticesOnFace,allVertices[[j]]],coeff[[j]],0]],{j,Length[allVertices]}];
Li=Table[If[j==k,1,0],{j,Length[vertices]},{k,allVertices}];
P2=Table[If[j==k,1,0],{j,Length[vertices]},{k,DeleteDuplicates[Join[allVertices,vertices]]}];
P2.L.small.Li
],{i,Length[faces]}]

findGeneratorsFromG[G_]:=findGeneratorsN[FindFace[PlanarGraph[AdjacencyGraph[If[#==-1,1,0]&/@#&/@G,VertexLabels->"Name"]]],Range[Length[G]],G]

search[generators_,current_,lastgen_,upperBound_,faces_]:=Block[{totals,invalid},
invalid[σ_,face_]:=σ==lastgen||AllTrue[Delete[σ.current,{#}&/@face],#>upperBound&];
totals[i_]:=If[invalid[generators[[i]],faces[[i]]],0,Length[current]-Length[faces[[i]]]-Length[Select[Delete[generators[[i]].current,{#}&/@faces[[i]]],#>upperBound&]]+search[generators,generators[[i]].current,generators[[i]],upperBound,faces]];
Total[totals/@Range[Length[generators]]]]

$RecursionLimit=16384

fractalDimensionFromG[G_,max_,delta_,root_]:=Block[{xs,ys,data,graph,faces,sigmas,nlm,params},
xs=Range[delta,max,delta];
graph=graphFromG[G];
faces=FindFace[graph];
sigmas=findGeneratorsFromG[G];
ys=search[sigmas,root,{},#,faces]+Length[root]&/@xs;
data={xs[[#]],ys[[#]]}&/@Range[Length[ys]];
nlm=NonlinearModelFit[data,a (x^b),{a,b},x];
params=nlm["BestFitParameters"];
Print[nlm["RSquared"]];
b/.params
]

pyramidG[n_]:=Table[Which[i==0&&j==0,1,i==0||j==0,-1,True,1-((2((Sin[((Abs[i-j]π)/n)])^2))/((Sin[(π/n)])^2))],{i,0,n},{j,0,n}]

pyramidFractalDimension[n_,max_,delta_]:=fractalDimensionFromG[N[pyramidG[n]],max,delta,N[Table[If[i==n+1,-1,((Sec[(π/n)]+Tan[(π/n)])/(Tan[(π/n)]))],{i,n+1}]]]

pyramidTuple[n_]:=N[Table[If[i==n+1,-1,((Sec[(π/n)]+Tan[(π/n)])/(Tan[(π/n)]))],{i,n+1}]]