print(`OlliRicci.txt: A Maple Package for Computing Ollivier Ricci Curvature of Graphs`):
print(`by Kaylee Weatherspoon`):
print(`for a list of procedures, type: Help();`):

Help:=proc():
if nargs=0 then 
print(`The functions are`):
print(`distributionx, costmatrix, wasserstein, wassermanualdist,OlliRicciEdge, EdgeCurvsList, OlliRicciEdgeMD,OlliRicci, OlliRicciW, OlliRicciWD, RicciCommDetect, labelcomms, Prune`):

elif nargs=1 and args[1]=distributionx then 
print(`distributionx is a helper function which inputs a graph as a Maple Graph object: Graph([1..n], {Eges}) and a vertex x and outputs [Neighbors of x, Vector] where the Vector is a column vector which is a uniform distribution on the neighbors of x. That is,  if x has n neighbors, the vector is a column vector of length n with all entries equal to 1/n`):
print(`Try distributionx(CycleGraph(5), 1)`):

elif nargs=1 and args[1]=costmatrix then
print(`costmatrix is a helper function which inputs a graph as a Maple Graph object: Graph([1..n], {Eges}) and two vertices x and y and computes the shortest paths between all neighbors of x and neighbors of y`):
print(`try costmatrix(CycleGraph(5), 1, 3)`):

elif nargs=1 and args[1]=wasserstein then 
print(`Given two vertices x,y, and a graph as a Maple Graph object: Graph([1..n], {Eges}),  wasserstein computes the Wasserstein/Earth Mover's distance between x and y in G, assuming a uniform distribution for x and y`):
print(`Try wasserstein(1, 3,CompleteGraph(5))`):

elif nargs=1 and args[1]=wassermanualdist then 
print(`The Inputs are (x, distx, y, disty, G), where x and y are vertices, and distx is of the form [distribution of x]. disty follows the same pattern, and G is a graph as a MapleGraph object: Graph([1..n], {Eges}). This serves the same purpose as wasserstein but allows the user to set the distribution to be non-uniform as long as the conditions of the LP are satisfied.`):
print(`Try wassermanualdist(1, [[2, 3, 4], [0.6, 0.3, 0.1]], 2, [[1, 3, 4], [0.2, 0.2, 0.6]], CompleteGraph(4))`):

elif nargs=1 and args[1]=OlliRicciEdge then
print(`Given G, a graph as a maple object Graph([1..n], {Eges}), two adjacent vertices x and y, and the distance or weight of the associated edge, OlliRicciEdge computes the Ricci Curvature of the edge xy.`):
print(`Try OlliRicciEdge(Graph([1,2,3,4,5], {{1,2}, {2,3}, {1,3}, {3,4}, {5, 1}}), 1, 3, 66)`):


elif nargs=1 and args[1]=OlliRicciEdgeMD then
print(`The input is (G, x, distx, y, disty, distance). Given G, which is a graph as a Maple object Graph([1..n], {Eges}), vertices x and y, their distributions given as [distribution of x] and the distance or weight of the edge, this outputs the Ollivier Ricci curvature of the edge xy`):
print(`Try OlliRicciEdgeMD(CompleteGraph(4), 1, [[2, 3, 4], [0.6, 0.3, 0.1]], 3, [[2, 1, 4], [0.2, 0.2, 0.6]], 2)`):

elif nargs=1 and args[1]=OlliRicci then 
print(`INPUT: Graph in Maple format:Graph([1..n], {Eges}). OUTPUT: The Ollivier Ricci curvature of G, which is the average of the Ollivier Ricci curvature over all edges.`):
print(`Try: OlliRicci(CompleteGraph(7))`):

elif nargs=1 and args[1]=OlliRicciW then 
print(`INPUT:Inputs a weighted graph in the form [n,{{u,v,weightuv}, {etc}}] OUTPUT: Ricci Curvature of Weighted Graph`):
 print(`Try:OlliRicciW([4, {{1, 2, 4}, {1, 3, 5}, {1, 4, 16}, {2, 3, 66}, {2, 4, 55}, {3, 4, 10}}])`):


elif nargs=1 and args[1]=OlliRicciWD then
print(`INPUT: (G1, D), where G1 is a weighted graph in the form [n,{{u,v,weightuv}, {etc}}] and D is a list of list of distributions, in order. The first entry in the list should be the list representing the distribution of the vertex called 1. Outputs the Ollivier Ricci Curvature of a weighted graph with manually chosen distributions, as long as there is a solution to the corresponding LP.`):
print(`Try OlliRicciWD([4, {{1, 2, 4}, {1, 3, 5}, {1, 4, 16}, {2, 3, 66}, {2, 4, 55}, {3, 4, 10}}], [[1/3, 1/3, 1/3], [1/3, 1/3, 1/3], [1/3, 1/3, 1/3], [1/3, 1/3, 1/3]])`):
print(`Try OlliRicciWD([4, {{1, 2, 4}, {1, 3, 5}, {1, 4, 16}, {2, 3, 66}, {2, 4, 55}, {3, 4, 10}}], [[0.6, 0.3, 0.1], [0.5, 0.25, 0.25], [0.3, 0.6, 0.1], [0.001, 0.899, 0.1]])`):


elif nargs=1 and args[1]=Prune then
print(`This inputs a graph G in Maple format and "prunes" the `):
print(`Try Prune(Graph([1,2,3,4,5], {{1,2},{1,3},{1,4},{3,4},{2,5},{1,5}}))`):


elif nargs=1 and args[1]=labelcomms then
print(`This is a helper function which inputs a graph G in Maple format and assigns a label to each vertex based on which connected component it is a member of.`):
print(`LabelComms(Graph([1,2,3,4,5], {{1,2},{2,3},{3,4},{4,5},{3,5}}))`):

elif nargs=1 and args[1]=RicciCommDetect then
print(`INPUT: (G, s) where G is a graph in Maple format and s is the size of the smallest community we want to detect. This uses OlliRicciEdgeD, LabelComms, and Prune to compute communities. It outputs a table which shows which  community each vertex belonds to.`):
print(`Try RicciCommDetect(Graph([1,2,3,4,5], {{1,2},{2,3},{3,4},{4,5},{3,5}}, 3))`):

elif nargs=1 and args[1]=EdgeCurvsList then
print(`INPUT: G, where G is a graph in Maple Format`):
print(`OUTPUT: A list of all edge curvatures and the corresponding edge for each.`):
print(`Try EdgeCurvsList(Graph([1, 2, 3, 4, 5], {{1, 2}, {1, 3}, {2, 3}, {3, 5}, {4, 5}}))`):



else 
print(`There is no help for`, args[1]):
fi:
end:



with(GraphTheory):
with(LinearAlgebra):
with(Optimization):



distributionx:=proc(G,x) local Neis, n, distvec, i:
    Neis:=Neighbors(G,x):
    n:=nops(Neis):
    distvec:= Vector(n);
#we assign a uniform distribution for convenience and simplicity
    distvec:=Vector([seq(1/n, i=1..n)]):
return([Neis, distvec]):
end:

#x is source, y is sink
costmatrix:=proc(G, x, y) local Nx, Ny, n, m, C, i, j, path:
    Nx:=Neighbors(G, x):
    Ny:=Neighbors(G,y):
    n:=nops(Nx):
    m:=nops(Ny):
    C:=Matrix(n,m):
    for i from 1 to n do 
        for j from 1 to m do 
            path:=ShortestPath(G, Nx[i], Ny[j]):
            if nops(path)=0 then 
                C[i,j]:=infinity:
            else 
                C[i,j]:=nops(path)-1
            fi:            
        od:
    od:
return(C);
end:

#wasserstein or earth mover's distance common choice for OR 
#a and b are distribution vectors, C is the cost matrix
wasserstein:=proc(x, y, G) local n,m, obj, distx, disty, Neisx, Neisy, C, flovars, constr, i, j, flovarslist, LPsol, wassersteindist:
    distx:=distributionx(G, x):
    disty:=distributionx(G, y):

    Neisx:=distx[1]:
    Neisy:=disty[1]:
    n:=nops(Neisx):
    m:=nops(Neisy):

    C:=costmatrix(G, x, y):

    if n=0 or m=0 then return("ERROR: isolated vertex"):
    fi:

    flovars:=Matrix(n,m,(i,j) -> Z[i,j]):

    obj:=add(add(C[i,j]*flovars[i,j], j=1..m),i=1..n):

    constr:={}:

 #outflow contstraints
    for i from 1 to n do 
        constr:=constr union {add(flovars[i,j], j=1..m)=distx[2][i]}:
    od:
 #inflow constraints
    for j from 1 to m do    
        constr:=constr union {add(flovars[i,j], i=1..n)=disty[2][j]}:
    od:
#nonnegativity
    for i from 1 to n do 
        for j from 1 to m do
            constr:=constr union {flovars[i,j]>=0}:
        od:
    od:
#formatting for LP solver
    flovarslist:=[seq(seq(flovars[i,j], j=1..m), i=1..n)]:

    LPsol:=LPSolve(obj, constr):
    wassersteindist:=LPsol[1]:
    RETURN(wassersteindist);

end:

#input distribution as [[neighbors], Vector(distribution)]:
wassermanualdist:=proc(x, distx, y, disty, G) local n,m, obj, Neisx, Neisy, C, flovars, constr, i, j, flovarslist, LPsol, wassersteindist:
   
    Neisx:=Neighbors(G,x):
    Neisy:=Neighbors(G,y):
    n:=nops(Neisx):
    m:=nops(Neisy):

    C:=costmatrix(G, x, y):

    if n=0 or m=0 then return("ERROR: isolated vertex"):
    fi:

    flovars:=Matrix(n,m,(i,j) -> Z[i,j]):

    obj:=add(add(C[i,j]*flovars[i,j], j=1..m),i=1..n):

    constr:={}:

 #outflow contstraints
    for i from 1 to n do 
        constr:=constr union {add(flovars[i,j], j=1..m)=distx[i]}:
    od:
 #inflow constraints
    for j from 1 to m do    
        constr:=constr union {add(flovars[i,j], i=1..n)=disty[j]}:
    od:
#nonnegativity
    for i from 1 to n do 
        for j from 1 to m do
            constr:=constr union {flovars[i,j]>=0}:
        od:
    od:
#formatting for LP solver
    flovarslist:=[seq(seq(flovars[i,j], j=1..m), i=1..n)]:

    LPsol:=LPSolve(obj, constr):
    wassersteindist:=LPsol[1]:
    RETURN(wassersteindist);

end:


#computing the Ollivier Ricci Curvature of an Edge
OlliRicciEdge:=proc(G, x, y, distance) local W:
    #distance is edge length/weight
    W:=wasserstein(x, y, G):

RETURN(1-(W/distance)):
end:


OlliRicciEdgeD:=proc(G, x, distx, y, disty, distance) local Neix, Neiy, C, W, DFXx, DFXy:
    #distance is edge length/weight
    W:=wassermanualdist(x, distx, y, disty, G):

RETURN(1-(W/distance)):
end:


EdgeCurvsList:=proc(G) local edge, beeeep:
beeeep:=[]:
for edge in Edges(G) do 
    beeeep:=[op(beeeep), [OlliRicciEdge(G, edge[1], edge[2], 1), {edge[1], edge[2]}]]:
od:
return(beeeep):
end:



#The Ollivier-Ricci Curvature of a Graph is the 
#Average of the Ollivier-Ricci Curvatures of the Edges

OlliRicci:=proc(G) local x, y, e, E, curvs, Size:
    E:=Edges(G):
    Size:=nops(E):
    curvs:=[]:
    for e in E do 
        x:=e[1]:
        y:=e[2]:
        curvs:=[op(curvs), OlliRicciEdge(G, x, y, 1)]:
    od:
RETURN(add(curvs)/Size);
end:




#Inputs a weighted graph in the form [n,{{u,v,weightuv}, {etc}}]
#Outputs Ollivier Ricci Curvature 

OlliRicciW:=proc(G1) local G, verts, wedges, edgecurvs, edges,i, f, e:
    verts:=[seq(i, i=1..G1[1])]:
    wedges:=G1[2]:
    edgecurvs:=[]:
#need to process G
    edges:={}:
    for f in wedges do 
        edges:= edges union {{f[1], f[2]}}:
    od:
    G:=Graph(verts, edges):
#build list of weighted curvatures 
    for e in wedges do 
        edgecurvs:=[op(edgecurvs), OlliRicciEdge(G, e[1], e[2], e[3])]:
    od:
    return(add(edgecurvs)/nops(wedges)):
end:


#allows manual distribution and manual edge weight
#D is list of list of distributions in order 
OlliRicciWD:=proc(G1, D) local G, verts, wedges, edgecurvs, edges,i, f, e:
    verts:=[seq(i, i=1..G1[1])]:
    wedges:=G1[2]:
    edgecurvs:=[]:
#need to process G 
    edges:={}:
    for f in wedges do 
        edges:= edges union {{f[1], f[2]}}:
    od:
    G:=Graph(verts, edges):
#build list of weighted curvatures 
    for e in wedges do 
        dist(e[1]):=[Neighbors(G, e[1]), Vector(D[e[1]])]:
        print(e, e[1], dist(e[1]));
        dist(e[2]):=[Neighbors(G, e[2]), Vector(D[e[2]])]:
        print(e, e[2], dist(e[2]));
        edgecurvs:=[op(edgecurvs), OlliRicciEdgeD(G, e[1], dist(e[1]), e[2], dist(e[2]), e[3])]:
    od:
    return(add(edgecurvs)/nops(wedges)):
end:



#Community Detection Algorithm from Sia/Jonckheree/Bogdan Paper

Prune:=proc(G) local Gcopy, ecurvs, x, y, edge, minEdge, minCurv, aff, u,v, Neis, e:
Gcopy:=CopyGraph(G):
ecurvs:=table():

for edge in Edges(Gcopy) do 
    x:=edge[1]:
    y:=edge[2]:
    ecurvs[[x,y]]:=OlliRicciEdge(Gcopy, x, y, 1):
od:

while true do
    minEdge:=NULL:
    minCurv:=0:

    for edge in Edges(Gcopy) do 
        x:=edge[1]:
        y:=edge[2]:
        if assigned(ecurvs[[x,y]]) and ecurvs[[x,y]]< minCurv then
            minCurv:=ecurvs[[x,y]]:
            minEdge:=[x,y]:
        fi:
    od:

    if minEdge =NULL then break: fi:

    x:=minEdge[1]:
    y:=minEdge[2]:
    RemoveEdge(Gcopy, x, y):

    ecurvs[[x,y]]:='undef':

    #optimizing: recomputing curvatures only for affected edges, not for all except removed
    for u in [x,y] do 
        Neis:=Neighbors(Gcopy, u):
        for v in Neis do 
            if v<u then 
                e:=[v,u]:
            else    
                e:=[u,v]:
            fi:
            ecurvs[[op(e)]]:=OlliRicciEdge(Gcopy, e[1], e[2], 1):
        od:
    od:
od:

RETURN(Gcopy):

end:


labelcomms:=proc(G) local conncomps, labels, i, conncomp, v:

conncomps:=ConnectedComponents(G):
labels:=table():

for i to nops(conncomps) do 
    conncomp:=conncomps[i]:
    for v in conncomp do
        labels[v]:=i:
    od:
od:

#table sending each node to a community
return(print(labels));

end:

#s is the size of the smallest community we care about
RicciCommDetect:=proc(G, s) local Gpruned, labels, conncomps, labels2, i, conncomp, v:

Gpruned:=Prune(G):

conncomps:=ConnectedComponents(Gpruned):
labels:=table():

i:=1:
for conncomp in conncomps do
    if nops(conncomp)>= s then 
        for v in conncomp do 
            labels[v]:= i:
        od:
    i:=i+1:
    fi:
od:
RETURN(print(labels)):
end: