On[Assert];

Nps = {1,1,5,1,1};
dim = 2*Length[Nps]-2;
n = dim/2;
NFP = Plus@@Nps;
Assert[Length[Nps] == dim/2+1];
Assert[Plus @@ Nps == NFP];
Assert[Nps == Reverse[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 == 8,
	If[NFP == 9,
		vertexCoordinates = {{0,0},{1,1},{2,2},{3,3},{4,4},{3,5},{2,6},{1,7},{0,8}};]];
d = dim/2;
cardI = NFP*dim/4;

liaux=Table[Select[ Flatten[Permutations /@ IntegerPartitions[If[i≤2,d-i,If[i==3||i==4||i==5,2,d-i+4]], {NFP}, Table[j, {j, 0, If[i≤2,d-i,If[i==3||i==4||i==5,2,d-i+4]]}]], 1], #[[1]] == 0 && If[i≠0,#[[2]] == 0,#[[2]]==1 ] && #[[3]]≤2&&  #[[4]]≤2 && #[[5]]≤2  && #[[6]]≤2 && #[[7]]≤2&&  #[[8]]≤3 && #[[9]]≤4 && 
If[i==3, #[[3]]==0,True] && If[i==4, #[[4]]==0 && #[[3]]==0,True] && If[i==5, #[[5]]==0 &&  #[[4]]==0 &&  #[[3]]==0 ,True]&& If[i==6, #[[6]]==0 &&#[[5]]==0 &&  #[[4]]==0 &&  #[[3]]==0 ,True] && If[i==7, #[[9]]==1  ,True]
&& #[[i + 1]] == 0 &] , {i, 0, 8}];

k=1;
lf={};
While[k≤Length[liaux[[1]]],
liaux2=ReplacePart[liaux,1->{liaux[[1]][[k]]}];
li=Tuples[liaux2];
Print["k: ", k];
 AppendTo[lf,Select[li, Total[#]==Total[#]==Table[If[i≤2,i,If[i≠8,Floor[(i+5)/4],4]],{i,0,8}]&]];
k=k+1]; 

lf=Flatten[lf,1]; 
li={}; 

(* Deleting duplicate by reversing the flow *)

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

(* 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], 
index = 1;
j=1;
While[j≤ NFP, 
r=1;
While[r≤ NFP, 
entry=lff[[p]][[j]][[r]];
k=1;
While[ k≤ entry,
e[[p]][[index]][[1]]=j;
e[[p]][[index]][[2]]=r;
++index;
++k
];
r=r+1];
j=j+1];
p=p+1];

(* Functions "deltai" and "deltat" for each multigraph of "lff": 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];

Index=2;
toPart = Sum[Nps[[p+1]]*(6*p*(p-1) + (5n - 3n^2)/2), {p,0,n}]/Index;
parts0 = IntegerPartitions[toPart, {cardI}];
parts=parts0*Index;
parts = Join @@ (Permutations/@parts);
parts = DeleteDuplicates[parts];
Monitor[
For[i = 1, i≤Length[lff], ++i,

toStore = {};

For[j = 1, j≤Length[parts], ++j,
	If[parts[[j]][[1]]==Index && parts[[j]][[cardI]]==Index && 0 == Det[A[[i]]-DiagonalMatrix[parts[[j]]]],AppendTo[toStore, parts[[j]]]];
];

count=i;
fileName = FileNameJoin[{NotebookDirectory[],"graphsb2_1b4_5", "graph"<>ToString[count]<>".txt"}];
If[FileExistsQ[fileName], DeleteFile[fileName]];
outfile = OpenWrite[fileName];
Export[outfile, toStore, "Table"];
Close[outfile];

];, ProgressIndicator[i, {1, Length[lff]}]];
Delete[parts];
Delete[toStore];

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

Print[];
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])];
Monitor[
While[k≤ Length[lff],
j=1;
    Parttoprk[[k]]={}; 
Partsmallrk[[k]]={};
c=k;
Auxgoodpart=Import[FileNameJoin[{NotebookDirectory[],"graphsb2_1b4_5","graph"<>ToString[c]<>".txt"}],"Table"];
Lpar=Length[Auxgoodpart];
While[j≤Lpar,
auxpart=Auxgoodpart[[j]];
            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]]];
j=j+1];
Export[FileNameJoin[{NotebookDirectory[],"graphsb2_1b4_5","Part"<>ToString[c]<>"toprk.txt"}],Parttoprk[[k]]];
Export[FileNameJoin[{NotebookDirectory[],"graphsb2_1b4_5","Part"<>ToString[c]<>"smallrk.txt"}],Partsmallrk[[k]]];
k=k+1];, ProgressIndicator[k, {1, Length[lff]}]];

(*Part 4 - 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[{},{i,1,Length[lff]}];
Monitor[
While[j≤Length[lff],
AuxIW={};
c=j;
AuxPart=Flatten[Import[FileNameJoin[{NotebookDirectory[],"graphsb2_1b4_5",ToString["Part"<>ToString[c]<>"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];, ProgressIndicator[j, {1, Length[lff]}]];

(* 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[{},{i,1,Length[lff]}];
Monitor[
While[j≤Length[lff],
AuxList={};
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; 
If[FIW[[j]]≠{},
Print[];
Print["Multigraph #"<>ToString[c]<>":"]; 
Print["The possible isotropy weights are:"];
Print[FIW[[j]]] ]];
j=j+1];, ProgressIndicator[j, {1, Length[lff]}]];

(* PART 5 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={};

Monitor[
While[j≤Length[lff],
c=j;
AuxPart=Flatten[Import[FileNameJoin[{NotebookDirectory[],"graphsb2_1b4_5",ToString["Part"<>ToString[c]<>"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}]]];
j=j+1];, ProgressIndicator[j, {1, Length[lff]}]];


(* 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_1b4_5","PNullSpaces.txt"}],PNull,"List"];

(* PART 5 b) ISOTROPY WEIGHTS *)

(* First we consider the span of the different null spaces 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 element of "lff" *)


Print[];
Print["Determining the possible isotropy weights..."];
p=1;
IW2=Table[0,{i,1,Length[lff]}];
While[p≤Length[PossibleNullSpaces2],

c=PossibleNullSpaces2[[p]][[1]];
Print["Multigraph #"<>ToString[c]];
LLength=Length[PossibleNullSpaces2[[p]][[2]]];
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];
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 multigraph *)

FIWaux=IW2;
p=1;
Print[];
While[p≤Length[PossibleNullSpaces2],
LLength=Length[PossibleNullSpaces2[[p]][[2]]];
NLLength=IntegerPart[LLength/10];
q=PossibleNullSpaces2[[p]][[1]];
c=q;
Print["Multigraph #"<>ToString[c]];
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_1b4_5","FIWaux"<>ToString[c]<>".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 ,2,If[j==3||j==4||j==5||j==6||j==7,3,If[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[],"graphsb2_1b4_5","PNullSpaces.txt"}],"List"]];

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

(* 2. 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_1b4_5","PNullSpaces.txt"}],"List"]];

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

(* 3. 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[],"graphsb2_1b4_5","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]]; 
c=q;
IWauxq=ToExpression[Import[FileNameJoin[{NotebookDirectory[],"graphsb2_1b4_5","NewShortFIWaux"<>ToString[c]<>".txt"}],"List"]];
k=1;
Print[];
Length2=Length[IWauxq];

Print["Multigraph #"<>ToString[c]<>":"];
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]];

If[Length[VVsmall]==3,
solsmall=Solve[Drop[NSLocsmall,-1]==0,Drop[VVsmall,2],PositiveIntegers];,
solsmall=Solve[NSLocsmall==0,Drop[VVsmall,-1],PositiveIntegers]];

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}]];
k=k+1];
DropFunction2[IWauxq,droplistsmall,q];
IWaux=NewList;
Export[FileNameJoin[{NotebookDirectory[],"graphsb2_1b4_5","FinalIWaux"<>ToString[c]<>".txt"}],IWaux,"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; FL=False;
IWFinal=Table[{},{i,1,Length[lff]}];
While[p≤Length[lff],
c=p;
If[MemberQ[p,PNull],IWFinal[[p]]=ToExpression[Import[FileNameJoin[{NotebookDirectory[],"graphsb2_1b4_5","FinalIWaux"<>ToString[c]<>".txt"}],"List"]]];
If[FIW[[p]]≠{}|| IWFinal[[p]]≠ {}, NS=Length[IWFinal[[p]]]+Length[FIW[[p]]];
Print["**************************************************************************"];
Print["All Possible Isotropy Weights for Multigraph #"<>ToString[c]<>":"];
Print[" "];
Print["There are at most "<>ToString[NS]<>" solutions:"]];
If[FIW[[p]]≠{}, 
Print[FIW[[p]]];
Print[]];
If[IWFinal[[p]]≠ {},k=1;While[k≤ Length[IWFinal[[p]]], Print[IWFinal[[p]][[k]]];Print[]; k=k+1]];
If[FIW[[p]]≠{}|| IWFinal[[p]]≠ {},Fl=True; Print[GraphPlot[lff[[p]],DirectedEdges->True,VertexLabeling->True,MultiedgeStyle->True,SelfLoopStyle->True,VertexCoordinateRules->vertexCoordinates,ImageSize->Large]]];
p=p+1];
If[FL==False, Print["There are no possible circle actions"]];