Initialize

Utils

Also dictionary from koller - util

Notation

On[Assert];
$AssertFunction=(Message[Assert::asrtf,#1];Print[#1];Abort[])&;
Needs["Notation`"];
InfixNotation[(NotationTemplateTag[∈]),MyElement];
InfixNotation[(NotationTemplateTag[∉]),NotMyElement];
AddInputAlias[(NotationTemplateTag[∈]),"in",#]&/@Notebooks[];
AddInputAlias[(NotationTemplateTag[∉]),"notin",#]&/@Notebooks[];

Mine

Tilde[a_,b_]:=Join[a,b];
MyElement[a_,b_List]:=MemberQ[b,a];
NotMyElement[a_,b_List]:=Not[MyElement[a,b]];
Tilde[x__]:=Join[x];
Subset[a_List,b_List]:=(a⋂b==a)
NotSubset[a_List,b_List]:=Not[Subset[a,b]];
Backslash[a_List,b_List]:=(Complement[a,b]);
BracketingBar[a_List]:=Length[a];
disjoint[A_,B_]:=Length[A⋂B]==0;


(* order preserving intersection *)
intersect[a_,b_]:=Select[a,MemberQ[b,#]&];
fl1[lst_]:=Flatten[lst,1];
pos1[lst_,item_]:=Position[lst,item][[1,1]];

(* custom tables *)
SetAttributes[sum,HoldAll];
SetAttributes[prod,HoldAll];
SetAttributes[table,HoldAll];
SetAttributes[anyTrue,HoldAll];
SetAttributes[allTrue,HoldAll];

Clear[allTrue,anyTrue];
table[var_∈lis_,expr_]:=ReleaseHold[Hold[expr]/.HoldPattern[var]->#]&/@lis;
prod[var_∈lis_,expr_]:=Times@@table[var∈lis,expr];
sum[var_∈lis_,expr_]:=Plus@@table[var∈lis,expr];
anyTrue[var_∈lis_,expr_]:=LengthWhile[lis,Not[TrueQ[ReleaseHold[Hold[expr]/.HoldPattern[var]->#]]]&]<Length[lis];
allTrue[var_∈lis_,expr_]:=LengthWhile[lis,TrueQ[ReleaseHold[Hold[expr]/.HoldPattern[var]->#]]&]==Length[lis];

(** new **)
(* gets all pairs of sets that intersect *)
intersectingPairs[bags_List]:=Module[{cands,result},
cands=Outer[List,bags,bags,1]//fl1;
Select[cands,#[[1]]≠#[[2]]&&!disjoint[#[[1]],#[[2]]]&&OrderedQ[#]&]
];

(* maximum spanning tree *)
makeSpanningBagTree[bags_]:=Module[{jtree,edgeCandidates,bestEdge,lens},
acyclicEdges[ed_]:=
Select[ed,AcyclicGraphQ[Graph[UndirectedEdge@@@(jtree∼{#})]]&];

findHeaviestEdge[ed_]:=(
lens=\[LeftBracketingBar]#[[1]]⋂#[[2]]\[RightBracketingBar]&/@ed;
If[lens=={},{},ed[[Last@Ordering@lens]]]
);

(* adds pair of bags with largest intersection to tree, removes cyclical candidates *)
growTree[]:=(
bestEdge=findHeaviestEdge[edgeCandidates];
jtree=jtree∼{bestEdge};
edgeCandidates=edgeCandidates∖{bestEdge};
edgeCandidates=acyclicEdges[edgeCandidates];
);
edgeCandidates=intersectingPairs[bags];
jtree={};
While[edgeCandidates≠{},growTree[]];
Sort[jtree]
];

(* selects smallest item given list of partial orders *)
(* selectSmallest[{1,2,3},{1,1,1},{1,1,2}] gives 1,2 *)
selectSmallest[items_,orders__]:=(
ordList=Transpose[{orders}];
Assert[Length[items]==Length[ordList]];
best=First@Sort@ordList;
poses=Position[ordList,_?(#==best&)];
Sort[Extract[items,poses]]
);

idx[list_,val_]:=If[val∈list,Position[list,val]//Flatten//First,-1];

(*** dictionaries ***)
(* extracts key/value from the downvalue *)
Clear[index,value,indices,items,values,indexQ,reverseLookup];
index[downvalue_,dict_]:=(downvalue[[1]]/.HoldPattern[dict[x_]]->x)//ReleaseHold;
value[downvalue_]:=downvalue[[-1]];
indices[dict_]:=Map[#[[1]]/.{HoldPattern[dict[x_]]->x}&,DownValues[dict]]//ReleaseHold;
values[dict_]:=Map[#[[-1]]&,DownValues[dict]];
items[dict_]:=Map[{index[#,dict],value[#]}&,DownValues[dict]];
indexQ[dict_,index_]:=If[MatchQ[dict[index],HoldPattern[dict[index]]],False,True];

(* take dictionary with multiple values per key, find keys matching given value *)
(* u. in jtCheck *)
reverseLookup[multimap_,val_]:=Select[indices[multimap],val∈multimap[#]&];

Junction Build

From clean-junction/clean-junction2.nb
In : g, adjMat (sparse array), V (vertex list in VertexList order)
   Out : jnodes, jedges,elimed
Entry point : buildJTree


Modified from clean-junction.nb to use vert2idx

vert2idx[verts_List]:=Flatten[Position[V,Alternatives@@verts]];
idx2vert[idx_List]:=V[[idx]];

touchingVerts[{}]:=V;
touchingVerts[v_List]:=(
vArr=SparseArray[vert2idx[v]->ConstantArray[1,\[LeftBracketingBar]v\[RightBracketingBar]],\[LeftBracketingBar]V\[RightBracketingBar]];
idx2vert[Cases[ArrayRules[(adjMat.vArr)(1-vArr)],HoldPattern[{x_}->_?Positive]->x]]
);
touchingVerts[v_Integer]:=touchingVerts[{v}];

connComps[v_]:=ConnectedComponents[Subgraph[g,v]];

(* approximates the number of cycles going through each point *)
computeAccumulation[mat_]:=(
r=Max@Eigenvalues@N@mat;
Diagonal[Inverse[IdentityMatrix[Length[mat]]-mat/(r*1.1)]]//Normal
);

gridCentralRow[k_]:=(Table[k(Floor[k/2])+j,{j,1,k}]);
gridCentralCol[k_]:=(Table[k j+Floor[k/2]+1,{j,0,k-1}]);


(* needs g, adjMat as sparseArray and V. Sets jnodes/jedges *)
buildJTree[method_:"minfill",tieBr_:"default"]:=(
Assert[Head[V]==List,Head[g]==Graph,Head[adjMat]==SparseArray];
elimed={}; (* elimination sequence *)
jnodes={};
tieBreaker=V;

(* how to pick among vertices that give equal fill *)
Which[
tieBr=="default",tieBreaker=V,
tieBr=="accumulation",tieBreaker=computeAccumulation[adjMat],
tieBr=="eccentricity",tieBreaker=VertexEccentricity[g,#]&/@V,
tieBr=="random",tieBreaker=RandomSample[V],
True,Assert["unknown tiebreaker"]
];

Which[
method=="minfill",While[\[LeftBracketingBar]elimed\[RightBracketingBar]<\[LeftBracketingBar]V\[RightBracketingBar],minFillEliminate[]],
method=="mincut",While[\[LeftBracketingBar]elimed\[RightBracketingBar]<\[LeftBracketingBar]V\[RightBracketingBar],minCutEliminate[]],

(* experimental methods for grids *)
method=="grid",
k=Floor[Sqrt[Length[V]]];
row=gridCentralRow[k];
col=gridCentralCol[k]∖row;
rowOrd=Ordering[VertexEccentricity[g,#]&/@row];
colOrd=Ordering[VertexEccentricity[g,#]&/@col];
elimVertex/@row[[rowOrd]];
elimVertex/@col[[colOrd]];
While[\[LeftBracketingBar]elimed\[RightBracketingBar]<\[LeftBracketingBar]V\[RightBracketingBar],minFillEliminate[]],

method=="grid2",
elimVertex/@{3,8,13,14,15,18,23,12,11};
While[\[LeftBracketingBar]elimed\[RightBracketingBar]<\[LeftBracketingBar]V\[RightBracketingBar],minFillEliminate[]],

method=="grid3",
elimVertex/@Ordering[VertexEccentricity[g,#]&/@V],

True,Assert["unknown method"]
];
jedges=makeSpanningBagTree[jnodes];
);

isBagRedundant[bag_]:=anyTrue[b∈jnodes,(b≠bag)∧(bag⊂b)];
elimVertex[v_]:=(
(*Print["Eliminating ",v," remain ",Length[V∖elimed]];*)
Assert[!(v∈elimed),"Duplicate elimination"];
comps=connComps[V∖elimed];
ourComp=Select[comps,v∈#&]//First;
boundary=touchingVerts[ourComp];
AppendTo[jnodes,boundary⋃{v}];
AppendTo[elimed,v];
jnodes=Select[jnodes,Not[isBagRedundant[#]]&];
);

(* returns negative number of edges in bag created for each member of verts *)
getFillValues[{}]:={};
getFillValues[verts_List]:=(
Assert[ValueQ[elimed]];
If[elimed=={},Return[ConstantArray[0,\[LeftBracketingBar]verts\[RightBracketingBar]]]];
elimArr=SparseArray[vert2idx[elimed]->ConstantArray[1,\[LeftBracketingBar]elimed\[RightBracketingBar]],\[LeftBracketingBar]V\[RightBracketingBar]];
getFillSingleVert[v_]:=(
vArr=SparseArray[vert2idx[{v}]->1,\[LeftBracketingBar]V\[RightBracketingBar]];
touchArr=adjMat.vArr;
cliqueArr=elimArr touchArr+vArr;
-cliqueArr.adjMat.cliqueArr/2
);
getFillSingleVert/@verts
);

minFillEliminate[]:=(
frontier=touchingVerts[elimed];
Assert[frontier≠{},"Either nothing left to eliminate or graph not connected"];
fillVals=getFillValues[frontier];
tieVals=tieBreaker[[frontier]];
elimVertex[First[selectSmallest[frontier,fillVals,tieVals]]]
);

minCutEliminate[]:=(
Assert[V∖elimed≠{},"nothing left to eliminate"];
comps=ConnectedComponents[Subgraph[g,V∖elimed]];
largestComp=First[Sort[comps,Length[#1]≥Length[#2]&]];
If[\[LeftBracketingBar]largestComp\[RightBracketingBar]<4,minFillEliminate[],minCutElimComp[largestComp]]
);

(* Splits components into 2 and eliminates the boundary in minfill order *)
Needs["GraphUtilities`"]
$ContextPath=DeleteCases[$ContextPath,"Combinatorica`"|"GraphUtilities`"];
minCutElimComp[comp_]:=(
children=GraphUtilities`MinCut[Subgraph[g,comp],2];
children=Sort[children,Length[#1]≤Length[#2]&];
compBoundary=touchingVerts[children[[1]]]⋂children[[2]];
remainBoundary=compBoundary;
While[remainBoundary≠{},
fillVals=getFillValues[remainBoundary];
tieVals=tieBreaker[[remainBoundary]];
cands=selectSmallest[remainBoundary,fillVals,tieVals];
elimVertex[First[cands]];
remainBoundary=remainBoundary∖elimed;
];
);

initBuildJTree[gr_]:=(
adjMat=AdjacencyMatrix[gr];
g=AdjacencyGraph[adjMat];
adjMat=AdjacencyMatrix[g];
V=VertexList[g];
);

(* checks that "jnodes/jedges" forms a junction tree of g *)
checkJunctionTree[g_,jnodes_,jedges_]:=(
Assert[Union[VertexList[g]]==Union[Flatten[jnodes]],"Variable preservation"];
Assert[allTrue[ed∈List@@@EdgeList[g],anyTrue[bag∈jnodes,ed⊂bag]],"Factor preservation"];
isSeparator[{i_,j_}]:=\[LeftBracketingBar]ConnectedComponents[Subgraph[g,V∖(i⋂j)]]\[RightBracketingBar]>1;
Assert[allTrue[edge∈jedges,isSeparator[edge]],"Running intersection"];
True
);

Domain

(* initializes: g,adjMat,V, cD,cX,cY, messages *)
Clear[initHardcore,findMaximalDomains];
Clear[g,V,adjMat,edgeDomains,vertexDomains,lField,nField,field,msg];
Clear[mdom,ch];
Clear[δ,τ,cavityGet];

graphStyle="SimpleLink";

(* finds domains not included in any other *)
findMaximalDomains[domains_List]:=(
covered[dom_]:=anyTrue[b∈domains,b≠dom∧dom⊂b];
Select[domains,!covered[#]&]
);

(* makes each factor a clique *)
makeInteractionGraph[domains_List]:=(
maxCliques=findMaximalDomains[domains];
edges=fl1[Subsets[#,{2}]&/@maxCliques];
Graph[UndirectedEdge@@@edges]
);

(* initialize interaction graph, local domains *)
initHardcore[graph_Graph,weight_]:=(
Assert[Head[graph]==Graph&&!DirectedGraphQ[graph]];
λ=weight;

(* interaction graph *)
originalEdges=List@@@EdgeList[graph];
edgeDomains=Sort/@originalEdges;
vertexDomains={#}&/@VertexList[graph];
domains=Sort[edgeDomains∼vertexDomains];
g=makeInteractionGraph[domains];
V=VertexList[g];
adjMat=AdjacencyMatrix[g];

(* local fields *)
Clear[cD,cX,cY];
initFieldEdge[edgeDom_]:=(
Assert[Head[edgeDom]==List];
cD[edgeDom]=edgeDom;
cX[edgeDom]=Subsets[edgeDom];;
cY[edgeDom]={1,1,1,0};
);
initFieldVertex[vertDom_]:=(
Assert[Head[vertDom]==List];
cD[vertDom]=vertDom;
cX[vertDom]=Subsets[vertDom];
cY[vertDom]={1,λ}
);
table[edgeD∈edgeDomains,initFieldEdge[edgeD]];
table[vertD∈vertexDomains,initFieldVertex[vertD]];
);

(*** Functions working with jedges/jnodes initialized ****)

(* assign domains to jnodes in sequence *)
assignDomains[]:=(
checkJunctionTree[g,jnodes,jedges];
unassigned=domains;
Clear[cl];
assignOneDomain[jnode_]:=(
coveredDomains=Select[unassigned,#⊂jnode&];
unassigned=unassigned∖coveredDomains;
cl[jnode]=coveredDomains;
);
assignOneDomain/@jnodes;
Assert[unassigned=={}];
);

(* checks that decomposition is a junction tree and every domains assigned to 1 bag *)
checkTreeDecomposition[]:=(
checkJunctionTree[g,jnodes,jedges];
Assert[allTrue[f∈domains,\[LeftBracketingBar]reverseLookup[cl,f]\[RightBracketingBar]==1],"factor preservation"];
);


mdom[k_]:=First[k]⋂Last[k];
ch[{i_List}]:=Cases[messages,{i,_}];
ch[{i_List,j_List}]:=Cases[messages,{j,Except[i]}];

initMessages[]:=(
checkJunctionTree[g,jnodes,jedges];
Assert[Sort[fl1[values[cl]]]==domains,"domains not assigned"];

messages=jedges∼(Reverse/@jedges);
initFieldMsg[msg_]:=(
{i,j}=msg;
atoms=i⋂j;
cD[msg]=atoms;
cX[msg]=Subsets[atoms];
cY[msg]=ConstantArray[1,\[LeftBracketingBar]Subsets[atoms]\[RightBracketingBar]];
);

table[msg∈messages,initFieldMsg[msg]];
);

(* message LHS *)
δ[msg_,x_]:=cavityGet[msg,x];
τ[factor_,x_]:=cavityGet[factor,x];
cavityGet[field_,x_]:=(
xx=x⋂cD[field];
i=idx[cX[field],xx];
If[i==-1,0,cY[field][[i]]]
);

(* message RHS *)
evalMsg[msg_,x_]:=(
{i,j}=msg;
sep=i⋂j;
lField[xx_]:=prod[dom∈cl[j],τ[dom,xx]];
nField[xx_]:=prod[m∈ch[msg],δ[m,xx]];
field[xx_]:=lField[xx]nField[xx];
Simplify[sum[y∈Subsets[j∖sep],field[x⋃y]]]
)/;msg∈messages;

evalClique[j_,x_]:=(
lField[xx_]:=prod[dom∈cl[j],τ[dom,xx]];
nField[xx_]:=prod[m∈ch[{j}],δ[m,xx]];
field[xx_]:=lField[xx]nField[xx];
field[x]
)/;j∈jnodes;

showMsg[msg_]:=Column[{#,evalMsg[msg,#]}&/@cX[msg]];
showClique[j_]:=Column[{#,evalClique[j,#]}&/@Subsets[j]];

updateMsg[msg_]:=(
newY=table[x∈cX[msg],evalMsg[msg,x]];
cY[msg]=newY;
);

treeUpdate[]:=(
checkTreeDecomposition[];
messageDeps=Table[{#\[DirectedEdge]p&/@ch[p]},{p,messages}]//Flatten;
messageGraph=Graph[messages,messageDeps,VertexLabels->Placed["Name",Tooltip]];

If[messageDeps≠{},(* TopologicalSort doesn't work on empty graphs *)
updateMsg/@TopologicalSort[messageGraph],
updateMsg/@messages
];
);

independencePoly[gr_]:=(
initHardcore[gr,x];
buildJTree["minfill"];
(* for planar graphs, use MinCut heurstic *)
(*buildJTree["mincut"];*)
assignDomains[];
initMessages[];
treeUpdate[];
evalClique[First[jnodes],#]&/@Subsets[First[jnodes]]//Total//Simplify
);

Visualize

(*isRedundant[bag_]:=AnyTrue[{b,bags},b≠bag∧bag⊂b];*)

(*** Graph Layout ***)
correctPR[pr_,scale_]:={pr[[1]]+{-1/scale,2/scale},pr[[2]]+{-2/scale,1/scale}};
getPlotRange[lst_]:=({xs,ys}=Transpose[lst];
(*two extra 1/scale offsets needed for exact match*){{Min[xs]-lineThickness/2,Max[xs]+lineThickness/2+1/masterScale},{Min[ys]-lineThickness/2-1/masterScale,Max[ys]+lineThickness/2}}
);

(*** Custom drawing ***)
drawGraphFragment[bag_,scale_:masterScale]:=(
bagGraph=Subgraph[g,bag];
allPoints=bag/.vcoords;
plotRange=getPlotRange[allPoints];
imSize=-scale*Subtract@@@plotRange;

Subgraph[g,bag,PlotRange->plotRange,GraphStyle->"SimpleLink",EdgeStyle->Transparent,ImageSize->imSize,VertexCoordinates->vcoords]
);
drawBag[bag_List,scale_:masterScale]:=(
bagGraph=Subgraph[g,bag,VertexCoordinates->vcoords,GraphStyle->graphStyle];
tileEnds=Select[Outer[List,bag,bag]//fl1,#[[1]]≤#[[2]]&];
tiles=FindShortestPath[CompleteGraph[Length[V]],#[[1]],#[[2]]]&/@tileEnds;

bagColor=Extract[colorList,Position[jnodes,bag]]//First;
tileOpts={bagColor,AbsoluteThickness[lineThickness*scale],CapForm["Round"],JoinForm["Round"],Opacity[.8]};
allPoints=Flatten[tiles/.vcoords,1];
plotRange=getPlotRange[allPoints];
imSize=-scale*Subtract@@@plotRange;

tileG=Graphics[tileOpts~Join~ (Line[#/.vcoords]&/@tiles),PlotRange->plotRange,ImageSize->imSize]
);

showBags[]:=(
(*Assert[Union[V]==Union[Flatten[jnodes]]];*)
Show[Sequence@@(drawBag/@jnodes),g,PlotRange->graphPR,ImageSize->graphIS]
);
Clear[showTree];
showTree[fixCorners_:True]:=Module[{},
bags=fl1[jedges];
Assert[\[LeftBracketingBar]jedges\[RightBracketingBar]>0];
bagCentroid[bag_]:=Mean[bag/.vcoords];

findExtremeBag[vec_]:=(
vertList=First/@vcoords;
coordList=Last/@vcoords;
goodBags=DeleteCases[bags,{}];
extremePos=First[Ordering[goodBags,1,bagCentroid[#1].vec>bagCentroid[#2].vec&]];
goodBags[[extremePos]]
);

extremeDirs={{1,1},{1,-1},{-1,1},{-1,-1}};
extremeBags=findExtremeBag/@extremeDirs;
extremePoses=bagCentroid/@extremeBags;

(*** Show bags as vertices of another GraphPlot ***)
If[fixCorners,
scale=masterScale/1.5;
GraphPlot[Rule@@@jedges,
VertexRenderingFunction->(Inset[Show[drawBag[#2,scale],drawGraphFragment[#2,scale]],#]&),
EdgeRenderingFunction->({Gray,Thick,Arrowheads[.0],Arrow[#1,.0]}&),
VertexCoordinateRules->Thread[Thread[extremeBags->extremePoses]],
PlotRange->graphPR,
ImageSize->1.5graphIS],
scale=masterScale/3;
GraphPlot[Rule@@@jedges,
VertexRenderingFunction->(Inset[Show[drawBag[#2,scale],drawGraphFragment[#2,scale]],#]&),
EdgeRenderingFunction->({Gray,Thick,Arrowheads[.0],Arrow[#1,.0]}&),
ImageSize->Large]
]
];

(* called from initGraph and initGraphNamed *)
initGraph0[gr_Graph]:=(
masterScale=50;
lineThickness=.5;
vertSize=1/2;
V=VertexList[gr];

vcoords=Thread[VertexList[gr]->(VertexCoordinates/.AbsoluteOptions[gr,VertexCoordinates])];

graphPR=getPlotRange[V/.vcoords];
graphPR=correctPR[graphPR,masterScale];
graphIS=-masterScale Subtract@@@graphPR;
graphStyle="DiagramBlack";
colorList=ColorData["Atoms","ColorList"];

makeVertex[{x_,y_},v_]:=(
plotRange=Transpose[{{x,y}-vertSize/2,{x,y}+vertSize/2}];
plotRange=correctPR[plotRange,masterScale];
imageSize=-masterScale Subtract@@@plotRange;
label={Black,Style[Text[v,{x,y}],Bold,14]};
label∼{Transparent,EdgeForm[Thick],Button[ Disk[{x,y},vertSize/2],Beep[]]}
);
shapeFunc[{x_,y_},v_,{w_,h_}]:=makeVertex[{x,y},v];
g=Graph[EdgeList[gr],VertexCoordinates->vcoords,PlotRange->graphPR,ImageSize->graphIS,VertexShapeFunction->shapeFunc]
);

initGraph[gr_Graph]:=initGraph0[gr];
initGraphNamed[graphName_]:=(
Check[GraphData[graphName],Assert[False,"graph not found ",graphName]];
initGraph0[GraphData[graphName]]
);

Usage

independencePoly[GridGraph[{3,3}]]==GraphData[{"Grid",{3,3}},"IndependencePolynomial"][x]

independencePoly[GridGraph[{6,6}]]

g=GridGraph[{3,3}];
initGraph[g]