(* ******************************************************************************************************)
(*                                                                                                      *)
(*                    This program is based on the algorithms presented in the paper                 *)
(*                                                                                                      *)
(*            New tools for classifying Hamiltonian circle actions with isolated fixed points        *)
(*                                                                                                      *)
(*                             by Leonor Godinho and Silvia Sabatini                                  *)
(*                                                                                                      *)
(*                                                                                                      *)
(*                                                                                                      *) 
(*  It determines all possible multisets of isotropy weights of  Hamiltonian circle actions with a  *)
(*  minimal number of fixed points that satisfy Property ((SuperPlus[(Subscript[P, 0])])) i.e. circle actions for which there  *) 
(*  exists a non-negative multigraph. All the references made are from the above paper.             *)
(*                                                                                                      *)
(* **************************************************************************************************** *)




(* Part I- Preparation for Non Negative Multigraphs *)


NotebookDirectory;
SetDirectory[NotebookDirectory[]];
Needs["GraphUtilities`"];
Needs["CCompilerDriver`"];

(* First we request the dimension of the manifold *)

dim=Input["This program determines all possible isotropy weights for Hamiltonian circle actions with a minimal number of fixed points and a non-negative multigraph, according to the paper New tools for classifying Hamiltonian circle actions with isolated fixed points  by Leonor Godinho and Silvia Sabatini. Please insert the dimension of the manifold:"];

(* Then we compute the number of fixed points, the value of the integral ∫c_1c_{n-1} and the number of edges (cardI) *)

d=dim/2;
NFP=d+1;
integral=d*NFP^2/2;
cardI=NFP*dim/4;

(* Possible graphs *)

(* First we request the vertex coordinates and we store them in the array "vertexcoordinates"(eg. in the 8-dim'l case with 5 fixed points these can be {0,-3},{1,-2},{2,0},{1,2},{0,3}) *)

k=1;
vertexcoordinates=Table[0,{i,1,NFP}];
While[k≤NFP,vertexcoordinates[[k]]=k-> Input["Please insert coordinates for vertex # "<>ToString[k]<>" in the form {...,...}"];k=k+1];

(* Then we construct the list "lf" of all possible NNMs - non-negative mutigraphs *)

li=Select[Tuples[Table[Flatten[Permutations/@IntegerPartitions[d-i,{d+1},Table[j,{j,0,d-i}]],1],{i,0,d}]],UpperTriangularize[#]==# &];
lf=Select[li,Total[#]==Table[i,{i,0,d}]&];

(* Finally, we delete unnecessary repetitions obtained by reversing the circle action, getting the final list of NNMs "lff" *)

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

(* We then print the final list of NNMs *) 

lgg=Table[GraphPlot[lff[[k]],DirectedEdges->True,VertexLabeling->True,MultiedgeStyle->True,SelfLoopStyle->True,VertexCoordinateRules->vertexcoordinates],{k,1,Length[lff]}];
Print["**************************************************************************"];
Print["All possible non-negative multigraphs (including those with cycles):"];
Print[lgg];

(* Auxiliary list of  multigraphs to determine those with multiple connected components whose indices are stored in the array "MultipleCompPos" *)

auxlff=Table[If[lff[[p]][[i]][[j]]≠0,1,0],{p,1,Length[lff]},{i,1,NFP},{j,1,NFP}];
AuxComponents=Table[WeakComponents[AdjacencyGraph[auxlff[[p]]]],{p,1,Length[lff]}];
MultipleCompPos=Cases[Table[If[Length[AuxComponents[[p]]]>1,p,0],{p,1,Length[lff]}],Except[0]];

(* Matrices A (see Section 4)*)

(* Edge function: e[[p]][[i]] is the i-th edge of the p-th multigraph of "lff"; e[[p]][[i]][[1]] is its initial point and e[[p]][[i]][[2]] is its terminal point (cf. Section 4.) *)
  
e=Array[0&,{Length[lff],cardI,2}];
p=1;
While[p≤ Length[lff], 
j=1;
While[j≤ NFP, 
r=1;
While[r≤ NFP, 
entry=lff[[p]][[j]][[r]]; 
If[entry≠ 0, 
k=1;
While[ k≤ entry, index=k+1/2(j-1)(2*NFP-j)+ Total[Table[lff[[p]][[j]][[m]],{m,1,r-1}]]; e[[p]][[index]][[1]]=j;e[[p]][[index]][[2]]=r; k=k+1]];
r=r+1];
j=j+1];
p=p+1];

(* Functions "deltai" and "deltat" for each multigraph of "lff" (vf. Section 4.): deltai[[p]][[h,m]]=delta(i(e_h),e_m) and  deltat[[p]][[h,m]]=delta(t(e_h),e_m), where e_h and e_m are respectively the k-th and m-th edges of the p-th multigraph of "lff" *)
  
deltai=Array[0&,{Length[lff],cardI,cardI}];
deltat=Array[0&,{Length[lff],cardI,cardI}];

(* Matrices A: A[[p]] is the matrix A for the p-th multigraph of "lff" *)

A=Array[0&,{Length[lff],cardI,cardI}];
p=1;
While[p≤ Length[lff],
h=1; 
While[h≤ cardI, 
m=1;
While[m≤ cardI, 
If[ e[[p]][[h]][[1]]==e[[p]][[m]][[1]] && e[[p]][[h]][[1]]≠ e[[p]][[m]][[2]] , deltai[[p]][[h]][[m]]=1];
If[ e[[p]][[h]][[1]]==e[[p]][[m]][[2]] && e[[p]][[h]][[1]]≠ e[[p]][[m]][[1]] , deltai[[p]][[h]][[m]]=-1];
If[ e[[p]][[h]][[2]]==e[[p]][[m]][[1]] && e[[p]][[h]][[2]]≠ e[[p]][[m]][[2]] , deltat[[p]][[h]][[m]]=1];
If[ e[[p]][[h]][[2]]==e[[p]][[m]][[2]] && e[[p]][[h]][[2]]≠ e[[p]][[m]][[1]] , deltat[[p]][[h]][[m]]=-1];
A[[p]][[h]][[m]]= deltai[[p]][[h]][[m]]-deltat[[p]][[h]][[m]];
      m=m+1];
h=h+1];
p=p+1];

(* 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["**************************************************************************"];
Print["List of  determinants to be used by the c++ programs:"];
NCArray=Table[0,{i,1,Length[lff]}];
auxedges=Table[0,{i,1,Length[lff]}];
k=1;
While[k≤ Length[lff],
Print["Multigraph #"<>ToString[k]<>":"];
dd=Diagonal[lff[[k]]];
nc=Total[dd];
If[nc≠0,
ap=Delete[Table[i,{i,1,NFP}],Position[dd,0]];
auxedges[[k]]=Partition[Flatten[Table[(2 d+2-ap[[i]]) (ap[[i]]-1)/2+mm,{i,1,Length[ap]},{mm,1,dd[[ap[[i]]]]}]],1];
NCArray[[k]]=nc;
Lae=Length[auxedges[[k]]];
If[MemberQ[MultipleCompPos,k],acomp=AuxComponents[[k]];
p=1;
blocks={};
While[p≤Length[acomp],i=1;
auxblocks={};
While[i≤cardI,If[MemberQ[acomp[[p]],e[[k]][[i]][[1]]],auxblocks=Append[auxblocks,i]];
i=i+1];
blocks=Append[blocks,auxblocks];
p=p+1];
dblocks=Table[Select[blocks[[i]],Not[MemberQ[Flatten[auxedges[[k]]],#]]&],{i,1,Length[blocks]}];
dblocks=Cases[dblocks,Except[{}]];
NMA=Table[A[[k]][[dblocks[[p]][[i]]]][[dblocks[[p]][[j]]]],{p,1,Length[dblocks]},{i,1,Length[dblocks[[p]]]},{j,1,Length[dblocks[[p]]]}];
p=1;
ndet=Table[0,{i,1,Length[NMA]}];
While[p≤Length[NMA],L=Array[l,Length[NMA[[p]]]];
ndet[[p]]=Simplify[Det[NMA[[p]]-DiagonalMatrix[L]]];
calcdet[L_]:=Do[ndet[[p]]=Replace[ndet[[p]],l[j]->L[[j]],{1,100}],{j,Length[NMA[[p]]]}];
calcdet[ToExpression[ToString[Table[Row[{"x",dblocks[[p]][[i]]-Length[Select[Flatten[auxedges[[k]]],#<dblocks[[p]][[i]]&]]}],{i,1,Length[dblocks[[p]]]}]]]];
p=p+1];
SumNdet=If[Length[ndet]>1,Sum[ndet[[r]]^2,{r,1,Length[ndet]}],ndet[[1]]];
Print[CForm[SumNdet]],delete=DeleteDuplicates[Join[Flatten[Table[{i,Flatten[auxedges[[k]]][[j]]},{i,1,cardI},{j,1,Lae}],1],Flatten[Table[{Flatten[auxedges[[k]]][[j]],i},{i,1,cardI},{j,1,Lae}],1]]];
NMA=Cases[Delete[A[[k]],delete],Except[{}]];
nne=cardI-nc;L=Array[l,nne];
ndet=Simplify[Det[NMA-DiagonalMatrix[L]]];
calcdet[L_]:=Do[ndet=Replace[ndet,l[j]->L[[j]],{1,100}],{j,nne}];
auxdet=ndet;
calcdet[ToExpression[ToString[Table[Row[{"x",i}],{i,1,nne}]]]];
Print[CForm[ndet]]];
Print[{"# of cycles="<>ToString[nc],"Connected Components="<>ToString[dblocks]}];Print[],

If[MemberQ[MultipleCompPos,k],
acomp=AuxComponents[[k]];
p=1;
blocks={};
While[p≤ Length[acomp], 
i=1;
auxblocks={};
While[i≤ cardI, 
If[MemberQ[acomp[[p]],e[[k]][[i]][[1]]],auxblocks=Append[auxblocks,i]];
i=i+1];
blocks=Append[blocks,auxblocks];
p=p+1];    
NMA=Table[A[[k]][[blocks[[p]][[i]]]][[blocks[[p]][[j]]]],{p,1,Length[blocks]},{i,1,Length[blocks[[p]]]},{j,1,Length[blocks[[p]]]}];

p=1; 
ndet=Table[0,{i,1,Length[NMA]}];
While[p≤ Length[NMA],
L=Array[l,Length[NMA[[p]]]];
ndet[[p]]=Simplify[Det[NMA[[p]]-DiagonalMatrix[L]]];
calcdet[L_]:= Do[ndet[[p]]=Replace[ndet[[p]],l[j]->L[[j]],{1,100}],{j,Length[NMA[[p]]]}];
calcdet[ToExpression[ToString[Table[Row[{"x",blocks[[p]][[i]]}],{i,1,Length[blocks[[p]]]}]]]];
p=p+1];

SumNdet=If[Length[ndet]>1,Sum[ndet[[r]]^2,{r,1,Length[ndet]}],ndet[[1]]];
Print[CForm[SumNdet]];
Print["Connected Components="<>ToString[blocks]],
     
L=Array[l,cardI];
ndet=Simplify[Det[A[[k]]-DiagonalMatrix[L]]];
calcdet[L_]:= Do[ndet=Replace[ndet,l[j]->L[[j]],{1,100}],{j,cardI}];  
calcdet[ToExpression[ToString[Table[Row[{"x",i}],{i,1,cardI}]]]];
Print[CForm[ndet]]];
Print[]];
k=k+1];

(* Part III - Sorting Partitions according to the the rank of the matrices A *)

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

Parttoprk=Table[0,{i,1,Length[lff]}];
Partsmallrk=Table[0,{i,1,Length[lff]}];

(* 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 the Partitions according to the rank of the corresponding matrices..."];
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])];

While[k≤ Length[lff],
Print[];
Print["Multigraph #"<>ToString[k]];
Print[];
j=1;
    Parttoprk[[k]]={}; 
Partsmallrk[[k]]={};
Auxgoodpart=Import[ToString["OutputPart"<>ToString[k]<>".txt"],"Table"];
Lpar=Length[Auxgoodpart];
NLpar=IntegerPart[Lpar/10];
While[j≤Lpar,
If[NCArray[[k]]==0,
auxpart=Auxgoodpart[[j]],
listedges=auxedges[[k]]/.{x_/;x>cardI-NCArray[[k]]}->{-1};
auxpart=Insert[Auxgoodpart[[j]],0,listedges]];
      Mat=A[[k]]-DiagonalMatrix[auxpart];
  lns=NullSpace[Mat];
  If[MatrixRank[Mat]==cardI-1, 
        s=lns[[1]];
        If[And@@(#1>0&)/@s||And@@(#1<0&)/@s,Parttoprk[[k]]=Append[Parttoprk[[k]],auxpart]],
   If[test2[lns], Partsmallrk[[k]]=Append[Partsmallrk[[k]],auxpart]]];
  If[NLpar≠0,If[Divisible[j,NLpar],If[j/NLpar≤ 10,Print[j/NLpar*10 "%"]]]];
j=j+1];
Export["Part"<>ToString[k]<>"toprk.txt",Parttoprk[[k]]];
Export["Part"<>ToString[k]<>"smallrk.txt",Partsmallrk[[k]]];
k=k+1];

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

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

(* 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[lff]}];
While[j≤Length[lff],
AuxIW={};
AuxPart=Flatten[Import[ToString["Part"<>ToString[j]<>"toprk.txt"],"List"]];
LAP=Length[AuxPart];
If[LAP≠0,
p=1;
Weight={};
While[p≤LAP,Mat=A[[j]]-DiagonalMatrix[ToExpression[AuxPart[[p]]]];
lns=Flatten[NullSpace[Mat]]; 
i=1;
AuxWeight={};
While[i≤ NFP, 
W={};
l=1;
While[l≤ cardI, 
If[ e[[j]][[l]][[1]]==i,W=Append[W,lns[[l]]]];
If[e[[j]][[l]][[2]]==i,W=Append[W,-lns[[l]]]];
l=l+1];
AuxWeight=Append[AuxWeight,W];
i=i+1];
Weight=Append[Weight,AuxWeight];
p=p+1],Weight={}];
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[lff]}];
While[j≤Length[lff],
AuxList={};
Print[];
Print["Multigraph #"<>ToString[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 - SMALLER RANK *)

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

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

(* 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;
PossibleNullSpaces={};
While[j≤Length[lff],
AuxPart=Flatten[Import[ToString["Part"<>ToString[j]<>"smallrk.txt"],"List"]];
LAP=Length[AuxPart];
AuxPossibleNullSpaces={};
If[LAP≠0,p=1;
While[p≤LAP,
flag=True;
AuxPartition=ToExpression[AuxPart[[p]]];
Mat=A[[j]]-DiagonalMatrix[AuxPartition];
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]];
p=p+1];
If[Length[AuxPossibleNullSpaces]≠0,PossibleNullSpaces=Append[PossibleNullSpaces,{j,AuxPossibleNullSpaces}]]];
Print["Multigraph #"<>ToString[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["PNullSpaces.txt",PNull,"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 #"<>ToString[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[[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[[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 #"<>ToString[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["FIWaux"<>ToString[q]<>".txt",FIWaux[[q]],"List"];
p=p+1];

(* Part V c) Performing several tets on the isotropy weights to reduce the number of possible cases *)

(* 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,j,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["PNullSpaces.txt","List"]];

p=1;
While[p≤Length[PNull],
k=1;
droplist4={};
q=PNull[[p]];
Print["Multigraph #"<>ToString[q]]; 
FIWauxq=ToExpression[Import["FIWaux"<>ToString[q]<>".txt","List"]];
TestRational[FIWauxq,cardI];
DropFunction[FIWauxq,droplist4];
aux3=NewList;
aux4=Simplify[aux3];
Export["NewFIWaux"<>ToString[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["PNullSpaces.txt","List"]];

p=1;
While[p≤Length[PNull],
q=PNull[[p]];
Print["Multigraph #"<>ToString[q]];
FIWauxq=ToExpression[Import["NewFIWaux"<>ToString[q]<>".txt","List"]];
LLength=Length[FIWauxq];
droplist={};
LLL=Tally[Table[e[[q]][[i]][[2]],{i,1,d}]]; 
l=1;
While[l≤Length[LLL],
If[LLL[[l,2]]==d&& d≠ 1,
LLL2=Flatten[Position[Table[e[[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[[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["NNewFIWaux"<>ToString[q]<>".txt",FIWauxIII,"List"];
p=p+1];

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

PNull=ToExpression[Import["PNullSpaces.txt","List"]];

p=1;
While[p≤Length[PNull],
q=PNull[[p]];
Print["Multigraph #"<>ToString[q]];
FIWauxq=ToExpression[Import["NNewFIWaux"<>ToString[q]<>".txt","List"]];
TestRational[FIWauxq,cardI];
DropFunction[FIWauxq,droplist4];
aux3=NewList;
aux4=Simplify[aux3];
Export["NNewFIWauxb"<>ToString[q]<>".txt",aux4,"List"];
p=p+1];

(* 4. In the cases wehere we have multiple edges we select possible cases according to Lemmas 8.1 and 8.2.
We save the resulting weights in the files "ShortFIWaux?.txt" where "?" is the number of the graph *)

PNull=ToExpression[Import["PNullSpaces.txt","List"]];

p=1;
While[p≤Length[PNull],
q=PNull[[p]];
Print["Multigraph #"<>ToString[q]];
FIWauxq=ToExpression[Import["NNewFIWauxb"<>ToString[q]<>".txt","List"]];
LLength=Length[FIWauxq];
droplist={};
LLL=Table[Tally[Table[e[[q]][[i]][[2]],{i,(j-1) (d+1/2-(j-1)/2)+1,j (d+1/2-j/2)}]],{j,1,NFP-2}];
fp=1;
While[fp≤NFP-2,
LLLfp=LLL[[fp]];
l=1;
While[l≤Length[LLLfp],
nme=LLLfp[[l,2]];
If[nme>1&&nme<d,
If[nme==d-1 && fp==2,
k=1;
While[k≤Length[FIWauxq],
G=PolynomialGCD@@Table[FIWauxq[[k]][[2]][[i]],{i,2,d}];
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]];
If[nme==d-2,
AAA=Select[e[[q]],#[[2]]==LLLfp[[l,1]]||#[[1]]==LLLfp[[l,1]]&];(* the edges that start or end at the terminal point of the multiple edges *)
BBB=Select[AAA,#[[1]]≠#[[2]]&];(* the edges of AAA that are not cycles *)
CCC=Select[BBB,#[[1]]≠fp&];(* the edges of AAA that are not cycles and are not part of this set multiple edges *)
AAA2=Select[e[[q]],#[[2]]==fp||#[[1]]==fp&];
BBB2=Select[AAA2,#[[1]]≠#[[2]]&];
CCC2=Select[BBB2,#[[2]]≠LLLfp[[l,1]]&];
If[CCC=={}||CCC2=={},
LLL2=Flatten[Position[Table[e[[q]][[(fp-1) d+i]][[2]],{i,1,d-(fp-1)}],LLLfp[[l,1]]]];
k=1;
While[k≤Length[FIWauxq],G=PolynomialGCD@@Table[FIWauxq[[k]][[fp]][[fp-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[{Abs[G]==1,FIWauxq[[k]][[fp]][[fp-1+LLL2[[1]]]]>0} ,Variables[G][[1]],Reals];
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]]];
r=2;
While[r≤nme,
LLL2=Flatten[Position[Table[e[[q]][[i]][[2]],{i,(fp-1) (d+1/2-(fp-1)/2)+1,fp (d+1/2-fp/2)}],LLLfp[[l,1]]]];
LLL3=Select[DeleteDuplicates[Table[DeleteDuplicates[Sort[Tuples[LLL2,r][[i]]]],{i,1,Length[Tuples[LLL2,r]]}]],Length[#]==r&];
k=1;
While[k≤LLength,
LLL4=Table[PolynomialGCD@@Table[FIWauxq[[k]][[fp]][[fp-1+LLL3[[vv]][[i]]]],{i,1,Length[LLL3[[vv]]]}],{vv,1,Length[LLL3]}];
rr=1;
While[rr≤Length[LLL4],
If[NumberQ[LLL4[[rr]]]&&LLL4[[rr]]>1,
AuxList=Table[fp-1+LLL3[[rr]][[i]],{i,1,Length[LLL3[[rr]]]}];
LLL5=Select[Table[i,{i,1,d}],Not[MemberQ[AuxList,#]]&];
LLL6=Table[GCD[LLL4[[rr]],FIWauxq[[k]][[fp]][[LLL5[[tt]]]]],{tt,1,Length[LLL5]}];
LLL7=DeleteCases[LLL6,1];
If[Length[LLL7]==0,
droplist=Append[droplist,k];
rr=Length[LLL4],

LLL8=Select[LLL7,#==LLL4[[rr]]||Not[NumberQ[#]]&];
If[Length[LLL8]==0,
droplist=Append[droplist,k];
rr=Length[LLL4],

 LLL9=Select[Table[FIWauxq[[k]][[LLLfp[[l,1]]]][[i]],{i,1,d}],Not[MemberQ[Table[-FIWauxq[[k]][[fp]][[fp-1+LLL3[[rr]][[i]]]],{i,1,Length[LLL3[[rr]]]}],#]]&];
If[Length[LLL9]==d-r,
LLL10=Table[GCD[LLL4[[rr]],LLL9[[i]]],{i,1,d-r}]; 
LLL11=DeleteCases[LLL10,1];
If[Length[LLL11]==0,
droplist=Append[droplist,k];
rr=Length[LLL4],

LLL12=Select[LLL11,#==LLL4[[rr]]||Not[NumberQ[#]]&];
If[Length[LLL12]==0,
droplist=Append[droplist,k];
rr=Length[LLL4]]]]]]];
rr=rr+1];
k=k+1];
r=r+1]];
l=l+1];
fp=fp+1];

DropFunction[FIWauxq,droplist];
FIWauxIV=NewList;
Export["ShortFIWaux"<>ToString[q]<>".txt",FIWauxIV,"List"];
p=p+1];

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

PNull=ToExpression[Import["PNullSpaces.txt","List"]];

p=1;
While[p≤Length[PNull],
k=1;
droplist4={};
q=PNull[[p]];
Print["Multigraph #"<>ToString[q]];
FIWauxq=ToExpression[Import["ShortFIWaux"<>ToString[q]<>".txt","List"]];
TestRational[FIWauxq,cardI];
DropFunction[FIWauxq,droplist4];
aux3=NewList;
aux4=Simplify[aux3];
Export["ShortFIWauxb"<>ToString[q]<>".txt",aux4,"List"];
p=p+1];

(* 6. Since we elliminated some common divisors we perform again the test on multiple edges to select possible cases according to Lemmas 8.1 and 8.2.
We save the resulting weights in the files "ShortFIWauxc?.txt" where "?" is the number of the graph *)

PNull=ToExpression[Import["PNullSpaces.txt","List"]];

p=1;
While[p≤Length[PNull],
q=PNull[[p]];
Print["Multigraph #"<>ToString[q]];
FIWauxq=ToExpression[Import["ShortFIWauxb"<>ToString[q]<>".txt","List"]];
LLength=Length[FIWauxq];
NLLength=IntegerPart[LLength/10];
droplist={};
LLL=Table[Tally[Table[e[[q]][[i]][[2]],{i,(j-1) (d+1/2-(j-1)/2)+1,j (d+1/2-j/2)}]],{j,1,NFP-2}];
fp=1;
While[fp≤NFP-2,
LLLfp=LLL[[fp]];
l=1;
While[l≤Length[LLLfp],
nme=LLLfp[[l,2]];
If[nme>1&&nme<d,

If[nme==d-1 && fp==2,
k=1;
While[k≤Length[FIWauxq],
G=PolynomialGCD@@Table[FIWauxq[[k]][[2]][[i]],{i,2,d}];
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]];

If[nme==d-2,
AAA=Select[e[[q]],#[[2]]==LLLfp[[l,1]]||#[[1]]==LLLfp[[l,1]]&];(* the edges that start or end at the terminal point of the multiple edges *)
BBB=Select[AAA,#[[1]]≠#[[2]]&];(* the edges of AAA that are not cycles *)
CCC=Select[BBB,#[[1]]≠fp&];(* the edges of AAA that are not cycles and are not part of this set multiple edges *)
AAA2=Select[e[[q]],#[[2]]==fp||#[[1]]==fp&];
BBB2=Select[AAA2,#[[1]]≠#[[2]]&];
CCC2=Select[BBB2,#[[2]]≠LLLfp[[l,1]]&];
If[CCC=={}||CCC2=={},
LLL2=Flatten[Position[Table[e[[q]][[(fp-1) d+i]][[2]],{i,1,d-(fp-1)}],LLLfp[[l,1]]]];
k=1;
While[k≤Length[FIWauxq],G=PolynomialGCD@@Table[FIWauxq[[k]][[fp]][[fp-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[{Abs[G]==1,FIWauxq[[k]][[fp]][[fp-1+LLL2[[1]]]]>0} ,Variables[G][[1]],Reals];
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]]];
r=2;
While[r≤nme,
LLL2=Flatten[Position[Table[e[[q]][[i]][[2]],{i,(fp-1) (d+1/2-(fp-1)/2)+1,fp (d+1/2-fp/2)}],LLLfp[[l,1]]]];
LLL3=Select[DeleteDuplicates[Table[DeleteDuplicates[Sort[Tuples[LLL2,r][[i]]]],{i,1,Length[Tuples[LLL2,r]]}]],Length[#]==r&];
k=1;
While[k≤LLength,LLL4=Table[PolynomialGCD@@Table[FIWauxq[[k]][[fp]][[fp-1+LLL3[[vv]][[i]]]],{i,1,Length[LLL3[[vv]]]}],{vv,1,Length[LLL3]}];
rr=1;
While[rr≤Length[LLL4],
If[NumberQ[LLL4[[rr]]]&&LLL4[[rr]]>1,
AuxList=Table[fp-1+LLL3[[rr]][[i]],{i,1,Length[LLL3[[rr]]]}];
LLL5=Select[Table[i,{i,1,d}],Not[MemberQ[AuxList,#]]&];
LLL6=Table[GCD[LLL4[[rr]],FIWauxq[[k]][[fp]][[LLL5[[tt]]]]],{tt,1,Length[LLL5]}];
LLL7=DeleteCases[LLL6,1];
If[Length[LLL7]==0,
droplist=Append[droplist,k];
rr=Length[LLL4],

LLL8=Select[LLL7,#==LLL4[[rr]]||Not[NumberQ[#]]&];
If[Length[LLL8]==0,
droplist=Append[droplist,k];
rr=Length[LLL4],

 LLL9=Select[Table[FIWauxq[[k]][[LLLfp[[l,1]]]][[i]],{i,1,d}],Not[MemberQ[Table[-FIWauxq[[k]][[fp]][[fp-1+LLL3[[rr]][[i]]]],{i,1,Length[LLL3[[rr]]]}],#]]&];
If[Length[LLL9]==d-r,
LLL10=Table[GCD[LLL4[[rr]],LLL9[[i]]],{i,1,d-r}]; 

LLL11=DeleteCases[LLL10,1];
If[Length[LLL11]==0,
droplist=Append[droplist,k];
rr=Length[LLL4],

LLL12=Select[LLL11,#==LLL4[[rr]]||Not[NumberQ[#]]&];
If[Length[LLL12]==0,
droplist=Append[droplist,k];
rr=Length[LLL4]]]]]]];
rr=rr+1];
k=k+1];
r=r+1]];
l=l+1];
fp=fp+1];

DropFunction[FIWauxq,droplist];
FIWauxIV=NewList;
Export["ShortFIWauxc"<>ToString[q]<>".txt",FIWauxIV,"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["PNullSpaces.txt","List"]];

p=1;
While[p≤Length[PNull],
q=PNull[[p]]; 
IWauxq=ToExpression[Import["ShortFIWauxc"<>ToString[q]<>".txt","List"]];
Print["Multigraph #"<>ToString[q]<>":"];
TestPositive[IWauxq,cardI];
DropFunction[IWauxq,droplistq];
IWListq=NewList;
Export["NewShortFIWaux"<>ToString[q]<>".txt",IWListq,"List"];
p=p+1];

(* 8. We can select the isotropy weights that may satisfy (8.2) *)


PNull=ToExpression[Import["PNullSpaces.txt","List"]];

Print[];
Print["Selecting the isotropy weights that satisfy (.....):"];
Print[];
p=1;
droplistsmall0={};
Bignumlist={};
While[p≤ Length[PNull], 
q=PNull[[p]]; 

IWauxq=ToExpression[Import["NewShortFIWaux"<>ToString[q]<>".txt","List"]];
k=1;
Length2=Length[IWauxq];
NLength2=IntegerPart[Length2/10];
Print["Multigraph #"<>ToString[q]];

While[k≤ Length2,

Array[b,cardI];
S0=Sum[IWauxq[[k]][[1]][[i]],{i,1,d}];
S1=Sum[IWauxq[[k]][[2]][[i]],{i,1,d}];
SNFPm1=Sum[IWauxq[[k]][[NFP-1]][[i]],{i,1,d}];
SNFP=Sum[IWauxq[[k]][[NFP]][[i]],{i,1,d}];
C1a=Simplify[(S1-S0)/IWauxq[[k]][[2]][[1]]];
If[d==4,If[q≠ 75,Div={1},Div={1,5}],Div=Select[Divisors[integral],#≤NFP &]];
C1b=Simplify[(SNFPm1-SNFP)/IWauxq[[k]][[NFP-1]][[d]]];

If[ToString[C1a]==ToString[C1b],

If[NumberQ[C1a]&&NumberQ[C1b]&&Not[MemberQ[Div,C1a]],droplistsmall0=Append[droplistsmall0,{q,k}]],

 Sol=Solve[C1a==C1b];

If[Length[Sol]==0, 
droplistsmall0=Append[droplistsmall0,{q,k}],

If[(NumberQ[C1a]&&Not[MemberQ[Div,C1a]])||(NumberQ[C1b]&& Not[MemberQ[Div,C1b]]),

(* One of C1a or C1b is a number but is not a divisor of "integral" smaller that NFP (if dim M=8 then it is not 1 nor 5) *)
 
droplistsmall0=Append[droplistsmall0,{q,k}],

VVsmall0=DeleteDuplicates[Flatten[Variables[{C1a,C1b}]]];

Listaux=VVsmall0/.Sol;
If[MemberQ[Table[MemberQ[NumberQ/@(Listaux[[i]]),False],{i,1,Length[Listaux]}],True]==False,

(* Then the elements of Listaux are all numbers *)

If[MemberQ[Table[MemberQ[Table[ IntegerQ[Listaux[[i]][[j]]],{j,1,Length[Listaux[[i]]]}],False],{i,1,Length[Listaux]}],False]==False,

droplistsmall0=Append[droplistsmall0,{q,k}],

vv=1;
While[vv≤ Length[Sol],
                   l=1;
aux2=IWauxq[[k]];
While[l≤NFP,
j=1;
aux={};
While[j≤dim/2,
aux=Append[aux,Simplify[IWauxq[[k]][[l]][[j]]/.Flatten[Sol[[vv]]]]];
j=j+1];
aux2=ReplacePart[aux2,l->aux];
l=l+1];
If[Length[Sol]==1, IWauxq=ReplacePart[IWauxq,k->aux2]];

Flag=True;
j=1;
loclistsmallaux={};
While[j≤dim/2-1,
Locsmall=0;
i=1;
While[i≤NFP,Locsmall=Locsmall+SymmetricPolynomial[j,aux2[[i]]]/SymmetricPolynomial[dim/2,aux2[[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[solsmall≠{},
If[Not[MemberQ[Table[MemberQ[Table[Positive[(VVsmall/.solsmall)[[i]][[j]]]&& IntegerQ[(VVsmall/.solsmall)[[i]][[j]]],{j,1,Length[(VVsmall/.solsmall)[[i]]]}],False],{i,1,Length[VVsmall/.solsmall]}],False]==False&&solsmall≠{{}}&&Length[Flatten[Variables/@Flatten[VVsmall/.solsmall]]]==0],
Flag=False,Bignumlist=Append[Bignumlist,NSLocsmall]]]];
Clear[solsmall];
vv=vv+1];
If[Flag, droplistsmall0=Append[droplistsmall0,{q,k}]]],

If[NumberQ[C1a]||NumberQ[C1b],
aux=Table[IWauxq[[k]]/.Sol[[i]],{i,1,Length[Sol]}];
l=1;
Flag=True;
While[l≤ Length[aux],
j=1;
loclistsmallaux={};
While[j≤dim/2-1,
Locsmall=0;
i=1;
While[i≤NFP,Locsmall=Locsmall+SymmetricPolynomial[j,aux[[l]][[i]]]/SymmetricPolynomial[dim/2,aux[[l]][[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[solsmall≠{},
If[Not[MemberQ[Table[MemberQ[Table[Positive[(VVsmall/.solsmall)[[i]][[j]]]&& IntegerQ[(VVsmall/.solsmall)[[i]][[j]]],{j,1,Length[(VVsmall/.solsmall)[[i]]]}],False],{i,1,Length[VVsmall/.solsmall]}],False]==False&&solsmall≠{{}}&&Length[Flatten[Variables/@Flatten[VVsmall/.solsmall]]]==0],
Flag=False,Bignumlist=Append[Bignumlist,NSLocsmall]]]];
Clear[solsmall];
l=l+1];
If[Flag, droplistsmall0=Append[droplistsmall0,{q,k}],If[Length[Sol]==1,ReplacePart[IWauxq,k->aux2]]],


TTT={};
auxa=Simplify[C1a/.Sol];
vv=1;
While[vv≤ Length[auxa],
ii=1;
While[ii≤Length[Div],
Saux=Flatten[Solve[auxa[[vv]]==Div[[ii]]],1];
If[Saux≠ {},TTT=Append[TTT,Saux]];
ii=ii+1];
vv=vv+1];

aux=Table[IWauxq[[k]]/.(TTT[[i]][[j]]),{i,1,Length[TTT]},{j,1,Length[TTT[[i]]]}];
If[Length[TTT]==1, IWauxq=ReplacePart[IWauxq,k->Flatten[aux,2]]];
kk=1;
numerators={};
Flag=True;
While[kk≤ Length[aux],
                   l=1;
While[l≤ Length[aux[[kk]]],
If[MemberQ[Flatten[Table[Simplify[Element[aux[[kk]][[l]][[i]][[j]],Integers],Table[Element[b[i],Integers],{i,1,cardI}]],{i,1,NFP},{j,1,d}]],False]==False,
j=1;
loclistsmallaux={};
While[j≤dim/2-1,
Locsmall=0;
i=1;
While[i≤NFP,Locsmall=Locsmall+SymmetricPolynomial[j,aux[[kk]][[l]][[i]]]/SymmetricPolynomial[dim/2,aux[[kk]][[l]][[i]]];
i=i+1];
SLocsmall=Together[Locsmall];
loclistsmallaux=Append[loclistsmallaux,SLocsmall];
j=j+1];
NSLocsmall=Numerator/@loclistsmallaux;
If[Not[MemberQ[numerators,NSLocsmall]]&& Not[MemberQ[Bignumlist,NSLocsmall]],
 numerators=Append[numerators,NSLocsmall];
VVsmall=DeleteDuplicates[Flatten[Variables/@NSLocsmall]];
solsmall=Solve[NSLocsmall==0];
If[solsmall≠ {},
If[Not[MemberQ[Table[MemberQ[Table[Block[{$MaxExtraPrecision=1000},Positive[(VVsmall/.solsmall)[[i]][[j]]]]&& Block[{$MaxExtraPrecision=1000},IntegerQ[(VVsmall/.solsmall)[[i]][[j]]]],{j,1,Length[(VVsmall/.solsmall)[[i]]]}],False],{i,1,Length[VVsmall/.solsmall]}],False]==False&&solsmall≠{{}}&&Length[Flatten[Variables/@Flatten[VVsmall/.solsmall]]]==0],
Flag=False,Bignumlist=Append[Bignumlist,NSLocsmall]]];Clear[solsmall]]];
l=l+1];
kk=kk+1];
If[Flag,droplistsmall0=Append[droplistsmall0,{q,k}]]]]]]];
If[NLength2≠0,If[Divisible[k,NLength2],If[k/NLength2≤10,Print[k/NLength2*10 "%"]]]];
k=k+1];
p=p+1];

p=1;
While[p≤Length[PNull],
q=PNull[[p]];
IWauxq=ToExpression[Import["NewShortFIWaux"<>ToString[q]<>".txt","List"]];
DropFunction2[IWauxq,droplistsmall0,q];
IWaux=NewList;
Export["C1FIWaux"<>ToString[q]<>".txt",IWaux,"List"];
p=p+1];

(* 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["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["C1FIWaux"<>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["C1FIWaux"<>ToString[q]<>".txt","List"]];
DropFunction2[IWauxq,droplistsmall,q];
IWaux=NewList;
Export["FinalIWaux"<>ToString[q]<>".txt",IWaux,"List"];
p=p+1];

(* 9. Test for possibility of integer solutions and save the possible ones in "FinalIW?.txt" *)

PNull=ToExpression[Import["PNullSpaces.txt","List"]];

p=1;
While[p≤Length[PNull],
q=PNull[[p]]; 
IWauxq=ToExpression[Import["FinalIWaux"<>ToString[q]<>".txt","List"]];
Print["Multigraph #"<>ToString[q]<>":"];
TestPositive[IWauxq,cardI];
DropFunction[IWauxq,droplistq];
IWListq=NewList;
Export["FinalIW"<>ToString[q]<>".txt",IWListq,"List"];
p=p+1];

(* 10. Delete duplicates: *)

PNull=ToExpression[Import["PNullSpaces.txt","List"]];

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

(* 11. Whenever there is a weight that is of the form b[i]/m we make b[i]=0 mod m, substituting b[i] by m*b[i] *)

p=1;
While[p≤Length[PNull],
q=PNull[[p]];
IWauxq=IWFinal[[q]];
k=1;
While[k≤ Length[IWauxq],
j=1;
While[j≤ NFP,
vv=1;
While[vv≤ d,
If[Variables[IWauxq[[k]][[j]][[vv]]]≠ {},Vaux=Variables[IWauxq[[k]][[j]][[vv]]];
If[Length[Vaux]==1,
LLLL=Denominator[Coefficient[IWauxq[[k]][[j]][[vv]],Variables[IWauxq[[k]][[j]][[vv]]]]][[1]];
If[LLLL≠1,V4=IWauxq[[k]]/.{Variables[IWauxq[[k]][[j]][[vv]]][[1]]-> LLLL*Variables[IWauxq[[k]][[j]][[vv]]][[1]]};
IWauxq=ReplacePart[IWauxq,k->V4]]]];
vv=vv+1];
j=j+1];
k=k+1];
IWFinal=ReplacePart[IWFinal, q-> IWauxq];
p=p+1];

(* 12. Then we delete duplicates again... *)

PNull=ToExpression[Import["PNullSpaces.txt","List"]];

p=1;
While[p≤Length[PNull],
q=PNull[[p]];
IWauxq=IWFinal[[q]];
FIWaux3=DeleteDuplicates[Table[Sort[Expand[IWauxq[[i,j]]]],{i,1,Length[IWauxq]},{j,1,NFP}]];
IWFinal=ReplacePart[IWFinal,q-> FIWaux3];
Export["ReallyLastIWaux"<>ToString[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;
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];