# Earth Mover Distance/Wasserstein Distance Of Two Lists

I am in need of calculating the Earth Mover Distance between two one-dimensional discrete distributions with finite support (a.k.a "two lists of identical length whose entries add up to 1"). I found out that EMD is one of the distance functions for ImageDistance[], but in that implementation, the command creates a histogram from the images and calculates their distance; I already have a histogram.

How can I access the EMD directly, to calculate it between two arrays of identical dimension and size representing two discrete distributions? (Right now I only need 1D, but why be modest?)

I have implemented a greedy algorithm (move each "dirt pile" to the nearest "hole" available), but would like to make sure I have the exact value.

It seems as though there is no built-in functionality according to the wisdom of this crowd, so I bit the bullet and dove in to find a solution that is more "from scratch". In particular, I utilized the "LinearOptimization" method to solve the equivalent transport problem. I'll share my code here, in case someone will find it useful and or can improve upon it down the road.

distance[vec1_List,vec2_List,p_]:=Module[{d},
(* The distance between two points in the support *)
d=1.Norm[vec1-vec2,p];
Return[d]
]

coordinates[index_Integer,dimension_List]:=Module[{i=index-1,dim=dimension,m,coord},
(* Convert the index of the flattened list into the index of the tensor, give its dimensions *)
coord={};
While[Length[dim]>0,
m=Mod[i,dim[[-1]]];
PrependTo[coord,m+1];
i=(i-m)/dim[[-1]];
dim=Delete[dim,-1];
];
Return[coord]
]

distanceMatrix[vec1_List,vec2_List,dimension_List,p_]:=Module[{m,dim=Length[vec1]},
(* This matrix contains the distances between any two points of the support *)
m=Table[distance[coordinates[i,dimension],coordinates[j,dimension],p],{i,1,dim},{j,1,dim}];
Return[m]
]

wassersteinLinProg[vec1_List,vec2_List,dim_List,p_]:=Module[{metric,res,sol},
If[Product[dim[[i]],{i,Length[dim]}]==Length[vec1]&&Length[vec1]==Length[vec2]&&Total[vec1]==Total[vec2],
metric=distanceMatrix[vec1,vec2,dim,p];
res=LinearOptimization[Total[Inactivate[metric*r],2],{Total[r,{2}]==vec2,Total[r]==vec1,r\[VectorGreaterEqual]0},r\[Element]Matrices[Dimensions[metric]]][[1]];
sol=r/.res;
Return[Total[metric*sol,2]]
,
Return[Infinity]
];
]

(* Usage *)
v1={1,0,0,0,0,0,0,0};
v2={0,0,0,0,0,0,0,1};
wassersteinLinProg[v1,v2,{8},2]
wassersteinLinProg[v1,v2,{4,2},2]
wassersteinLinProg[v1,v2,{2,2,2},1]


Returns

7.
3.16228
3.