Construction of MST complete interaction bubbles

Generate perfect 1-factorizations of the complete graph using backtracking, recursion, pruning, randomization, etc.

SetDirectory@NotebookDirectory[];

Estimating the number of graphs

Number of colorings of the complete graph on n vertices with n-1 colors

Table[(((n-1))^(Binomial[n,2])),{n,2,8,2}]

Choose all edges of the first color, then the second one, etc

Table[Product[Binomial[Binomial[n,2]-i(n/2),(n/2)],{i,0,n-2}],{n,2,8,2}]

Choose coloring of the first node, then the next one, etc

Table[Product[i!,{i,n-1}],{n,2,8,2}]

Bubble construction

Generating all matchings

Generating all the matchings (1-factors of the complete graph with n vertices) is equivalent to generating all possible Wick contractions of n fields, hence the name of this function.

ClearAll@Wick;
Wick[{}]={{}};
Wick[l_]:=
Flatten[Table[
Prepend[#,{l[[1]],l[[i]]}]&/@Wick@Drop[l[[2;;]],{i-1}]
,{i,2,Length@l}],1]

ClearAll@Factors;
Factors[n_Integer]:=Factors[n]=Wick@Range@n;

Building compatibility matrices

The compatibility matrix tells us if two 1-factors can be used in the same factorization. We use the asymmetric version to speed up the algorithm by choosing a specific order in which factors should be added.

ClearAll@CompatibilityMatrix;
CompatibilityMatrix[n_Integer]:=CompatibilityMatrix[n]=
FromDigits[#,2]&/@Table[
If[OrderedQ@{f1,f2},Boole[Intersection[f1,f2]==={}],0]
,{f1,Factors@n},{f2,Factors@n}];

We may want to recover the symmetric compatibility matrix for some purposes

ClearAll@SymmetricCompatibilityMatrix;
SymmetricCompatibilityMatrix[n_Integer]:=SymmetricCompatibilityMatrix[n]=
With[{m=IntegerDigits[#,2,(n-1)!!]&/@CompatibilityMatrix[n]},FromDigits[#,2]&/@(m+Transpose@m)]

We can actually cook up a compatibility matrix that incorporates the MST condition checking.

ClearAll@CompatibilityMatrix2;
CompatibilityMatrix2[n_Integer]:=CompatibilityMatrix2[n]=
FromDigits[#,2]&/@ParallelTable[
If[OrderedQ@{f1,f2},Boole[Intersection[f1,f2]==={}&&ConnectedGraphQ@Graph@Join[f1,f2]],0]
,{f1,Factors@n},{f2,Factors@n}];

Pre-computed matrices for n=12

(*Export["CompMatrix_r12.m",CompatibilityMatrix@12]*)

CompatibilityMatrix[12]=Import["CompMatrix_r12.m"];

(*Export["CompMatrix2_r12.m",CompatibilityMatrix2@12]*)

CompatibilityMatrix2[12]=Import["CompMatrix2_r12.m"];

Checking MST

MaximallySingleTraceQ[c_Integer,m_]:=And@@Flatten@Table[ConnectedGraphQ@Graph@Join[Position[m,c],Position[m,j]],{j,c-1}]

Step filtering (optional)

Using the asymmetric compatibility matrices, the 1-factors form a directed acyclic graph (DAG).
Paths in this graph are counted by entries in the matrix-power of the compatibility matrix: (Subscript[((C^k)), ij]) is the number of paths from factor i to factor j using k+1 factors in total.
The longest path in the DAG has n-1 factors, i.e. there is no path with n factors in total: (C^(n-1)) = (Subscript[0, (n-1)!! x (n-1)!!])
We can argue that each 1-factor has a specific “time” or level in which it should be incorporated to the graph: if we add it too early, it will impose too many constraints on the next ones, if we add it too late there won’t be enough compatible 1-factors to finish.
We compute a level mask by looking at which 1-factos are never again accessible after a given number of steps.

ClearAll@CompatibilityMatrixPower;
CompatibilityMatrixPower[n_Integer,0]:=UnitVector[(n-1)!!,1];
CompatibilityMatrixPower[n_Integer,e_Integer]:=CompatibilityMatrixPower[n,e]=CompatibilityMatrixPower[n,e-1].(IntegerDigits[#,2,(n-1)!!]&/@CompatibilityMatrix@n)//Unitize;

ClearAll@LevelMask;
LevelMask[n_Integer,c_Integer]:=
CompatibilityMatrixPower[n,c-1]-Sum[CompatibilityMatrixPower[n,e],{e,c,n-1}]//Positive//Boole//FromDigits[#,2]&

ClearAll@CompatibilityMatrixPower2;
CompatibilityMatrixPower2[n_Integer,0]:=UnitVector[(n-1)!!,1];
CompatibilityMatrixPower2[n_Integer,e_Integer]:=CompatibilityMatrixPower2[n,e]=CompatibilityMatrixPower2[n,e-1].(IntegerDigits[#,2,(n-1)!!]&/@CompatibilityMatrix2@n)//Unitize;

ClearAll@LevelMask2;
LevelMask2[n_Integer,c_Integer]:=LevelMask2[n,c]=
CompatibilityMatrixPower2[n,c-1]-Sum[CompatibilityMatrixPower2[n,e],{e,c,n-1}]//Positive//Boole//FromDigits[#,2]&

Isomorphism checking

ColorMatrix takes a list of 1-factors and builds an adjacency matrix A where (Subscript[A, ij]) is the color of the edge i<->j

ColorMatrix[n_Integer,p_List]:=
Block[{m=ConstantArray[0,{n,n}],c=0},
Normal@Sum[
c++;
c AdjacencyMatrix@Graph@Factors[n][[f]]
,{f,p}]
]

FromColorMatrix reconstructs a list of 1-factors from the given color matrix (the inverse operation of ColorMatrix).

FromColorMatrix[n_Integer,m_]:=Sort@Flatten@Table[FirstPosition[Factors@n,Sort@DeleteDuplicates[Sort/@Position[m,c]]],{c,n-1}]

CanonicalColoring finds a canonical representation of the adjacency matrix to check for bubble isomorphism (using only the color permutations in plist).
This is a very slow function, nauty would do much better!

ClearAll@CanonicalColoring;
CanonicalColoring[n_Integer,m_][plist_List]:=
Min@Table[
Block[{newm=m[[p,p]]},
newm=Replace[newm,First@#->Last@#&/@Transpose@{First@newm,Range@n-1},{2}];
FromDigits[Flatten@newm,n]
]
,{p,plist}]

FromCanonicalColoring recovers an adjacency matrix from its canonical representation.

FromCanonicalColoring[n_Integer,id_Integer]:=Partition[IntegerDigits[id,n,(n^2)],n];

ClearAll@PartialCanonicalColoring;
PartialCanonicalColoring[n_Integer,c_Integer,m_]:=Min@Table[FromDigits[Flatten@m[[p,p]],c+1],{p,Permutations@Range@n}]

These are subsets of permutations we may want to try.

ClearAll@Perms
Perms[n_Integer]:=Perms[n]=PermutationList[#,n]&/@GroupElements@GraphAutomorphismGroup@Graph@Partition[Range@n,2]

ClearAll@Perms2;
Perms2[n_Integer]:=Perms2[n]=Sort[Flatten/@Flatten[Permutations/@Factors@n,1]];

Bubble generation

GenerateBubbles is the basic recursive generator

ClearAll@GenerateBubbles;
GenerateBubbles[n_Integer,1]:={{1}};
GenerateBubbles[n_Integer,c_Integer]:=GenerateBubbles[n,c]=
Select[
Flatten[Table[
Block[{m=CompatibilityMatrix[n][[b]],p},
p=Select[Range[(n-1)!!]IntegerDigits[BitAnd@@m,2,(n-1)!!],Positive];
Append[b,#]&/@p
]
,{b,GenerateBubbles[n,c-1]}],1]
,MaximallySingleTraceQ[c,ColorMatrix[n,#]]&]

GenerateBubbles2 incorporates level masks to speed up the computations

ClearAll@GenerateBubbles2;
GenerateBubbles2[n_Integer,1]:={{1}};
GenerateBubbles2[n_Integer,c_Integer]:=GenerateBubbles2[n,c]=
Select[
Flatten[Table[
Block[{m=CompatibilityMatrix[n][[b]],p},
p=Select[Range[(n-1)!!]IntegerDigits[BitAnd[LevelMask[n,c],BitAnd@@m],2,(n-1)!!],Positive];
Append[b,#]&/@p
]
,{b,GenerateBubbles2[n,c-1]}],1]
,MaximallySingleTraceQ[c,ColorMatrix[n,#]]&]

GenerateBubbles3 uses the compatibility matrix to also account for the MST condition

ClearAll@GenerateBubbles3;
GenerateBubbles3[n_Integer,1]:={{1}};
GenerateBubbles3[n_Integer,c_Integer]:=GenerateBubbles3[n,c]=
Flatten[Table[
Block[{m=CompatibilityMatrix2[n][[b]],p},
p=Select[Range[(n-1)!!]IntegerDigits[BitAnd[LevelMask[n,c],BitAnd@@m],2,(n-1)!!],Positive];
Append[b,#]&/@p
]
,{b,GenerateBubbles3[n,c-1]}],1]

Results

n=6

There is a unique n=6 complete MST interaction bubble, up to isomorphisms.

Length/@Table[GenerateBubbles[6,c],{c,5}]//AbsoluteTiming

Length/@Table[GenerateBubbles2[6,c],{c,5}]//AbsoluteTiming

Length/@Table[GenerateBubbles3[6,c],{c,5}]//AbsoluteTiming

DeleteDuplicatesBy[GenerateBubbles2[6,5],CanonicalColoring[6,ColorMatrix[6,#]][Perms@6]&]

Length@%

n=8

There is a unique n=8 complete MST interaction bubble, up to isomorphisms.

Length/@Table[GenerateBubbles[8,c],{c,7}]//AbsoluteTiming

Length/@Table[GenerateBubbles2[8,c],{c,7}]//AbsoluteTiming

Length/@Table[GenerateBubbles3[8,c],{c,7}]//AbsoluteTiming

DeleteDuplicatesBy[GenerateBubbles2[8,7],CanonicalColoring[8,ColorMatrix[8,#]][Perms@8]&]

Length@%

n=10

There is a unique n=10 complete MST interaction bubble, up to isomorphisms.

Length/@Table[GenerateBubbles3[10,c],{c,9}]//AbsoluteTiming

DeleteDuplicatesBy[GenerateBubbles2[10,9],CanonicalColoring[10,ColorMatrix[10,#]][Perms@10]&]

Length@%

DeleteDuplicatesBy[%%,CanonicalColoring[10,ColorMatrix[10,#]][Perms2@10]&]

Length@%

n=12

With the right isomorphism checking function, we should find there are 5 inequivalent n=12 bubbles

First@GenerateBubbles2[12,2]

GenerateBubbles3[12,2]={{1,1070}};

Length@GenerateBubbles3[12,3]

Length@GenerateBubbles3[12,4]

Length@GenerateBubbles3[12,5]

Length@GenerateBubbles3[12,6]

Length@GenerateBubbles3[12,7]

Length@GenerateBubbles3[12,8]

Length@GenerateBubbles3[12,9]

Length@GenerateBubbles3[12,10]

Length@GenerateBubbles3[12,11]

Export["Bubbles_n12.m",GenerateBubbles3[12,11]]

DeleteDuplicatesBy[GenerateBubbles3[12,11],CanonicalColoring[12,ColorMatrix[12,#]][Perms@12]&]

Length@%

DeleteDuplicatesBy[%%,CanonicalColoring[12,ColorMatrix[12,#]][Perms2@12]&]

Length@%

Random bubble construction

Random generation of bubbles (we use the symmetric compatibility matrix to avoid skewing our choices)

ClearAll@GenerateRandomBubble;
GenerateRandomBubble[n_Integer,1]:={1};
GenerateRandomBubble[n_Integer,c_Integer]:=
Block[{b=GenerateRandomBubble[n,c-1],m,p,res},
m=SymmetricCompatibilityMatrix[n][[b]];
p=Select[Range[(n-1)!!]IntegerDigits[BitAnd@@m,2,(n-1)!!],Positive];
If[p=={},
GenerateRandomBubble[n,c],
res=Append[b,RandomChoice@p];
If[!MaximallySingleTraceQ[c,ColorMatrix[n,res]],
GenerateRandomBubble[n,c],
res
]
]
]

Random generation of bubbles using level masks, but otherwise same as above

ClearAll@GenerateRandomBubble2;
GenerateRandomBubble2[n_Integer,1]:={1};
GenerateRandomBubble2[n_Integer,c_Integer]:=
Block[{b=GenerateRandomBubble2[n,c-1],m,p,res},
If[c≥9,Print@c];
m=SymmetricCompatibilityMatrix[n][[b]];
p=Select[Range[(n-1)!!]IntegerDigits[BitAnd[LevelMask[n,c],BitAnd@@m],2,(n-1)!!],Positive];
If[p=={},
GenerateRandomBubble2[n,c],
res=Append[b,RandomChoice@p];
If[!MaximallySingleTraceQ[c,ColorMatrix[n,res]],
GenerateRandomBubble2[n,c],
res
]
]
]

Random generation of bubbles using level masks, but otherwise same as above

Do[
TotTime[c]=0;
NumGen[c]=0;
,{c,11}]

ClearAll@GenerateRandomBubbleTiming;
GenerateRandomBubbleTiming[n_Integer,1]:={1};
GenerateRandomBubbleTiming[n_Integer,c_Integer]:=
Block[{b,m,p,res,time},
{time,b}=AbsoluteTiming@GenerateRandomBubbleTiming[n,c-1];
TotTime[c-1]+=time;
NumGen[c-1]++;
m=SymmetricCompatibilityMatrix[n][[b]];
p=Select[Range[(n-1)!!]IntegerDigits[BitAnd[LevelMask[n,c],BitAnd@@m],2,(n-1)!!],Positive];
If[p=={},
GenerateRandomBubbleTiming[n,c],
res=Append[b,RandomChoice@p];
If[!MaximallySingleTraceQ[c,ColorMatrix[n,res]],
GenerateRandomBubbleTiming[n,c],
res
]
]
]

GenerateRandomBubble2[10,9]//AbsoluteTiming

GenerateRandomBubbleTiming[10,9]//AbsoluteTiming

GenerateRandomBubbleTiming[12,5]//AbsoluteTiming

GenerateRandomBubbleTiming[12,6]//AbsoluteTiming

GenerateRandomBubbleTiming[12,7]//AbsoluteTiming

GenerateRandomBubbleTiming[12,8]//AbsoluteTiming

GenerateRandomBubbleTiming[12,9]//AbsoluteTiming

GenerateRandomBubbleTiming[12,10]//AbsoluteTiming

GenerateRandomBubbleTiming[12,11]//AbsoluteTiming

Table[TotTime@c,{c,8}]

Table[NumGen@c,{c,8}]

%%/%

Plots for the problem statement

ClearAll@Colors;
Colors[3][i_]:={Red,Green,Blue}[[i]];
Colors[5][i_]:={Red,Green,Blue,Black,Cyan}[[i]];

Coloring[m_]:=
With[{n=Length@m},
CompleteGraph[n,EdgeStyle->Flatten@Table[i\[UndirectedEdge]j->Colors[n-1][m[[i,j]]],{i,n},{j,i-1}]]
]

Coloring@ColorMatrix[6,First@GenerateBubbles[6,5]]

Export["cMSTr5.pdf",%]

Coloring@ColorMatrix[4,First@GenerateBubbles[4,3]]

Export["tetrahedral.pdf",%]

Graph[{1\[UndirectedEdge]2,1\[UndirectedEdge]2,1\[UndirectedEdge]2}]

Export["simple.pdf",%]

Graph[{1\[UndirectedEdge]2,1\[UndirectedEdge]2,1\[UndirectedEdge]3,2\[UndirectedEdge]4,3\[UndirectedEdge]4,3\[UndirectedEdge]4}]

Export["pillow.pdf",%]