Code

Optimization

In[]:=
MinSumPermutation[mm_?MatrixQ]:=Module[{n,m,objective,constraints,vars,res,p},​​{n,m}=Dimensions[mm];objective=Total[Inactive[Times][p,mm],2];constraints={Table[Sum[Indexed[p,{i,j}],{j,m}]1,{i,n}],Table[Sum[Indexed[p,{i,j}],{i,n}]≤1,{j,m}],0p1};vars={Element[p,Matrices[{n,m},Integers]]};res=LinearOptimization[objective,constraints,vars];Flatten[SparseArray[p/.res]["NonzeroPositions"][[All,2]]]]

GeodesicsBundle

In[]:=
GeodesicsBundle[graph_,endpoints:{_Integer,_Integer},neighborhoodSize_:3,plottingFunction_:GeodesicsAll]:=Module[{neighborhoods,distanceMatrix,matchingEndpoints},​​neighborhoods=VertexList[NeighborhoodGraph[graph,#,neighborhoodSize]]&/@endpoints;​​distanceMatrix=Outer[GraphDistance[graph,#1,#2]&,neighborhoods〚1〛,neighborhoods〚2〛];​​matchingEndpoints=If[Dimensions[distanceMatrix]〚1〛≥Dimensions[distanceMatrix]〚2〛,​​Transpose[{neighborhoods〚1,MinSumPermutation[Transpose@distanceMatrix]〛,neighborhoods〚2〛}],Transpose[{neighborhoods〚1〛,neighborhoods〚2,MinSumPermutation[distanceMatrix]〛}]];​​plottingFunction[graph,matchingEndpoints]​​]
In[]:=
GeodesicsBundle[graph_,endpointFractions:{_Real,_Real},neighborhoodSize_:3]:=Module[{path,endpoints},​​path=Echo@FindShortestPath[graph,Sequence@@GraphAntipodes[graph]];​​endpoints=Echo@path〚Round[endpointFractionsLength[path]]/.01〛;​​GeodesicsBundle[graph,endpoints,neighborhoodSize]​​]
In[]:=
GeodesicsBundle[graph_,neighborhoodSize_:3]:=GeodesicsBundle[graph,{0.2,0.8},neighborhoodSize]

GeodesicsAll

In[]:=
ShortestPathsVertices[graph_,v1_,v2_]:=VertexList[graph]〚Position[MapThread[Plus,{GraphDistance[graph,v1],GraphDistance[graph,v2]}],GraphDistance[graph,v1,v2]]〚All,1〛〛
In[]:=
GeodesicsAll[g_Graph,vpairs:{{_,_}..}]:=HighlightGraph[g,Style[Subgraph[g,Catenate[ShortestPathsVertices[g,##]&@@@vpairs]],Thickness[.02],Red]]

Sphere

In[]:=
sphere=IndexGraph@MeshConnectivityGraph[DiscretizeGraphics[Sphere[],MaxCellMeasure{"Area".005}]];
GeodesicsBundle[sphere,{.3,.7},3]