Solving this question was quite challenging before I search journals and found out it's an NP-Complete (Source) and surprising it's has three answers in Code Golf Community (people compete usually on fewer lines-of-code solutions).
My first few attempts were using iterative approaches, simplifying the graph and solving it gradually. For example, if there is a 1-1 edge, then if you haven't a 2-node graph, it would never be connected to each other and similarly, nodes with capacity 8 should have 2 bridges on each side.
But soon or later you will reach a point in which you have nothing to simplify based on their capacity and/or number of edges. Sometimes you have to remove edges because it's obvious that edge will make the graph not connected (or inter-connected ;) or more specifically WeaklyConnectedGraphQ) that's where it becomes hard.
This solution will first simplify the graph (optional), then uses Google OR-Tools's constraint solver to find possible trees that cover all the vertices and satisfy the vertex capacity, then we check whether a combination is good or search for another.
Why Google OR-Tools when Mathematica has an entire guide page on optimization, well because in this particular case I couldn't use FindInstance as efficiently as OR-Tools which has a native library and has been used in the production environment.
Before we start, we need these 3 java libraries (all of them are from google):
Loading OR-tools via JLink + defining a class
<< JLink`;
AddToClassPath["C:\\protobuf-java-3.19.4.jar"];
AddToClassPath["C:\\ortools-win32-x86-64-9.3.10497.jar"];
AddToClassPath["C:\\ortools-java-9.3.10497.jar"];
LoadJavaClass["com.google.ortools.Loader"];
Loader`loadNativeLibraries[]
Hashiwokakero`Private`findCombination=ExternalEvaluate["Java","
import com.google.ortools.constraintsolver.Solver;
import com.google.ortools.constraintsolver.IntVar;
import com.google.ortools.constraintsolver.DecisionBuilder;
class find_combinations{
Solver solver;
IntVar[] edges_to_solve;
public find_combinations(int[][] edges_weight_min_max,int[] vertex_capacity,int[][] connections,int[][] intesections){
this.solver = new Solver(\"test\");
this.edges_to_solve = new IntVar[edges_weight_min_max.length];
// add min and max for each edge weight
for(int i =0;i < edges_weight_min_max.length;i++){
int min = edges_weight_min_max[i][0];
int max = edges_weight_min_max[i][1];
if (min==max){
this.edges_to_solve[i]=this.solver.makeIntConst(min);
}else{
this.edges_to_solve[i]=this.solver.makeIntVar(min, max);
}
}
// add vertex capacity constrain
for(int i =0;i < connections.length;i++){
IntVar[] temp;
if (connections[i][3]!=0){
temp=new IntVar[4];
}else if (connections[i][2]!=0){
temp=new IntVar[3];
}else if (connections[i][1]!=0){
temp=new IntVar[2];
}else{
temp=new IntVar[1];
}
for(int j =0;j < 4;j++){
if (connections[i][j]==0){
break;
}
temp[j]=this.edges_to_solve[connections[i][j]-1];
}
this.solver.addConstraint(this.solver.makeSumEquality(temp, vertex_capacity[i]));
}
// add intersection constrain
for(int i =0;i < intesections.length;i++){
if (intesections[i].length == 2){
this.solver.addConstraint(this.solver.makeEquality(this.solver.makeIntConst(0),this.solver.makeProd(this.edges_to_solve[intesections[i][0]-1],this.edges_to_solve[intesections[i][1]-1])));
}
}
DecisionBuilder db =
this.solver.makePhase(this.edges_to_solve, this.solver.INT_VAR_SIMPLE, this.solver.INT_VALUE_SIMPLE);
this.solver.newSearch(db);
}
public long[] get(){
if (this.solver.nextSolution()){
long[] result= new long[this.edges_to_solve.length];
for(int i =0;i < this.edges_to_solve.length;i++){
result[i]=this.edges_to_solve[i].value();
}
return result;
}
this.solver.endSearch();
return new long[0];
}
public void stop(){
this.solver.endSearch();
}
}
"];
Defining Hashiwokakero`simplify, Hashiwokakero`buildGraph and Hashiwokakero`solve
BeginPackage["Hashiwokakero`", {"JLink`"}];
ClearAll[simplify,buildGraph,solve];
Begin["`Private`"];
ClearAll[edgeSafeDelete,nodeWeightKIndices,selectIntersections,checkPossibleBridgesOverlaps,buildGraph];
edgeSafeDelete[graph_,{}]:=graph;
edgeSafeDelete[graph_,edges_]:=Block[{temp=Intersection[Join[#,Reverse/@#]&@EdgeList[graph],DeleteDuplicatesBy[edges,Sort]]},If[temp=!={},EdgeDelete[graph,temp],graph]]
nodeWeightKIndices[graph_,k_]:=Catenate@Position[AnnotationValue[graph,VertexWeight],k,{1}];
selectIntersections[{{x1_,y1_},{x2_,y2_}},listOfVectors_,edges_]:=Block[{data=Flatten[listOfVectors,{2,3}],temp},
temp=Transpose[{Sign[data[[1]]-x1]*Sign[x2-data[[1]]]+Sign[y1-data[[2]]]*Sign[data[[4]]-y1],Sign[x1-data[[1]]]*Sign[data[[3]]-x1]+Sign[data[[2]]-y1]*Sign[y2-data[[2]]]}];
edges[[Catenate@Position[temp,{2,_}|{_,2},{1}]]]
]
checkPossibleBridgesOverlaps[islands_,possibleBridges_]:=Table[
Catenate[Values@GroupBy[possibleBridges[[i]],Switch[islands[[#,2]]-islands[[i,2]],{_?Positive,0},1,{_?Negative,0},3,{0,_?Positive},2,{0,_?Negative},4]&,TakeSmallestBy[EuclideanDistance[islands[[i,2]],islands[[#,2]]]&,1]]],{i,Length@islands}]
buildGraph[coordinatesAndCapacity_]:=Block[{data=Transpose@Prepend[Transpose@coordinatesAndCapacity,Range@Length@coordinatesAndCapacity],data1PB,data1GraphRules,currentGraph},
data1PB=checkPossibleBridgesOverlaps[data,Catenate@Position[data,{_,{Except[#1],#2},_}|{_,{#1,Except[#2]},_},{1}]&@@@data[[All,2]]];
data1GraphRules=Catenate@MapIndexed[Thread[{First@#2,#1}]&,data1PB];
currentGraph=Graph[DeleteDuplicatesBy[data1GraphRules,Sort]
,VertexCoordinates->data[[All,2]]
,VertexLabels->Thread[Range[Length@data]->data[[All,3]]]
,VertexCapacity->data[[All,3]]
,VertexWeight->data[[All,3]]
];
Do[AnnotationValue[{currentGraph,edge},"EdgeWeight"]={0,1,2};,{edge,EdgeList[currentGraph]}];
currentGraph
]
checkPossibleBridgesOverlaps[islands_,possibleBridges_]:=Table[
Catenate[Values@GroupBy[possibleBridges[[i]],Switch[islands[[#,2]]-islands[[i,2]],{_?Positive,0},1,{_?Negative,0},3,{0,_?Positive},2,{0,_?Negative},4]&,TakeSmallestBy[EuclideanDistance[islands[[i,2]],islands[[#,2]]]&,1]]],{i,Length@islands}]
simplify[graph_]:=Block[{currentGraph=graph},
(*remove 1-1 bridge - one time*)
Block[{oneNodes=nodeWeightKIndices[currentGraph,1]},
currentGraph=EdgeDelete[currentGraph,EdgeList[currentGraph,(Alternatives@@oneNodes)\[UndirectedEdge](Alternatives@@oneNodes)]];
];
(*cap 2-2 edge weight - one time*)
Block[{oneNodes=nodeWeightKIndices[currentGraph,2]},
Do[AnnotationValue[{currentGraph,edge},"EdgeWeight"]={0,1};,{edge,EdgeList[currentGraph,(Alternatives@@oneNodes)\[UndirectedEdge](Alternatives@@oneNodes)]}];
];
(*assing {0,1} to all bridges to/from 1 capacity islands - one time*)
Block[{verticesWithOneCapacityIndex=nodeWeightKIndices[currentGraph,1]},
Do[AnnotationValue[{currentGraph,i},"EdgeWeight"]={0,1};,{i,EdgeList[currentGraph,(Alternatives@@verticesWithOneCapacityIndex)\[UndirectedEdge]_]}];
];
(*find all intersections *)
Block[{edges=EdgeList[currentGraph]
,edgesCoordinates=Partition[Extract[AnnotationValue[currentGraph,VertexCoordinates],List/@Catenate[List@@@EdgeList[currentGraph]]],2]},
Do[AnnotationValue[{currentGraph,edges[[i]]},"Intersections"]=selectIntersections[edgesCoordinates[[i]],edgesCoordinates,edges];,{i,Length@edges}];
];
(*filter ceratain edges like a node with capacity 5 and 3 bridges, surely has 3 bridges*)
Do[Block[{vertices=Select[nodeWeightKIndices[currentGraph,k],EdgeCount[currentGraph,#\[UndirectedEdge]_]==Ceiling[k/2]&],evenQ=EvenQ[k],edges},
edges=EdgeList[currentGraph,(Alternatives@@vertices)\[UndirectedEdge]_];
edges=Pick[edges,Length/@AnnotationValue[{currentGraph,edges},"EdgeWeight"],2|3];
AnnotationValue[{currentGraph,edges},"EdgeWeight"]=If[evenQ,2,Select[#>0&]/@AnnotationValue[{currentGraph,edges},"EdgeWeight"]];
currentGraph=edgeSafeDelete[currentGraph,Catenate@AnnotationValue[{currentGraph,edges},"Intersections"]];
];,{k,3,8}];
(*update all the intersections*)
Block[{edges=EdgeList[currentGraph]},
Do[AnnotationValue[{currentGraph,edge},"Intersections"]=Intersection[AnnotationValue[{currentGraph,edge},"Intersections"],Join[edges,Reverse/@edges]];,{edge,edges}];
];
currentGraph
]
solve[g_]:=Module[{graph=Hashiwokakero`simplify[g],edges,nodesCapacity,connections,intersections,edgesWeightMinMax,counter=1,combination},
edges=EdgeList[graph];
edgesWeightMinMax=MinMax/@AnnotationValue[{graph,edges},"EdgeWeight"];
nodesCapacity=AnnotationValue[graph,VertexCapacity];
connections=Table[Catenate@Position[edges,i\[UndirectedEdge]_|_\[UndirectedEdge]i,{1}],{i,VertexCount[graph]}];
connections=PadRight[connections,{VertexCount[graph],4},0];
intersections=Catenate@Table[Thread@{EdgeIndex[graph,edge],EdgeIndex[graph,#]&/@AnnotationValue[{graph,edge},"Intersections"]},{edge,edges}];
If[intersections=={},intersections={{}}];
intersections=DeleteDuplicatesBy[intersections,Sort];
JavaBlock[Module[{temp},
temp=JavaNew[findCombination,edgesWeightMinMax,nodesCapacity,connections,intersections];
combination=temp@get[];
If[combination==={},Return[{}]];
While[Not@WeaklyConnectedGraphQ@EdgeDelete[graph,Pick[edges,combination,0]],
counter+=1;
combination=temp@get[];
If[combination==={},Return[{}]];
];
temp@stop[];
]];
Print["Found answer in ",counter," iteration."];
EdgeDelete[EdgeAdd[graph,Pick[edges,combination,2]],Pick[edges,combination,0]]
]
End[];
EndPackage[];
How does it work?
Input format should be {{{x1,y1},capacity1}, ..., {{xN,yN},capacityN}}.
assuming your sample data:
data={{{0, 1}, 3}, {{1, 1}, 2}, {{1/3, 2/3}, 1}, {{2/3, 2/3}, 2}
,{{1, 1/3}, 1}, {{0, 0}, 4}, {{2/3, 0}, 3}};
First, build the graph (it also defines other properties):
buildGraph[data]
Output:
Then, apply solve to find a solution:
solve[buildGraph[data]]
(* Print: "Found answer in 1 iteration." *)
Output:
For more complex cases, you may want to simplify the graph first:
Under the hood
When you build the graph with buildGraph it will define the VertexCoordinates, VertexCapacity, VertexWeight and "EdgeWeight" + intersections for each edge under "Intersections" property.
Applying simplify, will:
- Remove
1-1 cases
- Limit edges connected to nodes with a capacity of 1 to a maximum of 1
- Limit
2-2 cases to the maximum 1 capacity
- For odd capacities, like nodes with capacity 3 and 2 edges, will raise the minimum edge capacity to 1
- For even capacities, like node with capacity 4 and 2 edges, will assign 2 for
"EdgeWeight" for each edge
For finding a solution, a list of each edge min, max weight, nodes capacities, connections for each node, and edges intersections will be sent and Google OR-Tools will find combinations of edge weight that satisfy node capacity and prevent intersection. It will send a combination to Wolfram-part and it checks whether the graph with the combinations is weakly connected or not. If it's not, it will ask for the next combination, if not, the search is ended.
You could manipulate the code to give you all the answers.
For the above complex case, it took 70 iterations to find the answer and that was under 1 second (~0.5 more specifically).
I've tried other approaches to tackle this challenge, but wolfram's efficiency in graph computation is not ready to compete with OR-Tools constraint solver at least in my testing on this challenge.
Bonus
One of the websites I found useful is menneske which could generate Hashi puzzles with different sizes and difficulties. It could be amusing but what if we were stuck? How can we bring that puzzle in Mathematica? well here is the javascript snippet:
result=[];
n=document.querySelectorAll('div.hashi table tr').length;
document.querySelector('div.hashi table').querySelectorAll('tr').forEach(function(row,row_index){
row.querySelectorAll('td').forEach(function(cell,col_index){
if(cell.innerText!==""){
result.push([[col_index,n-row_index],parseInt(cell.innerText)]);
}
})
})
JSON.stringify(result)
Run the script above in the console and use ImportString[...,"JSON"] to parse it and use it as input to the defined functions.
{a <-> b, a <-> f, a <-> d, b <-> e, c <-> f, d <-> g, f <-> g, f <-> g}you get another answer. $\endgroup${a <-> d, a <-> f, a <-> f, b <-> e, c <-> d, b <-> g, f <-> g, f <-> g}gives yet another. $\endgroup$