On[Assert];

(* Run this if restarting *)
Nps = {1, 2, 3, 2, 1};
   dim = 2*Length[Nps] - 2;
   n = dim/2;
   NFP = Plus @@ Nps;
   
   vertexCoordinates = {};
   For[i = 1, i <= Length[Nps], ++i,
     	For[j = 1, j <= Nps[[i]], ++j,
       
       	AppendTo[vertexCoordinates, {j - (Nps[[i]])/2, i}]; 
       ];
     ];
   If[dim == 6,
     	If[NFP == 4,
      		vertexCoordinates = {{0, 0}, {1, 1}, {1, 2}, {0, 3}};]];
   
   If[dim == 8,
     	If[NFP == 6,
      		vertexCoordinates = {{0, 0}, {1, 1}, {2, 2}, {2, 3}, {1, 4}, {0, 5}};]];
   d = dim/2;
   cardI = NFP*dim/4;

lff0=ToExpression[Import[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","graphsb2_2b4_2_61To74.txt"}],"List"]];
lff=ToExpression[Import[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","graphsb2_2b4_2_UpTo60.txt"}],"List"]];

Length[lff0]

lff0[[1]][[1]]

(* check if the reversing works fine to avoid Duplicates *)

listp={};
p=1;
While[p≤  Length[lff0], 
Baux=lff0[[p]]; 
Daux=Table[Baux[[NFP-(i-1)]][[NFP-(j-1)]],{i,1,NFP},{j,1,NFP}];If[Not[MemberQ[lff,Daux]],lff=Append[lff,lff0[[p]]];listp=Append[listp,p]];
p=p+1];

Print[Length[lff]]
Print[listp]

Export[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","graphsb2_2b4_2_Final.txt"}],lff,"List"];

GraphPlot[lff[[1]],DirectedEdges->True,VertexLabeling->True,MultiedgeStyle->True,SelfLoopStyle->True,VertexCoordinateRules->vertexCoordinates,ImageSize->Large]

GraphPlot[lff[[791234]],DirectedEdges->True,VertexLabeling->True,MultiedgeStyle->True,SelfLoopStyle->True,VertexCoordinateRules->vertexCoordinates,ImageSize->Large]

Length[lff]

lff=ToExpression[Import[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","graphsb2_2b4_2_Final.txt"}],"List"]];

(* Part II - List of determinants that are going to be used by the c++ programs *)

(* Attention: When the multigraph has multiple connected components it gives the sum of the squares of the determinants of the different corresponding matrices *) 

Print["**************************************************************************"];

lff=ToExpression[Import[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","graphsb2_2b4_2_Final.txt"}],"List"]];
p=1;
dets = {};
(*NCArray=Table[0,{i,1,Length[lff]}];
auxedges=Table[0,{i,1,Length[lff]}];*)
Monitor[
While[p≤ Length[lff], 
(* Edge function: e[[i]] is the i-th edge of the p-th multigraph of "lff"; e[[i]][[1]] is its initial point and e[[i]][[2]] is its terminal point (cf. Section 4.) *)
e=Array[0&,{cardI,2}];
index = 1;
j=1;
While[j≤ NFP, 
r=1;
While[r≤ NFP, 
entry=lff[[p]][[j]][[r]];
k=1;
While[ k≤ entry,
e[[index]][[1]]=j;
e[[index]][[2]]=r;
index=index+1;
k=k+1
];
r=r+1];
j=j+1];
(* Matrices A: A is the matrix A for the p-th multigraph of "lff" *)

A=Array[0&,{cardI,cardI}];
deltai=Array[0&,{cardI,cardI}];
deltat=Array[0&,{cardI,cardI}];
h=1; 
While[h≤ cardI, 
m=1;
While[m≤ cardI, 
If[ e[[h]][[1]]==e[[m]][[1]] && e[[h]][[1]]≠ e[[m]][[2]] , deltai[[h]][[m]]=1];
If[ e[[h]][[1]]==e[[m]][[2]] && e[[h]][[1]]≠ e[[m]][[1]] , deltai[[h]][[m]]=-1];
If[ e[[h]][[2]]==e[[m]][[1]] && e[[h]][[2]]≠ e[[m]][[2]] , deltat[[h]][[m]]=1];
If[ e[[h]][[2]]==e[[m]][[2]] && e[[h]][[2]]≠ e[[m]][[1]] , deltat[[h]][[m]]=-1];
A[[h]][[m]]= deltai[[h]][[m]]-deltat[[h]][[m]];
      m=m+1];
h=h+1];
          nc=0;
			L=Array[3&, cardI];
			ndet=Simplify[Det[A-DiagonalMatrix[L]]];
			If[ndet==0,AppendTo[dets, p]];
	p=p+1];
, ProgressIndicator[p, {0, Length[lff]}]];

(* Print[lgg] *)

(* Part II - List of matrices that are going to be used *)


Print["**************************************************************************"];
lff=ToExpression[Import[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","graphsb2_2b4_2_Final.txt"}],"List"]];

dets=ToExpression[Import[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","numberofgraphs.txt"}],"List"]]; 
p=1;
e=Array[0&,{ Length[dets],cardI,2}];
(*NCArray=Table[0,{i,1,Length[lff]}];
auxedges=Table[0,{i,1,Length[lff]}];*)
Monitor[
While[p≤ Length[dets], 
(* Edge function: e[[i]] is the i-th edge of the p-th multigraph of "lff"; e[[i]][[1]] is its initial point and e[[i]][[2]] is its terminal point (cf. Section 4.) *)
index = 1;
j=1;
While[j≤ NFP, 
r=1;
While[r≤ NFP, 
entry=lff[[dets[[p]]]][[j]][[r]];
k=1;
While[ k≤ entry,
e[[p]][[index]][[1]]=j;
e[[p]][[index]][[2]]=r;
index=index+1;
k=k+1
];
r=r+1];
j=j+1];
	p=p+1];, ProgressIndicator[p, {0, Length[dets]}]];

L=Array[3&, cardI]

(* Part III - Sorting Matrices  according to the their rank *)

Print[];
Print["**************************************************************************"];
Print["Working with NNMs"];
Print["**************************************************************************"];

dets=ToExpression[Import[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","numberofgraphs.txt"}],"List"]];
A=ToExpression[Import[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","Matrices.txt"}],"List"]];

Matricestoprk={};
Matricessmallrk={};

(* First we sort the partitions according to whether the rank of the matrices
 A-Diag(part) is (cardI-1) or smaller and store them in the files "Part?toprk.txt" and "Part?smallrk.txt"; when the rank is cardI-1 we discard the partitions that originate null spaces with no vectors with positive entries;when the rank is smaller than cardI-1 we discard those whose vectors have an entry equal to zero (they do not satisfy "test 2" below)*)

Print[];
Print["Sorting Matrices according to their rank..."];
k=1; 
test2=Function[x,(Val=False;s=1;
While[s≤cardI,
t=1;
If[Not[Val],
Val=True; 
While[t≤Length[x],
Val=(Val&&(x[[t]][[s]]==0));
t=t+1]];
s=s+1];
Not[Val])];
Monitor[
While[k≤ Length[dets],

 L=Array[3&, cardI];

      Mat=A[[k]]-DiagonalMatrix[L];
  lns=NullSpace[Mat];
  If[MatrixRank[Mat]==cardI-1, 
        s=lns[[1]];
        If[And@@(#1>0&)/@s||And@@(#1<0&)/@s,Matricestoprk=Append[Matricestoprk,k]],
   If[test2[lns], Matricessmallrk=Append[Matricessmallrk,k]]];

k=k+1];, ProgressIndicator[k, {0, Length[dets]}]];
Export[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","Matricestoprk.txt"}],Matricestoprk];
Export[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","Matricessmallrk.txt"}],Matricessmallrk];

Export[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","Matricestoprk.txt"}],Matricestoprk];

(*Part IV - Maximum Rank (RANK = cardI-1) *)

Print[];
Print["**************************************************************************"];
Print["Using matrices with rank equal to cardI-1:"]
A=ToExpression[Import[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","Matrices.txt"}],"List"]];
Matricestoprk=ToExpression[Import[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","Matricestoprk.txt"}],"List"]];
(* Isotropy Weights: Creates list of possible isotropy weights for NNMs and stores them in the array IW: for the p-th NNM, IW[[p]] stores its possible isotropy weights at the different fixed points *)

Print[];
Print["Determining the possible isotropy weights..."]
j=1;
IW=Table[0,{i,1,Length[Matricestoprk]}];
L=Array[3&, cardI];
Print[Length[Matricestoprk]];
While[j≤Length[Matricestoprk],
AuxIW={};
Weight={};
Mat=A[[Matricestoprk[[j]]]]-DiagonalMatrix[L];
lns=Flatten[NullSpace[Mat]]; 
i=1;
AuxWeight={};
While[i≤ NFP, 
W={};
l=1;
While[l≤ cardI, 
If[ e[[Matricestoprk[[j]]]][[l]][[1]]==i,W=Append[W,lns[[l]]]];
If[e[[Matricestoprk[[j]]]][[l]][[2]]==i,W=Append[W,-lns[[l]]]];
l=l+1];
AuxWeight=Append[AuxWeight,W];
i=i+1];
Weight=Append[Weight,AuxWeight];

AuxIW=Append[AuxIW,Weight];
IW[[j]]=Flatten[AuxIW,1];
j=j+1];

(* Using the ABBV localization formula we select those which satisfy ∫c_l = 0 for l=0,...,dim/2-1 and store them in the array "FIW" *)
Print[];
Print["Selecting those for which ∫c_l = 0 for l=0,...,dim/2-1..."];
j=1;
FIW=Table[0,{i,1,Length[Matricestoprk]}];
While[j≤Length[Matricestoprk],
AuxList={};
Print[];
Print["Multigraph #: ", dets[[Matricestoprk[[j]]]]];
k=1;
selectlist={};
While[k≤Length[IW[[j]]],
l=0;
Flag=True;
While[l≤dim/2-1,
Loc=0;
i=1;
While[i≤NFP,Loc=Loc+SymmetricPolynomial[l,IW[[j]][[k]][[i]]]/SymmetricPolynomial[dim/2,IW[[j]][[k]][[i]]];
i=i+1];
If[Loc≠ 0,l=dim/2;Flag=False,l=l+1]];
If[Flag,selectlist=Append[selectlist,IW[[j]][[k]]]];
k=k+1];
If[Length[selectlist]≠0,
Fselectlist=DeleteDuplicates[Table[Sort[selectlist[[i,j]]],{i,Length[selectlist]},{j,NFP}]];
FIW[[j]]=Fselectlist; 
Print[];
Print["The possible isotropy weights are:"];
Print[FIW[[j]]] ,
FIW[[j]]={};Print["There are no isotropy weights"]];
j=j+1];

(* PART V a) - Selection of Possible Null Spaces *)

Print[];
Print["**************************************************************************"];
Print["Using partitions that give matrices with rank smaller than cardI-1:"]

Matricessmallrk=ToExpression[Import[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","Matricessmallrk.txt"}],"List"]];
(* We build the list "PossisbleNullSpaces" that stores all null spaces in RowReduced form, discarding those which have a basis vector with a negative entry while all the other basis vectors have the same entry equal to zero, since the corresponding matrices cannot have a positive vector in their kernels *)

Print[];
Print["Determining the null spaces of the matrices in Row Reduced form when the rank is smaller than cardI-1..."];
j=1;
 L=Array[3&, cardI];
PossibleNullSpaces={};
While[j≤Length[Matricessmallrk],
AuxPossibleNullSpaces={};
flag=True;
Mat=A[[Matricessmallrk[[j]]]]-DiagonalMatrix[L];
lns=RowReduce[NullSpace[Mat]];
l=1;
While[l≤Length[lns],
t=1;
While[t≤cardI,If[lns[[l]][[t]]<0&&Select[Drop[lns[[All,t]],{l}],#≠0&]=={},flag=False];
If[flag,t=t+1,t=cardI+1]];
If[flag,l=l+1,l=Length[lns]+1]];
If[flag,AuxPossibleNullSpaces=Append[AuxPossibleNullSpaces,lns]];
If[Length[AuxPossibleNullSpaces]≠0,PossibleNullSpaces=Append[PossibleNullSpaces,{j,AuxPossibleNullSpaces}]];
Print["Multigraph #", dets[[Matricessmallrk[[j]]]]];
j=j+1];


(* We remove the elements of "PossisbleNullSpaces" for which the span of the above bases has an entry that is always negative, creating a new list of possible null spaces: "PossibleNullSpaces2" *)

Print[];
Print["Removing some Null Spaces whose bases have a negative entry..."];
p=1;
PossibleNullSpaces2={};
While[p≤Length[PossibleNullSpaces],
k=1;
AuxPossibleNullspaces={};
While[k≤Length[PossibleNullSpaces[[p]][[2]]],
flag2=True;
i=1;
While[i≤cardI,If[Select[PossibleNullSpaces[[p]][[2]][[k]][[All,i]],#>0&]=={},flag2=False];
If[flag2,i=i+1,i=cardI+1]];
If[flag2,AuxPossibleNullspaces=Append[AuxPossibleNullspaces,PossibleNullSpaces[[p]][[2]][[k]]]];
k=k+1];
If[Length[AuxPossibleNullspaces]≠0,PossibleNullSpaces2=Append[PossibleNullSpaces2,{PossibleNullSpaces[[p]][[1]],AuxPossibleNullspaces}]];
p=p+1];

p=1;
PNull={};
While[p≤Length[PossibleNullSpaces2], 
PNull=Append[PNull,PossibleNullSpaces2[[p]][[1]]];
p=p+1];
Export[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","PNullSpaces.txt"}],PNull,"List"];

Matricessmallrk=ToExpression[Import[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","Matricessmallrk.txt"}],"List"]];

(* PART V b) ISOTROPY WEIGHTS *)

(* First we consider the span of the different null spaces for NNMs to get the general expressions of the isotropy weights and store these weights in the array "IW2"; the value IW2[[p]] stores the possible isotropy weights at the different fixed points for the p-th NNM of "lff" *)


Print[];
Print["Determining the possible isotropy weights..."];
p=1;
IW2=Table[0,{i,1,Length[lff]}];
While[p≤Length[PossibleNullSpaces2],
Print[];
Print["Multigraph #:",dets[[Matricessmallrk[[PossibleNullSpaces[[p]][[1]]]]]]];
LLength=Length[PossibleNullSpaces2[[p]][[2]]];
(* NLLength=IntegerPart[LLength/10]; *)
k=1;
WeightList={};
While[k≤LLength,
Array[b,Length[PossibleNullSpaces2[[p]][[2]][[k]]]];
i=1;
AuxWeight={};
While[i≤NFP,
l=1;
Weight={};
While[l≤cardI,
 If[ e[[Matricessmallrk[[PossibleNullSpaces2[[p]][[1]]]]]][[l]][[2]]==i,
W=0;
j=1;
While[j≤Length[PossibleNullSpaces2[[p]][[2]][[k]]],W=W-b[j]*PossibleNullSpaces2[[p]][[2]][[k]][[j]][[l]];
j=j+1];
Weight=Append[Weight,W]];
If[ e[[Matricessmallrk[[PossibleNullSpaces2[[p]][[1]]]]]][[l]][[1]]==i,
W=0;
j=1;
While[j≤Length[PossibleNullSpaces2[[p]][[2]][[k]]],W=W+b[j]*PossibleNullSpaces2[[p]][[2]][[k]][[j]][[l]];
j=j+1];
Weight=Append[Weight,W]];
l=l+1];
  AuxWeight=Append[AuxWeight,Weight];
i=i+1];
WeightList=Append[WeightList,AuxWeight];
(* If[NLLength≠0,If[Divisible[k,NLLength],If[k/NLLength ≤10,Print[k/NLLength*10 "%"]]]]; *)
k=k+1]; 
IW2[[PossibleNullSpaces2[[p]][[1]]]]=WeightList;
p=p+1];

(* We then use the fact that the action is effective to elliminate common divisors at each fixed point *)

(* Attention: We Store the Isotropy weights obtained in the files "FIWaux?.txt", where "?" is the number of the graph *)

FIWaux=IW2;
p=1;
Print[];
While[p≤Length[PossibleNullSpaces2],
LLength=Length[PossibleNullSpaces2[[p]][[2]]];
NLLength=IntegerPart[LLength/10];
q=PossibleNullSpaces2[[p]][[1]];
Print["Multigraph #: ",dets[[Matricessmallrk[[q]]]]];
k=1;
While[k≤LLength,
Array[b,Length[PossibleNullSpaces2[[p]][[2]][[k]]]];
i=1;
While[i≤NFP,G=PolynomialGCD@@FIWaux[[q]][[k]][[i]];
If[ToString[G]≠ToString[1]&&Not[NumberQ[G]],
Laux=Solve[G==1,Variables[G][[1]]];
l=1;
aux2=FIWaux[[q]][[k]];
While[l≤NFP,
j=1;
aux={};
While[j≤dim/2,
aux=Append[aux,Simplify[FIWaux[[q]][[k]][[l]][[j]]/.Flatten[Laux]]];
j=j+1];
aux2=ReplacePart[aux2,l->aux];
l=l+1];
FIWaux[[q]]=ReplacePart[FIWaux[[q]],k->aux2]];
i=i+1];
If[NLLength≠0,If[Divisible[k,NLLength],If[k/NLLength ≤10,Print[k/NLLength*10 "%"]]]];k=k+1];
Export[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","FIWaux"<>ToString[dets[[Matricessmallrk[[q]]]]]<>".txt"}],FIWaux[[q]],"List"];
p=p+1];

(* Attention: this part contains several functions that are going to be used in other parts below *)

TestRational=Function[{x,ned},
(Length2=Length[x];
k=1;droplist4={};
While[k≤Length2,
Flag=True;
Array[b,ned];
i=1;
NFPaux=Length[x[[k]]];
While[i≤NFPaux,
j=1;
dimaux=Length[x[[k]][[i]]];
While[j≤dimaux,
v=x[[k]][[i]][[j]];
If[NumberQ[v]&&Element[v,Rationals]&&Not[IntegerQ[v]],Flag=False];
If[Flag,j=j+1,j=dimaux+1]];
If[Flag,i=i+1,i=NFPaux+1]];
If[Not[Flag],droplist4=Append[droplist4,k]];
k=k+1])];

TestPositive=Function[{x,ned},
(k=1;
droplistq={};
While[k≤Length[x],
Array[b,ned];
NFPaux=Length[x[[k]]];
dimaux=Length[x[[k]][[1]]];
Laux=Flatten[Table [x[[k]][[j]][[i]]>0 ,{j,1,NFPaux-1},{i,If[j==1,1,If[j==2 || j==3,2,If[j==4||j==5||j==6,3,If[j==7||j==8,4]]]],dimaux}]];
If[ToString[Reduce[Laux,DeleteDuplicates[Flatten[Table[Variables[x[[k]][[j]][[i]]],{j,1,NFPaux-1},{i,j,d}]]],Integers]]==ToString[False],droplistq=Append[droplistq,k]];
k=k+1])];

DropFunction=Function[{x,y},
(NewList={};
k=1;
While[k≤Length[x],
If[Not[MemberQ[y,k]],NewList=Append[NewList,x[[k]]]];
k=k+1])];

DropFunction2=Function[{x,y,q},
(NewList={};
k=1;
While[k≤Length[x],
If[Not[MemberQ[y,{q,k}]],NewList=Append[NewList,x[[k]]]];
k=k+1])];

(* 1. First we eliminate cases where we obtain rational weights after ellimination of common divisors. 
We store the possible weights in the files "NewFIWaux?.txt", where "?" is the number of the multigraph *)

(*PNull=ToExpression[Import[FileNameJoin[{NotebookDirectory[],"dim"<>ToString[dim]<>"nfp"<>ToString[NFP],"PNullSpaces.txt"}],"List"]];*)

p=1;
While[p≤Length[PNull],
k=1;
droplist4={};
q=PNull[[p]];
Print["Multigraph #: ",dets[[Matricessmallrk[[q]]]]]; 
FIWauxq=ToExpression[Import[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","FIWaux"<>ToString[dets[[Matricessmallrk[[q]]]]]<>".txt"}],"List"]];
TestRational[FIWauxq,cardI];
DropFunction[FIWauxq,droplist4];
aux3=NewList;
aux4=Simplify[aux3];
Export[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","NewFIWaux"<>ToString[dets[[Matricessmallrk[[q]]]]]<>".txt"}],aux4,"List"];
p=p+1];

(* 2. In the cases where we have multiple edges at the minimum we select possible cases according to Lemma 8.2 (1).
We save the resulting weights in the files "NNewFIWaux?.txt", where "?" is the number of the graph *)

(* PNull=ToExpression[Import[FileNameJoin[{NotebookDirectory[],"dim"<>ToString[dim]<>"nfp"<>ToString[NFP],"PNullSpaces.txt"}],"List"]]; *)

p=1;
While[p≤Length[PNull],
q=PNull[[p]];
Print["Multigraph #: ",dets[[Matricessmallrk[[q]]]]];
FIWauxq=ToExpression[Import[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","NewFIWaux"<>ToString[dets[[Matricessmallrk[[q]]]]]<>".txt"}],"List"]];
LLength=Length[FIWauxq];
droplist={};
LLL=Tally[Table[e[[Matricessmallrk[[q]]]][[i]][[2]],{i,1,d}]]; 
l=1;
While[l≤Length[LLL],
If[LLL[[l,2]]==d&& d≠ 1,
LLL2=Flatten[Position[Table[e[[Matricessmallrk[[q]]]][[i]][[2]],{i,1,d}],LLL[[l,1]]]];
k=1;
While[k≤LLength,
LLL3=Select[Tuples[LLL2,2],#[[1]]<#[[2]]&];
LLL4=Table[PolynomialGCD@@Table[{FIWauxq[[k]][[1]][[LLL3[[i]][[1]]]],FIWauxq[[k]][[1]][[LLL3[[i]][[2]]]]},{i,1,Length[LLL3]}][[i]],{i,1,Length[LLL3]}];
LLL5=DeleteCases[Cases[LLL4,_Integer],1];
If[Length[LLL5]>0,
droplist=Append[droplist,k],
cc=1;
While[cc≤Length[LLL4],
If[ToString[LLL4[[cc]]]≠ToString[1]&&Not[NumberQ[LLL4[[cc]]]],Laux=Solve[LLL4[[cc]]==1,Variables[LLL4[[cc]]][[1]]];
r=1;
aux2=FIWauxq[[k]];
While[r≤NFP,
j=1;
aux={};
While[j≤dim/2,
aux=Append[aux,Simplify[FIWauxq[[k]][[r]][[j]]/.Flatten[Laux]]];
j=j+1];
aux2=ReplacePart[aux2,r->aux];
r=r+1];
FIWauxq=Part[FIWauxq,k->aux2]];
cc=cc+1]];
k=k+1],
If[LLL[[l,2]]==d-1&& d-1≠ 1,
LLL2=Flatten[Position[Table[e[[Matricessmallrk[[q]]]][[i]][[2]],{i,1,d}],LLL[[l,1]]]];
k=1;
While[k≤Length[FIWauxq],
G=PolynomialGCD@@Table[FIWauxq[[k]][[1]][[LLL2[[i]]]],{i,1,Length[LLL2]}];
If[NumberQ[G]&&G>1,
droplist=Append[droplist,k],
If[ToString[G]≠ToString[1]&&Not[NumberQ[G]],
Laux=Solve[G==1,Variables[G][[1]]];
r=1;
aux2=FIWauxq[[k]];
While[r≤NFP,
j=1;
aux={};
While[j≤dim/2,
aux=Append[aux,Simplify[FIWauxq[[k]][[r]][[j]]/.Flatten[Laux]]];
j=j+1];
aux2=ReplacePart[aux2,r->aux];
r=r+1];
FIWauxq=ReplacePart[FIWauxq,k->aux2]]];
k=k+1]]];
l=l+1];

DropFunction[FIWauxq,droplist];
FIWauxIII=NewList;
Export[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","NNewFIWaux"<>ToString[dets[[Matricessmallrk[[q]]]]]<>".txt"}],FIWauxIII,"List"];
p=p+1];

(* 7. Now we select the cases with the right number and position of positive weights. The possible cases are stored in the files  "NewShortFIWaux?.txt" *)

PNull=ToExpression[Import[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","PNullSpaces.txt"}],"List"]]; 

p=1;
While[p≤Length[PNull],
q=PNull[[p]]; 
IWauxq=ToExpression[Import[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","NNewFIWaux"<>ToString[dets[[Matricessmallrk[[q]]]]]<>".txt"}],"List"]];
Print["Multigraph #: ",dets[[Matricessmallrk[[q]]]]];
TestPositive[IWauxq,cardI];
DropFunction[IWauxq,droplistq];
IWListq=NewList;
Export[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","NewShortFIWaux"<>ToString[dets[[Matricessmallrk[[q]]]]]<>".txt"}],IWListq,"List"];
p=p+1];

PNull=ToExpression[Import[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","PNullSpaces.txt"}],"List"]];

(* 9. Now we select the cases for which ∫c_j = 0 for j=0,...,dim/2-1 and store them in the files "FinalIWaux?.txt" *)

(* PNull=ToExpression[Import[FileNameJoin[{NotebookDirectory[],"dim"<>ToString[dim]<>"nfp"<>ToString[NFP],"PNullSpaces.txt"}],"List"]]; *)

droplistsmall={};
Bignumlist={};
Print[];
Print["Selecting the isotropy weights that give ∫c_j = 0 for j=0,...,dim/2-1..."];
p=1;
While[ p≤ Length[PNull], 
q=PNull[[p]]; 
IWauxq=ToExpression[Import[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","NewShortFIWaux"<>ToString[dets[[Matricessmallrk[[q]]]]]<>".txt"}],"List"]];
k=1;
Print[];
Length2=Length[IWauxq];
    NLength2=IntegerPart[Length2/10];
Print["Multigraph #: ",ToString[dets[[Matricessmallrk[[q]]]]]];
While[k≤ Length2,
Array[b,cardI];
loclistsmallaux={};
j=0;
While[j≤dim/2-1,
Locsmall=0;
i=1;
While[i≤NFP,Locsmall=Locsmall+SymmetricPolynomial[j,IWauxq[[k]][[i]]]/SymmetricPolynomial[dim/2,IWauxq[[k]][[i]]];
i=i+1];
SLocsmall=Together[Locsmall];
loclistsmallaux=Append[loclistsmallaux,SLocsmall];
j=j+1];
NSLocsmall=Numerator/@loclistsmallaux;
If[Not[MemberQ[Bignumlist,NSLocsmall]],

VVsmall=DeleteDuplicates[Flatten[Variables/@NSLocsmall]];
solsmall=Solve[NSLocsmall==0];
If[Length[solsmall]==0,droplistsmall=Append[droplistsmall,{q,k}];Bignumlist=Append[Bignumlist,NSLocsmall],

If[MemberQ[Table[Select[(VVsmall/.solsmall)[[i]],#==0||(Variables[#]=={}&&Not[Block[{$MaxExtraPrecision=1000},IntegerQ[#]]])||(Variables[#]=={}&&Not[Block[{$MaxExtraPrecision=1000},Positive[#]]]) &],{i,1,Length[VVsmall/.solsmall]}],{}]==False,droplistsmall=Append[droplistsmall,{q,k}];Bignumlist=Append[Bignumlist,NSLocsmall],

Listaux=Position[Table[Select[(VVsmall/.solsmall)[[i]],#==0||(Variables[#]=={}&&Not[IntegerQ[#]])||(Variables[#]=={}&&Not[Positive[#]]) &],{i,1,Length[VVsmall/.solsmall]}],{}];

If[Length[Listaux]==1,

aux2=IWauxq[[k]];
l=1;
While[l≤NFP,
j=1;
aux={};
While[j≤dim/2,
aux=Append[aux,Simplify[IWauxq[[k]][[l]][[j]]/.Flatten[solsmall[[Listaux[[1]]]]]]];
j=j+1];
aux2=ReplacePart[aux2,l->aux];
l=l+1];
IWauxq=ReplacePart[IWauxq,k->aux2]]];
Clear[solsmall]],
droplistsmall=Append[droplistsmall,{q,k}]];

If[NLength2≠0,If[Divisible[k,NLength2],If[k/NLength2≤10,Print[k/NLength2*10 "%"]]]];
k=k+1];
p=p+1];


IWFinal=Table[{},{i,1,Length[Matricessmallrk]}];
p=1;
While[p≤Length[PNull],
q=PNull[[p]];
IWauxq=ToExpression[Import[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","NewShortFIWaux"<>ToString[dets[[Matricessmallrk[[q]]]]]<>".txt"}],"List"]];
DropFunction2[IWauxq,droplistsmall,q];
IWaux=NewList;
Export[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","FinalIWaux"<>ToString[dets[[Matricessmallrk[[q]]]]]<>".txt"}],IWaux,"List"];
p=p+1];

(* 10. Delete duplicates: *)

(* PNull=ToExpression[Import[FileNameJoin[{NotebookDirectory[],"dim"<>ToString[dim]<>"nfp"<>ToString[NFP],"PNullSpaces.txt"}],"List"]]; *)

IWFinal=Table[{},{i,1,Length[PNull]}];
p=1;
While[p≤Length[PNull],
q=PNull[[p]];
IWauxq=ToExpression[Import[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","FinalIWaux"<>ToString[dets[[Matricessmallrk[[q]]]]]<>".txt"}],"List"]];
FIWaux3=DeleteDuplicates[Table[Sort[Expand[IWauxq[[i,j]]]],{i,1,Length[IWauxq]},{j,1,NFP}]];
IWFinal[[p]]=FIWaux3;
Export[FileNameJoin[{NotebookDirectory[],"graphsb2_2b4_2","LastIWaux"<>ToString[dets[[Matricessmallrk[[q]]]]]<>".txt"}],FIWaux3,"List"];
p=p+1];

(* Part VI - Conclusion - Finally we print all possible isotropy weights for each multigraph*)

Print[];
Print[];
Print["**************************************************************************"];
Print[];
Print["Conclusion:"];
Print[];
p=1;
FinalGraphs={};
While[p≤Length[PNull],
Print["**************************************************************************"];
Print["All Possible Isotropy Weights for Multigraph #"<>ToString[dets[[Matricessmallrk[[PNull[[p]]]]]]]<>":"];
Print[];
If[IWFinal[[p]]≠ {}, NS=Length[IWFinal[[p]]];
Print["There are at most "<>ToString[NS]<>" solutions:"]];

(*If[IWFinal[[p]]=={},Print["There are no possible isotropy weights" ]];*)
If[IWFinal[[p]]≠ {},k=1;While[k≤ Length[IWFinal[[p]]], Print[IWFinal[[p]][[k]]];Print[]; k=k+1]];

If[ IWFinal[[p]]≠ {},
lgg=GraphPlot[lff[[dets[[Matricessmallrk[[PNull[[p]]]]]]]],DirectedEdges->True,VertexLabeling->True,MultiedgeStyle->True,SelfLoopStyle->True,VertexCoordinateRules->{{(1/2),1},{0.2,2},{0.8,2},{-(1/2),3},{(1/2),3},{(3/2),3},{0.2,4},{0.8,4},{(1/2),5}},ImageSize->Large];
 Print[lgg];AppendTo[FinalGraphs,{p,lff[[dets[[Matricessmallrk[[PNull[[p]]]]]]]]}];Print[lff[[dets[[Matricessmallrk[[PNull[[p]]]]]]]]]];
p=p+1];
Print[FinalGraphs];

(* check if interchanging fixed points of index 2 works fine *)

FinalList={}; 
listp={};
listp2={};
p=1;
While[p≤  Length[FinalGraphs], 
Baux=FinalGraphs[[p]][[2]]; 
Daux=Baux;
If[Baux[[4]][[4]]==Baux[[4]][[6]]&&Baux[[4]][[6]]==Baux[[6]][[4]]&&Baux[[4]][[6]]==Baux[[6]][[6]],
  j=1;
While[j≤NFP,
If[j≠4&& j≠6, Daux[[4]][[j]] =Baux[[6]][[j]]; Daux[[6]][[j]] =Baux[[4]][[j]]]; j=j+1];
i=1;
While[i≤NFP,
If[i≠4&& i≠6,Daux[[i]][[4]] =Baux[[i]][[6]]; Daux[[i]][[6]] =Baux[[i]][[4]]]; i=i+1];,
If[Baux[[5]][[5]]==Baux[[5]][[6]]&&Baux[[6]][[5]]==Baux[[6]][[6]]&&Baux[[5]][[6]]==Baux[[6]][[5]],
  j=1;
While[j≤NFP,
If[j≠5&& j≠6, Daux[[5]][[j]] =Baux[[6]][[j]]; Daux[[6]][[j]] =Baux[[5]][[j]]]; j=j+1];
i=1;While[i≤NFP,
If[i≠5&& i≠6,Daux[[i]][[5]] =Baux[[i]][[6]]; Daux[[i]][[6]] =Baux[[i]][[5]]]; i=i+1], 
If[Baux[[1]][[1]]==Baux[[1]][[2]]&&Baux[[1]][[2]]==Baux[[2]][[1]]&&Baux[[2]][[1]]==Baux[[2]][[2]],
  j=1;
While[j≤NFP,
If[j≠1&& j≠1, Daux[[1]][[j]] =Baux[[2]][[j]]; Daux[[2]][[j]] =Baux[[1]][[j]]]; j=j+1];
i=1;While[i≤NFP,
If[i≠1&& i≠2,Daux[[i]][[1]] =Baux[[i]][[2]]; Daux[[i]][[2]] =Baux[[i]][[1]]]; i=i+1]] 
]];
If[Not[MemberQ[FinalList,Daux]],FinalList=Append[FinalList,Baux];listp=Append[listp,FinalGraphs[[p]][[1]]];listp2=Append[listp2,{FinalGraphs[[p]][[1]],Baux}]];
p=p+1];Print[Length[FinalList]];
Print[listp];Print[Length[FinalGraphs]];

(* check if interchanging fixed points of index 2 works fine *)

FinalList2={}; 
listq={};
listq2={};
p=1;
While[p≤  Length[FinalList], 
Baux=listp2[[p]][[2]]; 
Daux=Baux;
If[Baux[[5]][[5]]==Baux[[5]][[6]]&&Baux[[6]][[5]]==Baux[[6]][[6]]&&Baux[[5]][[6]]==Baux[[6]][[5]],
  j=1;
While[j≤NFP,
If[j≠5&& j≠6, Daux[[5]][[j]] =Baux[[6]][[j]]; Daux[[6]][[j]] =Baux[[5]][[j]]]; j=j+1];
i=1;While[i≤NFP,
If[i≠5&& i≠6,Daux[[i]][[5]] =Baux[[i]][[6]]; Daux[[i]][[6]] =Baux[[i]][[5]]]; i=i+1], 
If[Baux[[1]][[1]]==Baux[[1]][[2]]&&Baux[[1]][[2]]==Baux[[2]][[1]]&&Baux[[2]][[1]]==Baux[[2]][[2]],
  j=1;
While[j≤NFP,
If[j≠1&& j≠1, Daux[[1]][[j]] =Baux[[2]][[j]]; Daux[[2]][[j]] =Baux[[1]][[j]]]; j=j+1];
i=1;While[i≤NFP,
If[i≠1&& i≠2,Daux[[i]][[1]] =Baux[[i]][[2]]; Daux[[i]][[2]] =Baux[[i]][[1]]]; i=i+1]] 
];
If[Not[MemberQ[FinalList2,Daux]],FinalList2=Append[FinalList2,Baux];listq=Append[listq,listp2[[p]][[1]]];listq2=Append[listq2,{listp2[[p]][[1]],Baux}]];
p=p+1];Print[Length[FinalList2]];
Print[listq];

Print[Length[FinalList]];Print[Length[FinalGraphs]];

(* check if interchanging fixed points of index 2 works fine *)

FinalList3={}; 
listr={};
listr2={};
p=1;
While[p≤  Length[FinalList2], 
Baux=listq2[[p]][[2]]; 
Daux=Baux;

If[Baux[[7]][[8]]==Baux[[8]][[7]]&&Baux[[7]][[8]]==Baux[[7]][[7]]&&Baux[[7]][[8]]==Baux[[8]][[8]],
  j=1;
While[j≤NFP,
If[j≠7&& j≠8, Daux[[7]][[j]] =Baux[[8]][[j]]; Daux[[8]][[j]] =Baux[[7]][[j]]]; j=j+1];
i=1;While[i≤NFP,
If[i≠7&& i≠8,Daux[[i]][[7]] =Baux[[i]][[8]]; Daux[[i]][[8]] =Baux[[i]][[7]]]; i=i+1]
];
If[Not[MemberQ[FinalList3,Daux]],FinalList3=Append[FinalList3,Baux];listr=Append[listr,listq2[[p]][[1]]];listr2=Append[listr2,{listq2[[p]][[1]],Baux}]];
p=p+1];Print[Length[FinalList3]];
Print[listr];
Print[Length[FinalList2]];
Print[Length[FinalList]];Print[Length[FinalGraphs]];

(* check if interchanging fixed points of index 2 works fine *)

FinalList4={}; 
lists={};
lists2={};
p=1;
While[p≤  Length[FinalList3], 
Baux=listr2[[p]][[2]]; 
Daux=Baux;

If[Baux[[4]][[5]]==Baux[[5]][[4]]&&Baux[[4]][[5]]==Baux[[5]][[5]]&&Baux[[4]][[5]]==Baux[[4]][[4]],
  j=1;
While[j≤NFP,
If[j≠4&& j≠5, Daux[[4]][[j]] =Baux[[5]][[j]]; Daux[[5]][[j]] =Baux[[4]][[j]]]; j=j+1];
i=1;While[i≤NFP,
If[i≠4&& i≠5,Daux[[i]][[4]] =Baux[[i]][[5]]; Daux[[i]][[5]] =Baux[[i]][[4]]]; i=i+1]
];
If[Not[MemberQ[FinalList4,Daux]],FinalList4=Append[FinalList4,Baux];lists=Append[lists,listr2[[p]][[1]]];lists2=Append[lists2,{listr2[[p]][[1]],Baux}]];
p=p+1];Print[Length[FinalList4]];
Print[lists];
Print[Length[FinalList3]];
Print[Length[FinalList2]];
Print[Length[FinalList]];Print[Length[FinalGraphs]];

(* Part VI - Conclusion - Finally we print all possible isotropy weights for each multigraph*)

Print[];
Print[];
Print["**************************************************************************"];
Print[];
Print["Conclusion:"];
Print[];
t=1;
FinalGraphs2={};
While[t≤Length[lists],
p=lists[[t]];
Print["**************************************************************************"];
Print["All Possible Isotropy Weights for Multigraph #: ",dets[[Matricessmallrk[[PNull[[p]]]]]]];
Print[];
If[IWFinal[[p]]≠ {}, NS=Length[IWFinal[[p]]];
Print["There are at most "<>ToString[NS]<>" solutions:"]];

If[IWFinal[[p]]=={},Print["There are no possible isotropy weights" ]];
If[IWFinal[[p]]≠ {},k=1;While[k≤ Length[IWFinal[[p]]], Print[IWFinal[[p]][[k]]];Print[];

Locsmall=0;
Locsmall2c12c2=0;
Locsmallc22=0;
i=1;
While[i≤NFP,Locsmall=Locsmall+(SymmetricPolynomial[1,IWFinal[[p]][[k]][[i]]])^4/SymmetricPolynomial[dim/2,IWFinal[[p]][[k]][[i]]];
Locsmall2c12c2=Locsmall2c12c2+(SymmetricPolynomial[1,IWFinal[[p]][[k]][[i]]])^2*SymmetricPolynomial[2,IWFinal[[p]][[k]][[i]]]/SymmetricPolynomial[dim/2,IWFinal[[p]][[k]][[i]]];Locsmallc22=Locsmallc22+(SymmetricPolynomial[2,IWFinal[[p]][[k]][[i]]])^2/SymmetricPolynomial[dim/2,IWFinal[[p]][[k]][[i]]];
i=i+1];
loclistsmallaux=Together[Locsmall];
loclistsmallaux2=Together[Locsmall2c12c2];
loclistsmallaux3=Together[Locsmallc22];
Print["c1^4= ",loclistsmallaux];
Print["c1^2c2= ",loclistsmallaux2];
Print["c2^2= ",loclistsmallaux3];
Print[];
k=k+1]];

If[ IWFinal[[p]]≠ {},
lgg=GraphPlot[lff[[dets[[Matricessmallrk[[PNull[[p]]]]]]]],DirectedEdges->True,VertexLabeling->True,MultiedgeStyle->True,SelfLoopStyle->True,VertexCoordinateRules->{{(1/2),1},{0.2,2},{0.8,2},{-(1/2),3},{(1/2),3},{(3/2),3},{0.2,4},{0.8,4},{(1/2),5}},ImageSize->Large];
 Print[lgg];AppendTo[FinalGraphs2,{p,lff[[dets[[Matricessmallrk[[PNull[[p]]]]]]]]}];Print[lff[[dets[[Matricessmallrk[[PNull[[p]]]]]]]]]];
t=t+1];
Print[FinalGraphs2];

Simplify[-((2 (-51 b[1]-50 ((b[1])^2)+2 ((b[1])^3)+((b[1])^4)))/(b[1] (1+b[1])))]

Solve[-((2 (-1-51 b[1]-50 ((b[1])^2)+2 ((b[1])^3)+((b[1])^4)))/(b[1] (1+b[1])))==99, Integers]

-((2 (-1-51 b[2]-50 ((b[2])^2)+2 ((b[2])^3)+((b[2])^4)))/(b[2] (1+b[2])))

(* Part VI - Conclusion - Finally we print all possible isotropy weights for each multigraph*)

Print["**************************************************************************"];
Print["Conclusion:"];
Print["**************************************************************************"];
p=1;
While[p≤Length[lff],

If[FIW[[p]]≠{}|| IWFinal[[p]]≠ {}, NS=Length[IWFinal[[p]]]+Length[FIW[[p]]];
Print[];
Print["**************************************************************************"];
Print["All Possible Isotropy Weights for Multigraph #"<>ToString[p]<>":"];
Print["There are at most "<>ToString[NS]<>" solutions:"]];
If[FIW[[p]]≠{}, k=1;While[k≤ Length[FIW[[p]]], Print[FIW[[p]][[k]]];Print[]; 
Locsmall=0;
Locsmall2c12c2=0;
Locsmallc22=0;
i=1;
While[i≤NFP,
Locsmall=Locsmall+(SymmetricPolynomial[1,FIW[[p]][[k]][[i]]])^4/SymmetricPolynomial[dim/2,FIW[[p]][[k]][[i]]];
Locsmall2c12c2=Locsmall2c12c2+(SymmetricPolynomial[1,FIW[[p]][[k]][[i]]])^2*SymmetricPolynomial[2,FIW[[p]][[k]][[i]]]/SymmetricPolynomial[dim/2,FIW[[p]][[k]][[i]]];Locsmallc22=Locsmallc22+(SymmetricPolynomial[2,FIW[[p]][[k]][[i]]])^2/SymmetricPolynomial[dim/2,FIW[[p]][[k]][[i]]];
i=i+1];
loclistsmallaux=Together[Locsmall];
loclistsmallaux2=Together[Locsmall2c12c2];
loclistsmallaux3=Together[Locsmallc22];
Print["c1^4= ",loclistsmallaux];
Print["c1^2c2= ",loclistsmallaux2];
Print["c2^2= ",loclistsmallaux3];
Print[];
k=k+1]];

If[IWFinal[[p]]=={}&& FIW[[p]]=={},Print["There are no possible isotropy weights for multigraph # "<>ToString[p]]];

If[IWFinal[[p]]≠ {},k=1;While[k≤ Length[IWFinal[[p]]], Print[IWFinal[[p]][[k]]];Print[]; 
Locsmall=0;
Locsmall2c12c2=0;
Locsmallc22=0;
i=1;
While[i≤NFP,Locsmall=Locsmall+(SymmetricPolynomial[1,IWFinal[[p]][[k]][[i]]])^4/SymmetricPolynomial[dim/2,IWFinal[[p]][[k]][[i]]];
Locsmall2c12c2=Locsmall2c12c2+(SymmetricPolynomial[1,IWFinal[[p]][[k]][[i]]])^2*SymmetricPolynomial[2,IWFinal[[p]][[k]][[i]]]/SymmetricPolynomial[dim/2,IWFinal[[p]][[k]][[i]]];Locsmallc22=Locsmallc22+(SymmetricPolynomial[2,IWFinal[[p]][[k]][[i]]])^2/SymmetricPolynomial[dim/2,IWFinal[[p]][[k]][[i]]];
i=i+1];
loclistsmallaux=Together[Locsmall];
loclistsmallaux2=Together[Locsmall2c12c2];
loclistsmallaux3=Together[Locsmallc22];
Print["c1^4= ",loclistsmallaux];
Print["c1^2c2= ",loclistsmallaux2];
Print["c2^2= ",loclistsmallaux3];
Print[];
k=k+1]];
If[FIW[[p]]≠{}|| IWFinal[[p]]≠ {}, Print[lgg[[p]]];Print["**************************************************************************"];Print[]];
p=p+1];

IWFinal

(* Part VI - Conclusion - Finally we print all possible isotropy weights for each multigraph*)

Print[];
Print[];
Print["**************************************************************************"];
Print[];
Print["Conclusion:"];
Print[];
p=1;
While[p≤Length[lff],
Print["**************************************************************************"];
Print["All Possible Isotropy Weights for Multigraph #"<>ToString[p]<>":"];
Print[];


If[FIW[[p]]≠{}|| IWFinal[[p]]≠ {}, NS=Length[IWFinal[[p]]]+Length[FIW[[p]]];
Print["There are at most "<>ToString[NS]<>" solutions:"]];
If[FIW[[p]]≠{},
Print[FIW[[p]]];
Print[]];
If[IWFinal[[p]]=={}&& FIW[[p]]=={},Print["There are no possible isotropy weights" ]];
If[IWFinal[[p]]≠ {},k=1;While[k≤ Length[IWFinal[[p]]], Print[IWFinal[[p]][[k]]];Print[]; k=k+1]];
If[FIW[[p]]≠{}|| IWFinal[[p]]≠ {}, Print[lgg[[p]]]];
p=p+1];
(*******************************)

(* Appendix: Now we select the cases for which ∫c_1^4 = 512, 288 or 400  independetly of the weights and store them in the files "SecondReallyLastIWaux?.txt" *)

PNull=ToExpression[Import[FileNameJoin[{NotebookDirectory[],"dim"<>ToString[dim]<>"nfp"<>ToString[NFP],"PNullSpaces.txt"}],"List"]];

droplistsmall={};
Bignumlist={};
Print[];
Print["Selecting the isotropy weights that give ∫c_j = 0 for j=0,...,dim/2-1..."];
p=1;
While[ p≤ Length[PNull], 
q=PNull[[p]]; 
IWauxq=ToExpression[Import[FileNameJoin[{NotebookDirectory[],"dim"<>ToString[dim]<>"nfp"<>ToString[NFP],"NewShortFIWaux"<>ToString[q]<>".txt"}],"List"]];
k=1;
Print[];
Length2=Length[IWauxq];
    NLength2=IntegerPart[Length2/10];
Print["Multigraph #"<>ToString[q]<>":"];
While[k≤ Length2,
Array[b,cardI];
loclistsmallaux={};
j=0;
While[j≤dim/2-1,
Locsmall=0;
i=1;
While[i≤NFP,Locsmall=Locsmall+SymmetricPolynomial[j,IWauxq[[k]][[i]]]/SymmetricPolynomial[dim/2,IWauxq[[k]][[i]]];
i=i+1];
SLocsmall=Together[Locsmall];
loclistsmallaux=Append[loclistsmallaux,SLocsmall];
j=j+1];
NSLocsmall=Numerator/@loclistsmallaux;
If[Not[MemberQ[Bignumlist,NSLocsmall]],

VVsmall=DeleteDuplicates[Flatten[Variables/@NSLocsmall]];
solsmall=Solve[NSLocsmall==0];
If[Length[solsmall]==0,droplistsmall=Append[droplistsmall,{q,k}];Bignumlist=Append[Bignumlist,NSLocsmall],

If[MemberQ[Table[Select[(VVsmall/.solsmall)[[i]],#==0||(Variables[#]=={}&&Not[Block[{$MaxExtraPrecision=1000},IntegerQ[#]]])||(Variables[#]=={}&&Not[Block[{$MaxExtraPrecision=1000},Positive[#]]]) &],{i,1,Length[VVsmall/.solsmall]}],{}]==False,droplistsmall=Append[droplistsmall,{q,k}];Bignumlist=Append[Bignumlist,NSLocsmall],

Listaux=Position[Table[Select[(VVsmall/.solsmall)[[i]],#==0||(Variables[#]=={}&&Not[IntegerQ[#]])||(Variables[#]=={}&&Not[Positive[#]]) &],{i,1,Length[VVsmall/.solsmall]}],{}];

If[Length[Listaux]==1,

aux2=IWauxq[[k]];
l=1;
While[l≤NFP,
j=1;
aux={};
While[j≤dim/2,
aux=Append[aux,Simplify[IWauxq[[k]][[l]][[j]]/.Flatten[solsmall[[Listaux[[1]]]]]]];
j=j+1];
aux2=ReplacePart[aux2,l->aux];
l=l+1];
IWauxq=ReplacePart[IWauxq,k->aux2]]];
Clear[solsmall]],
droplistsmall=Append[droplistsmall,{q,k}]];

If[NLength2≠0,If[Divisible[k,NLength2],If[k/NLength2≤10,Print[k/NLength2*10 "%"]]]];
k=k+1];
p=p+1];


IWFinal=Table[{},{i,1,Length[lff]}];
p=1;
While[p≤Length[PNull],
q=PNull[[p]];
IWauxq=ToExpression[Import[FileNameJoin[{NotebookDirectory[],"dim"<>ToString[dim]<>"nfp"<>ToString[NFP],"NewShortFIWaux"<>ToString[q]<>".txt"}],"List"]];
DropFunction2[IWauxq,droplistsmall,q];
IWaux=NewList;
Export[FileNameJoin[{NotebookDirectory[],"dim"<>ToString[dim]<>"nfp"<>ToString[NFP],"FinalIWaux"<>ToString[q]<>".txt"}],IWaux,"List"];
p=p+1];

Table[N[4(2i+96+2Sqrt[i(i+96)])],{i,2,6}]

Table[N[(2(Sqrt[(i+96)]+Sqrt[i]))^2],{i,2,6}]

N[629/4]

157*4

Table[N[(2(Sqrt[(i+96)]-Sqrt[i]))^2],{i,2,6}]

FI={{1,1,2,2},{-1,1,1,2},{-1,1,1,2},{-2,-1,1,2},{-1,-1,1,1},{-2,-1,1,2},{-2,-1,-1,1},{-2,-1,-1,1},{-2,-2,-1,-1}}

k=1;

k=1;
 FI={{{1,2,2,3},{-3,2,1,4},{-1,-2,1,2},{-1 ,-2,-2,3},{1,2 ,-2,-3},{-1,-1,1,2},{-1,1,-2,-3},{-1,-1,-2,3},{-2 ,-1,-2,-3}} };
(* FI={{{1,1,2,2},{-1,1,1,2},{-1,1,1,2},{-2,-1,1,2},{-1,-1,1,1},{-2,-1,1,2},{-2,-1,-1,1},{-2,-1,-1,1},{-2,-2,-1,-1}}}; *)
Locsmall2c12c2=0;
Locsmall=0;
Locsmallc22=0;
i=1;
While[i≤NFP,
Locsmall=Locsmall+(SymmetricPolynomial[1,FI[[k]][[i]]])^4/SymmetricPolynomial[dim/2,FI[[k]][[i]]];
Locsmall2c12c2=Locsmall2c12c2+(SymmetricPolynomial[1,FI[[k]][[i]]])^2*SymmetricPolynomial[2,FI[[k]][[i]]]/SymmetricPolynomial[dim/2,FI[[k]][[i]]];Locsmallc22=Locsmallc22+(SymmetricPolynomial[2,FI[[k]][[i]]])^2/SymmetricPolynomial[dim/2,FI[[k]][[i]]];
i=i+1];
loclistsmallaux=Together[Locsmall];
loclistsmallaux2=Together[Locsmall2c12c2];
loclistsmallaux3=Together[Locsmallc22];
Print["c1^4= ",loclistsmallaux];
Print["c1^2c2= ",loclistsmallaux2];
Print["c2^2= ",loclistsmallaux3];
Print[];

NFP

dim

i

Locsmall