summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authortlacombe <lacombe1993@gmail.com>2020-01-13 15:13:35 +0100
committertlacombe <lacombe1993@gmail.com>2020-01-13 15:13:35 +0100
commite807654f30363ec92c0d24b702bfd081ec8aae5a (patch)
tree90c9eb2519a9dc16d9cbc2eef1a5cdbb5e7bf3f6
parent4bcdd64974900302f420fb08435275cc8faa794a (diff)
update variable names, going for order and internal_p
-rw-r--r--src/python/doc/wasserstein_distance_user.rst2
-rw-r--r--src/python/gudhi/wasserstein.py34
-rwxr-xr-xsrc/python/test/test_wasserstein_distance.py34
3 files changed, 35 insertions, 35 deletions
diff --git a/src/python/doc/wasserstein_distance_user.rst b/src/python/doc/wasserstein_distance_user.rst
index 8862a5ce..32999a0c 100644
--- a/src/python/doc/wasserstein_distance_user.rst
+++ b/src/python/doc/wasserstein_distance_user.rst
@@ -30,7 +30,7 @@ Note that persistence diagrams must be submitted as (n x 2) numpy arrays and mus
diag1 = np.array([[2.7, 3.7],[9.6, 14.],[34.2, 34.974]])
diag2 = np.array([[2.8, 4.45],[9.5, 14.1]])
- message = "Wasserstein distance value = " + '%.2f' % gudhi.wasserstein.wasserstein_distance(diag1, diag2, internal_p=2., q=1.)
+ message = "Wasserstein distance value = " + '%.2f' % gudhi.wasserstein.wasserstein_distance(diag1, diag2, order=1., internal_p=2.)
print(message)
The output is:
diff --git a/src/python/gudhi/wasserstein.py b/src/python/gudhi/wasserstein.py
index 2acf93d6..aef54f64 100644
--- a/src/python/gudhi/wasserstein.py
+++ b/src/python/gudhi/wasserstein.py
@@ -23,12 +23,12 @@ def _proj_on_diag(X):
return np.array([Z , Z]).T
-def _build_dist_matrix(X, Y, q=2., internal_p=2.):
+def _build_dist_matrix(X, Y, order=2., internal_p=2.):
'''
:param X: (n x 2) numpy.array encoding the (points of the) first diagram.
:param Y: (m x 2) numpy.array encoding the second diagram.
:param internal_p: Ground metric (i.e. norm l_q).
- :param q: exponent for the Wasserstein metric.
+ :param order: exponent for the Wasserstein metric.
:returns: (n+1) x (m+1) np.array encoding the cost matrix C.
For 1 <= i <= n, 1 <= j <= m, C[i,j] encodes the distance between X[i] and Y[j], while C[i, m+1] (resp. C[n+1, j]) encodes the distance (to the p) between X[i] (resp Y[j]) and its orthogonal proj onto the diagonal.
note also that C[n+1, m+1] = 0 (it costs nothing to move from the diagonal to the diagonal).
@@ -36,13 +36,13 @@ def _build_dist_matrix(X, Y, q=2., internal_p=2.):
Xdiag = _proj_on_diag(X)
Ydiag = _proj_on_diag(Y)
if np.isinf(internal_p):
- C = sc.cdist(X,Y, metric='chebyshev')**q
- Cxd = np.linalg.norm(X - Xdiag, ord=internal_p, axis=1)**q
- Cdy = np.linalg.norm(Y - Ydiag, ord=internal_p, axis=1)**q
+ C = sc.cdist(X,Y, metric='chebyshev')**order
+ Cxd = np.linalg.norm(X - Xdiag, ord=internal_p, axis=1)**order
+ Cdy = np.linalg.norm(Y - Ydiag, ord=internal_p, axis=1)**order
else:
- C = sc.cdist(X,Y, metric='minkowski', p=internal_p)**q
- Cxd = np.linalg.norm(X - Xdiag, ord=internal_p, axis=1)**q
- Cdy = np.linalg.norm(Y - Ydiag, ord=internal_p, axis=1)**q
+ C = sc.cdist(X,Y, metric='minkowski', p=internal_p)**order
+ Cxd = np.linalg.norm(X - Xdiag, ord=internal_p, axis=1)**order
+ Cdy = np.linalg.norm(Y - Ydiag, ord=internal_p, axis=1)**order
Cf = np.hstack((C, Cxd[:,None]))
Cdy = np.append(Cdy, 0)
@@ -51,23 +51,23 @@ def _build_dist_matrix(X, Y, q=2., internal_p=2.):
return Cf
-def _perstot(X, q, internal_p):
+def _perstot(X, order, internal_p):
'''
:param X: (n x 2) numpy.array (points of a given diagram).
:param internal_p: Ground metric on the (upper-half) plane (i.e. norm l_p in R^2); Default value is 2 (Euclidean norm).
- :param q: exponent for Wasserstein; Default value is 2.
+ :param order: exponent for Wasserstein; Default value is 2.
:returns: float, the total persistence of the diagram (that is, its distance to the empty diagram).
'''
Xdiag = _proj_on_diag(X)
- return (np.sum(np.linalg.norm(X - Xdiag, ord=internal_p, axis=1)**q))**(1./q)
+ return (np.sum(np.linalg.norm(X - Xdiag, ord=internal_p, axis=1)**order))**(1./order)
-def wasserstein_distance(X, Y, q=2., internal_p=2.):
+def wasserstein_distance(X, Y, order=2., internal_p=2.):
'''
:param X: (n x 2) numpy.array encoding the (finite points of the) first diagram. Must not contain essential points (i.e. with infinite coordinate).
:param Y: (m x 2) numpy.array encoding the second diagram.
:param internal_p: Ground metric on the (upper-half) plane (i.e. norm l_p in R^2); Default value is 2 (euclidean norm).
- :param q: exponent for Wasserstein; Default value is 2.
+ :param order: exponent for Wasserstein; Default value is 2.
:returns: the q-Wasserstein distance (1 <= q < infinity) with respect to the internal_p-norm as ground metric.
:rtype: float
'''
@@ -79,11 +79,11 @@ def wasserstein_distance(X, Y, q=2., internal_p=2.):
if Y.size == 0:
return 0.
else:
- return _perstot(Y, q, internal_p)
+ return _perstot(Y, order, internal_p)
elif Y.size == 0:
- return _perstot(X, q, internal_p)
+ return _perstot(X, order, internal_p)
- M = _build_dist_matrix(X, Y, q=q, internal_p=internal_p)
+ M = _build_dist_matrix(X, Y, order=order, internal_p=internal_p)
a = np.full(n+1, 1. / (n + m) ) # weight vector of the input diagram. Uniform here.
a[-1] = a[-1] * m # normalized so that we have a probability measure, required by POT
b = np.full(m+1, 1. / (n + m) ) # weight vector of the input diagram. Uniform here.
@@ -94,5 +94,5 @@ def wasserstein_distance(X, Y, q=2., internal_p=2.):
# The default numItermax=100000 is not sufficient for some examples with 5000 points, what is a good value?
ot_cost = (n+m) * ot.emd2(a, b, M, numItermax=2000000)
- return ot_cost ** (1./q)
+ return ot_cost ** (1./order)
diff --git a/src/python/test/test_wasserstein_distance.py b/src/python/test/test_wasserstein_distance.py
index 40112c1b..602e4bf1 100755
--- a/src/python/test/test_wasserstein_distance.py
+++ b/src/python/test/test_wasserstein_distance.py
@@ -23,26 +23,26 @@ def test_basic_wasserstein():
diag4 = np.array([[0, 3], [4, 8]])
emptydiag = np.array([[]])
- assert wasserstein_distance(emptydiag, emptydiag, internal_p=2., q=1.) == 0.
- assert wasserstein_distance(emptydiag, emptydiag, internal_p=np.inf, q=1.) == 0.
- assert wasserstein_distance(emptydiag, emptydiag, internal_p=np.inf, q=2.) == 0.
- assert wasserstein_distance(emptydiag, emptydiag, internal_p=2., q=2.) == 0.
+ assert wasserstein_distance(emptydiag, emptydiag, internal_p=2., order=1.) == 0.
+ assert wasserstein_distance(emptydiag, emptydiag, internal_p=np.inf, order=1.) == 0.
+ assert wasserstein_distance(emptydiag, emptydiag, internal_p=np.inf, order=2.) == 0.
+ assert wasserstein_distance(emptydiag, emptydiag, internal_p=2., order=2.) == 0.
- assert wasserstein_distance(diag3, emptydiag, internal_p=np.inf, q=1.) == 2.
- assert wasserstein_distance(diag3, emptydiag, internal_p=1., q=1.) == 4.
+ assert wasserstein_distance(diag3, emptydiag, internal_p=np.inf, order=1.) == 2.
+ assert wasserstein_distance(diag3, emptydiag, internal_p=1., order=1.) == 4.
- assert wasserstein_distance(diag4, emptydiag, internal_p=1., q=2.) == 5. # thank you Pythagorician triplets
- assert wasserstein_distance(diag4, emptydiag, internal_p=np.inf, q=2.) == 2.5
- assert wasserstein_distance(diag4, emptydiag, internal_p=2., q=2.) == 3.5355339059327378
+ assert wasserstein_distance(diag4, emptydiag, internal_p=1., order=2.) == 5. # thank you Pythagorician triplets
+ assert wasserstein_distance(diag4, emptydiag, internal_p=np.inf, order=2.) == 2.5
+ assert wasserstein_distance(diag4, emptydiag, internal_p=2., order=2.) == 3.5355339059327378
- assert wasserstein_distance(diag1, diag2, internal_p=2., q=1.) == 1.4453593023967701
- assert wasserstein_distance(diag1, diag2, internal_p=2.35, q=1.74) == 0.9772734057168739
+ assert wasserstein_distance(diag1, diag2, internal_p=2., order=1.) == 1.4453593023967701
+ assert wasserstein_distance(diag1, diag2, internal_p=2.35, order=1.74) == 0.9772734057168739
- assert wasserstein_distance(diag1, emptydiag, internal_p=2.35, q=1.7863) == 3.141592214572228
+ assert wasserstein_distance(diag1, emptydiag, internal_p=2.35, order=1.7863) == 3.141592214572228
- assert wasserstein_distance(diag3, diag4, internal_p=1., q=1.) == 3.
- assert wasserstein_distance(diag3, diag4, internal_p=np.inf, q=1.) == 3. # no diag matching here
- assert wasserstein_distance(diag3, diag4, internal_p=np.inf, q=2.) == np.sqrt(5)
- assert wasserstein_distance(diag3, diag4, internal_p=1., q=2.) == np.sqrt(5)
- assert wasserstein_distance(diag3, diag4, internal_p=4.5, q=2.) == np.sqrt(5)
+ assert wasserstein_distance(diag3, diag4, internal_p=1., order=1.) == 3.
+ assert wasserstein_distance(diag3, diag4, internal_p=np.inf, order=1.) == 3. # no diag matching here
+ assert wasserstein_distance(diag3, diag4, internal_p=np.inf, order=2.) == np.sqrt(5)
+ assert wasserstein_distance(diag3, diag4, internal_p=1., order=2.) == np.sqrt(5)
+ assert wasserstein_distance(diag3, diag4, internal_p=4.5, order=2.) == np.sqrt(5)