From 74ca2d77bca479785c581d2f48d87628d5c1444d Mon Sep 17 00:00:00 2001 From: Slasnista Date: Fri, 25 Aug 2017 15:12:20 +0200 Subject: refactoring examples according to new DA classes --- examples/da/plot_otda_classes.py | 142 +++++++++++++++++++++ examples/da/plot_otda_color_images.py | 151 ++++++++++++++++++++++ examples/da/plot_otda_d2.py | 163 ++++++++++++++++++++++++ examples/da/plot_otda_mapping.py | 119 +++++++++++++++++ examples/da/plot_otda_mapping_colors_images.py | 169 +++++++++++++++++++++++++ 5 files changed, 744 insertions(+) create mode 100644 examples/da/plot_otda_classes.py create mode 100644 examples/da/plot_otda_color_images.py create mode 100644 examples/da/plot_otda_d2.py create mode 100644 examples/da/plot_otda_mapping.py create mode 100644 examples/da/plot_otda_mapping_colors_images.py (limited to 'examples') diff --git a/examples/da/plot_otda_classes.py b/examples/da/plot_otda_classes.py new file mode 100644 index 0000000..1bfe2bb --- /dev/null +++ b/examples/da/plot_otda_classes.py @@ -0,0 +1,142 @@ +# -*- coding: utf-8 -*- +""" +======================== +OT for domain adaptation +======================== + +This example introduces a domain adaptation in a 2D setting and the 4 OTDA +approaches currently supported in POT. + +""" + +# Authors: Remi Flamary +# Stanilslas Chambon +# +# License: MIT License + +import matplotlib.pylab as pl +import ot + + +# number of source and target points to generate +ns = 150 +nt = 150 + +Xs, ys = ot.datasets.get_data_classif('3gauss', ns) +Xt, yt = ot.datasets.get_data_classif('3gauss2', nt) + +# Instantiate the different transport algorithms and fit them + +# EMD Transport +ot_emd = ot.da.EMDTransport() +ot_emd.fit(Xs=Xs, Xt=Xt) + +# Sinkhorn Transport +ot_sinkhorn = ot.da.SinkhornTransport(reg_e=1e-1) +ot_sinkhorn.fit(Xs=Xs, Xt=Xt) + +# Sinkhorn Transport with Group lasso regularization +ot_lpl1 = ot.da.SinkhornLpl1Transport(reg_e=1e-1, reg_cl=1e0) +ot_lpl1.fit(Xs=Xs, ys=ys, Xt=Xt) + +# Sinkhorn Transport with Group lasso regularization l1l2 +ot_l1l2 = ot.da.SinkhornL1l2Transport(reg_e=1e-1, reg_cl=2e0, max_iter=20, + verbose=True) +ot_l1l2.fit(Xs=Xs, ys=ys, Xt=Xt) + +# transport source samples onto target samples +transp_Xs_emd = ot_emd.transform(Xs=Xs) +transp_Xs_sinkhorn = ot_sinkhorn.transform(Xs=Xs) +transp_Xs_lpl1 = ot_lpl1.transform(Xs=Xs) +transp_Xs_l1l2 = ot_l1l2.transform(Xs=Xs) + +############################################################################## +# Fig 1 : plots source and target samples +############################################################################## + +pl.figure(1, figsize=(10, 5)) +pl.subplot(1, 2, 1) +pl.scatter(Xs[:, 0], Xs[:, 1], c=ys, marker='+', label='Source samples') +pl.xticks([]) +pl.yticks([]) +pl.legend(loc=0) +pl.title('Source samples') + +pl.subplot(1, 2, 2) +pl.scatter(Xt[:, 0], Xt[:, 1], c=yt, marker='o', label='Target samples') +pl.xticks([]) +pl.yticks([]) +pl.legend(loc=0) +pl.title('Target samples') +pl.tight_layout() + +############################################################################## +# Fig 2 : plot optimal couplings and transported samples +############################################################################## + +param_img = {'interpolation': 'nearest', 'cmap': 'spectral'} + +pl.figure(2, figsize=(15, 8)) +pl.subplot(2, 4, 1) +pl.imshow(ot_emd.coupling_, **param_img) +pl.xticks([]) +pl.yticks([]) +pl.title('Optimal coupling\nEMDTransport') + +pl.subplot(2, 4, 2) +pl.imshow(ot_sinkhorn.coupling_, **param_img) +pl.xticks([]) +pl.yticks([]) +pl.title('Optimal coupling\nSinkhornTransport') + +pl.subplot(2, 4, 3) +pl.imshow(ot_lpl1.coupling_, **param_img) +pl.xticks([]) +pl.yticks([]) +pl.title('Optimal coupling\nSinkhornLpl1Transport') + +pl.subplot(2, 4, 4) +pl.imshow(ot_l1l2.coupling_, **param_img) +pl.xticks([]) +pl.yticks([]) +pl.title('Optimal coupling\nSinkhornL1l2Transport') + +pl.subplot(2, 4, 5) +pl.scatter(Xt[:, 0], Xt[:, 1], c=yt, marker='o', + label='Target samples', alpha=0.3) +pl.scatter(transp_Xs_emd[:, 0], transp_Xs_emd[:, 1], c=ys, + marker='+', label='Transp samples', s=30) +pl.xticks([]) +pl.yticks([]) +pl.title('Transported samples\nEmdTransport') +pl.legend(loc="lower left") + +pl.subplot(2, 4, 6) +pl.scatter(Xt[:, 0], Xt[:, 1], c=yt, marker='o', + label='Target samples', alpha=0.3) +pl.scatter(transp_Xs_sinkhorn[:, 0], transp_Xs_sinkhorn[:, 1], c=ys, + marker='+', label='Transp samples', s=30) +pl.xticks([]) +pl.yticks([]) +pl.title('Transported samples\nSinkhornTransport') + +pl.subplot(2, 4, 7) +pl.scatter(Xt[:, 0], Xt[:, 1], c=yt, marker='o', + label='Target samples', alpha=0.3) +pl.scatter(transp_Xs_lpl1[:, 0], transp_Xs_lpl1[:, 1], c=ys, + marker='+', label='Transp samples', s=30) +pl.xticks([]) +pl.yticks([]) +pl.title('Transported samples\nSinkhornLpl1Transport') + +pl.subplot(2, 4, 8) +pl.scatter(Xt[:, 0], Xt[:, 1], c=yt, marker='o', + label='Target samples', alpha=0.3) +pl.scatter(transp_Xs_l1l2[:, 0], transp_Xs_l1l2[:, 1], c=ys, + marker='+', label='Transp samples', s=30) +pl.xticks([]) +pl.yticks([]) +pl.title('Transported samples\nSinkhornL1l2Transport') +pl.tight_layout() + +pl.show() diff --git a/examples/da/plot_otda_color_images.py b/examples/da/plot_otda_color_images.py new file mode 100644 index 0000000..a46ac29 --- /dev/null +++ b/examples/da/plot_otda_color_images.py @@ -0,0 +1,151 @@ +# -*- coding: utf-8 -*- +""" +======================================================== +OT for domain adaptation with image color adaptation [6] +======================================================== + +This example presents a way of transferring colors between two image +with Optimal Transport as introduced in [6] + +[6] Ferradans, S., Papadakis, N., Peyre, G., & Aujol, J. F. (2014). +Regularized discrete optimal transport. +SIAM Journal on Imaging Sciences, 7(3), 1853-1882. +""" + +# Authors: Remi Flamary +# Stanilslas Chambon +# +# License: MIT License + +import numpy as np +from scipy import ndimage +import matplotlib.pylab as pl + +import ot + + +def im2mat(I): + """Converts and image to matrix (one pixel per line)""" + return I.reshape((I.shape[0] * I.shape[1], I.shape[2])) + + +def mat2im(X, shape): + """Converts back a matrix to an image""" + return X.reshape(shape) + + +def minmax(I): + return np.clip(I, 0, 1) + + +# Loading images +I1 = ndimage.imread('../../data/ocean_day.jpg').astype(np.float64) / 256 +I2 = ndimage.imread('../../data/ocean_sunset.jpg').astype(np.float64) / 256 + +X1 = im2mat(I1) +X2 = im2mat(I2) + +# training samples +nb = 1000 +idx1 = np.random.randint(X1.shape[0], size=(nb,)) +idx2 = np.random.randint(X2.shape[0], size=(nb,)) + +Xs = X1[idx1, :] +Xt = X2[idx2, :] + +# EMDTransport +ot_emd = ot.da.EMDTransport() +ot_emd.fit(Xs=Xs, Xt=Xt) + +# SinkhornTransport +ot_sinkhorn = ot.da.SinkhornTransport(reg_e=1e-1) +ot_sinkhorn.fit(Xs=Xs, Xt=Xt) + +# prediction between images (using out of sample prediction as in [6]) +transp_Xs_emd = ot_emd.transform(Xs=X1) +transp_Xt_emd = ot_emd.inverse_transform(Xt=X2) + +transp_Xs_sinkhorn = ot_emd.transform(Xs=X1) +transp_Xt_sinkhorn = ot_emd.inverse_transform(Xt=X2) + +I1t = minmax(mat2im(transp_Xs_emd, I1.shape)) +I2t = minmax(mat2im(transp_Xt_emd, I2.shape)) + +I1te = minmax(mat2im(transp_Xs_sinkhorn, I1.shape)) +I2te = minmax(mat2im(transp_Xt_sinkhorn, I2.shape)) + +############################################################################## +# plot original image +############################################################################## + +pl.figure(1, figsize=(6.4, 3)) + +pl.subplot(1, 2, 1) +pl.imshow(I1) +pl.axis('off') +pl.title('Image 1') + +pl.subplot(1, 2, 2) +pl.imshow(I2) +pl.axis('off') +pl.title('Image 2') + +############################################################################## +# scatter plot of colors +############################################################################## + +pl.figure(2, figsize=(6.4, 3)) + +pl.subplot(1, 2, 1) +pl.scatter(Xs[:, 0], Xs[:, 2], c=Xs) +pl.axis([0, 1, 0, 1]) +pl.xlabel('Red') +pl.ylabel('Blue') +pl.title('Image 1') + +pl.subplot(1, 2, 2) +pl.scatter(Xt[:, 0], Xt[:, 2], c=Xt) +pl.axis([0, 1, 0, 1]) +pl.xlabel('Red') +pl.ylabel('Blue') +pl.title('Image 2') +pl.tight_layout() + +############################################################################## +# plot new images +############################################################################## + +pl.figure(3, figsize=(8, 4)) + +pl.subplot(2, 3, 1) +pl.imshow(I1) +pl.axis('off') +pl.title('Image 1') + +pl.subplot(2, 3, 2) +pl.imshow(I1t) +pl.axis('off') +pl.title('Image 1 Adapt') + +pl.subplot(2, 3, 3) +pl.imshow(I1te) +pl.axis('off') +pl.title('Image 1 Adapt (reg)') + +pl.subplot(2, 3, 4) +pl.imshow(I2) +pl.axis('off') +pl.title('Image 2') + +pl.subplot(2, 3, 5) +pl.imshow(I2t) +pl.axis('off') +pl.title('Image 2 Adapt') + +pl.subplot(2, 3, 6) +pl.imshow(I2te) +pl.axis('off') +pl.title('Image 2 Adapt (reg)') +pl.tight_layout() + +pl.show() diff --git a/examples/da/plot_otda_d2.py b/examples/da/plot_otda_d2.py new file mode 100644 index 0000000..78c0372 --- /dev/null +++ b/examples/da/plot_otda_d2.py @@ -0,0 +1,163 @@ +# -*- coding: utf-8 -*- +""" +============================== +OT for empirical distributions +============================== + +This example introduces a domain adaptation in a 2D setting. It explicits +the problem of domain adaptation and introduces some optimal transport +approaches to solve it. + +Quantities such as optimal couplings, greater coupling coefficients and +transported samples are represented in order to give a visual understanding +of what the transport methods are doing. +""" + +# Authors: Remi Flamary +# Stanilslas Chambon +# +# License: MIT License + +import matplotlib.pylab as pl +import ot + +# number of source and target points to generate +ns = 150 +nt = 150 + +Xs, ys = ot.datasets.get_data_classif('3gauss', ns) +Xt, yt = ot.datasets.get_data_classif('3gauss2', nt) + +# Cost matrix +M = ot.dist(Xs, Xt) + +# Instantiate the different transport algorithms and fit them + +# EMD Transport +ot_emd = ot.da.EMDTransport() +ot_emd.fit(Xs=Xs, Xt=Xt) + +# Sinkhorn Transport +ot_sinkhorn = ot.da.SinkhornTransport(reg_e=1e-1) +ot_sinkhorn.fit(Xs=Xs, Xt=Xt) + +# Sinkhorn Transport with Group lasso regularization +ot_lpl1 = ot.da.SinkhornLpl1Transport(reg_e=1e-1, reg_cl=1e0) +ot_lpl1.fit(Xs=Xs, ys=ys, Xt=Xt) + +# transport source samples onto target samples +transp_Xs_emd = ot_emd.transform(Xs=Xs) +transp_Xs_sinkhorn = ot_sinkhorn.transform(Xs=Xs) +transp_Xs_lpl1 = ot_lpl1.transform(Xs=Xs) + +############################################################################## +# Fig 1 : plots source and target samples + matrix of pairwise distance +############################################################################## + +pl.figure(1, figsize=(10, 10)) +pl.subplot(2, 2, 1) +pl.scatter(Xs[:, 0], Xs[:, 1], c=ys, marker='+', label='Source samples') +pl.xticks([]) +pl.yticks([]) +pl.legend(loc=0) +pl.title('Source samples') + +pl.subplot(2, 2, 2) +pl.scatter(Xt[:, 0], Xt[:, 1], c=yt, marker='o', label='Target samples') +pl.xticks([]) +pl.yticks([]) +pl.legend(loc=0) +pl.title('Target samples') + +pl.subplot(2, 2, 3) +pl.imshow(M, interpolation='nearest') +pl.xticks([]) +pl.yticks([]) +pl.title('Matrix of pairwise distances') +pl.tight_layout() + +############################################################################## +# Fig 2 : plots optimal couplings for the different methods +############################################################################## + +pl.figure(2, figsize=(10, 6)) + +pl.subplot(2, 3, 1) +pl.imshow(ot_emd.coupling_, interpolation='nearest') +pl.xticks([]) +pl.yticks([]) +pl.title('Optimal coupling\nEMDTransport') + +pl.subplot(2, 3, 2) +pl.imshow(ot_sinkhorn.coupling_, interpolation='nearest') +pl.xticks([]) +pl.yticks([]) +pl.title('Optimal coupling\nSinkhornTransport') + +pl.subplot(2, 3, 3) +pl.imshow(ot_lpl1.coupling_, interpolation='nearest') +pl.xticks([]) +pl.yticks([]) +pl.title('Optimal coupling\nSinkhornLpl1Transport') + +pl.subplot(2, 3, 4) +ot.plot.plot2D_samples_mat(Xs, Xt, ot_emd.coupling_, c=[.5, .5, 1]) +pl.scatter(Xs[:, 0], Xs[:, 1], c=ys, marker='+', label='Source samples') +pl.scatter(Xt[:, 0], Xt[:, 1], c=yt, marker='o', label='Target samples') +pl.xticks([]) +pl.yticks([]) +pl.title('Main coupling coefficients\nEMDTransport') + +pl.subplot(2, 3, 5) +ot.plot.plot2D_samples_mat(Xs, Xt, ot_sinkhorn.coupling_, c=[.5, .5, 1]) +pl.scatter(Xs[:, 0], Xs[:, 1], c=ys, marker='+', label='Source samples') +pl.scatter(Xt[:, 0], Xt[:, 1], c=yt, marker='o', label='Target samples') +pl.xticks([]) +pl.yticks([]) +pl.title('Main coupling coefficients\nSinkhornTransport') + +pl.subplot(2, 3, 6) +ot.plot.plot2D_samples_mat(Xs, Xt, ot_lpl1.coupling_, c=[.5, .5, 1]) +pl.scatter(Xs[:, 0], Xs[:, 1], c=ys, marker='+', label='Source samples') +pl.scatter(Xt[:, 0], Xt[:, 1], c=yt, marker='o', label='Target samples') +pl.xticks([]) +pl.yticks([]) +pl.title('Main coupling coefficients\nSinkhornLpl1Transport') +pl.tight_layout() + +############################################################################## +# Fig 3 : plot transported samples +############################################################################## + +# display transported samples +pl.figure(4, figsize=(10, 4)) +pl.subplot(1, 3, 1) +pl.scatter(Xt[:, 0], Xt[:, 1], c=yt, marker='o', + label='Target samples', alpha=0.5) +pl.scatter(transp_Xs_emd[:, 0], transp_Xs_emd[:, 1], c=ys, + marker='+', label='Transp samples', s=30) +pl.title('Transported samples\nEmdTransport') +pl.legend(loc=0) +pl.xticks([]) +pl.yticks([]) + +pl.subplot(1, 3, 2) +pl.scatter(Xt[:, 0], Xt[:, 1], c=yt, marker='o', + label='Target samples', alpha=0.5) +pl.scatter(transp_Xs_sinkhorn[:, 0], transp_Xs_sinkhorn[:, 1], c=ys, + marker='+', label='Transp samples', s=30) +pl.title('Transported samples\nSinkhornTransport') +pl.xticks([]) +pl.yticks([]) + +pl.subplot(1, 3, 3) +pl.scatter(Xt[:, 0], Xt[:, 1], c=yt, marker='o', + label='Target samples', alpha=0.5) +pl.scatter(transp_Xs_lpl1[:, 0], transp_Xs_lpl1[:, 1], c=ys, + marker='+', label='Transp samples', s=30) +pl.title('Transported samples\nSinkhornLpl1Transport') +pl.xticks([]) +pl.yticks([]) + +pl.tight_layout() +pl.show() diff --git a/examples/da/plot_otda_mapping.py b/examples/da/plot_otda_mapping.py new file mode 100644 index 0000000..ed234f5 --- /dev/null +++ b/examples/da/plot_otda_mapping.py @@ -0,0 +1,119 @@ +# -*- coding: utf-8 -*- +""" +=============================================== +OT mapping estimation for domain adaptation [8] +=============================================== + +This example presents how to use MappingTransport to estimate at the same +time both the coupling transport and approximate the transport map with either +a linear or a kernelized mapping as introduced in [8] + +[8] M. Perrot, N. Courty, R. Flamary, A. Habrard, + "Mapping estimation for discrete optimal transport", + Neural Information Processing Systems (NIPS), 2016. +""" + +# Authors: Remi Flamary +# Stanilslas Chambon +# +# License: MIT License + +import numpy as np +import matplotlib.pylab as pl +import ot + + +np.random.seed(0) + +############################################################################## +# generate +############################################################################## + +n = 100 # nb samples in source and target datasets +theta = 2 * np.pi / 20 +nz = 0.1 +Xs, ys = ot.datasets.get_data_classif('gaussrot', n, nz=nz) +Xs_new, _ = ot.datasets.get_data_classif('gaussrot', n, nz=nz) +Xt, yt = ot.datasets.get_data_classif('gaussrot', n, theta=theta, nz=nz) + +# one of the target mode changes its variance (no linear mapping) +Xt[yt == 2] *= 3 +Xt = Xt + 4 + + +# MappingTransport with linear kernel +ot_mapping_linear = ot.da.MappingTransport( + kernel="linear", mu=1e0, eta=1e-8, bias=True, + max_iter=20, verbose=True) + +ot_mapping_linear.fit( + Xs=Xs, Xt=Xt) + +# for original source samples, transform applies barycentric mapping +transp_Xs_linear = ot_mapping_linear.transform(Xs=Xs) + +# for out of source samples, transform applies the linear mapping +transp_Xs_linear_new = ot_mapping_linear.transform(Xs=Xs_new) + + +# MappingTransport with gaussian kernel +ot_mapping_gaussian = ot.da.MappingTransport( + kernel="gaussian", eta=1e-5, mu=1e-1, bias=True, sigma=1, + max_iter=10, verbose=True) +ot_mapping_gaussian.fit(Xs=Xs, Xt=Xt) + +# for original source samples, transform applies barycentric mapping +transp_Xs_gaussian = ot_mapping_gaussian.transform(Xs=Xs) + +# for out of source samples, transform applies the gaussian mapping +transp_Xs_gaussian_new = ot_mapping_gaussian.transform(Xs=Xs_new) + + +############################################################################## +# plot data +############################################################################## + +pl.figure(1, (10, 5)) +pl.clf() +pl.scatter(Xs[:, 0], Xs[:, 1], c=ys, marker='+', label='Source samples') +pl.scatter(Xt[:, 0], Xt[:, 1], c=yt, marker='o', label='Target samples') +pl.legend(loc=0) +pl.title('Source and target distributions') + +############################################################################## +# plot transported samples +############################################################################## + +pl.figure(2) +pl.clf() +pl.subplot(2, 2, 1) +pl.scatter(Xt[:, 0], Xt[:, 1], c=yt, marker='o', + label='Target samples', alpha=.2) +pl.scatter(transp_Xs_linear[:, 0], transp_Xs_linear[:, 1], c=ys, marker='+', + label='Mapped source samples') +pl.title("Bary. mapping (linear)") +pl.legend(loc=0) + +pl.subplot(2, 2, 2) +pl.scatter(Xt[:, 0], Xt[:, 1], c=yt, marker='o', + label='Target samples', alpha=.2) +pl.scatter(transp_Xs_linear_new[:, 0], transp_Xs_linear_new[:, 1], + c=ys, marker='+', label='Learned mapping') +pl.title("Estim. mapping (linear)") + +pl.subplot(2, 2, 3) +pl.scatter(Xt[:, 0], Xt[:, 1], c=yt, marker='o', + label='Target samples', alpha=.2) +pl.scatter(transp_Xs_gaussian[:, 0], transp_Xs_gaussian[:, 1], c=ys, + marker='+', label='barycentric mapping') +pl.title("Bary. mapping (kernel)") + +pl.subplot(2, 2, 4) +pl.scatter(Xt[:, 0], Xt[:, 1], c=yt, marker='o', + label='Target samples', alpha=.2) +pl.scatter(transp_Xs_gaussian_new[:, 0], transp_Xs_gaussian_new[:, 1], c=ys, + marker='+', label='Learned mapping') +pl.title("Estim. mapping (kernel)") +pl.tight_layout() + +pl.show() diff --git a/examples/da/plot_otda_mapping_colors_images.py b/examples/da/plot_otda_mapping_colors_images.py new file mode 100644 index 0000000..56b5a6f --- /dev/null +++ b/examples/da/plot_otda_mapping_colors_images.py @@ -0,0 +1,169 @@ +# -*- coding: utf-8 -*- +""" +==================================================================================== +OT for domain adaptation with image color adaptation [6] with mapping estimation [8] +==================================================================================== + +[6] Ferradans, S., Papadakis, N., Peyre, G., & Aujol, J. F. (2014). Regularized + discrete optimal transport. SIAM Journal on Imaging Sciences, 7(3), + 1853-1882. +[8] M. Perrot, N. Courty, R. Flamary, A. Habrard, "Mapping estimation for + discrete optimal transport", Neural Information Processing Systems (NIPS), + 2016. + +""" + +# Authors: Remi Flamary +# Stanilslas Chambon +# +# License: MIT License + +import numpy as np +from scipy import ndimage +import matplotlib.pylab as pl +import ot + + +def im2mat(I): + """Converts and image to matrix (one pixel per line)""" + return I.reshape((I.shape[0] * I.shape[1], I.shape[2])) + + +def mat2im(X, shape): + """Converts back a matrix to an image""" + return X.reshape(shape) + + +def minmax(I): + return np.clip(I, 0, 1) + + +############################################################################## +# Generate data +############################################################################## + +# Loading images +# I1 = ndimage.imread('../../data/ocean_day.jpg').astype(np.float64) / 256 +# I2 = ndimage.imread('../../data/ocean_sunset.jpg').astype(np.float64) / 256 + +I1 = ndimage.imread('data/ocean_day.jpg').astype(np.float64) / 256 +I2 = ndimage.imread('data/ocean_sunset.jpg').astype(np.float64) / 256 + + +X1 = im2mat(I1) +X2 = im2mat(I2) + +# training samples +nb = 1000 +idx1 = np.random.randint(X1.shape[0], size=(nb,)) +idx2 = np.random.randint(X2.shape[0], size=(nb,)) + +Xs = X1[idx1, :] +Xt = X2[idx2, :] + + +############################################################################## +# Domain adaptation for pixel distribution transfer +############################################################################## + +# EMDTransport +ot_emd = ot.da.EMDTransport() +ot_emd.fit(Xs=Xs, Xt=Xt) +transp_Xs_emd = ot_emd.transform(Xs=X1) +Image_emd = minmax(mat2im(transp_Xs_emd, I1.shape)) + +# SinkhornTransport +ot_sinkhorn = ot.da.SinkhornTransport(reg_e=1e-1) +ot_sinkhorn.fit(Xs=Xs, Xt=Xt) +transp_Xs_sinkhorn = ot_emd.transform(Xs=X1) +Image_sinkhorn = minmax(mat2im(transp_Xs_sinkhorn, I1.shape)) + +ot_mapping_linear = ot.da.MappingTransport( + mu=1e0, eta=1e-8, bias=True, max_iter=20, verbose=True) +ot_mapping_linear.fit(Xs=Xs, Xt=Xt) + +X1tl = ot_mapping_linear.transform(X1) +Image_mapping_linear = minmax(mat2im(X1tl, I1.shape)) + +ot_mapping_gaussian = ot.da.MappingTransport( + mu=1e0, eta=1e-2, sigma=1, bias=False, max_iter=10, verbose=True) +ot_mapping_gaussian.fit(Xs=Xs, Xt=Xt) + +X1tn = ot_mapping_gaussian.transform(X1) # use the estimated mapping +Image_mapping_gaussian = minmax(mat2im(X1tn, I1.shape)) + +############################################################################## +# plot original images +############################################################################## + +pl.figure(1, figsize=(6.4, 3)) +pl.subplot(1, 2, 1) +pl.imshow(I1) +pl.axis('off') +pl.title('Image 1') + +pl.subplot(1, 2, 2) +pl.imshow(I2) +pl.axis('off') +pl.title('Image 2') +pl.tight_layout() + +############################################################################## +# plot pixel values distribution +############################################################################## + +pl.figure(2, figsize=(6.4, 5)) + +pl.subplot(1, 2, 1) +pl.scatter(Xs[:, 0], Xs[:, 2], c=Xs) +pl.axis([0, 1, 0, 1]) +pl.xlabel('Red') +pl.ylabel('Blue') +pl.title('Image 1') + +pl.subplot(1, 2, 2) +pl.scatter(Xt[:, 0], Xt[:, 2], c=Xt) +pl.axis([0, 1, 0, 1]) +pl.xlabel('Red') +pl.ylabel('Blue') +pl.title('Image 2') +pl.tight_layout() + +############################################################################## +# plot transformed images +############################################################################## + +pl.figure(2, figsize=(10, 5)) + +pl.subplot(2, 3, 1) +pl.imshow(I1) +pl.axis('off') +pl.title('Im. 1') + +pl.subplot(2, 3, 4) +pl.imshow(I2) +pl.axis('off') +pl.title('Im. 2') + +pl.subplot(2, 3, 2) +pl.imshow(Image_emd) +pl.axis('off') +pl.title('EmdTransport') + +pl.subplot(2, 3, 5) +pl.imshow(Image_sinkhorn) +pl.axis('off') +pl.title('SinkhornTransport') + +pl.subplot(2, 3, 3) +pl.imshow(Image_mapping_linear) +pl.axis('off') +pl.title('MappingTransport (linear)') + +pl.subplot(2, 3, 6) +pl.imshow(Image_mapping_gaussian) +pl.axis('off') +pl.title('MappingTransport (gaussian)') +pl.tight_layout() + +pl.show() -- cgit v1.2.3 From f80693b545ac40817cb95eb31260fe7598514b96 Mon Sep 17 00:00:00 2001 From: Slasnista Date: Fri, 25 Aug 2017 15:15:52 +0200 Subject: small corrections for examples --- examples/da/plot_otda_classes.py | 2 +- examples/da/plot_otda_color_images.py | 2 +- examples/da/plot_otda_d2.py | 2 +- examples/da/plot_otda_mapping.py | 2 +- examples/da/plot_otda_mapping_colors_images.py | 6 +++--- 5 files changed, 7 insertions(+), 7 deletions(-) (limited to 'examples') diff --git a/examples/da/plot_otda_classes.py b/examples/da/plot_otda_classes.py index 1bfe2bb..e5c82fb 100644 --- a/examples/da/plot_otda_classes.py +++ b/examples/da/plot_otda_classes.py @@ -10,7 +10,7 @@ approaches currently supported in POT. """ # Authors: Remi Flamary -# Stanilslas Chambon +# Stanislas Chambon # # License: MIT License diff --git a/examples/da/plot_otda_color_images.py b/examples/da/plot_otda_color_images.py index a46ac29..bca7350 100644 --- a/examples/da/plot_otda_color_images.py +++ b/examples/da/plot_otda_color_images.py @@ -13,7 +13,7 @@ SIAM Journal on Imaging Sciences, 7(3), 1853-1882. """ # Authors: Remi Flamary -# Stanilslas Chambon +# Stanislas Chambon # # License: MIT License diff --git a/examples/da/plot_otda_d2.py b/examples/da/plot_otda_d2.py index 78c0372..1d2192f 100644 --- a/examples/da/plot_otda_d2.py +++ b/examples/da/plot_otda_d2.py @@ -14,7 +14,7 @@ of what the transport methods are doing. """ # Authors: Remi Flamary -# Stanilslas Chambon +# Stanislas Chambon # # License: MIT License diff --git a/examples/da/plot_otda_mapping.py b/examples/da/plot_otda_mapping.py index ed234f5..6d83507 100644 --- a/examples/da/plot_otda_mapping.py +++ b/examples/da/plot_otda_mapping.py @@ -14,7 +14,7 @@ a linear or a kernelized mapping as introduced in [8] """ # Authors: Remi Flamary -# Stanilslas Chambon +# Stanislas Chambon # # License: MIT License diff --git a/examples/da/plot_otda_mapping_colors_images.py b/examples/da/plot_otda_mapping_colors_images.py index 56b5a6f..05d9046 100644 --- a/examples/da/plot_otda_mapping_colors_images.py +++ b/examples/da/plot_otda_mapping_colors_images.py @@ -14,7 +14,7 @@ OT for domain adaptation with image color adaptation [6] with mapping estimation """ # Authors: Remi Flamary -# Stanilslas Chambon +# Stanislas Chambon # # License: MIT License @@ -82,14 +82,14 @@ ot_mapping_linear = ot.da.MappingTransport( mu=1e0, eta=1e-8, bias=True, max_iter=20, verbose=True) ot_mapping_linear.fit(Xs=Xs, Xt=Xt) -X1tl = ot_mapping_linear.transform(X1) +X1tl = ot_mapping_linear.transform(Xs=X1) Image_mapping_linear = minmax(mat2im(X1tl, I1.shape)) ot_mapping_gaussian = ot.da.MappingTransport( mu=1e0, eta=1e-2, sigma=1, bias=False, max_iter=10, verbose=True) ot_mapping_gaussian.fit(Xs=Xs, Xt=Xt) -X1tn = ot_mapping_gaussian.transform(X1) # use the estimated mapping +X1tn = ot_mapping_gaussian.transform(Xs=X1) # use the estimated mapping Image_mapping_gaussian = minmax(mat2im(X1tn, I1.shape)) ############################################################################## -- cgit v1.2.3 From 892d7ce10912de3d07692b5083578932a32ab61f Mon Sep 17 00:00:00 2001 From: Slasnista Date: Fri, 25 Aug 2017 15:18:05 +0200 Subject: set properly path of data --- examples/da/plot_otda_mapping_colors_images.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) (limited to 'examples') diff --git a/examples/da/plot_otda_mapping_colors_images.py b/examples/da/plot_otda_mapping_colors_images.py index 05d9046..4209020 100644 --- a/examples/da/plot_otda_mapping_colors_images.py +++ b/examples/da/plot_otda_mapping_colors_images.py @@ -43,11 +43,8 @@ def minmax(I): ############################################################################## # Loading images -# I1 = ndimage.imread('../../data/ocean_day.jpg').astype(np.float64) / 256 -# I2 = ndimage.imread('../../data/ocean_sunset.jpg').astype(np.float64) / 256 - -I1 = ndimage.imread('data/ocean_day.jpg').astype(np.float64) / 256 -I2 = ndimage.imread('data/ocean_sunset.jpg').astype(np.float64) / 256 +I1 = ndimage.imread('../../data/ocean_day.jpg').astype(np.float64) / 256 +I2 = ndimage.imread('../../data/ocean_sunset.jpg').astype(np.float64) / 256 X1 = im2mat(I1) -- cgit v1.2.3 From a8fa91bec26caa93329e61a104e0ad6afdf37363 Mon Sep 17 00:00:00 2001 From: Slasnista Date: Mon, 28 Aug 2017 11:03:28 +0200 Subject: handling input arguments in fit, transform... methods + remove old examples --- examples/plot_OTDA_2D.py | 126 ----------------- examples/plot_OTDA_classes.py | 117 ---------------- examples/plot_OTDA_color_images.py | 152 -------------------- examples/plot_OTDA_mapping.py | 124 ----------------- examples/plot_OTDA_mapping_color_images.py | 169 ---------------------- ot/da.py | 217 ++++++++++++++++------------- 6 files changed, 121 insertions(+), 784 deletions(-) delete mode 100644 examples/plot_OTDA_2D.py delete mode 100644 examples/plot_OTDA_classes.py delete mode 100644 examples/plot_OTDA_color_images.py delete mode 100644 examples/plot_OTDA_mapping.py delete mode 100644 examples/plot_OTDA_mapping_color_images.py (limited to 'examples') diff --git a/examples/plot_OTDA_2D.py b/examples/plot_OTDA_2D.py deleted file mode 100644 index f2108c6..0000000 --- a/examples/plot_OTDA_2D.py +++ /dev/null @@ -1,126 +0,0 @@ -# -*- coding: utf-8 -*- -""" -============================== -OT for empirical distributions -============================== - -""" - -# Author: Remi Flamary -# -# License: MIT License - -import numpy as np -import matplotlib.pylab as pl -import ot - - -#%% parameters - -n = 150 # nb bins - -xs, ys = ot.datasets.get_data_classif('3gauss', n) -xt, yt = ot.datasets.get_data_classif('3gauss2', n) - -a, b = ot.unif(n), ot.unif(n) -# loss matrix -M = ot.dist(xs, xt) -# M/=M.max() - -#%% plot samples - -pl.figure(1) -pl.subplot(2, 2, 1) -pl.scatter(xs[:, 0], xs[:, 1], c=ys, marker='+', label='Source samples') -pl.legend(loc=0) -pl.title('Source distributions') - -pl.subplot(2, 2, 2) -pl.scatter(xt[:, 0], xt[:, 1], c=yt, marker='o', label='Target samples') -pl.legend(loc=0) -pl.title('target distributions') - -pl.figure(2) -pl.imshow(M, interpolation='nearest') -pl.title('Cost matrix M') - - -#%% OT estimation - -# EMD -G0 = ot.emd(a, b, M) - -# sinkhorn -lambd = 1e-1 -Gs = ot.sinkhorn(a, b, M, lambd) - - -# Group lasso regularization -reg = 1e-1 -eta = 1e0 -Gg = ot.da.sinkhorn_lpl1_mm(a, ys.astype(np.int), b, M, reg, eta) - - -#%% visu matrices - -pl.figure(3) - -pl.subplot(2, 3, 1) -pl.imshow(G0, interpolation='nearest') -pl.title('OT matrix ') - -pl.subplot(2, 3, 2) -pl.imshow(Gs, interpolation='nearest') -pl.title('OT matrix Sinkhorn') - -pl.subplot(2, 3, 3) -pl.imshow(Gg, interpolation='nearest') -pl.title('OT matrix Group lasso') - -pl.subplot(2, 3, 4) -ot.plot.plot2D_samples_mat(xs, xt, G0, c=[.5, .5, 1]) -pl.scatter(xs[:, 0], xs[:, 1], c=ys, marker='+', label='Source samples') -pl.scatter(xt[:, 0], xt[:, 1], c=yt, marker='o', label='Target samples') - - -pl.subplot(2, 3, 5) -ot.plot.plot2D_samples_mat(xs, xt, Gs, c=[.5, .5, 1]) -pl.scatter(xs[:, 0], xs[:, 1], c=ys, marker='+', label='Source samples') -pl.scatter(xt[:, 0], xt[:, 1], c=yt, marker='o', label='Target samples') - -pl.subplot(2, 3, 6) -ot.plot.plot2D_samples_mat(xs, xt, Gg, c=[.5, .5, 1]) -pl.scatter(xs[:, 0], xs[:, 1], c=ys, marker='+', label='Source samples') -pl.scatter(xt[:, 0], xt[:, 1], c=yt, marker='o', label='Target samples') -pl.tight_layout() - -#%% sample interpolation - -xst0 = n * G0.dot(xt) -xsts = n * Gs.dot(xt) -xstg = n * Gg.dot(xt) - -pl.figure(4, figsize=(8, 3)) -pl.subplot(1, 3, 1) -pl.scatter(xt[:, 0], xt[:, 1], c=yt, marker='o', - label='Target samples', alpha=0.5) -pl.scatter(xst0[:, 0], xst0[:, 1], c=ys, - marker='+', label='Transp samples', s=30) -pl.title('Interp samples') -pl.legend(loc=0) - -pl.subplot(1, 3, 2) -pl.scatter(xt[:, 0], xt[:, 1], c=yt, marker='o', - label='Target samples', alpha=0.5) -pl.scatter(xsts[:, 0], xsts[:, 1], c=ys, - marker='+', label='Transp samples', s=30) -pl.title('Interp samples Sinkhorn') - -pl.subplot(1, 3, 3) -pl.scatter(xt[:, 0], xt[:, 1], c=yt, marker='o', - label='Target samples', alpha=0.5) -pl.scatter(xstg[:, 0], xstg[:, 1], c=ys, - marker='+', label='Transp samples', s=30) -pl.title('Interp samples Grouplasso') -pl.tight_layout() -pl.show() diff --git a/examples/plot_OTDA_classes.py b/examples/plot_OTDA_classes.py deleted file mode 100644 index 53e4bae..0000000 --- a/examples/plot_OTDA_classes.py +++ /dev/null @@ -1,117 +0,0 @@ -# -*- coding: utf-8 -*- -""" -======================== -OT for domain adaptation -======================== - -""" - -# Author: Remi Flamary -# -# License: MIT License - -import matplotlib.pylab as pl -import ot - - -#%% parameters - -n = 150 # nb samples in source and target datasets - -xs, ys = ot.datasets.get_data_classif('3gauss', n) -xt, yt = ot.datasets.get_data_classif('3gauss2', n) - - -#%% plot samples - -pl.figure(1, figsize=(6.4, 3)) - -pl.subplot(1, 2, 1) -pl.scatter(xs[:, 0], xs[:, 1], c=ys, marker='+', label='Source samples') -pl.legend(loc=0) -pl.title('Source distributions') - -pl.subplot(1, 2, 2) -pl.scatter(xt[:, 0], xt[:, 1], c=yt, marker='o', label='Target samples') -pl.legend(loc=0) -pl.title('target distributions') - - -#%% OT estimation - -# LP problem -da_emd = ot.da.OTDA() # init class -da_emd.fit(xs, xt) # fit distributions -xst0 = da_emd.interp() # interpolation of source samples - -# sinkhorn regularization -lambd = 1e-1 -da_entrop = ot.da.OTDA_sinkhorn() -da_entrop.fit(xs, xt, reg=lambd) -xsts = da_entrop.interp() - -# non-convex Group lasso regularization -reg = 1e-1 -eta = 1e0 -da_lpl1 = ot.da.OTDA_lpl1() -da_lpl1.fit(xs, ys, xt, reg=reg, eta=eta) -xstg = da_lpl1.interp() - -# True Group lasso regularization -reg = 1e-1 -eta = 2e0 -da_l1l2 = ot.da.OTDA_l1l2() -da_l1l2.fit(xs, ys, xt, reg=reg, eta=eta, numItermax=20, verbose=True) -xstgl = da_l1l2.interp() - -#%% plot interpolated source samples - -param_img = {'interpolation': 'nearest', 'cmap': 'spectral'} - -pl.figure(2, figsize=(8, 4.5)) -pl.subplot(2, 4, 1) -pl.imshow(da_emd.G, **param_img) -pl.title('OT matrix') - -pl.subplot(2, 4, 2) -pl.imshow(da_entrop.G, **param_img) -pl.title('OT matrix\nsinkhorn') - -pl.subplot(2, 4, 3) -pl.imshow(da_lpl1.G, **param_img) -pl.title('OT matrix\nnon-convex Group Lasso') - -pl.subplot(2, 4, 4) -pl.imshow(da_l1l2.G, **param_img) -pl.title('OT matrix\nGroup Lasso') - -pl.subplot(2, 4, 5) -pl.scatter(xt[:, 0], xt[:, 1], c=yt, marker='o', - label='Target samples', alpha=0.3) -pl.scatter(xst0[:, 0], xst0[:, 1], c=ys, - marker='+', label='Transp samples', s=30) -pl.title('Interp samples') -pl.legend(loc=0) - -pl.subplot(2, 4, 6) -pl.scatter(xt[:, 0], xt[:, 1], c=yt, marker='o', - label='Target samples', alpha=0.3) -pl.scatter(xsts[:, 0], xsts[:, 1], c=ys, - marker='+', label='Transp samples', s=30) -pl.title('Interp samples\nSinkhorn') - -pl.subplot(2, 4, 7) -pl.scatter(xt[:, 0], xt[:, 1], c=yt, marker='o', - label='Target samples', alpha=0.3) -pl.scatter(xstg[:, 0], xstg[:, 1], c=ys, - marker='+', label='Transp samples', s=30) -pl.title('Interp samples\nnon-convex Group Lasso') - -pl.subplot(2, 4, 8) -pl.scatter(xt[:, 0], xt[:, 1], c=yt, marker='o', - label='Target samples', alpha=0.3) -pl.scatter(xstgl[:, 0], xstgl[:, 1], c=ys, - marker='+', label='Transp samples', s=30) -pl.title('Interp samples\nGroup Lasso') -pl.tight_layout() -pl.show() diff --git a/examples/plot_OTDA_color_images.py b/examples/plot_OTDA_color_images.py deleted file mode 100644 index c5ff873..0000000 --- a/examples/plot_OTDA_color_images.py +++ /dev/null @@ -1,152 +0,0 @@ -# -*- coding: utf-8 -*- -""" -======================================================== -OT for domain adaptation with image color adaptation [6] -======================================================== - -[6] Ferradans, S., Papadakis, N., Peyre, G., & Aujol, J. F. (2014). -Regularized discrete optimal transport. -SIAM Journal on Imaging Sciences, 7(3), 1853-1882. -""" - -# Author: Remi Flamary -# -# License: MIT License - -import numpy as np -from scipy import ndimage -import matplotlib.pylab as pl -import ot - - -#%% Loading images - -I1 = ndimage.imread('../data/ocean_day.jpg').astype(np.float64) / 256 -I2 = ndimage.imread('../data/ocean_sunset.jpg').astype(np.float64) / 256 - -#%% Plot images - -pl.figure(1, figsize=(6.4, 3)) - -pl.subplot(1, 2, 1) -pl.imshow(I1) -pl.axis('off') -pl.title('Image 1') - -pl.subplot(1, 2, 2) -pl.imshow(I2) -pl.axis('off') -pl.title('Image 2') - -pl.show() - -#%% Image conversion and dataset generation - - -def im2mat(I): - """Converts and image to matrix (one pixel per line)""" - return I.reshape((I.shape[0] * I.shape[1], I.shape[2])) - - -def mat2im(X, shape): - """Converts back a matrix to an image""" - return X.reshape(shape) - - -X1 = im2mat(I1) -X2 = im2mat(I2) - -# training samples -nb = 1000 -idx1 = np.random.randint(X1.shape[0], size=(nb,)) -idx2 = np.random.randint(X2.shape[0], size=(nb,)) - -xs = X1[idx1, :] -xt = X2[idx2, :] - -#%% Plot image distributions - - -pl.figure(2, figsize=(6.4, 3)) - -pl.subplot(1, 2, 1) -pl.scatter(xs[:, 0], xs[:, 2], c=xs) -pl.axis([0, 1, 0, 1]) -pl.xlabel('Red') -pl.ylabel('Blue') -pl.title('Image 1') - -pl.subplot(1, 2, 2) -pl.scatter(xt[:, 0], xt[:, 2], c=xt) -pl.axis([0, 1, 0, 1]) -pl.xlabel('Red') -pl.ylabel('Blue') -pl.title('Image 2') -pl.tight_layout() - -#%% domain adaptation between images - -# LP problem -da_emd = ot.da.OTDA() # init class -da_emd.fit(xs, xt) # fit distributions - -# sinkhorn regularization -lambd = 1e-1 -da_entrop = ot.da.OTDA_sinkhorn() -da_entrop.fit(xs, xt, reg=lambd) - -#%% prediction between images (using out of sample prediction as in [6]) - -X1t = da_emd.predict(X1) -X2t = da_emd.predict(X2, -1) - -X1te = da_entrop.predict(X1) -X2te = da_entrop.predict(X2, -1) - - -def minmax(I): - return np.clip(I, 0, 1) - - -I1t = minmax(mat2im(X1t, I1.shape)) -I2t = minmax(mat2im(X2t, I2.shape)) - -I1te = minmax(mat2im(X1te, I1.shape)) -I2te = minmax(mat2im(X2te, I2.shape)) - -#%% plot all images - -pl.figure(2, figsize=(8, 4)) - -pl.subplot(2, 3, 1) -pl.imshow(I1) -pl.axis('off') -pl.title('Image 1') - -pl.subplot(2, 3, 2) -pl.imshow(I1t) -pl.axis('off') -pl.title('Image 1 Adapt') - -pl.subplot(2, 3, 3) -pl.imshow(I1te) -pl.axis('off') -pl.title('Image 1 Adapt (reg)') - -pl.subplot(2, 3, 4) -pl.imshow(I2) -pl.axis('off') -pl.title('Image 2') - -pl.subplot(2, 3, 5) -pl.imshow(I2t) -pl.axis('off') -pl.title('Image 2 Adapt') - -pl.subplot(2, 3, 6) -pl.imshow(I2te) -pl.axis('off') -pl.title('Image 2 Adapt (reg)') -pl.tight_layout() - -pl.show() diff --git a/examples/plot_OTDA_mapping.py b/examples/plot_OTDA_mapping.py deleted file mode 100644 index a0d7f8b..0000000 --- a/examples/plot_OTDA_mapping.py +++ /dev/null @@ -1,124 +0,0 @@ -# -*- coding: utf-8 -*- -""" -=============================================== -OT mapping estimation for domain adaptation [8] -=============================================== - -[8] M. Perrot, N. Courty, R. Flamary, A. Habrard, - "Mapping estimation for discrete optimal transport", - Neural Information Processing Systems (NIPS), 2016. -""" - -# Author: Remi Flamary -# -# License: MIT License - -import numpy as np -import matplotlib.pylab as pl -import ot - - -#%% dataset generation - -np.random.seed(0) # makes example reproducible - -n = 100 # nb samples in source and target datasets -theta = 2 * np.pi / 20 -nz = 0.1 -xs, ys = ot.datasets.get_data_classif('gaussrot', n, nz=nz) -xt, yt = ot.datasets.get_data_classif('gaussrot', n, theta=theta, nz=nz) - -# one of the target mode changes its variance (no linear mapping) -xt[yt == 2] *= 3 -xt = xt + 4 - - -#%% plot samples - -pl.figure(1, (6.4, 3)) -pl.clf() -pl.scatter(xs[:, 0], xs[:, 1], c=ys, marker='+', label='Source samples') -pl.scatter(xt[:, 0], xt[:, 1], c=yt, marker='o', label='Target samples') -pl.legend(loc=0) -pl.title('Source and target distributions') - - -#%% OT linear mapping estimation - -eta = 1e-8 # quadratic regularization for regression -mu = 1e0 # weight of the OT linear term -bias = True # estimate a bias - -ot_mapping = ot.da.OTDA_mapping_linear() -ot_mapping.fit(xs, xt, mu=mu, eta=eta, bias=bias, numItermax=20, verbose=True) - -xst = ot_mapping.predict(xs) # use the estimated mapping -xst0 = ot_mapping.interp() # use barycentric mapping - - -pl.figure(2) -pl.clf() -pl.subplot(2, 2, 1) -pl.scatter(xt[:, 0], xt[:, 1], c=yt, marker='o', - label='Target samples', alpha=.3) -pl.scatter(xst0[:, 0], xst0[:, 1], c=ys, - marker='+', label='barycentric mapping') -pl.title("barycentric mapping") - -pl.subplot(2, 2, 2) -pl.scatter(xt[:, 0], xt[:, 1], c=yt, marker='o', - label='Target samples', alpha=.3) -pl.scatter(xst[:, 0], xst[:, 1], c=ys, marker='+', label='Learned mapping') -pl.title("Learned mapping") -pl.tight_layout() - -#%% Kernel mapping estimation - -eta = 1e-5 # quadratic regularization for regression -mu = 1e-1 # weight of the OT linear term -bias = True # estimate a bias -sigma = 1 # sigma bandwidth fot gaussian kernel - - -ot_mapping_kernel = ot.da.OTDA_mapping_kernel() -ot_mapping_kernel.fit( - xs, xt, mu=mu, eta=eta, sigma=sigma, bias=bias, numItermax=10, verbose=True) - -xst_kernel = ot_mapping_kernel.predict(xs) # use the estimated mapping -xst0_kernel = ot_mapping_kernel.interp() # use barycentric mapping - - -#%% Plotting the mapped samples - -pl.figure(2) -pl.clf() -pl.subplot(2, 2, 1) -pl.scatter(xt[:, 0], xt[:, 1], c=yt, marker='o', - label='Target samples', alpha=.2) -pl.scatter(xst0[:, 0], xst0[:, 1], c=ys, marker='+', - label='Mapped source samples') -pl.title("Bary. mapping (linear)") -pl.legend(loc=0) - -pl.subplot(2, 2, 2) -pl.scatter(xt[:, 0], xt[:, 1], c=yt, marker='o', - label='Target samples', alpha=.2) -pl.scatter(xst[:, 0], xst[:, 1], c=ys, marker='+', label='Learned mapping') -pl.title("Estim. mapping (linear)") - -pl.subplot(2, 2, 3) -pl.scatter(xt[:, 0], xt[:, 1], c=yt, marker='o', - label='Target samples', alpha=.2) -pl.scatter(xst0_kernel[:, 0], xst0_kernel[:, 1], c=ys, - marker='+', label='barycentric mapping') -pl.title("Bary. mapping (kernel)") - -pl.subplot(2, 2, 4) -pl.scatter(xt[:, 0], xt[:, 1], c=yt, marker='o', - label='Target samples', alpha=.2) -pl.scatter(xst_kernel[:, 0], xst_kernel[:, 1], c=ys, - marker='+', label='Learned mapping') -pl.title("Estim. mapping (kernel)") -pl.tight_layout() - -pl.show() diff --git a/examples/plot_OTDA_mapping_color_images.py b/examples/plot_OTDA_mapping_color_images.py deleted file mode 100644 index 8064b25..0000000 --- a/examples/plot_OTDA_mapping_color_images.py +++ /dev/null @@ -1,169 +0,0 @@ -# -*- coding: utf-8 -*- -""" -==================================================================================== -OT for domain adaptation with image color adaptation [6] with mapping estimation [8] -==================================================================================== - -[6] Ferradans, S., Papadakis, N., Peyre, G., & Aujol, J. F. (2014). Regularized - discrete optimal transport. SIAM Journal on Imaging Sciences, 7(3), 1853-1882. -[8] M. Perrot, N. Courty, R. Flamary, A. Habrard, "Mapping estimation for - discrete optimal transport", Neural Information Processing Systems (NIPS), 2016. - -""" - -# Author: Remi Flamary -# -# License: MIT License - -import numpy as np -from scipy import ndimage -import matplotlib.pylab as pl -import ot - - -#%% Loading images - -I1 = ndimage.imread('../data/ocean_day.jpg').astype(np.float64) / 256 -I2 = ndimage.imread('../data/ocean_sunset.jpg').astype(np.float64) / 256 - -#%% Plot images - -pl.figure(1, figsize=(6.4, 3)) -pl.subplot(1, 2, 1) -pl.imshow(I1) -pl.axis('off') -pl.title('Image 1') - -pl.subplot(1, 2, 2) -pl.imshow(I2) -pl.axis('off') -pl.title('Image 2') -pl.tight_layout() - - -#%% Image conversion and dataset generation - -def im2mat(I): - """Converts and image to matrix (one pixel per line)""" - return I.reshape((I.shape[0] * I.shape[1], I.shape[2])) - - -def mat2im(X, shape): - """Converts back a matrix to an image""" - return X.reshape(shape) - - -X1 = im2mat(I1) -X2 = im2mat(I2) - -# training samples -nb = 1000 -idx1 = np.random.randint(X1.shape[0], size=(nb,)) -idx2 = np.random.randint(X2.shape[0], size=(nb,)) - -xs = X1[idx1, :] -xt = X2[idx2, :] - -#%% Plot image distributions - - -pl.figure(2, figsize=(6.4, 5)) - -pl.subplot(1, 2, 1) -pl.scatter(xs[:, 0], xs[:, 2], c=xs) -pl.axis([0, 1, 0, 1]) -pl.xlabel('Red') -pl.ylabel('Blue') -pl.title('Image 1') - -pl.subplot(1, 2, 2) -pl.scatter(xt[:, 0], xt[:, 2], c=xt) -pl.axis([0, 1, 0, 1]) -pl.xlabel('Red') -pl.ylabel('Blue') -pl.title('Image 2') -pl.tight_layout() - - -#%% domain adaptation between images - -def minmax(I): - return np.clip(I, 0, 1) - - -# LP problem -da_emd = ot.da.OTDA() # init class -da_emd.fit(xs, xt) # fit distributions - -X1t = da_emd.predict(X1) # out of sample -I1t = minmax(mat2im(X1t, I1.shape)) - -# sinkhorn regularization -lambd = 1e-1 -da_entrop = ot.da.OTDA_sinkhorn() -da_entrop.fit(xs, xt, reg=lambd) - -X1te = da_entrop.predict(X1) -I1te = minmax(mat2im(X1te, I1.shape)) - -# linear mapping estimation -eta = 1e-8 # quadratic regularization for regression -mu = 1e0 # weight of the OT linear term -bias = True # estimate a bias - -ot_mapping = ot.da.OTDA_mapping_linear() -ot_mapping.fit(xs, xt, mu=mu, eta=eta, bias=bias, numItermax=20, verbose=True) - -X1tl = ot_mapping.predict(X1) # use the estimated mapping -I1tl = minmax(mat2im(X1tl, I1.shape)) - -# nonlinear mapping estimation -eta = 1e-2 # quadratic regularization for regression -mu = 1e0 # weight of the OT linear term -bias = False # estimate a bias -sigma = 1 # sigma bandwidth fot gaussian kernel - - -ot_mapping_kernel = ot.da.OTDA_mapping_kernel() -ot_mapping_kernel.fit( - xs, xt, mu=mu, eta=eta, sigma=sigma, bias=bias, numItermax=10, verbose=True) - -X1tn = ot_mapping_kernel.predict(X1) # use the estimated mapping -I1tn = minmax(mat2im(X1tn, I1.shape)) - -#%% plot images - -pl.figure(2, figsize=(8, 4)) - -pl.subplot(2, 3, 1) -pl.imshow(I1) -pl.axis('off') -pl.title('Im. 1') - -pl.subplot(2, 3, 2) -pl.imshow(I2) -pl.axis('off') -pl.title('Im. 2') - -pl.subplot(2, 3, 3) -pl.imshow(I1t) -pl.axis('off') -pl.title('Im. 1 Interp LP') - -pl.subplot(2, 3, 4) -pl.imshow(I1te) -pl.axis('off') -pl.title('Im. 1 Interp Entrop') - -pl.subplot(2, 3, 5) -pl.imshow(I1tl) -pl.axis('off') -pl.title('Im. 1 Linear mapping') - -pl.subplot(2, 3, 6) -pl.imshow(I1tn) -pl.axis('off') -pl.title('Im. 1 nonlinear mapping') -pl.tight_layout() - -pl.show() diff --git a/ot/da.py b/ot/da.py index 8c62669..369b6a2 100644 --- a/ot/da.py +++ b/ot/da.py @@ -976,36 +976,41 @@ class BaseTransport(BaseEstimator): Returns self. """ - # pairwise distance - self.cost_ = dist(Xs, Xt, metric=self.metric) + if Xs is not None and Xt is not None: + # pairwise distance + self.cost_ = dist(Xs, Xt, metric=self.metric) - if (ys is not None) and (yt is not None): + if (ys is not None) and (yt is not None): - if self.limit_max != np.infty: - self.limit_max = self.limit_max * np.max(self.cost_) + if self.limit_max != np.infty: + self.limit_max = self.limit_max * np.max(self.cost_) - # assumes labeled source samples occupy the first rows - # and labeled target samples occupy the first columns - classes = np.unique(ys) - for c in classes: - idx_s = np.where((ys != c) & (ys != -1)) - idx_t = np.where(yt == c) + # assumes labeled source samples occupy the first rows + # and labeled target samples occupy the first columns + classes = np.unique(ys) + for c in classes: + idx_s = np.where((ys != c) & (ys != -1)) + idx_t = np.where(yt == c) - # all the coefficients corresponding to a source sample - # and a target sample : - # with different labels get a infinite - for j in idx_t[0]: - self.cost_[idx_s[0], j] = self.limit_max + # all the coefficients corresponding to a source sample + # and a target sample : + # with different labels get a infinite + for j in idx_t[0]: + self.cost_[idx_s[0], j] = self.limit_max - # distribution estimation - self.mu_s = self.distribution_estimation(Xs) - self.mu_t = self.distribution_estimation(Xt) + # distribution estimation + self.mu_s = self.distribution_estimation(Xs) + self.mu_t = self.distribution_estimation(Xt) - # store arrays of samples - self.Xs = Xs - self.Xt = Xt + # store arrays of samples + self.Xs = Xs + self.Xt = Xt - return self + return self + else: + print("POT-Warning") + print("Please provide both Xs and Xt arguments when calling") + print("fit method") def fit_transform(self, Xs=None, ys=None, Xt=None, yt=None): """Build a coupling matrix from source and target sets of samples @@ -1053,42 +1058,47 @@ class BaseTransport(BaseEstimator): The transport source samples. """ - if np.array_equal(self.Xs, Xs): - # perform standard barycentric mapping - transp = self.coupling_ / np.sum(self.coupling_, 1)[:, None] + if Xs is not None: + if np.array_equal(self.Xs, Xs): + # perform standard barycentric mapping + transp = self.coupling_ / np.sum(self.coupling_, 1)[:, None] - # set nans to 0 - transp[~ np.isfinite(transp)] = 0 + # set nans to 0 + transp[~ np.isfinite(transp)] = 0 - # compute transported samples - transp_Xs = np.dot(transp, self.Xt) - else: - # perform out of sample mapping - indices = np.arange(Xs.shape[0]) - batch_ind = [ - indices[i:i + batch_size] - for i in range(0, len(indices), batch_size)] + # compute transported samples + transp_Xs = np.dot(transp, self.Xt) + else: + # perform out of sample mapping + indices = np.arange(Xs.shape[0]) + batch_ind = [ + indices[i:i + batch_size] + for i in range(0, len(indices), batch_size)] - transp_Xs = [] - for bi in batch_ind: + transp_Xs = [] + for bi in batch_ind: - # get the nearest neighbor in the source domain - D0 = dist(Xs[bi], self.Xs) - idx = np.argmin(D0, axis=1) + # get the nearest neighbor in the source domain + D0 = dist(Xs[bi], self.Xs) + idx = np.argmin(D0, axis=1) - # transport the source samples - transp = self.coupling_ / np.sum(self.coupling_, 1)[:, None] - transp[~ np.isfinite(transp)] = 0 - transp_Xs_ = np.dot(transp, self.Xt) + # transport the source samples + transp = self.coupling_ / np.sum( + self.coupling_, 1)[:, None] + transp[~ np.isfinite(transp)] = 0 + transp_Xs_ = np.dot(transp, self.Xt) - # define the transported points - transp_Xs_ = transp_Xs_[idx, :] + Xs[bi] - self.Xs[idx, :] + # define the transported points + transp_Xs_ = transp_Xs_[idx, :] + Xs[bi] - self.Xs[idx, :] - transp_Xs.append(transp_Xs_) + transp_Xs.append(transp_Xs_) - transp_Xs = np.concatenate(transp_Xs, axis=0) + transp_Xs = np.concatenate(transp_Xs, axis=0) - return transp_Xs + return transp_Xs + else: + print("POT-Warning") + print("Please provide Xs argument when calling transform method") def inverse_transform(self, Xs=None, ys=None, Xt=None, yt=None, batch_size=128): @@ -1113,41 +1123,46 @@ class BaseTransport(BaseEstimator): The transported target samples. """ - if np.array_equal(self.Xt, Xt): - # perform standard barycentric mapping - transp_ = self.coupling_.T / np.sum(self.coupling_, 0)[:, None] + if Xt is not None: + if np.array_equal(self.Xt, Xt): + # perform standard barycentric mapping + transp_ = self.coupling_.T / np.sum(self.coupling_, 0)[:, None] - # set nans to 0 - transp_[~ np.isfinite(transp_)] = 0 + # set nans to 0 + transp_[~ np.isfinite(transp_)] = 0 - # compute transported samples - transp_Xt = np.dot(transp_, self.Xs) - else: - # perform out of sample mapping - indices = np.arange(Xt.shape[0]) - batch_ind = [ - indices[i:i + batch_size] - for i in range(0, len(indices), batch_size)] + # compute transported samples + transp_Xt = np.dot(transp_, self.Xs) + else: + # perform out of sample mapping + indices = np.arange(Xt.shape[0]) + batch_ind = [ + indices[i:i + batch_size] + for i in range(0, len(indices), batch_size)] - transp_Xt = [] - for bi in batch_ind: + transp_Xt = [] + for bi in batch_ind: - D0 = dist(Xt[bi], self.Xt) - idx = np.argmin(D0, axis=1) + D0 = dist(Xt[bi], self.Xt) + idx = np.argmin(D0, axis=1) - # transport the target samples - transp_ = self.coupling_.T / np.sum(self.coupling_, 0)[:, None] - transp_[~ np.isfinite(transp_)] = 0 - transp_Xt_ = np.dot(transp_, self.Xs) + # transport the target samples + transp_ = self.coupling_.T / np.sum( + self.coupling_, 0)[:, None] + transp_[~ np.isfinite(transp_)] = 0 + transp_Xt_ = np.dot(transp_, self.Xs) - # define the transported points - transp_Xt_ = transp_Xt_[idx, :] + Xt[bi] - self.Xt[idx, :] + # define the transported points + transp_Xt_ = transp_Xt_[idx, :] + Xt[bi] - self.Xt[idx, :] - transp_Xt.append(transp_Xt_) + transp_Xt.append(transp_Xt_) - transp_Xt = np.concatenate(transp_Xt, axis=0) + transp_Xt = np.concatenate(transp_Xt, axis=0) - return transp_Xt + return transp_Xt + else: + print("POT-Warning") + print("Please provide Xt argument when calling inverse_transform") class SinkhornTransport(BaseTransport): @@ -1413,15 +1428,20 @@ class SinkhornLpl1Transport(BaseTransport): Returns self. """ - super(SinkhornLpl1Transport, self).fit(Xs, ys, Xt, yt) + if Xs is not None and Xt is not None and ys is not None: - self.coupling_ = sinkhorn_lpl1_mm( - a=self.mu_s, labels_a=ys, b=self.mu_t, M=self.cost_, - reg=self.reg_e, eta=self.reg_cl, numItermax=self.max_iter, - numInnerItermax=self.max_inner_iter, stopInnerThr=self.tol, - verbose=self.verbose) + super(SinkhornLpl1Transport, self).fit(Xs, ys, Xt, yt) - return self + self.coupling_ = sinkhorn_lpl1_mm( + a=self.mu_s, labels_a=ys, b=self.mu_t, M=self.cost_, + reg=self.reg_e, eta=self.reg_cl, numItermax=self.max_iter, + numInnerItermax=self.max_inner_iter, stopInnerThr=self.tol, + verbose=self.verbose) + + return self + else: + print("POT-Warning") + print("Please provide both Xs, Xt, ys arguments to fit method") class SinkhornL1l2Transport(BaseTransport): @@ -1517,22 +1537,27 @@ class SinkhornL1l2Transport(BaseTransport): Returns self. """ - super(SinkhornL1l2Transport, self).fit(Xs, ys, Xt, yt) + if Xs is not None and Xt is not None and ys is not None: - returned_ = sinkhorn_l1l2_gl( - a=self.mu_s, labels_a=ys, b=self.mu_t, M=self.cost_, - reg=self.reg_e, eta=self.reg_cl, numItermax=self.max_iter, - numInnerItermax=self.max_inner_iter, stopInnerThr=self.tol, - verbose=self.verbose, log=self.log) + super(SinkhornL1l2Transport, self).fit(Xs, ys, Xt, yt) - # deal with the value of log - if self.log: - self.coupling_, self.log_ = returned_ - else: - self.coupling_ = returned_ - self.log_ = dict() + returned_ = sinkhorn_l1l2_gl( + a=self.mu_s, labels_a=ys, b=self.mu_t, M=self.cost_, + reg=self.reg_e, eta=self.reg_cl, numItermax=self.max_iter, + numInnerItermax=self.max_inner_iter, stopInnerThr=self.tol, + verbose=self.verbose, log=self.log) - return self + # deal with the value of log + if self.log: + self.coupling_, self.log_ = returned_ + else: + self.coupling_ = returned_ + self.log_ = dict() + + return self + else: + print("POT-Warning") + print("Please, provide both Xs, Xt and ys argument to fit method") class MappingTransport(BaseEstimator): -- cgit v1.2.3 From 7ab9037f1e4a08439083d1bc5705be5ed2e9e10a Mon Sep 17 00:00:00 2001 From: Nicolas Courty Date: Mon, 28 Aug 2017 14:41:09 +0200 Subject: Gromov-Wasserstein distance --- README.md | 4 +- examples/plot_gromov.py | 98 ++++++++++ ot/__init__.py | 6 +- ot/gromov.py | 482 ++++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 588 insertions(+), 2 deletions(-) create mode 100644 examples/plot_gromov.py create mode 100644 ot/gromov.py (limited to 'examples') diff --git a/README.md b/README.md index 7a65106..a42f17b 100644 --- a/README.md +++ b/README.md @@ -16,7 +16,7 @@ It provides the following solvers: * Conditional gradient [6] and Generalized conditional gradient for regularized OT [7]. * Joint OT matrix and mapping estimation [8]. * Wasserstein Discriminant Analysis [11] (requires autograd + pymanopt). - +* Gromov-Wasserstein distances [12] Some demonstrations (both in Python and Jupyter Notebook format) are available in the examples folder. @@ -182,3 +182,5 @@ You can also post bug reports and feature requests in Github issues. Make sure t [10] Chizat, L., Peyré, G., Schmitzer, B., & Vialard, F. X. (2016). [Scaling algorithms for unbalanced transport problems](https://arxiv.org/pdf/1607.05816.pdf). arXiv preprint arXiv:1607.05816. [11] Flamary, R., Cuturi, M., Courty, N., & Rakotomamonjy, A. (2016). [Wasserstein Discriminant Analysis](https://arxiv.org/pdf/1608.08063.pdf). arXiv preprint arXiv:1608.08063. + +[12] Peyré, Gabriel, Marco Cuturi, and Justin Solomon, [Gromov-Wasserstein averaging of kernel and distance matrices](http://proceedings.mlr.press/v48/peyre16.html) International Conference on Machine Learning (ICML). 2016. diff --git a/examples/plot_gromov.py b/examples/plot_gromov.py new file mode 100644 index 0000000..11e5336 --- /dev/null +++ b/examples/plot_gromov.py @@ -0,0 +1,98 @@ +# -*- coding: utf-8 -*- +""" +==================== +Gromov-Wasserstein example +==================== + +This example is designed to show how to use the Gromov-Wassertsein distance +computation in POT. + + +""" + +# Author: Erwan Vautier +# Nicolas Courty +# +# License: MIT License + +import scipy as sp +import numpy as np + +import ot +import matplotlib.pylab as pl +from mpl_toolkits.mplot3d import Axes3D + + + +""" +Sample two Gaussian distributions (2D and 3D) +==================== + +The Gromov-Wasserstein distance allows to compute distances with samples that do not belong to the same metric space. For +demonstration purpose, we sample two Gaussian distributions in 2- and 3-dimensional spaces. + +""" +n=30 # nb samples + +mu_s=np.array([0,0]) +cov_s=np.array([[1,0],[0,1]]) + +mu_t=np.array([4,4,4]) +cov_t=np.array([[1,0,0],[0,1,0],[0,0,1]]) + + + +xs=ot.datasets.get_2D_samples_gauss(n,mu_s,cov_s) +P=sp.linalg.sqrtm(cov_t) +xt= np.random.randn(n,3).dot(P)+mu_t + + + +""" +Plotting the distributions +==================== +""" +fig=pl.figure() +ax1=fig.add_subplot(121) +ax1.plot(xs[:,0],xs[:,1],'+b',label='Source samples') +ax2=fig.add_subplot(122,projection='3d') +ax2.scatter(xt[:,0],xt[:,1],xt[:,2],color='r') +pl.show() + + +""" +Compute distance kernels, normalize them and then display +==================== +""" + +C1=sp.spatial.distance.cdist(xs,xs) +C2=sp.spatial.distance.cdist(xt,xt) + +C1/=C1.max() +C2/=C2.max() + +pl.figure() +pl.subplot(121) +pl.imshow(C1) +pl.subplot(122) +pl.imshow(C2) +pl.show() + +""" +Compute Gromov-Wasserstein plans and distance +==================== +""" + +p=ot.unif(n) +q=ot.unif(n) + +gw=ot.gromov_wasserstein(C1,C2,p,q,'square_loss',epsilon=5e-4) +gw_dist=ot.gromov_wasserstein2(C1,C2,p,q,'square_loss',epsilon=5e-4) + +print('Gromov-Wasserstein distances between the distribution: '+str(gw_dist)) + +pl.figure() +pl.imshow(gw,cmap='jet') +pl.colorbar() +pl.show() + diff --git a/ot/__init__.py b/ot/__init__.py index 6d4c4c6..a295e1b 100644 --- a/ot/__init__.py +++ b/ot/__init__.py @@ -5,6 +5,7 @@ """ # Author: Remi Flamary +# Nicolas Courty # # License: MIT License @@ -17,11 +18,13 @@ from . import utils from . import datasets from . import plot from . import da +from . import gromov # OT functions from .lp import emd, emd2 from .bregman import sinkhorn, sinkhorn2, barycenter from .da import sinkhorn_lpl1_mm +from .gromov import gromov_wasserstein, gromov_wasserstein2 # utils functions from .utils import dist, unif, tic, toc, toq @@ -30,4 +33,5 @@ __version__ = "0.3.1" __all__ = ["emd", "emd2", "sinkhorn", "sinkhorn2", "utils", 'datasets', 'bregman', 'lp', 'plot', 'tic', 'toc', 'toq', - 'dist', 'unif', 'barycenter', 'sinkhorn_lpl1_mm', 'da', 'optim'] + 'dist', 'unif', 'barycenter', 'sinkhorn_lpl1_mm', 'da', 'optim', + 'gromov_wasserstein','gromov_wasserstein2'] diff --git a/ot/gromov.py b/ot/gromov.py new file mode 100644 index 0000000..f3c62c9 --- /dev/null +++ b/ot/gromov.py @@ -0,0 +1,482 @@ + +# -*- coding: utf-8 -*- +""" +Gromov-Wasserstein transport method +=================================== + + +""" + +# Author: Erwan Vautier +# Nicolas Courty +# +# License: MIT License + +import numpy as np + +from .bregman import sinkhorn +from .utils import dist + +def square_loss(a,b): + """ + Returns the value of L(a,b)=(1/2)*|a-b|^2 + """ + + return (1/2)*(a-b)**2 + +def kl_loss(a,b): + """ + Returns the value of L(a,b)=a*log(a/b)-a+b + """ + + return a*np.log(a/b)-a+b + +def tensor_square_loss(C1,C2,T): + """ + Returns the value of \mathcal{L}(C1,C2) \otimes T with the square loss + + function as the loss function of Gromow-Wasserstein discrepancy. + + Where : + + C1 : Metric cost matrix in the source space + C2 : Metric cost matrix in the target space + T : A coupling between those two spaces + + The square-loss function L(a,b)=(1/2)*|a-b|^2 is read as : + L(a,b) = f1(a)+f2(b)-h1(a)*h2(b) with : + f1(a)=(a^2)/2 + f2(b)=(b^2)/2 + h1(a)=a + h2(b)=b + + Parameters + ---------- + C1 : np.ndarray(ns,ns) + Metric cost matrix in the source space + C2 : np.ndarray(nt,nt) + Metric costfr matrix in the target space + T : np.ndarray(ns,nt) + Coupling between source and target spaces + + + Returns + ------- + tens : (ns*nt) ndarray + \mathcal{L}(C1,C2) \otimes T tensor-matrix multiplication result + + + """ + + + C1=np.asarray(C1,dtype=np.float64) + C2=np.asarray(C2,dtype=np.float64) + T=np.asarray(T,dtype=np.float64) + + + def f1(a): + return (a**2)/2 + + def f2(b): + return (b**2)/2 + + def h1(a): + return a + + def h2(b): + return b + + tens=-np.dot(h1(C1),T).dot(h2(C2).T) + tens=tens-tens.min() + + + return np.array(tens) + + +def tensor_kl_loss(C1,C2,T): + """ + Returns the value of \mathcal{L}(C1,C2) \otimes T with the square loss + + function as the loss function of Gromow-Wasserstein discrepancy. + + Where : + + C1 : Metric cost matrix in the source space + C2 : Metric cost matrix in the target space + T : A coupling between those two spaces + + The square-loss function L(a,b)=(1/2)*|a-b|^2 is read as : + L(a,b) = f1(a)+f2(b)-h1(a)*h2(b) with : + f1(a)=a*log(a)-a + f2(b)=b + h1(a)=a + h2(b)=log(b) + + Parameters + ---------- + C1 : np.ndarray(ns,ns) + Metric cost matrix in the source space + C2 : np.ndarray(nt,nt) + Metric costfr matrix in the target space + T : np.ndarray(ns,nt) + Coupling between source and target spaces + + + Returns + ------- + tens : (ns*nt) ndarray + \mathcal{L}(C1,C2) \otimes T tensor-matrix multiplication result + + References + ---------- + + .. [12] Peyré, Gabriel, Marco Cuturi, and Justin Solomon, "Gromov-Wasserstein averaging of kernel and distance matrices." International Conference on Machine Learning (ICML). 2016. + + """ + + + C1=np.asarray(C1,dtype=np.float64) + C2=np.asarray(C2,dtype=np.float64) + T=np.asarray(T,dtype=np.float64) + + + def f1(a): + return a*np.log(a+1e-15)-a + + def f2(b): + return b + + def h1(a): + return a + + def h2(b): + return np.log(b+1e-15) + + tens=-np.dot(h1(C1),T).dot(h2(C2).T) + tens=tens-tens.min() + + + + return np.array(tens) + +def update_square_loss(p,lambdas,T,Cs): + """ + Updates C according to the L2 Loss kernel with the S Ts couplings calculated at each iteration + + + Parameters + ---------- + p : np.ndarray(N,) + weights in the targeted barycenter + lambdas : list of the S spaces' weights + T : list of S np.ndarray(ns,N) + the S Ts couplings calculated at each iteration + Cs : Cs : list of S np.ndarray(ns,ns) + Metric cost matrices + + Returns + ---------- + C updated + + + """ + tmpsum=np.sum([lambdas[s]*np.dot(T[s].T,Cs[s]).dot(T[s]) for s in range(len(T))]) + ppt=np.dot(p,p.T) + + return(np.divide(tmpsum,ppt)) + +def update_kl_loss(p,lambdas,T,Cs): + """ + Updates C according to the KL Loss kernel with the S Ts couplings calculated at each iteration + + + Parameters + ---------- + p : np.ndarray(N,) + weights in the targeted barycenter + lambdas : list of the S spaces' weights + T : list of S np.ndarray(ns,N) + the S Ts couplings calculated at each iteration + Cs : Cs : list of S np.ndarray(ns,ns) + Metric cost matrices + + Returns + ---------- + C updated + + + """ + tmpsum=np.sum([lambdas[s]*np.dot(T[s].T,Cs[s]).dot(T[s]) for s in range(len(T))]) + ppt=np.dot(p,p.T) + + return(np.exp(np.divide(tmpsum,ppt))) + + +def gromov_wasserstein(C1,C2,p,q,loss_fun,epsilon,numItermax = 1000, stopThr=1e-9,verbose=False, log=False): + """ + Returns the gromov-wasserstein coupling between the two measured similarity matrices + + (C1,p) and (C2,q) + + The function solves the following optimization problem: + + .. math:: + \GW = arg\min_T \sum_{i,j,k,l} L(C1_{i,k},C2_{j,l})*T_{i,j}*T_{k,l}-\epsilon(H(T)) + + s.t. \GW 1 = p + + \GW^T 1= q + + \GW\geq 0 + + Where : + + C1 : Metric cost matrix in the source space + C2 : Metric cost matrix in the target space + p : distribution in the source space + q : distribution in the target space + L : loss function to account for the misfit between the similarity matrices + H : entropy + + + Parameters + ---------- + C1 : np.ndarray(ns,ns) + Metric cost matrix in the source space + C2 : np.ndarray(nt,nt) + Metric costfr matrix in the target space + p : np.ndarray(ns,) + distribution in the source space + q : np.ndarray(nt) + distribution in the target space + loss_fun : loss function used for the solver either 'square_loss' or 'kl_loss' + epsilon : float + Regularization term >0 + numItermax : int, optional + Max number of iterations + stopThr : float, optional + Stop threshold on error (>0) + verbose : bool, optional + Print information along iterations + log : bool, optional + record log if True + forcing : np.ndarray(N,2) + list of forced couplings (where N is the number of forcing) + + Returns + ------- + T : coupling between the two spaces that minimizes : + \sum_{i,j,k,l} L(C1_{i,k},C2_{j,l})*T_{i,j}*T_{k,l}-\epsilon(H(T)) + + """ + + + C1=np.asarray(C1,dtype=np.float64) + C2=np.asarray(C2,dtype=np.float64) + + T=np.dot(p,q.T) #Initialization + + cpt = 0 + err=1 + + while (err>stopThr and cpt0 + numItermax : int, optional + Max number of iterations + stopThr : float, optional + Stop threshold on error (>0) + verbose : bool, optional + Print information along iterations + log : bool, optional + record log if True + forcing : np.ndarray(N,2) + list of forced couplings (where N is the number of forcing) + + Returns + ------- + T : coupling between the two spaces that minimizes : + \sum_{i,j,k,l} L(C1_{i,k},C2_{j,l})*T_{i,j}*T_{k,l}-\epsilon(H(T)) + + """ + + if log: + gw,logv=gromov_wasserstein(C1,C2,p,q,loss_fun,epsilon,numItermax,stopThr,verbose,log) + else: + gw=gromov_wasserstein(C1,C2,p,q,loss_fun,epsilon,numItermax,stopThr,verbose,log) + + if loss_fun=='square_loss': + gw_dist=np.sum(gw*tensor_square_loss(C1,C2,gw)) + + + elif loss_fun=='kl_loss': + gw_dist=np.sum(gw*tensor_kl_loss(C1,C2,gw)) + + if log: + return gw_dist, logv + else: + return gw_dist + + +def gromov_barycenters(N,Cs,ps,p,lambdas,loss_fun,epsilon,numItermax = 1000, stopThr=1e-9, verbose=False, log=False): + """ + Returns the gromov-wasserstein barycenters of S measured similarity matrices + + (Cs)_{s=1}^{s=S} + + The function solves the following optimization problem: + + .. math:: + C = argmin_C\in R^NxN \sum_s \lambda_s GW(C,Cs,p,ps) + + + Where : + + Cs : metric cost matrix + ps : distribution + + Parameters + ---------- + N : Integer + Size of the targeted barycenter + Cs : list of S np.ndarray(ns,ns) + Metric cost matrices + ps : list of S np.ndarray(ns,) + sample weights in the S spaces + p : np.ndarray(N,) + weights in the targeted barycenter + lambdas : list of the S spaces' weights + L : tensor-matrix multiplication function based on specific loss function + update : function(p,lambdas,T,Cs) that updates C according to a specific Kernel + with the S Ts couplings calculated at each iteration + epsilon : float + Regularization term >0 + numItermax : int, optional + Max number of iterations + stopThr : float, optional + Stop threshol on error (>0) + verbose : bool, optional + Print information along iterations + log : bool, optional + record log if True + + Returns + ------- + C : Similarity matrix in the barycenter space (permutated arbitrarily) + + """ + + + S=len(Cs) + + Cs=[np.asarray(Cs[s],dtype=np.float64) for s in range(S)] + lambdas=np.asarray(lambdas,dtype=np.float64) + + T=[0 for s in range(S)] + + #Initialization of C : random SPD matrix + xalea=np.random.randn(N,2) + C=dist(xalea,xalea) + C/=C.max() + + cpt=0 + err=1 + + error=[] + + while(err>stopThr and cpt Date: Mon, 28 Aug 2017 15:04:04 +0200 Subject: gromov:flake8 and other --- examples/plot_gromov.py | 63 ++++---- ot/gromov.py | 382 ++++++++++++++++++++++++------------------------ 2 files changed, 216 insertions(+), 229 deletions(-) (limited to 'examples') diff --git a/examples/plot_gromov.py b/examples/plot_gromov.py index 11e5336..a33fde1 100644 --- a/examples/plot_gromov.py +++ b/examples/plot_gromov.py @@ -3,11 +3,8 @@ ==================== Gromov-Wasserstein example ==================== - -This example is designed to show how to use the Gromov-Wassertsein distance -computation in POT. - - +This example is designed to show how to use the Gromov-Wassertsein distance +computation in POT. """ # Author: Erwan Vautier @@ -20,43 +17,38 @@ import numpy as np import ot import matplotlib.pylab as pl -from mpl_toolkits.mplot3d import Axes3D - """ Sample two Gaussian distributions (2D and 3D) ==================== - -The Gromov-Wasserstein distance allows to compute distances with samples that do not belong to the same metric space. For -demonstration purpose, we sample two Gaussian distributions in 2- and 3-dimensional spaces. - +The Gromov-Wasserstein distance allows to compute distances with samples that do not belong to the same metric space. +For demonstration purpose, we sample two Gaussian distributions in 2- and 3-dimensional spaces. """ -n=30 # nb samples -mu_s=np.array([0,0]) -cov_s=np.array([[1,0],[0,1]]) +n = 30 # nb samples -mu_t=np.array([4,4,4]) -cov_t=np.array([[1,0,0],[0,1,0],[0,0,1]]) +mu_s = np.array([0, 0]) +cov_s = np.array([[1, 0], [0, 1]]) +mu_t = np.array([4, 4, 4]) +cov_t = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) -xs=ot.datasets.get_2D_samples_gauss(n,mu_s,cov_s) -P=sp.linalg.sqrtm(cov_t) -xt= np.random.randn(n,3).dot(P)+mu_t - +xs = ot.datasets.get_2D_samples_gauss(n, mu_s, cov_s) +P = sp.linalg.sqrtm(cov_t) +xt = np.random.randn(n, 3).dot(P) + mu_t """ Plotting the distributions ==================== """ -fig=pl.figure() -ax1=fig.add_subplot(121) -ax1.plot(xs[:,0],xs[:,1],'+b',label='Source samples') -ax2=fig.add_subplot(122,projection='3d') -ax2.scatter(xt[:,0],xt[:,1],xt[:,2],color='r') +fig = pl.figure() +ax1 = fig.add_subplot(121) +ax1.plot(xs[:, 0], xs[:, 1], '+b', label='Source samples') +ax2 = fig.add_subplot(122, projection='3d') +ax2.scatter(xt[:, 0], xt[:, 1], xt[:, 2], color='r') pl.show() @@ -65,11 +57,11 @@ Compute distance kernels, normalize them and then display ==================== """ -C1=sp.spatial.distance.cdist(xs,xs) -C2=sp.spatial.distance.cdist(xt,xt) +C1 = sp.spatial.distance.cdist(xs, xs) +C2 = sp.spatial.distance.cdist(xt, xt) -C1/=C1.max() -C2/=C2.max() +C1 /= C1.max() +C2 /= C2.max() pl.figure() pl.subplot(121) @@ -83,16 +75,15 @@ Compute Gromov-Wasserstein plans and distance ==================== """ -p=ot.unif(n) -q=ot.unif(n) +p = ot.unif(n) +q = ot.unif(n) -gw=ot.gromov_wasserstein(C1,C2,p,q,'square_loss',epsilon=5e-4) -gw_dist=ot.gromov_wasserstein2(C1,C2,p,q,'square_loss',epsilon=5e-4) +gw = ot.gromov_wasserstein(C1, C2, p, q, 'square_loss', epsilon=5e-4) +gw_dist = ot.gromov_wasserstein2(C1, C2, p, q, 'square_loss', epsilon=5e-4) -print('Gromov-Wasserstein distances between the distribution: '+str(gw_dist)) +print('Gromov-Wasserstein distances between the distribution: ' + str(gw_dist)) pl.figure() -pl.imshow(gw,cmap='jet') +pl.imshow(gw, cmap='jet') pl.colorbar() pl.show() - diff --git a/ot/gromov.py b/ot/gromov.py index f3c62c9..7cf3b42 100644 --- a/ot/gromov.py +++ b/ot/gromov.py @@ -17,39 +17,41 @@ import numpy as np from .bregman import sinkhorn from .utils import dist -def square_loss(a,b): + +def square_loss(a, b): """ Returns the value of L(a,b)=(1/2)*|a-b|^2 """ - - return (1/2)*(a-b)**2 -def kl_loss(a,b): + return (1 / 2) * (a - b)**2 + + +def kl_loss(a, b): """ Returns the value of L(a,b)=a*log(a/b)-a+b """ - - return a*np.log(a/b)-a+b -def tensor_square_loss(C1,C2,T): + return a * np.log(a / b) - a + b + + +def tensor_square_loss(C1, C2, T): """ - Returns the value of \mathcal{L}(C1,C2) \otimes T with the square loss - + Returns the value of \mathcal{L}(C1,C2) \otimes T with the square loss function as the loss function of Gromow-Wasserstein discrepancy. - + Where : - + C1 : Metric cost matrix in the source space C2 : Metric cost matrix in the target space T : A coupling between those two spaces - + The square-loss function L(a,b)=(1/2)*|a-b|^2 is read as : L(a,b) = f1(a)+f2(b)-h1(a)*h2(b) with : f1(a)=(a^2)/2 f2(b)=(b^2)/2 h1(a)=a h2(b)=b - + Parameters ---------- C1 : np.ndarray(ns,ns) @@ -64,54 +66,50 @@ def tensor_square_loss(C1,C2,T): ------- tens : (ns*nt) ndarray \mathcal{L}(C1,C2) \otimes T tensor-matrix multiplication result - - + + """ - - C1=np.asarray(C1,dtype=np.float64) - C2=np.asarray(C2,dtype=np.float64) - T=np.asarray(T,dtype=np.float64) - - + C1 = np.asarray(C1, dtype=np.float64) + C2 = np.asarray(C2, dtype=np.float64) + T = np.asarray(T, dtype=np.float64) + def f1(a): - return (a**2)/2 - + return (a**2) / 2 + def f2(b): - return (b**2)/2 - + return (b**2) / 2 + def h1(a): return a - + def h2(b): return b - - tens=-np.dot(h1(C1),T).dot(h2(C2).T) - tens=tens-tens.min() - + tens = -np.dot(h1(C1), T).dot(h2(C2).T) + tens = tens - tens.min() + return np.array(tens) - - -def tensor_kl_loss(C1,C2,T): + + +def tensor_kl_loss(C1, C2, T): """ - Returns the value of \mathcal{L}(C1,C2) \otimes T with the square loss - + Returns the value of \mathcal{L}(C1,C2) \otimes T with the square loss function as the loss function of Gromow-Wasserstein discrepancy. - + Where : - + C1 : Metric cost matrix in the source space C2 : Metric cost matrix in the target space T : A coupling between those two spaces - + The square-loss function L(a,b)=(1/2)*|a-b|^2 is read as : L(a,b) = f1(a)+f2(b)-h1(a)*h2(b) with : f1(a)=a*log(a)-a f2(b)=b h1(a)=a h2(b)=log(b) - + Parameters ---------- C1 : np.ndarray(ns,ns) @@ -126,44 +124,41 @@ def tensor_kl_loss(C1,C2,T): ------- tens : (ns*nt) ndarray \mathcal{L}(C1,C2) \otimes T tensor-matrix multiplication result - + References ---------- - + .. [12] Peyré, Gabriel, Marco Cuturi, and Justin Solomon, "Gromov-Wasserstein averaging of kernel and distance matrices." International Conference on Machine Learning (ICML). 2016. - + """ - - C1=np.asarray(C1,dtype=np.float64) - C2=np.asarray(C2,dtype=np.float64) - T=np.asarray(T,dtype=np.float64) - - + C1 = np.asarray(C1, dtype=np.float64) + C2 = np.asarray(C2, dtype=np.float64) + T = np.asarray(T, dtype=np.float64) + def f1(a): - return a*np.log(a+1e-15)-a - + return a * np.log(a + 1e-15) - a + def f2(b): return b - + def h1(a): return a - + def h2(b): - return np.log(b+1e-15) - - tens=-np.dot(h1(C1),T).dot(h2(C2).T) - tens=tens-tens.min() + return np.log(b + 1e-15) + + tens = -np.dot(h1(C1), T).dot(h2(C2).T) + tens = tens - tens.min() - - return np.array(tens) -def update_square_loss(p,lambdas,T,Cs): + +def update_square_loss(p, lambdas, T, Cs): """ Updates C according to the L2 Loss kernel with the S Ts couplings calculated at each iteration - - + + Parameters ---------- p : np.ndarray(N,) @@ -173,23 +168,25 @@ def update_square_loss(p,lambdas,T,Cs): the S Ts couplings calculated at each iteration Cs : Cs : list of S np.ndarray(ns,ns) Metric cost matrices - - Returns + + Returns ---------- C updated - - + + """ - tmpsum=np.sum([lambdas[s]*np.dot(T[s].T,Cs[s]).dot(T[s]) for s in range(len(T))]) - ppt=np.dot(p,p.T) - - return(np.divide(tmpsum,ppt)) + tmpsum = np.sum([lambdas[s] * np.dot(T[s].T, Cs[s]).dot(T[s]) + for s in range(len(T))]) + ppt = np.dot(p, p.T) -def update_kl_loss(p,lambdas,T,Cs): + return(np.divide(tmpsum, ppt)) + + +def update_kl_loss(p, lambdas, T, Cs): """ Updates C according to the KL Loss kernel with the S Ts couplings calculated at each iteration - - + + Parameters ---------- p : np.ndarray(N,) @@ -199,25 +196,26 @@ def update_kl_loss(p,lambdas,T,Cs): the S Ts couplings calculated at each iteration Cs : Cs : list of S np.ndarray(ns,ns) Metric cost matrices - - Returns + + Returns ---------- C updated - - + + """ - tmpsum=np.sum([lambdas[s]*np.dot(T[s].T,Cs[s]).dot(T[s]) for s in range(len(T))]) - ppt=np.dot(p,p.T) - - return(np.exp(np.divide(tmpsum,ppt))) + tmpsum = np.sum([lambdas[s] * np.dot(T[s].T, Cs[s]).dot(T[s]) + for s in range(len(T))]) + ppt = np.dot(p, p.T) + + return(np.exp(np.divide(tmpsum, ppt))) -def gromov_wasserstein(C1,C2,p,q,loss_fun,epsilon,numItermax = 1000, stopThr=1e-9,verbose=False, log=False): +def gromov_wasserstein(C1, C2, p, q, loss_fun, epsilon, numItermax=1000, stopThr=1e-9, verbose=False, log=False): """ Returns the gromov-wasserstein coupling between the two measured similarity matrices - + (C1,p) and (C2,q) - + The function solves the following optimization problem: .. math:: @@ -228,17 +226,17 @@ def gromov_wasserstein(C1,C2,p,q,loss_fun,epsilon,numItermax = 1000, stopThr=1e- \GW^T 1= q \GW\geq 0 - + Where : - + C1 : Metric cost matrix in the source space C2 : Metric cost matrix in the target space p : distribution in the source space - q : distribution in the target space + q : distribution in the target space L : loss function to account for the misfit between the similarity matrices H : entropy - - + + Parameters ---------- C1 : np.ndarray(ns,ns) @@ -259,80 +257,80 @@ def gromov_wasserstein(C1,C2,p,q,loss_fun,epsilon,numItermax = 1000, stopThr=1e- verbose : bool, optional Print information along iterations log : bool, optional - record log if True + record log if True forcing : np.ndarray(N,2) list of forced couplings (where N is the number of forcing) - + Returns ------- T : coupling between the two spaces that minimizes : \sum_{i,j,k,l} L(C1_{i,k},C2_{j,l})*T_{i,j}*T_{k,l}-\epsilon(H(T)) - + """ - - C1=np.asarray(C1,dtype=np.float64) - C2=np.asarray(C2,dtype=np.float64) + C1 = np.asarray(C1, dtype=np.float64) + C2 = np.asarray(C2, dtype=np.float64) + + T = np.dot(p, q.T) # Initialization - T=np.dot(p,q.T) #Initialization - cpt = 0 - err=1 - - while (err>stopThr and cpt stopThr and cpt < numItermax): + + Tprev = T + + if loss_fun == 'square_loss': + tens = tensor_square_loss(C1, C2, T) + + elif loss_fun == 'kl_loss': + tens = tensor_kl_loss(C1, C2, T) + + T = sinkhorn(p, q, tens, epsilon) + + if cpt % 10 == 0: # we can speed up the process by checking for the error only all the 10th iterations - err=np.linalg.norm(T-Tprev) - + err = np.linalg.norm(T - Tprev) + if log: log['err'].append(err) if verbose: - if cpt%200 ==0: - print('{:5s}|{:12s}'.format('It.','Err')+'\n'+'-'*19) - print('{:5d}|{:8e}|'.format(cpt,err)) - - cpt=cpt+1 + if cpt % 200 == 0: + print('{:5s}|{:12s}'.format( + 'It.', 'Err') + '\n' + '-' * 19) + print('{:5d}|{:8e}|'.format(cpt, err)) + + cpt = cpt + 1 if log: - return T,log + return T, log else: return T -def gromov_wasserstein2(C1,C2,p,q,loss_fun,epsilon,numItermax = 1000, stopThr=1e-9,verbose=False, log=False): +def gromov_wasserstein2(C1, C2, p, q, loss_fun, epsilon, numItermax=1000, stopThr=1e-9, verbose=False, log=False): """ Returns the gromov-wasserstein discrepancy between the two measured similarity matrices - + (C1,p) and (C2,q) - + The function solves the following optimization problem: .. math:: \GW_Dist = \min_T \sum_{i,j,k,l} L(C1_{i,k},C2_{j,l})*T_{i,j}*T_{k,l}-\epsilon(H(T)) - + Where : - + C1 : Metric cost matrix in the source space C2 : Metric cost matrix in the target space p : distribution in the source space - q : distribution in the target space + q : distribution in the target space L : loss function to account for the misfit between the similarity matrices H : entropy - - + + Parameters ---------- C1 : np.ndarray(ns,ns) @@ -353,55 +351,56 @@ def gromov_wasserstein2(C1,C2,p,q,loss_fun,epsilon,numItermax = 1000, stopThr=1e verbose : bool, optional Print information along iterations log : bool, optional - record log if True + record log if True forcing : np.ndarray(N,2) list of forced couplings (where N is the number of forcing) - + Returns ------- T : coupling between the two spaces that minimizes : \sum_{i,j,k,l} L(C1_{i,k},C2_{j,l})*T_{i,j}*T_{k,l}-\epsilon(H(T)) - + """ - + if log: - gw,logv=gromov_wasserstein(C1,C2,p,q,loss_fun,epsilon,numItermax,stopThr,verbose,log) + gw, logv = gromov_wasserstein( + C1, C2, p, q, loss_fun, epsilon, numItermax, stopThr, verbose, log) else: - gw=gromov_wasserstein(C1,C2,p,q,loss_fun,epsilon,numItermax,stopThr,verbose,log) + gw = gromov_wasserstein(C1, C2, p, q, loss_fun, + epsilon, numItermax, stopThr, verbose, log) + + if loss_fun == 'square_loss': + gw_dist = np.sum(gw * tensor_square_loss(C1, C2, gw)) - if loss_fun=='square_loss': - gw_dist=np.sum(gw*tensor_square_loss(C1,C2,gw)) - - - elif loss_fun=='kl_loss': - gw_dist=np.sum(gw*tensor_kl_loss(C1,C2,gw)) + elif loss_fun == 'kl_loss': + gw_dist = np.sum(gw * tensor_kl_loss(C1, C2, gw)) if log: return gw_dist, logv - else: + else: return gw_dist -def gromov_barycenters(N,Cs,ps,p,lambdas,loss_fun,epsilon,numItermax = 1000, stopThr=1e-9, verbose=False, log=False): +def gromov_barycenters(N, Cs, ps, p, lambdas, loss_fun, epsilon, numItermax=1000, stopThr=1e-9, verbose=False, log=False): """ Returns the gromov-wasserstein barycenters of S measured similarity matrices - + (Cs)_{s=1}^{s=S} - + The function solves the following optimization problem: .. math:: C = argmin_C\in R^NxN \sum_s \lambda_s GW(C,Cs,p,ps) - + Where : - + Cs : metric cost matrix ps : distribution - + Parameters ---------- - N : Integer + N : Integer Size of the targeted barycenter Cs : list of S np.ndarray(ns,ns) Metric cost matrices @@ -422,61 +421,58 @@ def gromov_barycenters(N,Cs,ps,p,lambdas,loss_fun,epsilon,numItermax = 1000, sto verbose : bool, optional Print information along iterations log : bool, optional - record log if True - + record log if True + Returns ------- C : Similarity matrix in the barycenter space (permutated arbitrarily) - + """ - - - S=len(Cs) - - Cs=[np.asarray(Cs[s],dtype=np.float64) for s in range(S)] - lambdas=np.asarray(lambdas,dtype=np.float64) - - T=[0 for s in range(S)] - - #Initialization of C : random SPD matrix - xalea=np.random.randn(N,2) - C=dist(xalea,xalea) - C/=C.max() - - cpt=0 - err=1 - - error=[] - - while(err>stopThr and cpt stopThr and cpt < numItermax): + + Cprev = C + + T = [gromov_wasserstein(Cs[s], C, ps[s], p, loss_fun, epsilon, + numItermax, 1e-5, verbose, log) for s in range(S)] + + if loss_fun == 'square_loss': + C = update_square_loss(p, lambdas, T, Cs) + + elif loss_fun == 'kl_loss': + C = update_kl_loss(p, lambdas, T, Cs) + + if cpt % 10 == 0: # we can speed up the process by checking for the error only all the 10th iterations - err=np.linalg.norm(C-Cprev) + err = np.linalg.norm(C - Cprev) error.append(err) - + if log: log['err'].append(err) if verbose: - if cpt%200 ==0: - print('{:5s}|{:12s}'.format('It.','Err')+'\n'+'-'*19) - print('{:5d}|{:8e}|'.format(cpt,err)) - - cpt=cpt+1 - - return C - + if cpt % 200 == 0: + print('{:5s}|{:12s}'.format( + 'It.', 'Err') + '\n' + '-' * 19) + print('{:5d}|{:8e}|'.format(cpt, err)) + cpt = cpt + 1 - \ No newline at end of file + return C -- cgit v1.2.3 From a29e22db4772ebc4a8266c917e2e662f624c6baa Mon Sep 17 00:00:00 2001 From: Slasnista Date: Tue, 29 Aug 2017 09:05:01 +0200 Subject: addressed AG comments + adding random seed --- examples/da/plot_otda_classes.py | 2 ++ examples/da/plot_otda_color_images.py | 3 ++- examples/da/plot_otda_d2.py | 14 ++++++++------ examples/da/plot_otda_mapping.py | 14 +++++++------- examples/da/plot_otda_mapping_colors_images.py | 2 ++ 5 files changed, 21 insertions(+), 14 deletions(-) (limited to 'examples') diff --git a/examples/da/plot_otda_classes.py b/examples/da/plot_otda_classes.py index e5c82fb..6870fa4 100644 --- a/examples/da/plot_otda_classes.py +++ b/examples/da/plot_otda_classes.py @@ -15,8 +15,10 @@ approaches currently supported in POT. # License: MIT License import matplotlib.pylab as pl +import numpy as np import ot +np.random.seed(42) # number of source and target points to generate ns = 150 diff --git a/examples/da/plot_otda_color_images.py b/examples/da/plot_otda_color_images.py index bca7350..805d0b0 100644 --- a/examples/da/plot_otda_color_images.py +++ b/examples/da/plot_otda_color_images.py @@ -20,9 +20,10 @@ SIAM Journal on Imaging Sciences, 7(3), 1853-1882. import numpy as np from scipy import ndimage import matplotlib.pylab as pl - import ot +np.random.seed(42) + def im2mat(I): """Converts and image to matrix (one pixel per line)""" diff --git a/examples/da/plot_otda_d2.py b/examples/da/plot_otda_d2.py index 1d2192f..8833eb2 100644 --- a/examples/da/plot_otda_d2.py +++ b/examples/da/plot_otda_d2.py @@ -19,17 +19,19 @@ of what the transport methods are doing. # License: MIT License import matplotlib.pylab as pl +import numpy as np import ot -# number of source and target points to generate -ns = 150 -nt = 150 +np.random.seed(42) -Xs, ys = ot.datasets.get_data_classif('3gauss', ns) -Xt, yt = ot.datasets.get_data_classif('3gauss2', nt) +n_samples_source = 150 +n_samples_target = 150 + +Xs, ys = ot.datasets.get_data_classif('3gauss', n_samples_source) +Xt, yt = ot.datasets.get_data_classif('3gauss2', n_samples_target) # Cost matrix -M = ot.dist(Xs, Xt) +M = ot.dist(Xs, Xt, metric='sqeuclidean') # Instantiate the different transport algorithms and fit them diff --git a/examples/da/plot_otda_mapping.py b/examples/da/plot_otda_mapping.py index 6d83507..aea7f09 100644 --- a/examples/da/plot_otda_mapping.py +++ b/examples/da/plot_otda_mapping.py @@ -23,7 +23,7 @@ import matplotlib.pylab as pl import ot -np.random.seed(0) +np.random.seed(42) ############################################################################## # generate @@ -31,10 +31,11 @@ np.random.seed(0) n = 100 # nb samples in source and target datasets theta = 2 * np.pi / 20 -nz = 0.1 -Xs, ys = ot.datasets.get_data_classif('gaussrot', n, nz=nz) -Xs_new, _ = ot.datasets.get_data_classif('gaussrot', n, nz=nz) -Xt, yt = ot.datasets.get_data_classif('gaussrot', n, theta=theta, nz=nz) +noise_level = 0.1 +Xs, ys = ot.datasets.get_data_classif('gaussrot', n, nz=noise_level) +Xs_new, _ = ot.datasets.get_data_classif('gaussrot', n, nz=noise_level) +Xt, yt = ot.datasets.get_data_classif( + 'gaussrot', n, theta=theta, nz=noise_level) # one of the target mode changes its variance (no linear mapping) Xt[yt == 2] *= 3 @@ -46,8 +47,7 @@ ot_mapping_linear = ot.da.MappingTransport( kernel="linear", mu=1e0, eta=1e-8, bias=True, max_iter=20, verbose=True) -ot_mapping_linear.fit( - Xs=Xs, Xt=Xt) +ot_mapping_linear.fit(Xs=Xs, Xt=Xt) # for original source samples, transform applies barycentric mapping transp_Xs_linear = ot_mapping_linear.transform(Xs=Xs) diff --git a/examples/da/plot_otda_mapping_colors_images.py b/examples/da/plot_otda_mapping_colors_images.py index 4209020..6c024ea 100644 --- a/examples/da/plot_otda_mapping_colors_images.py +++ b/examples/da/plot_otda_mapping_colors_images.py @@ -23,6 +23,8 @@ from scipy import ndimage import matplotlib.pylab as pl import ot +np.random.seed(42) + def im2mat(I): """Converts and image to matrix (one pixel per line)""" -- cgit v1.2.3 From 65de6fc9add57b95b8968e1e75fe1af342f81d01 Mon Sep 17 00:00:00 2001 From: Slasnista Date: Tue, 29 Aug 2017 13:34:34 +0200 Subject: pass on examples | introduced RandomState --- examples/da/plot_otda_classes.py | 20 +++++++++++++------- examples/da/plot_otda_color_images.py | 19 ++++++++++++++++--- examples/da/plot_otda_d2.py | 12 ++++++++++-- examples/da/plot_otda_mapping.py | 21 ++++++++++++++------- examples/da/plot_otda_mapping_colors_images.py | 9 ++++++--- 5 files changed, 59 insertions(+), 22 deletions(-) (limited to 'examples') diff --git a/examples/da/plot_otda_classes.py b/examples/da/plot_otda_classes.py index 6870fa4..ec57a37 100644 --- a/examples/da/plot_otda_classes.py +++ b/examples/da/plot_otda_classes.py @@ -15,19 +15,23 @@ approaches currently supported in POT. # License: MIT License import matplotlib.pylab as pl -import numpy as np import ot -np.random.seed(42) -# number of source and target points to generate -ns = 150 -nt = 150 +############################################################################## +# generate data +############################################################################## + +n_source_samples = 150 +n_target_samples = 150 + +Xs, ys = ot.datasets.get_data_classif('3gauss', n_source_samples) +Xt, yt = ot.datasets.get_data_classif('3gauss2', n_target_samples) -Xs, ys = ot.datasets.get_data_classif('3gauss', ns) -Xt, yt = ot.datasets.get_data_classif('3gauss2', nt) +############################################################################## # Instantiate the different transport algorithms and fit them +############################################################################## # EMD Transport ot_emd = ot.da.EMDTransport() @@ -52,6 +56,7 @@ transp_Xs_sinkhorn = ot_sinkhorn.transform(Xs=Xs) transp_Xs_lpl1 = ot_lpl1.transform(Xs=Xs) transp_Xs_l1l2 = ot_l1l2.transform(Xs=Xs) + ############################################################################## # Fig 1 : plots source and target samples ############################################################################## @@ -72,6 +77,7 @@ pl.legend(loc=0) pl.title('Target samples') pl.tight_layout() + ############################################################################## # Fig 2 : plot optimal couplings and transported samples ############################################################################## diff --git a/examples/da/plot_otda_color_images.py b/examples/da/plot_otda_color_images.py index 805d0b0..3984afb 100644 --- a/examples/da/plot_otda_color_images.py +++ b/examples/da/plot_otda_color_images.py @@ -22,7 +22,8 @@ from scipy import ndimage import matplotlib.pylab as pl import ot -np.random.seed(42) + +r = np.random.RandomState(42) def im2mat(I): @@ -39,6 +40,10 @@ def minmax(I): return np.clip(I, 0, 1) +############################################################################## +# generate data +############################################################################## + # Loading images I1 = ndimage.imread('../../data/ocean_day.jpg').astype(np.float64) / 256 I2 = ndimage.imread('../../data/ocean_sunset.jpg').astype(np.float64) / 256 @@ -48,12 +53,17 @@ X2 = im2mat(I2) # training samples nb = 1000 -idx1 = np.random.randint(X1.shape[0], size=(nb,)) -idx2 = np.random.randint(X2.shape[0], size=(nb,)) +idx1 = r.randint(X1.shape[0], size=(nb,)) +idx2 = r.randint(X2.shape[0], size=(nb,)) Xs = X1[idx1, :] Xt = X2[idx2, :] + +############################################################################## +# Instantiate the different transport algorithms and fit them +############################################################################## + # EMDTransport ot_emd = ot.da.EMDTransport() ot_emd.fit(Xs=Xs, Xt=Xt) @@ -75,6 +85,7 @@ I2t = minmax(mat2im(transp_Xt_emd, I2.shape)) I1te = minmax(mat2im(transp_Xs_sinkhorn, I1.shape)) I2te = minmax(mat2im(transp_Xt_sinkhorn, I2.shape)) + ############################################################################## # plot original image ############################################################################## @@ -91,6 +102,7 @@ pl.imshow(I2) pl.axis('off') pl.title('Image 2') + ############################################################################## # scatter plot of colors ############################################################################## @@ -112,6 +124,7 @@ pl.ylabel('Blue') pl.title('Image 2') pl.tight_layout() + ############################################################################## # plot new images ############################################################################## diff --git a/examples/da/plot_otda_d2.py b/examples/da/plot_otda_d2.py index 8833eb2..3daa0a6 100644 --- a/examples/da/plot_otda_d2.py +++ b/examples/da/plot_otda_d2.py @@ -19,10 +19,12 @@ of what the transport methods are doing. # License: MIT License import matplotlib.pylab as pl -import numpy as np import ot -np.random.seed(42) + +############################################################################## +# generate data +############################################################################## n_samples_source = 150 n_samples_target = 150 @@ -33,7 +35,10 @@ Xt, yt = ot.datasets.get_data_classif('3gauss2', n_samples_target) # Cost matrix M = ot.dist(Xs, Xt, metric='sqeuclidean') + +############################################################################## # Instantiate the different transport algorithms and fit them +############################################################################## # EMD Transport ot_emd = ot.da.EMDTransport() @@ -52,6 +57,7 @@ transp_Xs_emd = ot_emd.transform(Xs=Xs) transp_Xs_sinkhorn = ot_sinkhorn.transform(Xs=Xs) transp_Xs_lpl1 = ot_lpl1.transform(Xs=Xs) + ############################################################################## # Fig 1 : plots source and target samples + matrix of pairwise distance ############################################################################## @@ -78,6 +84,7 @@ pl.yticks([]) pl.title('Matrix of pairwise distances') pl.tight_layout() + ############################################################################## # Fig 2 : plots optimal couplings for the different methods ############################################################################## @@ -127,6 +134,7 @@ pl.yticks([]) pl.title('Main coupling coefficients\nSinkhornLpl1Transport') pl.tight_layout() + ############################################################################## # Fig 3 : plot transported samples ############################################################################## diff --git a/examples/da/plot_otda_mapping.py b/examples/da/plot_otda_mapping.py index aea7f09..09d2cb4 100644 --- a/examples/da/plot_otda_mapping.py +++ b/examples/da/plot_otda_mapping.py @@ -23,25 +23,31 @@ import matplotlib.pylab as pl import ot -np.random.seed(42) - ############################################################################## -# generate +# generate data ############################################################################## -n = 100 # nb samples in source and target datasets +n_source_samples = 100 +n_target_samples = 100 theta = 2 * np.pi / 20 noise_level = 0.1 -Xs, ys = ot.datasets.get_data_classif('gaussrot', n, nz=noise_level) -Xs_new, _ = ot.datasets.get_data_classif('gaussrot', n, nz=noise_level) + +Xs, ys = ot.datasets.get_data_classif( + 'gaussrot', n_source_samples, nz=noise_level) +Xs_new, _ = ot.datasets.get_data_classif( + 'gaussrot', n_source_samples, nz=noise_level) Xt, yt = ot.datasets.get_data_classif( - 'gaussrot', n, theta=theta, nz=noise_level) + 'gaussrot', n_target_samples, theta=theta, nz=noise_level) # one of the target mode changes its variance (no linear mapping) Xt[yt == 2] *= 3 Xt = Xt + 4 +############################################################################## +# Instantiate the different transport algorithms and fit them +############################################################################## + # MappingTransport with linear kernel ot_mapping_linear = ot.da.MappingTransport( kernel="linear", mu=1e0, eta=1e-8, bias=True, @@ -80,6 +86,7 @@ pl.scatter(Xt[:, 0], Xt[:, 1], c=yt, marker='o', label='Target samples') pl.legend(loc=0) pl.title('Source and target distributions') + ############################################################################## # plot transported samples ############################################################################## diff --git a/examples/da/plot_otda_mapping_colors_images.py b/examples/da/plot_otda_mapping_colors_images.py index 6c024ea..a628b05 100644 --- a/examples/da/plot_otda_mapping_colors_images.py +++ b/examples/da/plot_otda_mapping_colors_images.py @@ -23,7 +23,7 @@ from scipy import ndimage import matplotlib.pylab as pl import ot -np.random.seed(42) +r = np.random.RandomState(42) def im2mat(I): @@ -54,8 +54,8 @@ X2 = im2mat(I2) # training samples nb = 1000 -idx1 = np.random.randint(X1.shape[0], size=(nb,)) -idx2 = np.random.randint(X2.shape[0], size=(nb,)) +idx1 = r.randint(X1.shape[0], size=(nb,)) +idx2 = r.randint(X2.shape[0], size=(nb,)) Xs = X1[idx1, :] Xt = X2[idx2, :] @@ -91,6 +91,7 @@ ot_mapping_gaussian.fit(Xs=Xs, Xt=Xt) X1tn = ot_mapping_gaussian.transform(Xs=X1) # use the estimated mapping Image_mapping_gaussian = minmax(mat2im(X1tn, I1.shape)) + ############################################################################## # plot original images ############################################################################## @@ -107,6 +108,7 @@ pl.axis('off') pl.title('Image 2') pl.tight_layout() + ############################################################################## # plot pixel values distribution ############################################################################## @@ -128,6 +130,7 @@ pl.ylabel('Blue') pl.title('Image 2') pl.tight_layout() + ############################################################################## # plot transformed images ############################################################################## -- cgit v1.2.3 From 3007f1da1094f93fa4216386666085cf60316b04 Mon Sep 17 00:00:00 2001 From: Nicolas Courty Date: Thu, 31 Aug 2017 16:44:18 +0200 Subject: Minor corrections suggested by @agramfort + new barycenter example + test function --- README.md | 2 +- data/carre.png | Bin 0 -> 168 bytes data/coeur.png | Bin 0 -> 225 bytes data/rond.png | Bin 0 -> 230 bytes data/triangle.png | Bin 0 -> 254 bytes examples/plot_gromov.py | 14 +-- examples/plot_gromov_barycenter.py | 240 +++++++++++++++++++++++++++++++++++++ ot/gromov.py | 36 +++--- test/test_gromov.py | 38 ++++++ 9 files changed, 302 insertions(+), 28 deletions(-) create mode 100755 data/carre.png create mode 100755 data/coeur.png create mode 100755 data/rond.png create mode 100755 data/triangle.png create mode 100755 examples/plot_gromov_barycenter.py create mode 100644 test/test_gromov.py (limited to 'examples') diff --git a/README.md b/README.md index a42f17b..53672ae 100644 --- a/README.md +++ b/README.md @@ -183,4 +183,4 @@ You can also post bug reports and feature requests in Github issues. Make sure t [11] Flamary, R., Cuturi, M., Courty, N., & Rakotomamonjy, A. (2016). [Wasserstein Discriminant Analysis](https://arxiv.org/pdf/1608.08063.pdf). arXiv preprint arXiv:1608.08063. -[12] Peyré, Gabriel, Marco Cuturi, and Justin Solomon, [Gromov-Wasserstein averaging of kernel and distance matrices](http://proceedings.mlr.press/v48/peyre16.html) International Conference on Machine Learning (ICML). 2016. +[12] Gabriel Peyré, Marco Cuturi, and Justin Solomon, [Gromov-Wasserstein averaging of kernel and distance matrices](http://proceedings.mlr.press/v48/peyre16.html) International Conference on Machine Learning (ICML). 2016. diff --git a/data/carre.png b/data/carre.png new file mode 100755 index 0000000..45ff0ef Binary files /dev/null and b/data/carre.png differ diff --git a/data/coeur.png b/data/coeur.png new file mode 100755 index 0000000..3f511a6 Binary files /dev/null and b/data/coeur.png differ diff --git a/data/rond.png b/data/rond.png new file mode 100755 index 0000000..1c1a068 Binary files /dev/null and b/data/rond.png differ diff --git a/data/triangle.png b/data/triangle.png new file mode 100755 index 0000000..ca36d09 Binary files /dev/null and b/data/triangle.png differ diff --git a/examples/plot_gromov.py b/examples/plot_gromov.py index a33fde1..9bbdbde 100644 --- a/examples/plot_gromov.py +++ b/examples/plot_gromov.py @@ -1,8 +1,8 @@ # -*- coding: utf-8 -*- """ -==================== +========================== Gromov-Wasserstein example -==================== +========================== This example is designed to show how to use the Gromov-Wassertsein distance computation in POT. """ @@ -14,14 +14,14 @@ computation in POT. import scipy as sp import numpy as np +import matplotlib.pylab as pl import ot -import matplotlib.pylab as pl """ Sample two Gaussian distributions (2D and 3D) -==================== +============================================= The Gromov-Wasserstein distance allows to compute distances with samples that do not belong to the same metric space. For demonstration purpose, we sample two Gaussian distributions in 2- and 3-dimensional spaces. """ @@ -42,7 +42,7 @@ xt = np.random.randn(n, 3).dot(P) + mu_t """ Plotting the distributions -==================== +========================== """ fig = pl.figure() ax1 = fig.add_subplot(121) @@ -54,7 +54,7 @@ pl.show() """ Compute distance kernels, normalize them and then display -==================== +========================================================= """ C1 = sp.spatial.distance.cdist(xs, xs) @@ -72,7 +72,7 @@ pl.show() """ Compute Gromov-Wasserstein plans and distance -==================== +============================================= """ p = ot.unif(n) diff --git a/examples/plot_gromov_barycenter.py b/examples/plot_gromov_barycenter.py new file mode 100755 index 0000000..6a72b3b --- /dev/null +++ b/examples/plot_gromov_barycenter.py @@ -0,0 +1,240 @@ +# -*- coding: utf-8 -*- +""" +===================================== +Gromov-Wasserstein Barycenter example +===================================== +This example is designed to show how to use the Gromov-Wassertsein distance +computation in POT. +""" + +# Author: Erwan Vautier +# Nicolas Courty +# +# License: MIT License + + +import numpy as np +import scipy as sp + +import scipy.ndimage as spi +import matplotlib.pylab as pl +from sklearn import manifold +from sklearn.decomposition import PCA + +import ot + +""" + +Smacof MDS +========== +This function allows to find an embedding of points given a dissimilarity matrix +that will be given by the output of the algorithm +""" + + +def smacof_mds(C, dim, maxIter=3000, eps=1e-9): + """ + Returns an interpolated point cloud following the dissimilarity matrix C using SMACOF + multidimensional scaling (MDS) in specific dimensionned target space + + Parameters + ---------- + C : np.ndarray(ns,ns) + dissimilarity matrix + dim : Integer + dimension of the targeted space + maxIter : Maximum number of iterations of the SMACOF algorithm for a single run + + eps : relative tolerance w.r.t stress to declare converge + + + Returns + ------- + npos : R**dim ndarray + Embedded coordinates of the interpolated point cloud (defined with one isometry) + + + """ + + seed = np.random.RandomState(seed=3) + + mds = manifold.MDS( + dim, + max_iter=3000, + eps=1e-9, + dissimilarity='precomputed', + n_init=1) + pos = mds.fit(C).embedding_ + + nmds = manifold.MDS( + 2, + max_iter=3000, + eps=1e-9, + dissimilarity="precomputed", + random_state=seed, + n_init=1) + npos = nmds.fit_transform(C, init=pos) + + return npos + + +""" +Data preparation +================ +The four distributions are constructed from 4 simple images +""" + + +def im2mat(I): + """Converts and image to matrix (one pixel per line)""" + return I.reshape((I.shape[0] * I.shape[1], I.shape[2])) + + +carre = spi.imread('../data/carre.png').astype(np.float64) / 256 +rond = spi.imread('../data/rond.png').astype(np.float64) / 256 +triangle = spi.imread('../data/triangle.png').astype(np.float64) / 256 +fleche = spi.imread('../data/coeur.png').astype(np.float64) / 256 + +shapes = [carre, rond, triangle, fleche] + +S = 4 +xs = [[] for i in range(S)] + + +for nb in range(4): + for i in range(8): + for j in range(8): + if shapes[nb][i, j] < 0.95: + xs[nb].append([j, 8 - i]) + +xs = np.array([np.array(xs[0]), np.array(xs[1]), + np.array(xs[2]), np.array(xs[3])]) + + +""" +Barycenter computation +====================== +The four distributions are constructed from 4 simple images +""" +ns = [len(xs[s]) for s in range(S)] +N = 30 + +"""Compute all distances matrices for the four shapes""" +Cs = [sp.spatial.distance.cdist(xs[s], xs[s]) for s in range(S)] +Cs = [cs / cs.max() for cs in Cs] + +ps = [ot.unif(ns[s]) for s in range(S)] +p = ot.unif(N) + + +lambdast = [[float(i) / 3, float(3 - i) / 3] for i in [1, 2]] + +Ct01 = [0 for i in range(2)] +for i in range(2): + Ct01[i] = ot.gromov.gromov_barycenters(N, [Cs[0], Cs[1]], [ + ps[0], ps[1]], p, lambdast[i], 'square_loss', 5e-4, numItermax=100, stopThr=1e-3) + +Ct02 = [0 for i in range(2)] +for i in range(2): + Ct02[i] = ot.gromov.gromov_barycenters(N, [Cs[0], Cs[2]], [ + ps[0], ps[2]], p, lambdast[i], 'square_loss', 5e-4, numItermax=100, stopThr=1e-3) + +Ct13 = [0 for i in range(2)] +for i in range(2): + Ct13[i] = ot.gromov.gromov_barycenters(N, [Cs[1], Cs[3]], [ + ps[1], ps[3]], p, lambdast[i], 'square_loss', 5e-4, numItermax=100, stopThr=1e-3) + +Ct23 = [0 for i in range(2)] +for i in range(2): + Ct23[i] = ot.gromov.gromov_barycenters(N, [Cs[2], Cs[3]], [ + ps[2], ps[3]], p, lambdast[i], 'square_loss', 5e-4, numItermax=100, stopThr=1e-3) + +""" +Visualization +============= +""" + +"""The PCA helps in getting consistency between the rotations""" + +clf = PCA(n_components=2) +npos = [0, 0, 0, 0] +npos = [smacof_mds(Cs[s], 2) for s in range(S)] + +npost01 = [0, 0] +npost01 = [smacof_mds(Ct01[s], 2) for s in range(2)] +npost01 = [clf.fit_transform(npost01[s]) for s in range(2)] + +npost02 = [0, 0] +npost02 = [smacof_mds(Ct02[s], 2) for s in range(2)] +npost02 = [clf.fit_transform(npost02[s]) for s in range(2)] + +npost13 = [0, 0] +npost13 = [smacof_mds(Ct13[s], 2) for s in range(2)] +npost13 = [clf.fit_transform(npost13[s]) for s in range(2)] + +npost23 = [0, 0] +npost23 = [smacof_mds(Ct23[s], 2) for s in range(2)] +npost23 = [clf.fit_transform(npost23[s]) for s in range(2)] + + +fig = pl.figure(figsize=(10, 10)) + +ax1 = pl.subplot2grid((4, 4), (0, 0)) +pl.xlim((-1, 1)) +pl.ylim((-1, 1)) +ax1.scatter(npos[0][:, 0], npos[0][:, 1], color='r') + +ax2 = pl.subplot2grid((4, 4), (0, 1)) +pl.xlim((-1, 1)) +pl.ylim((-1, 1)) +ax2.scatter(npost01[1][:, 0], npost01[1][:, 1], color='b') + +ax3 = pl.subplot2grid((4, 4), (0, 2)) +pl.xlim((-1, 1)) +pl.ylim((-1, 1)) +ax3.scatter(npost01[0][:, 0], npost01[0][:, 1], color='b') + +ax4 = pl.subplot2grid((4, 4), (0, 3)) +pl.xlim((-1, 1)) +pl.ylim((-1, 1)) +ax4.scatter(npos[1][:, 0], npos[1][:, 1], color='r') + +ax5 = pl.subplot2grid((4, 4), (1, 0)) +pl.xlim((-1, 1)) +pl.ylim((-1, 1)) +ax5.scatter(npost02[1][:, 0], npost02[1][:, 1], color='b') + +ax6 = pl.subplot2grid((4, 4), (1, 3)) +pl.xlim((-1, 1)) +pl.ylim((-1, 1)) +ax6.scatter(npost13[1][:, 0], npost13[1][:, 1], color='b') + +ax7 = pl.subplot2grid((4, 4), (2, 0)) +pl.xlim((-1, 1)) +pl.ylim((-1, 1)) +ax7.scatter(npost02[0][:, 0], npost02[0][:, 1], color='b') + +ax8 = pl.subplot2grid((4, 4), (2, 3)) +pl.xlim((-1, 1)) +pl.ylim((-1, 1)) +ax8.scatter(npost13[0][:, 0], npost13[0][:, 1], color='b') + +ax9 = pl.subplot2grid((4, 4), (3, 0)) +pl.xlim((-1, 1)) +pl.ylim((-1, 1)) +ax9.scatter(npos[2][:, 0], npos[2][:, 1], color='r') + +ax10 = pl.subplot2grid((4, 4), (3, 1)) +pl.xlim((-1, 1)) +pl.ylim((-1, 1)) +ax10.scatter(npost23[1][:, 0], npost23[1][:, 1], color='b') + +ax11 = pl.subplot2grid((4, 4), (3, 2)) +pl.xlim((-1, 1)) +pl.ylim((-1, 1)) +ax11.scatter(npost23[0][:, 0], npost23[0][:, 1], color='b') + +ax12 = pl.subplot2grid((4, 4), (3, 3)) +pl.xlim((-1, 1)) +pl.ylim((-1, 1)) +ax12.scatter(npos[3][:, 0], npos[3][:, 1], color='r') diff --git a/ot/gromov.py b/ot/gromov.py index 7cf3b42..421ed3f 100644 --- a/ot/gromov.py +++ b/ot/gromov.py @@ -23,7 +23,7 @@ def square_loss(a, b): Returns the value of L(a,b)=(1/2)*|a-b|^2 """ - return (1 / 2) * (a - b)**2 + return 0.5 * (a - b)**2 def kl_loss(a, b): @@ -54,9 +54,9 @@ def tensor_square_loss(C1, C2, T): Parameters ---------- - C1 : np.ndarray(ns,ns) + C1 : ndarray, shape (ns, ns) Metric cost matrix in the source space - C2 : np.ndarray(nt,nt) + C2 : ndarray, shape (nt, nt) Metric costfr matrix in the target space T : np.ndarray(ns,nt) Coupling between source and target spaces @@ -87,7 +87,7 @@ def tensor_square_loss(C1, C2, T): return b tens = -np.dot(h1(C1), T).dot(h2(C2).T) - tens = tens - tens.min() + tens -= tens.min() return np.array(tens) @@ -112,9 +112,9 @@ def tensor_kl_loss(C1, C2, T): Parameters ---------- - C1 : np.ndarray(ns,ns) + C1 : ndarray, shape (ns, ns) Metric cost matrix in the source space - C2 : np.ndarray(nt,nt) + C2 : ndarray, shape (nt, nt) Metric costfr matrix in the target space T : np.ndarray(ns,nt) Coupling between source and target spaces @@ -149,7 +149,7 @@ def tensor_kl_loss(C1, C2, T): return np.log(b + 1e-15) tens = -np.dot(h1(C1), T).dot(h2(C2).T) - tens = tens - tens.min() + tens -= tens.min() return np.array(tens) @@ -175,9 +175,8 @@ def update_square_loss(p, lambdas, T, Cs): """ - tmpsum = np.sum([lambdas[s] * np.dot(T[s].T, Cs[s]).dot(T[s]) - for s in range(len(T))]) - ppt = np.dot(p, p.T) + tmpsum = sum([lambdas[s] * np.dot(T[s].T, Cs[s]).dot(T[s]) for s in range(len(T))]) + ppt = np.outer(p, p) return(np.divide(tmpsum, ppt)) @@ -203,9 +202,8 @@ def update_kl_loss(p, lambdas, T, Cs): """ - tmpsum = np.sum([lambdas[s] * np.dot(T[s].T, Cs[s]).dot(T[s]) - for s in range(len(T))]) - ppt = np.dot(p, p.T) + tmpsum = sum([lambdas[s] * np.dot(T[s].T, Cs[s]).dot(T[s]) for s in range(len(T))]) + ppt = np.outer(p, p) return(np.exp(np.divide(tmpsum, ppt))) @@ -239,9 +237,9 @@ def gromov_wasserstein(C1, C2, p, q, loss_fun, epsilon, numItermax=1000, stopThr Parameters ---------- - C1 : np.ndarray(ns,ns) + C1 : ndarray, shape (ns, ns) Metric cost matrix in the source space - C2 : np.ndarray(nt,nt) + C2 : ndarray, shape (nt, nt) Metric costfr matrix in the target space p : np.ndarray(ns,) distribution in the source space @@ -271,7 +269,7 @@ def gromov_wasserstein(C1, C2, p, q, loss_fun, epsilon, numItermax=1000, stopThr C1 = np.asarray(C1, dtype=np.float64) C2 = np.asarray(C2, dtype=np.float64) - T = np.dot(p, q.T) # Initialization + T = np.outer(p, q) # Initialization cpt = 0 err = 1 @@ -333,9 +331,9 @@ def gromov_wasserstein2(C1, C2, p, q, loss_fun, epsilon, numItermax=1000, stopTh Parameters ---------- - C1 : np.ndarray(ns,ns) + C1 : ndarray, shape (ns, ns) Metric cost matrix in the source space - C2 : np.ndarray(nt,nt) + C2 : ndarray, shape (nt, nt) Metric costfr matrix in the target space p : np.ndarray(ns,) distribution in the source space @@ -434,8 +432,6 @@ def gromov_barycenters(N, Cs, ps, p, lambdas, loss_fun, epsilon, numItermax=1000 Cs = [np.asarray(Cs[s], dtype=np.float64) for s in range(S)] lambdas = np.asarray(lambdas, dtype=np.float64) - T = [0 for s in range(S)] - # Initialization of C : random SPD matrix xalea = np.random.randn(N, 2) C = dist(xalea, xalea) diff --git a/test/test_gromov.py b/test/test_gromov.py new file mode 100644 index 0000000..75eeaab --- /dev/null +++ b/test/test_gromov.py @@ -0,0 +1,38 @@ +"""Tests for module gromov """ + +# Author: Erwan Vautier +# Nicolas Courty +# +# License: MIT License + +import numpy as np +import ot + + +def test_gromov(): + n = 50 # nb samples + + mu_s = np.array([0, 0]) + cov_s = np.array([[1, 0], [0, 1]]) + + xs = ot.datasets.get_2D_samples_gauss(n, mu_s, cov_s) + + xt = [xs[n - (i + 1)] for i in range(n)] + xt = np.array(xt) + + p = ot.unif(n) + q = ot.unif(n) + + C1 = ot.dist(xs, xs) + C2 = ot.dist(xt, xt) + + C1 /= C1.max() + C2 /= C2.max() + + G = ot.gromov_wasserstein(C1, C2, p, q, 'square_loss', epsilon=5e-4) + + # check constratints + np.testing.assert_allclose( + p, G.sum(1), atol=1e-04) # cf convergence gromov + np.testing.assert_allclose( + q, G.sum(0), atol=1e-04) # cf convergence gromov -- cgit v1.2.3 From bc68cc3e8b23ad7d542518ba8ffa665094d57663 Mon Sep 17 00:00:00 2001 From: Nicolas Courty Date: Thu, 31 Aug 2017 17:17:30 +0200 Subject: minor corrections --- examples/plot_gromov_barycenter.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) (limited to 'examples') diff --git a/examples/plot_gromov_barycenter.py b/examples/plot_gromov_barycenter.py index 6a72b3b..da52768 100755 --- a/examples/plot_gromov_barycenter.py +++ b/examples/plot_gromov_barycenter.py @@ -32,18 +32,19 @@ that will be given by the output of the algorithm """ -def smacof_mds(C, dim, maxIter=3000, eps=1e-9): +def smacof_mds(C, dim, max_iter=3000, eps=1e-9): """ Returns an interpolated point cloud following the dissimilarity matrix C using SMACOF multidimensional scaling (MDS) in specific dimensionned target space Parameters ---------- - C : np.ndarray(ns,ns) + C : ndarray, shape (ns, ns) dissimilarity matrix - dim : Integer + dim : int dimension of the targeted space - maxIter : Maximum number of iterations of the SMACOF algorithm for a single run + max_iter : int + Maximum number of iterations of the SMACOF algorithm for a single run eps : relative tolerance w.r.t stress to declare converge @@ -60,7 +61,7 @@ def smacof_mds(C, dim, maxIter=3000, eps=1e-9): mds = manifold.MDS( dim, - max_iter=3000, + max_iter=max_iter, eps=1e-9, dissimilarity='precomputed', n_init=1) @@ -68,7 +69,7 @@ def smacof_mds(C, dim, maxIter=3000, eps=1e-9): nmds = manifold.MDS( 2, - max_iter=3000, + max_iter=max_iter, eps=1e-9, dissimilarity="precomputed", random_state=seed, -- cgit v1.2.3