Dependencies

Globals:
masterScale,lineThickness,vertSize, vcoords, graphIS, graphPR
g, V, jnodes, jedges

temps: elim,bags,elimPartialOrder

Code

Utils

(* custom ∈ *)
On[Assert];
$AssertFunction=(Message[Assert::asrtf,#1];Print[#1];Abort[])&;

Needs["Notation`"];
InfixNotation[(NotationTemplateTag[∈]),MyElement];
InfixNotation[(NotationTemplateTag[∉]),NotMyElement];
(* add custom alias to all open notebooks *)
AddInputAlias[(NotationTemplateTag[∈]),"in",#]&/@Notebooks[];
AddInputAlias[(NotationTemplateTag[∉]),"notin",#]&/@Notebooks[];

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]];

connComps[v_]:=ConnectedComponents[Subgraph[g,v]]
touchingVerts[v_List]:=VertexList[NeighborhoodGraph[g,v]]∖v;
touchingVerts[v_Integer]:=touchingVerts[{v}];

(* 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]->#]&/@vals;
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 **)

(* orders list according to sequence of partial orders *)
orderListPartial[list_List,orders_List]:=(
combOrder=Transpose[orders[[All,list]]];
list[[Ordering[combOrder]]]
);

orderList[list_List,orders__]:=(
orderListPartial[list,{orders}]
);


(* 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 *)
makeSpanningTree[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]
];

Initialize

(* called from initGraph and initGraphNamed *)
initGraph0[gr_Graph]:=(
elimed={};
jnodes={};
V=Union[VertexList[gr]];
elimPartialOrder=V;

masterScale=50;
lineThickness=.5;
vertSize=1/2;

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]]
);

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]
]
]

Eliminate

isRedundant[bag_]:=anyTrue[b∈jnodes,(b≠bag)∧(bag⊂b)];
elim[v_]:=(
Print["Eliminating ",v];
If[v∈elimed,Return[]];
comps=connComps[V∖elimed];
ourComp=Select[comps,v∈#&]//First;
boundary=touchingVerts[ourComp];
AppendTo[jnodes,boundary⋃{v}];
elimed=elimed⋃{v};
jnodes=Select[jnodes,Not[isRedundant[#]]&];
);

smallestFill[]:=(
If[elimed=={},Return[First@Ordering[VertexEccentricity[g,#]&/@V]]];

cands=touchingVerts[elimed];
(* next one looks like minimum cardinality search *)
candLen[vert_]:=Length[touchingVerts[elimed⋃{vert}]];

(*candLen[vert_]:=(
comps=connComps[V∖elimed];
ourComp=Select[comps,vert∈#&]//First;
boundary=touchingVerts[ourComp];
vboundary=boundary⋂touchingVerts[vert];
existingEdges=EdgeList[Subgraph[g,{vert}⋃vboundary]];
((\[LeftBracketingBar]vboundary+1\[RightBracketingBar])^2)-\[LeftBracketingBar]existingEdges\[RightBracketingBar]
);*)

lens=candLen/@cands;

minCands=Extract[cands,Position[lens,_?(#==Min[lens]&)]];
(*Print["Redundancy: ",Count[lens,Min[lens]]," Boundary: ",Min[lens]];*)
(* minimum bandwidth not great *)
(*First[orderList[minCands,bandOrder]]*)
First[orderList[minCands,elimPartialOrder]]
)/;(V≠elimed);

(* returns jnodes, jedges, jcl *)
buildJTree[graph_]:=(
initGraph[graph];

While[elimed≠V,elim[smallestFill[]]];
jedges=makeSpanningTree[jnodes]
);

Usage

buildJTree[GridGraph[{3,3}]];
showTree[False]

SetDirectory[NotebookDirectory[]];
buildJTree[Import["./libtw/graphs/barley.dgf"]];
showBags[]

GraphPlot[Rule@@@jedges,VertexLabeling->True,ImageSize->Large]

ListPlot[Sort[Length/@jnodes]]

buildJTree[GridGraph[{3,2}]];
showTree[False]