(*Mathematica code to calculate Bakry Emery curvature*)
(*Copy code to Wolfram Cloud or Mathematica, then Shift + Enter to execute*)
(*--------------------------------------------------------------------------------------------------*)
(*Random weighted directed Graph Generation*)
RandAdj[VV_,EE_] := SparseArray[RandomInteger[{1, VV}, {EE, 2}] -> RandomReal[1,EE],{VV,VV}];
AdjToLaplace[A_] := A - DiagonalMatrix[A.SparseArray[Table[1,Length[A]]]];
RandL[VV_,EE_]:=AdjToLaplace[RandAdj[VV,EE]];
SeedRandom[123];
L=RandL[1000,20000];
(*RandL gives a directed graph with 1000 nodes and around 20000 edges*)
(*--------------------------------------------------------------------------------------------------*)
(*Ricci calculation given Laplacian L*)
(*L must be a square matrix with non negative off diagonal and row sums zero*)
L2=L.L;
A=SparseArray[L-DiagonalMatrix[Diagonal[L]]];
VV=Length[L];
b2all=Table[Flatten[L2[[k]]["NonzeroPositions"]],{k,VV}];
b1=Table[Flatten[A[[k]]["NonzeroPositions"]],{k,VV}];
b2=Table[Complement[b2all[[k]],b1[[k]],{k}],{k,VV}];
fg=SparseArray[Table[{i,i,i}->1,{i,VV}]];
T[ten_]:=TensorTranspose[ten,{1,3,2}];
Gam=(L.fg - fg.L - T[fg.L]);
Gam2=(L.Gam - Gam.L - T[Gam.L])/2;
inv := Compile[{{x, _Real}},If[x!=0,1/x,0]];
Inv[x_]:=Map[inv,x,{2}] ;
Q[k_,b1_,b2_,G2_]:=If[b2=={},G2[[b1,b1]],G2[[b1,b1]]-G2[[b1,b2]].Inv[G2[[b2,b2]]] . G2[[b2,b1]]];
Ricc[k_]:=If[b1[[k]]=={},0,Min[Eigenvalues[{Q[k,b1[[k]],b2[[k]],Gam2[[k]]],Gam[[k]][[b1[[k]],b1[[k]]]]}]]];
Ric=Table[Ricc[k],{k,VV}];
(*Ric gives the list of the Bakry Emery curvature of every node*)
(*--------------------------------------------------------------------------------------------------*)
(*Outputs average Ricci curvature*)
Mean[Ric]