diff options
| author | Antoine Rolet <antoine.rolet@gmail.com> | 2017-09-05 15:30:50 +0900 |
|---|---|---|
| committer | Antoine Rolet <antoine.rolet@gmail.com> | 2017-09-05 15:30:50 +0900 |
| commit | 13dfb3ddbbd8926b4751b82dd41c5570253b1f07 (patch) | |
| tree | b28098e98640c64483a599103e2fdb5df46d2c79 | |
| parent | 185eb3e2ef34b5ce6b8f90a28a5bcc78432b7fd3 (diff) | |
| parent | 16697047eff9326a0ecb483317c13a854a3d3a71 (diff) | |
Merge remote-tracking branch 'upstream/master'
54 files changed, 4781 insertions, 1984 deletions
@@ -97,3 +97,9 @@ ENV/ # Rope project settings .ropeproject + +# Mac stuff +.DS_Store + +# coverage output folder +cov_html/ diff --git a/.travis.yml b/.travis.yml index 6126b3a..ec2b3d2 100644 --- a/.travis.yml +++ b/.travis.yml @@ -13,9 +13,18 @@ matrix: python: 2.7 before_install: - ./.travis/before_install.sh +before_script: # configure a headless display to test plot generation + - "export DISPLAY=:99.0" + - "sh -e /etc/init.d/xvfb start" + - sleep 3 # give xvfb some time to start # command to install dependencies install: - pip install -r requirements.txt - - python setup.py install -# command to run tests -script: python test/test_load_module.py -v + - pip install flake8 pytest + - pip install . +# command to run tests + check syntax style +script: + - python setup.py develop + - flake8 examples/ ot/ test/ + - python -m py.test -v test/ + # - py.test ot test diff --git a/CODE_OF_CONDUCT.md b/CODE_OF_CONDUCT.md new file mode 100644 index 0000000..9c1c621 --- /dev/null +++ b/CODE_OF_CONDUCT.md @@ -0,0 +1,74 @@ +# Contributor Covenant Code of Conduct + +## Our Pledge + +In the interest of fostering an open and welcoming environment, we as +contributors and maintainers pledge to making participation in our project and +our community a harassment-free experience for everyone, regardless of age, body +size, disability, ethnicity, gender identity and expression, level of experience, +nationality, personal appearance, race, religion, or sexual identity and +orientation. + +## Our Standards + +Examples of behavior that contributes to creating a positive environment +include: + +* Using welcoming and inclusive language +* Being respectful of differing viewpoints and experiences +* Gracefully accepting constructive criticism +* Focusing on what is best for the community +* Showing empathy towards other community members + +Examples of unacceptable behavior by participants include: + +* The use of sexualized language or imagery and unwelcome sexual attention or +advances +* Trolling, insulting/derogatory comments, and personal or political attacks +* Public or private harassment +* Publishing others' private information, such as a physical or electronic + address, without explicit permission +* Other conduct which could reasonably be considered inappropriate in a + professional setting + +## Our Responsibilities + +Project maintainers are responsible for clarifying the standards of acceptable +behavior and are expected to take appropriate and fair corrective action in +response to any instances of unacceptable behavior. + +Project maintainers have the right and responsibility to remove, edit, or +reject comments, commits, code, wiki edits, issues, and other contributions +that are not aligned to this Code of Conduct, or to ban temporarily or +permanently any contributor for other behaviors that they deem inappropriate, +threatening, offensive, or harmful. + +## Scope + +This Code of Conduct applies both within project spaces and in public spaces +when an individual is representing the project or its community. Examples of +representing a project or community include using an official project e-mail +address, posting via an official social media account, or acting as an appointed +representative at an online or offline event. Representation of a project may be +further defined and clarified by project maintainers. + +## Enforcement + +Instances of abusive, harassing, or otherwise unacceptable behavior may be +reported by contacting the project team. All +complaints will be reviewed and investigated and will result in a response that +is deemed necessary and appropriate to the circumstances. The project team is +obligated to maintain confidentiality with regard to the reporter of an incident. +Further details of specific enforcement policies may be posted separately. + +Project maintainers who do not follow or enforce the Code of Conduct in good +faith may face temporary or permanent repercussions as determined by other +members of the project's leadership. + +## Attribution + +This Code of Conduct is adapted from the [Contributor Covenant][homepage], version 1.4, +available at [http://contributor-covenant.org/version/1/4][version] + +[homepage]: http://contributor-covenant.org +[version]: http://contributor-covenant.org/version/1/4/ diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md new file mode 100644 index 0000000..84ef29a --- /dev/null +++ b/CONTRIBUTING.md @@ -0,0 +1,204 @@ + +Contributing to POT +=================== + + +First off, thank you for considering contributing to POT. + +How to contribute +----------------- + +The preferred workflow for contributing to POT is to fork the +[main repository](https://github.com/rflamary/POT) on +GitHub, clone, and develop on a branch. Steps: + +1. Fork the [project repository](https://github.com/rflamary/POT) + by clicking on the 'Fork' button near the top right of the page. This creates + a copy of the code under your GitHub user account. For more details on + how to fork a repository see [this guide](https://help.github.com/articles/fork-a-repo/). + +2. Clone your fork of the scikit-learn repo from your GitHub account to your local disk: + + ```bash + $ git clone git@github.com:YourLogin/POT.git + $ cd POT + ``` + +3. Create a ``feature`` branch to hold your development changes: + + ```bash + $ git checkout -b my-feature + ``` + + Always use a ``feature`` branch. It's good practice to never work on the ``master`` branch! + +4. Develop the feature on your feature branch. Add changed files using ``git add`` and then ``git commit`` files: + + ```bash + $ git add modified_files + $ git commit + ``` + + to record your changes in Git, then push the changes to your GitHub account with: + + ```bash + $ git push -u origin my-feature + ``` + +5. Follow [these instructions](https://help.github.com/articles/creating-a-pull-request-from-a-fork) +to create a pull request from your fork. This will send an email to the committers. + +(If any of the above seems like magic to you, please look up the +[Git documentation](https://git-scm.com/documentation) on the web, or ask a friend or another contributor for help.) + +Pull Request Checklist +---------------------- + +We recommended that your contribution complies with the +following rules before you submit a pull request: + +- Follow the PEP8 Guidelines. + +- If your pull request addresses an issue, please use the pull request title + to describe the issue and mention the issue number in the pull request description. This will make sure a link back to the original issue is + created. + +- All public methods should have informative docstrings with sample + usage presented as doctests when appropriate. + +- Please prefix the title of your pull request with `[MRG]` (Ready for + Merge), if the contribution is complete and ready for a detailed review. + Two core developers will review your code and change the prefix of the pull + request to `[MRG + 1]` and `[MRG + 2]` on approval, making it eligible + for merging. An incomplete contribution -- where you expect to do more work before + receiving a full review -- should be prefixed `[WIP]` (to indicate a work + in progress) and changed to `[MRG]` when it matures. WIPs may be useful + to: indicate you are working on something to avoid duplicated work, + request broad review of functionality or API, or seek collaborators. + WIPs often benefit from the inclusion of a + [task list](https://github.com/blog/1375-task-lists-in-gfm-issues-pulls-comments) + in the PR description. + + +- When adding additional functionality, provide at least one + example script in the ``examples/`` folder. Have a look at other + examples for reference. Examples should demonstrate why the new + functionality is useful in practice and, if possible, compare it + to other methods available in scikit-learn. + +- Documentation and high-coverage tests are necessary for enhancements to be + accepted. Bug-fixes or new features should be provided with + [non-regression tests](https://en.wikipedia.org/wiki/Non-regression_testing). + These tests verify the correct behavior of the fix or feature. In this + manner, further modifications on the code base are granted to be consistent + with the desired behavior. + For the Bug-fixes case, at the time of the PR, this tests should fail for + the code base in master and pass for the PR code. + +- At least one paragraph of narrative documentation with links to + references in the literature (with PDF links when possible) and + the example. + +You can also check for common programming errors with the following +tools: + + +- No pyflakes warnings, check with: + + ```bash + $ pip install pyflakes + $ pyflakes path/to/module.py + ``` + +- No PEP8 warnings, check with: + + ```bash + $ pip install pep8 + $ pep8 path/to/module.py + ``` + +- AutoPEP8 can help you fix some of the easy redundant errors: + + ```bash + $ pip install autopep8 + $ autopep8 path/to/pep8.py + ``` + +Bonus points for contributions that include a performance analysis with +a benchmark script and profiling output (please report on the mailing +list or on the GitHub issue). + +Filing bugs +----------- +We use Github issues to track all bugs and feature requests; feel free to +open an issue if you have found a bug or wish to see a feature implemented. + +It is recommended to check that your issue complies with the +following rules before submitting: + +- Verify that your issue is not being currently addressed by other + [issues](https://github.com/rflamary/POT/issues?q=) + or [pull requests](https://github.com/rflamary/POT/pulls?q=). + +- Please ensure all code snippets and error messages are formatted in + appropriate code blocks. + See [Creating and highlighting code blocks](https://help.github.com/articles/creating-and-highlighting-code-blocks). + +- Please include your operating system type and version number, as well + as your Python, scikit-learn, numpy, and scipy versions. This information + can be found by running the following code snippet: + + ```python + import platform; print(platform.platform()) + import sys; print("Python", sys.version) + import numpy; print("NumPy", numpy.__version__) + import scipy; print("SciPy", scipy.__version__) + import ot; print("POT", ot.__version__) + ``` + +- Please be specific about what estimators and/or functions are involved + and the shape of the data, as appropriate; please include a + [reproducible](http://stackoverflow.com/help/mcve) code snippet + or link to a [gist](https://gist.github.com). If an exception is raised, + please provide the traceback. + +New contributor tips +-------------------- + +A great way to start contributing to scikit-learn is to pick an item +from the list of [Easy issues](https://github.com/scikit-learn/scikit-learn/issues?labels=Easy) +in the issue tracker. Resolving these issues allow you to start +contributing to the project without much prior knowledge. Your +assistance in this area will be greatly appreciated by the more +experienced developers as it helps free up their time to concentrate on +other issues. + +Documentation +------------- + +We are glad to accept any sort of documentation: function docstrings, +reStructuredText documents (like this one), tutorials, etc. +reStructuredText documents live in the source code repository under the +doc/ directory. + +You can edit the documentation using any text editor and then generate +the HTML output by typing ``make html`` from the doc/ directory. +Alternatively, ``make`` can be used to quickly generate the +documentation without the example gallery. The resulting HTML files will +be placed in ``_build/html/`` and are viewable in a web browser. See the +``README`` file in the ``doc/`` directory for more information. + +For building the documentation, you will need +[sphinx](http://sphinx.pocoo.org/), +[matplotlib](http://matplotlib.org/), and +[pillow](http://pillow.readthedocs.io/en/latest/). + +When you are writing documentation, it is important to keep a good +compromise between mathematical and algorithmic details, and give +intuition to the reader on what the algorithm does. It is best to always +start with a small paragraph with a hand-waving explanation of what the +method does to the data and a figure (coming from an example) +illustrating it. + + +This Contrubution guide is strongly inpired by the one of the [scikit-learn](https://github.com/scikit-learn/scikit-learn) team. diff --git a/MANIFEST.in b/MANIFEST.in index 2fd6a10..e0acb7a 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -1,5 +1,6 @@ graft ot/lp/ include README.md +include LICENSE include ot/lp/core.h include ot/lp/EMD.h include ot/lp/EMD_wrapper.cpp @@ -31,19 +31,28 @@ sremove : tr '\n' '\0' < files.txt | sudo xargs -0 rm -f -- rm files.txt -clean : +clean : FORCE $(PYTHON) setup.py clean + +pep8 : + flake8 examples/ ot/ test/ + +test : FORCE pep8 + python -m py.test -v test/ --cov=ot --cov-report html:cov_html -test: - pytest +pytest : FORCE + python -m py.test -v test/ --cov=ot -uploadpypi: +uploadpypi : #python setup.py register python setup.py sdist upload -r pypi -rdoc: +rdoc : pandoc --from=markdown --to=rst --output=docs/source/readme.rst README.md notebook : ipython notebook --matplotlib=inline --notebook-dir=notebooks/ + + +FORCE : @@ -22,34 +22,36 @@ Some demonstrations (both in Python and Jupyter Notebook format) are available i ## Installation -The Library has been tested on Linux and MacOSX. It requires a C++ compiler for using the EMD solver and rely on the following Python modules: +The library has been tested on Linux, MacOSX and Windows. It requires a C++ compiler for using the EMD solver and relies on the following Python modules: - Numpy (>=1.11) - Scipy (>=0.17) - Cython (>=0.23) - Matplotlib (>=1.5) +#### Pip installation -Under debian based linux the dependencies can be installed with +You can install the toolbox through PyPI with: ``` -sudo apt-get install python-numpy python-scipy python-matplotlib cython +pip install POT ``` - -To install the library, you can install it locally (after downloading it) on you machine using +or get the very latest version by downloading it and then running: ``` python setup.py install --user # for user install (no root) ``` -The toolbox is also available on PyPI with a possibly slightly older version. You can install it with: +#### Anaconda installation with conda-forge + +If you use the Anaconda python distribution, POT is available in [conda-forge](https://conda-forge.org). To install it and the required dependencies: ``` -pip install POT +conda install -c conda-forge pot ``` +#### Post installation check After a correct installation, you should be able to import the module without errors: ```python import ot ``` - Note that for easier access the module is name ot instead of pot. @@ -130,9 +132,12 @@ The contributors to this library are: * [Rémi Flamary](http://remi.flamary.com/) * [Nicolas Courty](http://people.irisa.fr/Nicolas.Courty/) +* [Alexandre Gramfort](http://alexandre.gramfort.net/) * [Laetitia Chapel](http://people.irisa.fr/Laetitia.Chapel/) * [Michael Perrot](http://perso.univ-st-etienne.fr/pem82055/) (Mapping estimation) * [Léo Gautheron](https://github.com/aje) (GPU implementation) +* [Nathalie Gayraud](https://www.linkedin.com/in/nathalie-t-h-gayraud/?ppe=1) +* [Stanislas Chambon](https://slasnista.github.io/) This toolbox benefit a lot from open source research and we would like to thank the following persons for providing some code (in various languages): @@ -141,6 +146,21 @@ This toolbox benefit a lot from open source research and we would like to thank * [Antoine Rolet](https://arolet.github.io/) ( Mex file for EMD ) * [Marco Cuturi](http://marcocuturi.net/) (Sinkhorn Knopp in Matlab/Cuda) + +## Contributions and code of conduct + +Every contribution is welcome and should respect the [contribution guidelines](CONTRIBUTING.md). Each member of the project is expected to follow the [code of conduct](CODE_OF_CONDUCT.md). + +## Support + +You can ask questions and join the development discussion: + +* On the [POT Slack channel](https://pot-toolbox.slack.com) +* On the POT [mailing list](https://mail.python.org/mm3/mailman3/lists/pot.python.org/) + + +You can also post bug reports and feature requests in Github issues. Make sure to read our [guidelines](CONTRIBUTING.md) first. + ## References [1] Bonneel, N., Van De Panne, M., Paris, S., & Heidrich, W. (2011, December). [Displacement interpolation using Lagrangian mass transport](https://people.csail.mit.edu/sparis/publi/2011/sigasia/Bonneel_11_Displacement_Interpolation.pdf). In ACM Transactions on Graphics (TOG) (Vol. 30, No. 6, p. 158). ACM. diff --git a/docs/source/readme.rst b/docs/source/readme.rst index 611001b..c1e0017 100644 --- a/docs/source/readme.rst +++ b/docs/source/readme.rst @@ -28,8 +28,8 @@ available in the examples folder. Installation ------------ -The Library has been tested on Linux and MacOSX. It requires a C++ -compiler for using the EMD solver and rely on the following Python +The library has been tested on Linux, MacOSX and Windows. It requires a +C++ compiler for using the EMD solver and relies on the following Python modules: - Numpy (>=1.11) @@ -37,25 +37,34 @@ modules: - Cython (>=0.23) - Matplotlib (>=1.5) -Under debian based linux the dependencies can be installed with +Pip installation +^^^^^^^^^^^^^^^^ + +You can install the toolbox through PyPI with: :: - sudo apt-get install python-numpy python-scipy python-matplotlib cython + pip install POT -To install the library, you can install it locally (after downloading -it) on you machine using +or get the very latest version by downloading it and then running: :: python setup.py install --user # for user install (no root) -The toolbox is also available on PyPI with a possibly slightly older -version. You can install it with: +Anaconda installation with conda-forge +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +If you use the Anaconda python distribution, POT is available in +`conda-forge <https://conda-forge.org>`__. To install it and the +required dependencies: :: - pip install POT + conda install -c conda-forge pot + +Post installation check +^^^^^^^^^^^^^^^^^^^^^^^ After a correct installation, you should be able to import the module without errors: @@ -109,6 +118,7 @@ Short examples # a,b are 1D histograms (sum to 1 and positive) # M is the ground cost matrix Wd=ot.emd2(a,b,M) # exact linear program + Wd_reg=ot.sinkhorn2(a,b,M,reg) # entropic regularized OT # if b is a matrix compute all distances to a and return a vector - Compute OT matrix @@ -117,8 +127,8 @@ Short examples # a,b are 1D histograms (sum to 1 and positive) # M is the ground cost matrix - Totp=ot.emd(a,b,M) # exact linear program - Totp_reg=ot.sinkhorn(a,b,M,reg) # entropic regularized OT + T=ot.emd(a,b,M) # exact linear program + T_reg=ot.sinkhorn(a,b,M,reg) # entropic regularized OT - Compute Wasserstein barycenter @@ -172,6 +182,7 @@ The contributors to this library are: - `Rémi Flamary <http://remi.flamary.com/>`__ - `Nicolas Courty <http://people.irisa.fr/Nicolas.Courty/>`__ +- `Alexandre Gramfort <http://alexandre.gramfort.net/>`__ - `Laetitia Chapel <http://people.irisa.fr/Laetitia.Chapel/>`__ - `Michael Perrot <http://perso.univ-st-etienne.fr/pem82055/>`__ (Mapping estimation) @@ -189,6 +200,25 @@ languages): - `Marco Cuturi <http://marcocuturi.net/>`__ (Sinkhorn Knopp in Matlab/Cuda) +Contributions and code of conduct +--------------------------------- + +Every contribution is welcome and should respect the `contribution +guidelines <CONTRIBUTING.md>`__. Each member of the project is expected +to follow the `code of conduct <CODE_OF_CONDUCT.md>`__. + +Support +------- + +You can ask questions and join the development discussion: + +- On the `POT Slack channel <https://pot-toolbox.slack.com>`__ +- On the POT `mailing + list <https://mail.python.org/mm3/mailman3/lists/pot.python.org/>`__ + +You can also post bug reports and feature requests in Github issues. +Make sure to read our `guidelines <CONTRIBUTING.md>`__ first. + References ---------- diff --git a/examples/da/plot_otda_classes.py b/examples/da/plot_otda_classes.py new file mode 100644 index 0000000..ec57a37 --- /dev/null +++ b/examples/da/plot_otda_classes.py @@ -0,0 +1,150 @@ +# -*- 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 <remi.flamary@unice.fr> +# Stanislas Chambon <stan.chambon@gmail.com> +# +# License: MIT License + +import matplotlib.pylab as pl +import ot + + +############################################################################## +# 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) + + +############################################################################## +# 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..3984afb --- /dev/null +++ b/examples/da/plot_otda_color_images.py @@ -0,0 +1,165 @@ +# -*- 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 <remi.flamary@unice.fr> +# Stanislas Chambon <stan.chambon@gmail.com> +# +# License: MIT License + +import numpy as np +from scipy import ndimage +import matplotlib.pylab as pl +import ot + + +r = np.random.RandomState(42) + + +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 + +X1 = im2mat(I1) +X2 = im2mat(I2) + +# training samples +nb = 1000 +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) + +# 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..3daa0a6 --- /dev/null +++ b/examples/da/plot_otda_d2.py @@ -0,0 +1,173 @@ +# -*- 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 <remi.flamary@unice.fr> +# Stanislas Chambon <stan.chambon@gmail.com> +# +# License: MIT License + +import matplotlib.pylab as pl +import ot + + +############################################################################## +# generate data +############################################################################## + +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, metric='sqeuclidean') + + +############################################################################## +# 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..09d2cb4 --- /dev/null +++ b/examples/da/plot_otda_mapping.py @@ -0,0 +1,126 @@ +# -*- 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 <remi.flamary@unice.fr> +# Stanislas Chambon <stan.chambon@gmail.com> +# +# License: MIT License + +import numpy as np +import matplotlib.pylab as pl +import ot + + +############################################################################## +# generate data +############################################################################## + +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_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_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, + 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..a628b05 --- /dev/null +++ b/examples/da/plot_otda_mapping_colors_images.py @@ -0,0 +1,171 @@ +# -*- 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 <remi.flamary@unice.fr> +# Stanislas Chambon <stan.chambon@gmail.com> +# +# License: MIT License + +import numpy as np +from scipy import ndimage +import matplotlib.pylab as pl +import ot + +r = np.random.RandomState(42) + + +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 + + +X1 = im2mat(I1) +X2 = im2mat(I2) + +# training samples +nb = 1000 +idx1 = r.randint(X1.shape[0], size=(nb,)) +idx2 = r.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(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(Xs=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() diff --git a/examples/plot_OTDA_2D.py b/examples/plot_OTDA_2D.py deleted file mode 100644 index a1fb804..0000000 --- a/examples/plot_OTDA_2D.py +++ /dev/null @@ -1,120 +0,0 @@ -# -*- coding: utf-8 -*- -""" -============================== -OT for empirical distributions -============================== - -""" - -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') - -#%% sample interpolation - -xst0=n*G0.dot(xt) -xsts=n*Gs.dot(xt) -xstg=n*Gg.dot(xt) - -pl.figure(4) -pl.subplot(2,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(2,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(2,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')
\ No newline at end of file diff --git a/examples/plot_OTDA_classes.py b/examples/plot_OTDA_classes.py deleted file mode 100644 index 089b45b..0000000 --- a/examples/plot_OTDA_classes.py +++ /dev/null @@ -1,112 +0,0 @@ -# -*- coding: utf-8 -*- -""" -======================== -OT for domain adaptation -======================== - -""" - -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) - -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') - - -#%% 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 -pl.figure(4,(15,8)) - -param_img={'interpolation':'nearest','cmap':'jet'} - -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 sinkhorn') - -pl.subplot(2,4,3) -pl.imshow(da_lpl1.G,**param_img) -pl.title('OT matrix non-convex Group Lasso') - -pl.subplot(2,4,4) -pl.imshow(da_l1l2.G,**param_img) -pl.title('OT matrix Group 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 Sinkhorn') - -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 non-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 Group Lasso')
\ No newline at end of file diff --git a/examples/plot_OTDA_color_images.py b/examples/plot_OTDA_color_images.py deleted file mode 100644 index 68eee44..0000000 --- a/examples/plot_OTDA_color_images.py +++ /dev/null @@ -1,145 +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. -""" - -import numpy as np -import scipy.ndimage as spi -import matplotlib.pylab as pl -import ot - - -#%% Loading images - -I1=spi.imread('../data/ocean_day.jpg').astype(np.float64)/256 -I2=spi.imread('../data/ocean_sunset.jpg').astype(np.float64)/256 - -#%% Plot images - -pl.figure(1) - -pl.subplot(1,2,1) -pl.imshow(I1) -pl.title('Image 1') - -pl.subplot(1,2,2) -pl.imshow(I2) -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,(10,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.imshow(I2) -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.show() - - - -#%% 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.minimum(np.maximum(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,(10,8)) - -pl.subplot(2,3,1) - -pl.imshow(I1) -pl.title('Image 1') - -pl.subplot(2,3,2) -pl.imshow(I1t) -pl.title('Image 1 Adapt') - - -pl.subplot(2,3,3) -pl.imshow(I1te) -pl.title('Image 1 Adapt (reg)') - -pl.subplot(2,3,4) - -pl.imshow(I2) -pl.title('Image 2') - -pl.subplot(2,3,5) -pl.imshow(I2t) -pl.title('Image 2 Adapt') - - -pl.subplot(2,3,6) -pl.imshow(I2te) -pl.title('Image 2 Adapt (reg)') - -pl.show() diff --git a/examples/plot_OTDA_mapping.py b/examples/plot_OTDA_mapping.py deleted file mode 100644 index 78b57e7..0000000 --- a/examples/plot_OTDA_mapping.py +++ /dev/null @@ -1,110 +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. -""" - -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,(8,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') - - - -#%% 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,(10,7)) -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") - - - -#%% 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,(10,7)) -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)") diff --git a/examples/plot_OTDA_mapping_color_images.py b/examples/plot_OTDA_mapping_color_images.py deleted file mode 100644 index f07dc6c..0000000 --- a/examples/plot_OTDA_mapping_color_images.py +++ /dev/null @@ -1,158 +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. - -""" - -import numpy as np -import scipy.ndimage as spi -import matplotlib.pylab as pl -import ot - - -#%% Loading images - -I1=spi.imread('../data/ocean_day.jpg').astype(np.float64)/256 -I2=spi.imread('../data/ocean_sunset.jpg').astype(np.float64)/256 - -#%% Plot images - -pl.figure(1) - -pl.subplot(1,2,1) -pl.imshow(I1) -pl.title('Image 1') - -pl.subplot(1,2,2) -pl.imshow(I2) -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,(10,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.imshow(I2) -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.show() - - - -#%% domain adaptation between images -def minmax(I): - return np.minimum(np.maximum(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,(10,8)) - -pl.subplot(2,3,1) - -pl.imshow(I1) -pl.title('Im. 1') - -pl.subplot(2,3,2) - -pl.imshow(I2) -pl.title('Im. 2') - - -pl.subplot(2,3,3) -pl.imshow(I1t) -pl.title('Im. 1 Interp LP') - -pl.subplot(2,3,4) -pl.imshow(I1te) -pl.title('Im. 1 Interp Entrop') - - -pl.subplot(2,3,5) -pl.imshow(I1tl) -pl.title('Im. 1 Linear mapping') - -pl.subplot(2,3,6) -pl.imshow(I1tn) -pl.title('Im. 1 nonlinear mapping') - -pl.show() diff --git a/examples/plot_OT_1D.py b/examples/plot_OT_1D.py index 6661aa3..0f3a26a 100644 --- a/examples/plot_OT_1D.py +++ b/examples/plot_OT_1D.py @@ -4,53 +4,57 @@ 1D optimal transport ==================== -@author: rflamary """ +# Author: Remi Flamary <remi.flamary@unice.fr> +# +# License: MIT License + import numpy as np import matplotlib.pylab as pl import ot from ot.datasets import get_1D_gauss as gauss - #%% parameters -n=100 # nb bins +n = 100 # nb bins # bin positions -x=np.arange(n,dtype=np.float64) +x = np.arange(n, dtype=np.float64) # Gaussian distributions -a=gauss(n,m=20,s=5) # m= mean, s= std -b=gauss(n,m=60,s=10) +a = gauss(n, m=20, s=5) # m= mean, s= std +b = gauss(n, m=60, s=10) # loss matrix -M=ot.dist(x.reshape((n,1)),x.reshape((n,1))) -M/=M.max() +M = ot.dist(x.reshape((n, 1)), x.reshape((n, 1))) +M /= M.max() #%% plot the distributions -pl.figure(1) -pl.plot(x,a,'b',label='Source distribution') -pl.plot(x,b,'r',label='Target distribution') +pl.figure(1, figsize=(6.4, 3)) +pl.plot(x, a, 'b', label='Source distribution') +pl.plot(x, b, 'r', label='Target distribution') pl.legend() #%% plot distributions and loss matrix -pl.figure(2) -ot.plot.plot1D_mat(a,b,M,'Cost matrix M') +pl.figure(2, figsize=(5, 5)) +ot.plot.plot1D_mat(a, b, M, 'Cost matrix M') #%% EMD -G0=ot.emd(a,b,M) +G0 = ot.emd(a, b, M) -pl.figure(3) -ot.plot.plot1D_mat(a,b,G0,'OT matrix G0') +pl.figure(3, figsize=(5, 5)) +ot.plot.plot1D_mat(a, b, G0, 'OT matrix G0') #%% Sinkhorn -lambd=1e-3 -Gs=ot.sinkhorn(a,b,M,lambd,verbose=True) +lambd = 1e-3 +Gs = ot.sinkhorn(a, b, M, lambd, verbose=True) + +pl.figure(4, figsize=(5, 5)) +ot.plot.plot1D_mat(a, b, Gs, 'OT matrix Sinkhorn') -pl.figure(4) -ot.plot.plot1D_mat(a,b,Gs,'OT matrix Sinkhorn') +pl.show() diff --git a/examples/plot_OT_2D_samples.py b/examples/plot_OT_2D_samples.py index edfb781..023e645 100644 --- a/examples/plot_OT_2D_samples.py +++ b/examples/plot_OT_2D_samples.py @@ -4,57 +4,60 @@ 2D Optimal transport between empirical distributions ==================================================== -@author: rflamary """ +# Author: Remi Flamary <remi.flamary@unice.fr> +# +# License: MIT License + import numpy as np import matplotlib.pylab as pl import ot #%% parameters and data generation -n=50 # nb samples +n = 50 # nb samples -mu_s=np.array([0,0]) -cov_s=np.array([[1,0],[0,1]]) +mu_s = np.array([0, 0]) +cov_s = np.array([[1, 0], [0, 1]]) -mu_t=np.array([4,4]) -cov_t=np.array([[1,-.8],[-.8,1]]) +mu_t = np.array([4, 4]) +cov_t = np.array([[1, -.8], [-.8, 1]]) -xs=ot.datasets.get_2D_samples_gauss(n,mu_s,cov_s) -xt=ot.datasets.get_2D_samples_gauss(n,mu_t,cov_t) +xs = ot.datasets.get_2D_samples_gauss(n, mu_s, cov_s) +xt = ot.datasets.get_2D_samples_gauss(n, mu_t, cov_t) -a,b = ot.unif(n),ot.unif(n) # uniform distribution on samples +a, b = np.ones((n,)) / n, np.ones((n,)) / n # uniform distribution on samples # loss matrix -M=ot.dist(xs,xt) -M/=M.max() +M = ot.dist(xs, xt) +M /= M.max() #%% plot samples pl.figure(1) -pl.plot(xs[:,0],xs[:,1],'+b',label='Source samples') -pl.plot(xt[:,0],xt[:,1],'xr',label='Target samples') +pl.plot(xs[:, 0], xs[:, 1], '+b', label='Source samples') +pl.plot(xt[:, 0], xt[:, 1], 'xr', label='Target samples') pl.legend(loc=0) -pl.title('Source and traget distributions') +pl.title('Source and target distributions') pl.figure(2) -pl.imshow(M,interpolation='nearest') +pl.imshow(M, interpolation='nearest') pl.title('Cost matrix M') #%% EMD -G0=ot.emd(a,b,M) +G0 = ot.emd(a, b, M) pl.figure(3) -pl.imshow(G0,interpolation='nearest') +pl.imshow(G0, interpolation='nearest') pl.title('OT matrix G0') pl.figure(4) -ot.plot.plot2D_samples_mat(xs,xt,G0,c=[.5,.5,1]) -pl.plot(xs[:,0],xs[:,1],'+b',label='Source samples') -pl.plot(xt[:,0],xt[:,1],'xr',label='Target samples') +ot.plot.plot2D_samples_mat(xs, xt, G0, c=[.5, .5, 1]) +pl.plot(xs[:, 0], xs[:, 1], '+b', label='Source samples') +pl.plot(xt[:, 0], xt[:, 1], 'xr', label='Target samples') pl.legend(loc=0) pl.title('OT matrix with samples') @@ -62,17 +65,19 @@ pl.title('OT matrix with samples') #%% sinkhorn # reg term -lambd=5e-4 +lambd = 5e-4 -Gs=ot.sinkhorn(a,b,M,lambd) +Gs = ot.sinkhorn(a, b, M, lambd) pl.figure(5) -pl.imshow(Gs,interpolation='nearest') +pl.imshow(Gs, interpolation='nearest') pl.title('OT matrix sinkhorn') pl.figure(6) -ot.plot.plot2D_samples_mat(xs,xt,Gs,color=[.5,.5,1]) -pl.plot(xs[:,0],xs[:,1],'+b',label='Source samples') -pl.plot(xt[:,0],xt[:,1],'xr',label='Target samples') +ot.plot.plot2D_samples_mat(xs, xt, Gs, color=[.5, .5, 1]) +pl.plot(xs[:, 0], xs[:, 1], '+b', label='Source samples') +pl.plot(xt[:, 0], xt[:, 1], 'xr', label='Target samples') pl.legend(loc=0) pl.title('OT matrix Sinkhorn with samples') + +pl.show() diff --git a/examples/plot_OT_L1_vs_L2.py b/examples/plot_OT_L1_vs_L2.py index 9bb92fe..dfc9462 100644 --- a/examples/plot_OT_L1_vs_L2.py +++ b/examples/plot_OT_L1_vs_L2.py @@ -4,13 +4,16 @@ 2D Optimal transport for different metrics ========================================== -Stole the figure idea from Fig. 1 and 2 in +Stole the figure idea from Fig. 1 and 2 in https://arxiv.org/pdf/1706.07650.pdf -@author: rflamary """ +# Author: Remi Flamary <remi.flamary@unice.fr> +# +# License: MIT License + import numpy as np import matplotlib.pylab as pl import ot @@ -20,89 +23,93 @@ import ot for data in range(2): if data: - n=20 # nb samples - xs=np.zeros((n,2)) - xs[:,0]=np.arange(n)+1 - xs[:,1]=(np.arange(n)+1)*-0.001 # to make it strictly convex... - - xt=np.zeros((n,2)) - xt[:,1]=np.arange(n)+1 + n = 20 # nb samples + xs = np.zeros((n, 2)) + xs[:, 0] = np.arange(n) + 1 + xs[:, 1] = (np.arange(n) + 1) * -0.001 # to make it strictly convex... + + xt = np.zeros((n, 2)) + xt[:, 1] = np.arange(n) + 1 else: - - n=50 # nb samples - xtot=np.zeros((n+1,2)) - xtot[:,0]=np.cos((np.arange(n+1)+1.0)*0.9/(n+2)*2*np.pi) - xtot[:,1]=np.sin((np.arange(n+1)+1.0)*0.9/(n+2)*2*np.pi) - - xs=xtot[:n,:] - xt=xtot[1:,:] - - - - a,b = ot.unif(n),ot.unif(n) # uniform distribution on samples - + + n = 50 # nb samples + xtot = np.zeros((n + 1, 2)) + xtot[:, 0] = np.cos( + (np.arange(n + 1) + 1.0) * 0.9 / (n + 2) * 2 * np.pi) + xtot[:, 1] = np.sin( + (np.arange(n + 1) + 1.0) * 0.9 / (n + 2) * 2 * np.pi) + + xs = xtot[:n, :] + xt = xtot[1:, :] + + a, b = ot.unif(n), ot.unif(n) # uniform distribution on samples + # loss matrix - M1=ot.dist(xs,xt,metric='euclidean') - M1/=M1.max() - + M1 = ot.dist(xs, xt, metric='euclidean') + M1 /= M1.max() + # loss matrix - M2=ot.dist(xs,xt,metric='sqeuclidean') - M2/=M2.max() - + M2 = ot.dist(xs, xt, metric='sqeuclidean') + M2 /= M2.max() + # loss matrix - Mp=np.sqrt(ot.dist(xs,xt,metric='euclidean')) - Mp/=Mp.max() - + Mp = np.sqrt(ot.dist(xs, xt, metric='euclidean')) + Mp /= Mp.max() + #%% plot samples - - pl.figure(1+3*data) + + pl.figure(1 + 3 * data, figsize=(7, 3)) pl.clf() - pl.plot(xs[:,0],xs[:,1],'+b',label='Source samples') - pl.plot(xt[:,0],xt[:,1],'xr',label='Target samples') + pl.plot(xs[:, 0], xs[:, 1], '+b', label='Source samples') + pl.plot(xt[:, 0], xt[:, 1], 'xr', label='Target samples') pl.axis('equal') pl.title('Source and traget distributions') - - pl.figure(2+3*data,(15,5)) - pl.subplot(1,3,1) - pl.imshow(M1,interpolation='nearest') - pl.title('Eucidean cost') - pl.subplot(1,3,2) - pl.imshow(M2,interpolation='nearest') + + pl.figure(2 + 3 * data, figsize=(7, 3)) + + pl.subplot(1, 3, 1) + pl.imshow(M1, interpolation='nearest') + pl.title('Euclidean cost') + + pl.subplot(1, 3, 2) + pl.imshow(M2, interpolation='nearest') pl.title('Squared Euclidean cost') - - pl.subplot(1,3,3) - pl.imshow(Mp,interpolation='nearest') + + pl.subplot(1, 3, 3) + pl.imshow(Mp, interpolation='nearest') pl.title('Sqrt Euclidean cost') + pl.tight_layout() + #%% EMD - - G1=ot.emd(a,b,M1) - G2=ot.emd(a,b,M2) - Gp=ot.emd(a,b,Mp) - - pl.figure(3+3*data,(15,5)) - - pl.subplot(1,3,1) - ot.plot.plot2D_samples_mat(xs,xt,G1,c=[.5,.5,1]) - pl.plot(xs[:,0],xs[:,1],'+b',label='Source samples') - pl.plot(xt[:,0],xt[:,1],'xr',label='Target samples') + G1 = ot.emd(a, b, M1) + G2 = ot.emd(a, b, M2) + Gp = ot.emd(a, b, Mp) + + pl.figure(3 + 3 * data, figsize=(7, 3)) + + pl.subplot(1, 3, 1) + ot.plot.plot2D_samples_mat(xs, xt, G1, c=[.5, .5, 1]) + pl.plot(xs[:, 0], xs[:, 1], '+b', label='Source samples') + pl.plot(xt[:, 0], xt[:, 1], 'xr', label='Target samples') pl.axis('equal') - #pl.legend(loc=0) + # pl.legend(loc=0) pl.title('OT Euclidean') - - pl.subplot(1,3,2) - - ot.plot.plot2D_samples_mat(xs,xt,G2,c=[.5,.5,1]) - pl.plot(xs[:,0],xs[:,1],'+b',label='Source samples') - pl.plot(xt[:,0],xt[:,1],'xr',label='Target samples') + + pl.subplot(1, 3, 2) + ot.plot.plot2D_samples_mat(xs, xt, G2, c=[.5, .5, 1]) + pl.plot(xs[:, 0], xs[:, 1], '+b', label='Source samples') + pl.plot(xt[:, 0], xt[:, 1], 'xr', label='Target samples') pl.axis('equal') - #pl.legend(loc=0) + # pl.legend(loc=0) pl.title('OT squared Euclidean') - - pl.subplot(1,3,3) - - ot.plot.plot2D_samples_mat(xs,xt,Gp,c=[.5,.5,1]) - pl.plot(xs[:,0],xs[:,1],'+b',label='Source samples') - pl.plot(xt[:,0],xt[:,1],'xr',label='Target samples') + + pl.subplot(1, 3, 3) + ot.plot.plot2D_samples_mat(xs, xt, Gp, c=[.5, .5, 1]) + pl.plot(xs[:, 0], xs[:, 1], '+b', label='Source samples') + pl.plot(xt[:, 0], xt[:, 1], 'xr', label='Target samples') pl.axis('equal') - #pl.legend(loc=0) + # pl.legend(loc=0) pl.title('OT sqrt Euclidean') + pl.tight_layout() + +pl.show() diff --git a/examples/plot_WDA.py b/examples/plot_WDA.py index 5fa2ab1..42789f2 100644 --- a/examples/plot_WDA.py +++ b/examples/plot_WDA.py @@ -4,93 +4,97 @@ Wasserstein Discriminant Analysis ================================= -@author: rflamary """ +# Author: Remi Flamary <remi.flamary@unice.fr> +# +# License: MIT License + import numpy as np import matplotlib.pylab as pl -import ot -from ot.datasets import get_1D_gauss as gauss + from ot.dr import wda, fda #%% parameters -n=1000 # nb samples in source and target datasets -nz=0.2 +n = 1000 # nb samples in source and target datasets +nz = 0.2 # generate circle dataset -t=np.random.rand(n)*2*np.pi -ys=np.floor((np.arange(n)*1.0/n*3))+1 -xs=np.concatenate((np.cos(t).reshape((-1,1)),np.sin(t).reshape((-1,1))),1) -xs=xs*ys.reshape(-1,1)+nz*np.random.randn(n,2) +t = np.random.rand(n) * 2 * np.pi +ys = np.floor((np.arange(n) * 1.0 / n * 3)) + 1 +xs = np.concatenate( + (np.cos(t).reshape((-1, 1)), np.sin(t).reshape((-1, 1))), 1) +xs = xs * ys.reshape(-1, 1) + nz * np.random.randn(n, 2) -t=np.random.rand(n)*2*np.pi -yt=np.floor((np.arange(n)*1.0/n*3))+1 -xt=np.concatenate((np.cos(t).reshape((-1,1)),np.sin(t).reshape((-1,1))),1) -xt=xt*yt.reshape(-1,1)+nz*np.random.randn(n,2) +t = np.random.rand(n) * 2 * np.pi +yt = np.floor((np.arange(n) * 1.0 / n * 3)) + 1 +xt = np.concatenate( + (np.cos(t).reshape((-1, 1)), np.sin(t).reshape((-1, 1))), 1) +xt = xt * yt.reshape(-1, 1) + nz * np.random.randn(n, 2) -nbnoise=8 +nbnoise = 8 -xs=np.hstack((xs,np.random.randn(n,nbnoise))) -xt=np.hstack((xt,np.random.randn(n,nbnoise))) +xs = np.hstack((xs, np.random.randn(n, nbnoise))) +xt = np.hstack((xt, np.random.randn(n, nbnoise))) #%% plot samples -pl.figure(1,(10,5)) +pl.figure(1, figsize=(6.4, 3.5)) -pl.subplot(1,2,1) -pl.scatter(xt[:,0],xt[:,1],c=ys,marker='+',label='Source samples') +pl.subplot(1, 2, 1) +pl.scatter(xt[:, 0], xt[:, 1], c=ys, marker='+', label='Source samples') pl.legend(loc=0) pl.title('Discriminant dimensions') -pl.subplot(1,2,2) -pl.scatter(xt[:,2],xt[:,3],c=ys,marker='+',label='Source samples') +pl.subplot(1, 2, 2) +pl.scatter(xt[:, 2], xt[:, 3], c=ys, marker='+', label='Source samples') pl.legend(loc=0) pl.title('Other dimensions') -pl.show() +pl.tight_layout() #%% Compute FDA -p=2 +p = 2 -Pfda,projfda = fda(xs,ys,p) +Pfda, projfda = fda(xs, ys, p) #%% Compute WDA -p=2 -reg=1e0 -k=10 -maxiter=100 +p = 2 +reg = 1e0 +k = 10 +maxiter = 100 -Pwda,projwda = wda(xs,ys,p,reg,k,maxiter=maxiter) +Pwda, projwda = wda(xs, ys, p, reg, k, maxiter=maxiter) #%% plot samples -xsp=projfda(xs) -xtp=projfda(xt) +xsp = projfda(xs) +xtp = projfda(xt) -xspw=projwda(xs) -xtpw=projwda(xt) +xspw = projwda(xs) +xtpw = projwda(xt) -pl.figure(1,(10,10)) +pl.figure(2) -pl.subplot(2,2,1) -pl.scatter(xsp[:,0],xsp[:,1],c=ys,marker='+',label='Projected samples') +pl.subplot(2, 2, 1) +pl.scatter(xsp[:, 0], xsp[:, 1], c=ys, marker='+', label='Projected samples') pl.legend(loc=0) pl.title('Projected training samples FDA') - -pl.subplot(2,2,2) -pl.scatter(xtp[:,0],xtp[:,1],c=ys,marker='+',label='Projected samples') +pl.subplot(2, 2, 2) +pl.scatter(xtp[:, 0], xtp[:, 1], c=ys, marker='+', label='Projected samples') pl.legend(loc=0) pl.title('Projected test samples FDA') - -pl.subplot(2,2,3) -pl.scatter(xspw[:,0],xspw[:,1],c=ys,marker='+',label='Projected samples') +pl.subplot(2, 2, 3) +pl.scatter(xspw[:, 0], xspw[:, 1], c=ys, marker='+', label='Projected samples') pl.legend(loc=0) pl.title('Projected training samples WDA') - -pl.subplot(2,2,4) -pl.scatter(xtpw[:,0],xtpw[:,1],c=ys,marker='+',label='Projected samples') +pl.subplot(2, 2, 4) +pl.scatter(xtpw[:, 0], xtpw[:, 1], c=ys, marker='+', label='Projected samples') pl.legend(loc=0) pl.title('Projected test samples WDA') +pl.tight_layout() + +pl.show() diff --git a/examples/plot_barycenter_1D.py b/examples/plot_barycenter_1D.py index 30eecbf..875f44c 100644 --- a/examples/plot_barycenter_1D.py +++ b/examples/plot_barycenter_1D.py @@ -4,135 +4,134 @@ 1D Wasserstein barycenter demo ============================== - -@author: rflamary """ +# Author: Remi Flamary <remi.flamary@unice.fr> +# +# License: MIT License + import numpy as np import matplotlib.pylab as pl import ot -from mpl_toolkits.mplot3d import Axes3D #necessary for 3d plot even if not used +# necessary for 3d plot even if not used +from mpl_toolkits.mplot3d import Axes3D # noqa from matplotlib.collections import PolyCollection #%% parameters -n=100 # nb bins +n = 100 # nb bins # bin positions -x=np.arange(n,dtype=np.float64) +x = np.arange(n, dtype=np.float64) # Gaussian distributions -a1=ot.datasets.get_1D_gauss(n,m=20,s=5) # m= mean, s= std -a2=ot.datasets.get_1D_gauss(n,m=60,s=8) +a1 = ot.datasets.get_1D_gauss(n, m=20, s=5) # m= mean, s= std +a2 = ot.datasets.get_1D_gauss(n, m=60, s=8) # creating matrix A containing all distributions -A=np.vstack((a1,a2)).T -nbd=A.shape[1] +A = np.vstack((a1, a2)).T +n_distributions = A.shape[1] # loss matrix + normalization -M=ot.utils.dist0(n) -M/=M.max() +M = ot.utils.dist0(n) +M /= M.max() #%% plot the distributions -pl.figure(1) -for i in range(nbd): - pl.plot(x,A[:,i]) +pl.figure(1, figsize=(6.4, 3)) +for i in range(n_distributions): + pl.plot(x, A[:, i]) pl.title('Distributions') +pl.tight_layout() #%% barycenter computation -alpha=0.2 # 0<=alpha<=1 -weights=np.array([1-alpha,alpha]) +alpha = 0.2 # 0<=alpha<=1 +weights = np.array([1 - alpha, alpha]) # l2bary -bary_l2=A.dot(weights) +bary_l2 = A.dot(weights) # wasserstein -reg=1e-3 -bary_wass=ot.bregman.barycenter(A,M,reg,weights) +reg = 1e-3 +bary_wass = ot.bregman.barycenter(A, M, reg, weights) pl.figure(2) pl.clf() -pl.subplot(2,1,1) -for i in range(nbd): - pl.plot(x,A[:,i]) +pl.subplot(2, 1, 1) +for i in range(n_distributions): + pl.plot(x, A[:, i]) pl.title('Distributions') -pl.subplot(2,1,2) -pl.plot(x,bary_l2,'r',label='l2') -pl.plot(x,bary_wass,'g',label='Wasserstein') +pl.subplot(2, 1, 2) +pl.plot(x, bary_l2, 'r', label='l2') +pl.plot(x, bary_wass, 'g', label='Wasserstein') pl.legend() pl.title('Barycenters') - +pl.tight_layout() #%% barycenter interpolation -nbalpha=11 -alphalist=np.linspace(0,1,nbalpha) +n_alpha = 11 +alpha_list = np.linspace(0, 1, n_alpha) -B_l2=np.zeros((n,nbalpha)) +B_l2 = np.zeros((n, n_alpha)) -B_wass=np.copy(B_l2) +B_wass = np.copy(B_l2) -for i in range(0,nbalpha): - alpha=alphalist[i] - weights=np.array([1-alpha,alpha]) - B_l2[:,i]=A.dot(weights) - B_wass[:,i]=ot.bregman.barycenter(A,M,reg,weights) +for i in range(0, n_alpha): + alpha = alpha_list[i] + weights = np.array([1 - alpha, alpha]) + B_l2[:, i] = A.dot(weights) + B_wass[:, i] = ot.bregman.barycenter(A, M, reg, weights) #%% plot interpolation -pl.figure(3,(10,5)) +pl.figure(3) -#pl.subplot(1,2,1) -cmap=pl.cm.get_cmap('viridis') +cmap = pl.cm.get_cmap('viridis') verts = [] -zs = alphalist -for i,z in enumerate(zs): - ys = B_l2[:,i] +zs = alpha_list +for i, z in enumerate(zs): + ys = B_l2[:, i] verts.append(list(zip(x, ys))) ax = pl.gcf().gca(projection='3d') -poly = PolyCollection(verts,facecolors=[cmap(a) for a in alphalist]) +poly = PolyCollection(verts, facecolors=[cmap(a) for a in alpha_list]) poly.set_alpha(0.7) ax.add_collection3d(poly, zs=zs, zdir='y') - ax.set_xlabel('x') ax.set_xlim3d(0, n) ax.set_ylabel('$\\alpha$') -ax.set_ylim3d(0,1) +ax.set_ylim3d(0, 1) ax.set_zlabel('') -ax.set_zlim3d(0, B_l2.max()*1.01) +ax.set_zlim3d(0, B_l2.max() * 1.01) pl.title('Barycenter interpolation with l2') +pl.tight_layout() -pl.show() - -pl.figure(4,(10,5)) - -#pl.subplot(1,2,1) -cmap=pl.cm.get_cmap('viridis') +pl.figure(4) +cmap = pl.cm.get_cmap('viridis') verts = [] -zs = alphalist -for i,z in enumerate(zs): - ys = B_wass[:,i] +zs = alpha_list +for i, z in enumerate(zs): + ys = B_wass[:, i] verts.append(list(zip(x, ys))) ax = pl.gcf().gca(projection='3d') -poly = PolyCollection(verts,facecolors=[cmap(a) for a in alphalist]) +poly = PolyCollection(verts, facecolors=[cmap(a) for a in alpha_list]) poly.set_alpha(0.7) ax.add_collection3d(poly, zs=zs, zdir='y') - ax.set_xlabel('x') ax.set_xlim3d(0, n) ax.set_ylabel('$\\alpha$') -ax.set_ylim3d(0,1) +ax.set_ylim3d(0, 1) ax.set_zlabel('') -ax.set_zlim3d(0, B_l2.max()*1.01) +ax.set_zlim3d(0, B_l2.max() * 1.01) pl.title('Barycenter interpolation with Wasserstein') +pl.tight_layout() -pl.show()
\ No newline at end of file +pl.show() diff --git a/examples/plot_compute_emd.py b/examples/plot_compute_emd.py index f2cdc35..893eecf 100644 --- a/examples/plot_compute_emd.py +++ b/examples/plot_compute_emd.py @@ -4,9 +4,12 @@ 1D optimal transport ==================== -@author: rflamary """ +# Author: Remi Flamary <remi.flamary@unice.fr> +# +# License: MIT License + import numpy as np import matplotlib.pylab as pl import ot @@ -15,60 +18,63 @@ from ot.datasets import get_1D_gauss as gauss #%% parameters -n=100 # nb bins -n_target=50 # nb target distributions +n = 100 # nb bins +n_target = 50 # nb target distributions # bin positions -x=np.arange(n,dtype=np.float64) +x = np.arange(n, dtype=np.float64) -lst_m=np.linspace(20,90,n_target) +lst_m = np.linspace(20, 90, n_target) # Gaussian distributions -a=gauss(n,m=20,s=5) # m= mean, s= std +a = gauss(n, m=20, s=5) # m= mean, s= std -B=np.zeros((n,n_target)) +B = np.zeros((n, n_target)) -for i,m in enumerate(lst_m): - B[:,i]=gauss(n,m=m,s=5) +for i, m in enumerate(lst_m): + B[:, i] = gauss(n, m=m, s=5) # loss matrix and normalization -M=ot.dist(x.reshape((n,1)),x.reshape((n,1)),'euclidean') -M/=M.max() -M2=ot.dist(x.reshape((n,1)),x.reshape((n,1)),'sqeuclidean') -M2/=M2.max() +M = ot.dist(x.reshape((n, 1)), x.reshape((n, 1)), 'euclidean') +M /= M.max() +M2 = ot.dist(x.reshape((n, 1)), x.reshape((n, 1)), 'sqeuclidean') +M2 /= M2.max() #%% plot the distributions pl.figure(1) -pl.subplot(2,1,1) -pl.plot(x,a,'b',label='Source distribution') +pl.subplot(2, 1, 1) +pl.plot(x, a, 'b', label='Source distribution') pl.title('Source distribution') -pl.subplot(2,1,2) -pl.plot(x,B,label='Target distributions') +pl.subplot(2, 1, 2) +pl.plot(x, B, label='Target distributions') pl.title('Target distributions') +pl.tight_layout() #%% Compute and plot distributions and loss matrix -d_emd=ot.emd2(a,B,M) # direct computation of EMD -d_emd2=ot.emd2(a,B,M2) # direct computation of EMD with loss M3 +d_emd = ot.emd2(a, B, M) # direct computation of EMD +d_emd2 = ot.emd2(a, B, M2) # direct computation of EMD with loss M3 pl.figure(2) -pl.plot(d_emd,label='Euclidean EMD') -pl.plot(d_emd2,label='Squared Euclidean EMD') +pl.plot(d_emd, label='Euclidean EMD') +pl.plot(d_emd2, label='Squared Euclidean EMD') pl.title('EMD distances') pl.legend() #%% -reg=1e-2 -d_sinkhorn=ot.sinkhorn2(a,B,M,reg) -d_sinkhorn2=ot.sinkhorn2(a,B,M2,reg) +reg = 1e-2 +d_sinkhorn = ot.sinkhorn2(a, B, M, reg) +d_sinkhorn2 = ot.sinkhorn2(a, B, M2, reg) pl.figure(2) pl.clf() -pl.plot(d_emd,label='Euclidean EMD') -pl.plot(d_emd2,label='Squared Euclidean EMD') -pl.plot(d_sinkhorn,'+',label='Euclidean Sinkhorn') -pl.plot(d_sinkhorn2,'+',label='Squared Euclidean Sinkhorn') +pl.plot(d_emd, label='Euclidean EMD') +pl.plot(d_emd2, label='Squared Euclidean EMD') +pl.plot(d_sinkhorn, '+', label='Euclidean Sinkhorn') +pl.plot(d_sinkhorn2, '+', label='Squared Euclidean Sinkhorn') pl.title('EMD distances') -pl.legend()
\ No newline at end of file +pl.legend() + +pl.show() diff --git a/examples/plot_optim_OTreg.py b/examples/plot_optim_OTreg.py index 8abb426..276b250 100644 --- a/examples/plot_optim_OTreg.py +++ b/examples/plot_optim_OTreg.py @@ -12,63 +12,80 @@ import matplotlib.pylab as pl import ot - #%% parameters -n=100 # nb bins +n = 100 # nb bins # bin positions -x=np.arange(n,dtype=np.float64) +x = np.arange(n, dtype=np.float64) # Gaussian distributions -a=ot.datasets.get_1D_gauss(n,m=20,s=5) # m= mean, s= std -b=ot.datasets.get_1D_gauss(n,m=60,s=10) +a = ot.datasets.get_1D_gauss(n, m=20, s=5) # m= mean, s= std +b = ot.datasets.get_1D_gauss(n, m=60, s=10) # loss matrix -M=ot.dist(x.reshape((n,1)),x.reshape((n,1))) -M/=M.max() +M = ot.dist(x.reshape((n, 1)), x.reshape((n, 1))) +M /= M.max() #%% EMD -G0=ot.emd(a,b,M) +G0 = ot.emd(a, b, M) -pl.figure(3) -ot.plot.plot1D_mat(a,b,G0,'OT matrix G0') +pl.figure(3, figsize=(5, 5)) +ot.plot.plot1D_mat(a, b, G0, 'OT matrix G0') #%% Example with Frobenius norm regularization -def f(G): return 0.5*np.sum(G**2) -def df(G): return G -reg=1e-1 +def f(G): + return 0.5 * np.sum(G**2) + + +def df(G): + return G -Gl2=ot.optim.cg(a,b,M,reg,f,df,verbose=True) + +reg = 1e-1 + +Gl2 = ot.optim.cg(a, b, M, reg, f, df, verbose=True) pl.figure(3) -ot.plot.plot1D_mat(a,b,Gl2,'OT matrix Frob. reg') +ot.plot.plot1D_mat(a, b, Gl2, 'OT matrix Frob. reg') #%% Example with entropic regularization -def f(G): return np.sum(G*np.log(G)) -def df(G): return np.log(G)+1 -reg=1e-3 +def f(G): + return np.sum(G * np.log(G)) + -Ge=ot.optim.cg(a,b,M,reg,f,df,verbose=True) +def df(G): + return np.log(G) + 1. -pl.figure(4) -ot.plot.plot1D_mat(a,b,Ge,'OT matrix Entrop. reg') + +reg = 1e-3 + +Ge = ot.optim.cg(a, b, M, reg, f, df, verbose=True) + +pl.figure(4, figsize=(5, 5)) +ot.plot.plot1D_mat(a, b, Ge, 'OT matrix Entrop. reg') #%% Example with Frobenius norm + entropic regularization with gcg -def f(G): return 0.5*np.sum(G**2) -def df(G): return G -reg1=1e-3 -reg2=1e-1 +def f(G): + return 0.5 * np.sum(G**2) + + +def df(G): + return G + + +reg1 = 1e-3 +reg2 = 1e-1 -Gel2=ot.optim.gcg(a,b,M,reg1,reg2,f,df,verbose=True) +Gel2 = ot.optim.gcg(a, b, M, reg1, reg2, f, df, verbose=True) -pl.figure(5) -ot.plot.plot1D_mat(a,b,Gel2,'OT entropic + matrix Frob. reg') -pl.show()
\ No newline at end of file +pl.figure(5, figsize=(5, 5)) +ot.plot.plot1D_mat(a, b, Gel2, 'OT entropic + matrix Frob. reg') +pl.show() diff --git a/ot/__init__.py b/ot/__init__.py index a79a5ce..6d4c4c6 100644 --- a/ot/__init__.py +++ b/ot/__init__.py @@ -4,6 +4,10 @@ """ +# Author: Remi Flamary <remi.flamary@unice.fr> +# +# License: MIT License + # All submodules and packages from . import lp @@ -24,6 +28,6 @@ from .utils import dist, unif, tic, toc, toq __version__ = "0.3.1" -__all__ = ["emd", "emd2", "sinkhorn","sinkhorn2", "utils", 'datasets', +__all__ = ["emd", "emd2", "sinkhorn", "sinkhorn2", "utils", 'datasets', 'bregman', 'lp', 'plot', 'tic', 'toc', 'toq', 'dist', 'unif', 'barycenter', 'sinkhorn_lpl1_mm', 'da', 'optim'] diff --git a/ot/bregman.py b/ot/bregman.py index 0d68602..d63c51d 100644 --- a/ot/bregman.py +++ b/ot/bregman.py @@ -3,9 +3,15 @@ Bregman projections for regularized OT """ +# Author: Remi Flamary <remi.flamary@unice.fr> +# Nicolas Courty <ncourty@irisa.fr> +# +# License: MIT License + import numpy as np -def sinkhorn(a,b, M, reg,method='sinkhorn', numItermax = 1000, stopThr=1e-9, verbose=False, log=False,**kwargs): + +def sinkhorn(a, b, M, reg, method='sinkhorn', numItermax=1000, stopThr=1e-9, verbose=False, log=False, **kwargs): u""" Solve the entropic regularization optimal transport problem and return the OT matrix @@ -92,22 +98,29 @@ def sinkhorn(a,b, M, reg,method='sinkhorn', numItermax = 1000, stopThr=1e-9, ver """ - if method.lower()=='sinkhorn': - sink= lambda: sinkhorn_knopp(a,b, M, reg,numItermax=numItermax, - stopThr=stopThr, verbose=verbose, log=log,**kwargs) - elif method.lower()=='sinkhorn_stabilized': - sink= lambda: sinkhorn_stabilized(a,b, M, reg,numItermax=numItermax, - stopThr=stopThr, verbose=verbose, log=log, **kwargs) - elif method.lower()=='sinkhorn_epsilon_scaling': - sink= lambda: sinkhorn_epsilon_scaling(a,b, M, reg,numItermax=numItermax, - stopThr=stopThr, verbose=verbose, log=log, **kwargs) + if method.lower() == 'sinkhorn': + def sink(): + return sinkhorn_knopp(a, b, M, reg, numItermax=numItermax, + stopThr=stopThr, verbose=verbose, log=log, **kwargs) + elif method.lower() == 'sinkhorn_stabilized': + def sink(): + return sinkhorn_stabilized(a, b, M, reg, numItermax=numItermax, + stopThr=stopThr, verbose=verbose, log=log, **kwargs) + elif method.lower() == 'sinkhorn_epsilon_scaling': + def sink(): + return sinkhorn_epsilon_scaling( + a, b, M, reg, numItermax=numItermax, + stopThr=stopThr, verbose=verbose, log=log, **kwargs) else: print('Warning : unknown method using classic Sinkhorn Knopp') - sink= lambda: sinkhorn_knopp(a,b, M, reg, **kwargs) + + def sink(): + return sinkhorn_knopp(a, b, M, reg, **kwargs) return sink() -def sinkhorn2(a,b, M, reg,method='sinkhorn', numItermax = 1000, stopThr=1e-9, verbose=False, log=False,**kwargs): + +def sinkhorn2(a, b, M, reg, method='sinkhorn', numItermax=1000, stopThr=1e-9, verbose=False, log=False, **kwargs): u""" Solve the entropic regularization optimal transport problem and return the loss @@ -170,7 +183,7 @@ def sinkhorn2(a,b, M, reg,method='sinkhorn', numItermax = 1000, stopThr=1e-9, ve >>> M=[[0.,1.],[1.,0.]] >>> ot.sinkhorn2(a,b,M,1) array([ 0.26894142]) - + References @@ -194,27 +207,33 @@ def sinkhorn2(a,b, M, reg,method='sinkhorn', numItermax = 1000, stopThr=1e-9, ve """ - if method.lower()=='sinkhorn': - sink= lambda: sinkhorn_knopp(a,b, M, reg,numItermax=numItermax, - stopThr=stopThr, verbose=verbose, log=log,**kwargs) - elif method.lower()=='sinkhorn_stabilized': - sink= lambda: sinkhorn_stabilized(a,b, M, reg,numItermax=numItermax, - stopThr=stopThr, verbose=verbose, log=log, **kwargs) - elif method.lower()=='sinkhorn_epsilon_scaling': - sink= lambda: sinkhorn_epsilon_scaling(a,b, M, reg,numItermax=numItermax, - stopThr=stopThr, verbose=verbose, log=log, **kwargs) + if method.lower() == 'sinkhorn': + def sink(): + return sinkhorn_knopp(a, b, M, reg, numItermax=numItermax, + stopThr=stopThr, verbose=verbose, log=log, **kwargs) + elif method.lower() == 'sinkhorn_stabilized': + def sink(): + return sinkhorn_stabilized(a, b, M, reg, numItermax=numItermax, + stopThr=stopThr, verbose=verbose, log=log, **kwargs) + elif method.lower() == 'sinkhorn_epsilon_scaling': + def sink(): + return sinkhorn_epsilon_scaling( + a, b, M, reg, numItermax=numItermax, + stopThr=stopThr, verbose=verbose, log=log, **kwargs) else: print('Warning : unknown method using classic Sinkhorn Knopp') - sink= lambda: sinkhorn_knopp(a,b, M, reg, **kwargs) - - b=np.asarray(b,dtype=np.float64) - if len(b.shape)<2: - b=b.reshape((-1,1)) + + def sink(): + return sinkhorn_knopp(a, b, M, reg, **kwargs) + + b = np.asarray(b, dtype=np.float64) + if len(b.shape) < 2: + b = b.reshape((-1, 1)) return sink() -def sinkhorn_knopp(a,b, M, reg, numItermax = 1000, stopThr=1e-9, verbose=False, log=False,**kwargs): +def sinkhorn_knopp(a, b, M, reg, numItermax=1000, stopThr=1e-9, verbose=False, log=False, **kwargs): """ Solve the entropic regularization optimal transport problem and return the OT matrix @@ -290,100 +309,101 @@ def sinkhorn_knopp(a,b, M, reg, numItermax = 1000, stopThr=1e-9, verbose=False, """ - a=np.asarray(a,dtype=np.float64) - b=np.asarray(b,dtype=np.float64) - M=np.asarray(M,dtype=np.float64) - - - if len(a)==0: - a=np.ones((M.shape[0],),dtype=np.float64)/M.shape[0] - if len(b)==0: - b=np.ones((M.shape[1],),dtype=np.float64)/M.shape[1] + a = np.asarray(a, dtype=np.float64) + b = np.asarray(b, dtype=np.float64) + M = np.asarray(M, dtype=np.float64) + if len(a) == 0: + a = np.ones((M.shape[0],), dtype=np.float64) / M.shape[0] + if len(b) == 0: + b = np.ones((M.shape[1],), dtype=np.float64) / M.shape[1] # init data Nini = len(a) Nfin = len(b) - if len(b.shape)>1: - nbb=b.shape[1] + if len(b.shape) > 1: + nbb = b.shape[1] else: - nbb=0 - + nbb = 0 if log: - log={'err':[]} + log = {'err': []} - # we assume that no distances are null except those of the diagonal of distances + # we assume that no distances are null except those of the diagonal of + # distances if nbb: - u = np.ones((Nini,nbb))/Nini - v = np.ones((Nfin,nbb))/Nfin + u = np.ones((Nini, nbb)) / Nini + v = np.ones((Nfin, nbb)) / Nfin else: - u = np.ones(Nini)/Nini - v = np.ones(Nfin)/Nfin + u = np.ones(Nini) / Nini + v = np.ones(Nfin) / Nfin + # print(reg) - #print(reg) + K = np.exp(-M / reg) + # print(np.min(K)) - K = np.exp(-M/reg) - #print(np.min(K)) - - Kp = (1/a).reshape(-1, 1) * K + Kp = (1 / a).reshape(-1, 1) * K cpt = 0 - err=1 - while (err>stopThr and cpt<numItermax): + err = 1 + while (err > stopThr and cpt < numItermax): uprev = u vprev = v KtransposeU = np.dot(K.T, u) v = np.divide(b, KtransposeU) - u = 1./np.dot(Kp,v) + u = 1. / np.dot(Kp, v) - if (np.any(KtransposeU==0) or - np.any(np.isnan(u)) or np.any(np.isnan(v)) or - np.any(np.isinf(u)) or np.any(np.isinf(v))): + if (np.any(KtransposeU == 0) or + np.any(np.isnan(u)) or np.any(np.isnan(v)) or + np.any(np.isinf(u)) or np.any(np.isinf(v))): # we have reached the machine precision # come back to previous solution and quit loop print('Warning: numerical errors at iteration', cpt) u = uprev v = vprev break - if cpt%10==0: - # we can speed up the process by checking for the error only all the 10th iterations + if cpt % 10 == 0: + # we can speed up the process by checking for the error only all + # the 10th iterations if nbb: - err = np.sum((u-uprev)**2)/np.sum((u)**2)+np.sum((v-vprev)**2)/np.sum((v)**2) + err = np.sum((u - uprev)**2) / np.sum((u)**2) + \ + np.sum((v - vprev)**2) / np.sum((v)**2) else: transp = u.reshape(-1, 1) * (K * v) - err = np.linalg.norm((np.sum(transp,axis=0)-b))**2 + err = np.linalg.norm((np.sum(transp, axis=0) - b))**2 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: - log['u']=u - log['v']=v + log['u'] = u + log['v'] = v - if nbb: #return only loss - res=np.zeros((nbb)) + if nbb: # return only loss + res = np.zeros((nbb)) for i in range(nbb): - res[i]=np.sum(u[:,i].reshape((-1,1))*K*v[:,i].reshape((1,-1))*M) + res[i] = np.sum( + u[:, i].reshape((-1, 1)) * K * v[:, i].reshape((1, -1)) * M) if log: - return res,log + return res, log else: return res - else: # return OT matrix + else: # return OT matrix if log: - return u.reshape((-1,1))*K*v.reshape((1,-1)),log + return u.reshape((-1, 1)) * K * v.reshape((1, -1)), log else: - return u.reshape((-1,1))*K*v.reshape((1,-1)) + return u.reshape((-1, 1)) * K * v.reshape((1, -1)) -def sinkhorn_stabilized(a,b, M, reg, numItermax = 1000,tau=1e3, stopThr=1e-9,warmstart=None, verbose=False,print_period=20, log=False,**kwargs): +def sinkhorn_stabilized(a, b, M, reg, numItermax=1000, tau=1e3, stopThr=1e-9, warmstart=None, verbose=False, print_period=20, log=False, **kwargs): """ Solve the entropic regularization OT problem with log stabilization @@ -468,21 +488,21 @@ def sinkhorn_stabilized(a,b, M, reg, numItermax = 1000,tau=1e3, stopThr=1e-9,war """ - a=np.asarray(a,dtype=np.float64) - b=np.asarray(b,dtype=np.float64) - M=np.asarray(M,dtype=np.float64) + a = np.asarray(a, dtype=np.float64) + b = np.asarray(b, dtype=np.float64) + M = np.asarray(M, dtype=np.float64) - if len(a)==0: - a=np.ones((M.shape[0],),dtype=np.float64)/M.shape[0] - if len(b)==0: - b=np.ones((M.shape[1],),dtype=np.float64)/M.shape[1] + if len(a) == 0: + a = np.ones((M.shape[0],), dtype=np.float64) / M.shape[0] + if len(b) == 0: + b = np.ones((M.shape[1],), dtype=np.float64) / M.shape[1] # test if multiple target - if len(b.shape)>1: - nbb=b.shape[1] - a=a[:,np.newaxis] + if len(b.shape) > 1: + nbb = b.shape[1] + a = a[:, np.newaxis] else: - nbb=0 + nbb = 0 # init data na = len(a) @@ -490,81 +510,80 @@ def sinkhorn_stabilized(a,b, M, reg, numItermax = 1000,tau=1e3, stopThr=1e-9,war cpt = 0 if log: - log={'err':[]} + log = {'err': []} - # we assume that no distances are null except those of the diagonal of distances + # we assume that no distances are null except those of the diagonal of + # distances if warmstart is None: - alpha,beta=np.zeros(na),np.zeros(nb) + alpha, beta = np.zeros(na), np.zeros(nb) else: - alpha,beta=warmstart + alpha, beta = warmstart if nbb: - u,v = np.ones((na,nbb))/na,np.ones((nb,nbb))/nb + u, v = np.ones((na, nbb)) / na, np.ones((nb, nbb)) / nb else: - u,v = np.ones(na)/na,np.ones(nb)/nb + u, v = np.ones(na) / na, np.ones(nb) / nb - def get_K(alpha,beta): + def get_K(alpha, beta): """log space computation""" - return np.exp(-(M-alpha.reshape((na,1))-beta.reshape((1,nb)))/reg) + return np.exp(-(M - alpha.reshape((na, 1)) - beta.reshape((1, nb))) / reg) - def get_Gamma(alpha,beta,u,v): + def get_Gamma(alpha, beta, u, v): """log space gamma computation""" - return np.exp(-(M-alpha.reshape((na,1))-beta.reshape((1,nb)))/reg+np.log(u.reshape((na,1)))+np.log(v.reshape((1,nb)))) + return np.exp(-(M - alpha.reshape((na, 1)) - beta.reshape((1, nb))) / reg + np.log(u.reshape((na, 1))) + np.log(v.reshape((1, nb)))) - #print(np.min(K)) + # print(np.min(K)) - K=get_K(alpha,beta) + K = get_K(alpha, beta) transp = K - loop=1 + loop = 1 cpt = 0 - err=1 + err = 1 while loop: - - uprev = u vprev = v # sinkhorn update - v = b/(np.dot(K.T,u)+1e-16) - u = a/(np.dot(K,v)+1e-16) - + v = b / (np.dot(K.T, u) + 1e-16) + u = a / (np.dot(K, v) + 1e-16) # remove numerical problems and store them in K - if np.abs(u).max()>tau or np.abs(v).max()>tau: + if np.abs(u).max() > tau or np.abs(v).max() > tau: if nbb: - alpha,beta=alpha+reg*np.max(np.log(u),1),beta+reg*np.max(np.log(v)) + alpha, beta = alpha + reg * \ + np.max(np.log(u), 1), beta + reg * np.max(np.log(v)) else: - alpha,beta=alpha+reg*np.log(u),beta+reg*np.log(v) + alpha, beta = alpha + reg * np.log(u), beta + reg * np.log(v) if nbb: - u,v = np.ones((na,nbb))/na,np.ones((nb,nbb))/nb + u, v = np.ones((na, nbb)) / na, np.ones((nb, nbb)) / nb else: - u,v = np.ones(na)/na,np.ones(nb)/nb - K=get_K(alpha,beta) - + u, v = np.ones(na) / na, np.ones(nb) / nb + K = get_K(alpha, beta) - if cpt%print_period==0: - # we can speed up the process by checking for the error only all the 10th iterations + if cpt % print_period == 0: + # we can speed up the process by checking for the error only all + # the 10th iterations if nbb: - err = np.sum((u-uprev)**2)/np.sum((u)**2)+np.sum((v-vprev)**2)/np.sum((v)**2) + err = np.sum((u - uprev)**2) / np.sum((u)**2) + \ + np.sum((v - vprev)**2) / np.sum((v)**2) else: - transp = get_Gamma(alpha,beta,u,v) - err = np.linalg.norm((np.sum(transp,axis=0)-b))**2 + transp = get_Gamma(alpha, beta, u, v) + err = np.linalg.norm((np.sum(transp, axis=0) - b))**2 if log: log['err'].append(err) if verbose: - if cpt%(print_period*20) ==0: - print('{:5s}|{:12s}'.format('It.','Err')+'\n'+'-'*19) - print('{:5d}|{:8e}|'.format(cpt,err)) + if cpt % (print_period * 20) == 0: + print( + '{:5s}|{:12s}'.format('It.', 'Err') + '\n' + '-' * 19) + print('{:5d}|{:8e}|'.format(cpt, err)) + if err <= stopThr: + loop = False - if err<=stopThr: - loop=False - - if cpt>=numItermax: - loop=False - + if cpt >= numItermax: + loop = False if np.any(np.isnan(u)) or np.any(np.isnan(v)): # we have reached the machine precision @@ -574,34 +593,34 @@ def sinkhorn_stabilized(a,b, M, reg, numItermax = 1000,tau=1e3, stopThr=1e-9,war v = vprev break - cpt = cpt +1 + cpt = cpt + 1 - - #print('err=',err,' cpt=',cpt) + # print('err=',err,' cpt=',cpt) if log: - log['logu']=alpha/reg+np.log(u) - log['logv']=beta/reg+np.log(v) - log['alpha']=alpha+reg*np.log(u) - log['beta']=beta+reg*np.log(v) - log['warmstart']=(log['alpha'],log['beta']) + log['logu'] = alpha / reg + np.log(u) + log['logv'] = beta / reg + np.log(v) + log['alpha'] = alpha + reg * np.log(u) + log['beta'] = beta + reg * np.log(v) + log['warmstart'] = (log['alpha'], log['beta']) if nbb: - res=np.zeros((nbb)) + res = np.zeros((nbb)) for i in range(nbb): - res[i]=np.sum(get_Gamma(alpha,beta,u[:,i],v[:,i])*M) - return res,log + res[i] = np.sum(get_Gamma(alpha, beta, u[:, i], v[:, i]) * M) + return res, log else: - return get_Gamma(alpha,beta,u,v),log + return get_Gamma(alpha, beta, u, v), log else: if nbb: - res=np.zeros((nbb)) + res = np.zeros((nbb)) for i in range(nbb): - res[i]=np.sum(get_Gamma(alpha,beta,u[:,i],v[:,i])*M) + res[i] = np.sum(get_Gamma(alpha, beta, u[:, i], v[:, i]) * M) return res else: - return get_Gamma(alpha,beta,u,v) + return get_Gamma(alpha, beta, u, v) + -def sinkhorn_epsilon_scaling(a,b, M, reg, numItermax = 100, epsilon0=1e4, numInnerItermax = 100,tau=1e3, stopThr=1e-9,warmstart=None, verbose=False,print_period=10, log=False,**kwargs): +def sinkhorn_epsilon_scaling(a, b, M, reg, numItermax=100, epsilon0=1e4, numInnerItermax=100, tau=1e3, stopThr=1e-9, warmstart=None, verbose=False, print_period=10, log=False, **kwargs): """ Solve the entropic regularization optimal transport problem with log stabilization and epsilon scaling. @@ -690,14 +709,14 @@ def sinkhorn_epsilon_scaling(a,b, M, reg, numItermax = 100, epsilon0=1e4, numInn """ - a=np.asarray(a,dtype=np.float64) - b=np.asarray(b,dtype=np.float64) - M=np.asarray(M,dtype=np.float64) + a = np.asarray(a, dtype=np.float64) + b = np.asarray(b, dtype=np.float64) + M = np.asarray(M, dtype=np.float64) - if len(a)==0: - a=np.ones((M.shape[0],),dtype=np.float64)/M.shape[0] - if len(b)==0: - b=np.ones((M.shape[1],),dtype=np.float64)/M.shape[1] + if len(a) == 0: + a = np.ones((M.shape[0],), dtype=np.float64) / M.shape[0] + if len(b) == 0: + b = np.ones((M.shape[1],), dtype=np.float64) / M.shape[1] # init data na = len(a) @@ -705,88 +724,94 @@ def sinkhorn_epsilon_scaling(a,b, M, reg, numItermax = 100, epsilon0=1e4, numInn # nrelative umerical precision with 64 bits numItermin = 35 - numItermax=max(numItermin,numItermax) # ensure that last velue is exact - + numItermax = max(numItermin, numItermax) # ensure that last velue is exact cpt = 0 if log: - log={'err':[]} + log = {'err': []} - # we assume that no distances are null except those of the diagonal of distances + # we assume that no distances are null except those of the diagonal of + # distances if warmstart is None: - alpha,beta=np.zeros(na),np.zeros(nb) + alpha, beta = np.zeros(na), np.zeros(nb) else: - alpha,beta=warmstart - + alpha, beta = warmstart - def get_K(alpha,beta): + def get_K(alpha, beta): """log space computation""" - return np.exp(-(M-alpha.reshape((na,1))-beta.reshape((1,nb)))/reg) + return np.exp(-(M - alpha.reshape((na, 1)) - beta.reshape((1, nb))) / reg) - #print(np.min(K)) - def get_reg(n): # exponential decreasing - return (epsilon0-reg)*np.exp(-n)+reg + # print(np.min(K)) + def get_reg(n): # exponential decreasing + return (epsilon0 - reg) * np.exp(-n) + reg - loop=1 + loop = 1 cpt = 0 - err=1 + err = 1 while loop: - regi=get_reg(cpt) + regi = get_reg(cpt) - G,logi=sinkhorn_stabilized(a,b, M, regi, numItermax = numInnerItermax, stopThr=1e-9,warmstart=(alpha,beta), verbose=False,print_period=20,tau=tau, log=True) + G, logi = sinkhorn_stabilized(a, b, M, regi, numItermax=numInnerItermax, stopThr=1e-9, warmstart=( + alpha, beta), verbose=False, print_period=20, tau=tau, log=True) - alpha=logi['alpha'] - beta=logi['beta'] + alpha = logi['alpha'] + beta = logi['beta'] - if cpt>=numItermax: - loop=False + if cpt >= numItermax: + loop = False - if cpt%(print_period)==0: # spsion nearly converged - # we can speed up the process by checking for the error only all the 10th iterations + if cpt % (print_period) == 0: # spsion nearly converged + # we can speed up the process by checking for the error only all + # the 10th iterations transp = G - err = np.linalg.norm((np.sum(transp,axis=0)-b))**2+np.linalg.norm((np.sum(transp,axis=1)-a))**2 + err = np.linalg.norm( + (np.sum(transp, axis=0) - b))**2 + np.linalg.norm((np.sum(transp, axis=1) - a))**2 if log: log['err'].append(err) if verbose: - if cpt%(print_period*10) ==0: - print('{:5s}|{:12s}'.format('It.','Err')+'\n'+'-'*19) - print('{:5d}|{:8e}|'.format(cpt,err)) + if cpt % (print_period * 10) == 0: + print( + '{:5s}|{:12s}'.format('It.', 'Err') + '\n' + '-' * 19) + print('{:5d}|{:8e}|'.format(cpt, err)) - if err<=stopThr and cpt>numItermin: - loop=False + if err <= stopThr and cpt > numItermin: + loop = False - cpt = cpt +1 - #print('err=',err,' cpt=',cpt) + cpt = cpt + 1 + # print('err=',err,' cpt=',cpt) if log: - log['alpha']=alpha - log['beta']=beta - log['warmstart']=(log['alpha'],log['beta']) - return G,log + log['alpha'] = alpha + log['beta'] = beta + log['warmstart'] = (log['alpha'], log['beta']) + return G, log else: return G -def geometricBar(weights,alldistribT): +def geometricBar(weights, alldistribT): """return the weighted geometric mean of distributions""" - assert(len(weights)==alldistribT.shape[1]) - return np.exp(np.dot(np.log(alldistribT),weights.T)) + assert(len(weights) == alldistribT.shape[1]) + return np.exp(np.dot(np.log(alldistribT), weights.T)) + def geometricMean(alldistribT): """return the geometric mean of distributions""" - return np.exp(np.mean(np.log(alldistribT),axis=1)) + return np.exp(np.mean(np.log(alldistribT), axis=1)) -def projR(gamma,p): + +def projR(gamma, p): """return the KL projection on the row constrints """ - return np.multiply(gamma.T,p/np.maximum(np.sum(gamma,axis=1),1e-10)).T + return np.multiply(gamma.T, p / np.maximum(np.sum(gamma, axis=1), 1e-10)).T + -def projC(gamma,q): +def projC(gamma, q): """return the KL projection on the column constrints """ - return np.multiply(gamma,q/np.maximum(np.sum(gamma,axis=0),1e-10)) + return np.multiply(gamma, q / np.maximum(np.sum(gamma, axis=0), 1e-10)) -def barycenter(A,M,reg, weights=None, numItermax = 1000, stopThr=1e-4,verbose=False,log=False): +def barycenter(A, M, reg, weights=None, numItermax=1000, stopThr=1e-4, verbose=False, log=False): """Compute the entropic regularized wasserstein barycenter of distributions A The function solves the following optimization problem: @@ -837,49 +862,49 @@ def barycenter(A,M,reg, weights=None, numItermax = 1000, stopThr=1e-4,verbose=Fa """ - if weights is None: - weights=np.ones(A.shape[1])/A.shape[1] + weights = np.ones(A.shape[1]) / A.shape[1] else: - assert(len(weights)==A.shape[1]) + assert(len(weights) == A.shape[1]) if log: - log={'err':[]} + log = {'err': []} - #M = M/np.median(M) # suggested by G. Peyre - K = np.exp(-M/reg) + # M = M/np.median(M) # suggested by G. Peyre + K = np.exp(-M / reg) cpt = 0 - err=1 + err = 1 - UKv=np.dot(K,np.divide(A.T,np.sum(K,axis=0)).T) - u = (geometricMean(UKv)/UKv.T).T + UKv = np.dot(K, np.divide(A.T, np.sum(K, axis=0)).T) + u = (geometricMean(UKv) / UKv.T).T - while (err>stopThr and cpt<numItermax): - cpt = cpt +1 - UKv=u*np.dot(K,np.divide(A,np.dot(K,u))) - u = (u.T*geometricBar(weights,UKv)).T/UKv + while (err > stopThr and cpt < numItermax): + cpt = cpt + 1 + UKv = u * np.dot(K, np.divide(A, np.dot(K, u))) + u = (u.T * geometricBar(weights, UKv)).T / UKv - if cpt%10==1: - err=np.sum(np.std(UKv,axis=1)) + if cpt % 10 == 1: + err = np.sum(np.std(UKv, axis=1)) # log and verbose print 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)) + if cpt % 200 == 0: + print( + '{:5s}|{:12s}'.format('It.', 'Err') + '\n' + '-' * 19) + print('{:5d}|{:8e}|'.format(cpt, err)) if log: - log['niter']=cpt - return geometricBar(weights,UKv),log + log['niter'] = cpt + return geometricBar(weights, UKv), log else: - return geometricBar(weights,UKv) + return geometricBar(weights, UKv) -def unmix(a,D,M,M0,h0,reg,reg0,alpha,numItermax = 1000, stopThr=1e-3,verbose=False,log=False): +def unmix(a, D, M, M0, h0, reg, reg0, alpha, numItermax=1000, stopThr=1e-3, verbose=False, log=False): """ Compute the unmixing of an observation with a given dictionary using Wasserstein distance @@ -942,44 +967,45 @@ def unmix(a,D,M,M0,h0,reg,reg0,alpha,numItermax = 1000, stopThr=1e-3,verbose=Fal """ - #M = M/np.median(M) - K = np.exp(-M/reg) + # M = M/np.median(M) + K = np.exp(-M / reg) - #M0 = M0/np.median(M0) - K0 = np.exp(-M0/reg0) + # M0 = M0/np.median(M0) + K0 = np.exp(-M0 / reg0) old = h0 - err=1 - cpt=0 - #log = {'niter':0, 'all_err':[]} + err = 1 + cpt = 0 + # log = {'niter':0, 'all_err':[]} if log: - log={'err':[]} - - - while (err>stopThr and cpt<numItermax): - K = projC(K,a) - K0 = projC(K0,h0) - new = np.sum(K0,axis=1) - inv_new = np.dot(D,new) # we recombine the current selection from dictionnary - other = np.sum(K,axis=1) - delta = np.exp(alpha*np.log(other)+(1-alpha)*np.log(inv_new)) # geometric interpolation - K = projR(K,delta) - K0 = np.dot(np.diag(np.dot(D.T,delta/inv_new)),K0) - - err=np.linalg.norm(np.sum(K0,axis=1)-old) + log = {'err': []} + + while (err > stopThr and cpt < numItermax): + K = projC(K, a) + K0 = projC(K0, h0) + new = np.sum(K0, axis=1) + # we recombine the current selection from dictionnary + inv_new = np.dot(D, new) + other = np.sum(K, axis=1) + # geometric interpolation + delta = np.exp(alpha * np.log(other) + (1 - alpha) * np.log(inv_new)) + K = projR(K, delta) + K0 = np.dot(np.diag(np.dot(D.T, delta / inv_new)), K0) + + err = np.linalg.norm(np.sum(K0, axis=1) - old) old = new 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)) + if cpt % 200 == 0: + print('{:5s}|{:12s}'.format('It.', 'Err') + '\n' + '-' * 19) + print('{:5d}|{:8e}|'.format(cpt, err)) - cpt = cpt+1 + cpt = cpt + 1 if log: - log['niter']=cpt - return np.sum(K0,axis=1),log + log['niter'] = cpt + return np.sum(K0, axis=1), log else: - return np.sum(K0,axis=1) + return np.sum(K0, axis=1) @@ -3,22 +3,34 @@ Domain adaptation with optimal transport """ +# Author: Remi Flamary <remi.flamary@unice.fr> +# Nicolas Courty <ncourty@irisa.fr> +# Michael Perrot <michael.perrot@univ-st-etienne.fr> +# +# License: MIT License + import numpy as np + from .bregman import sinkhorn from .lp import emd -from .utils import unif,dist,kernel +from .utils import unif, dist, kernel, cost_normalization +from .utils import check_params, deprecated, BaseEstimator from .optim import cg from .optim import gcg -def sinkhorn_lpl1_mm(a,labels_a, b, M, reg, eta=0.1,numItermax = 10,numInnerItermax = 200,stopInnerThr=1e-9,verbose=False,log=False): +def sinkhorn_lpl1_mm(a, labels_a, b, M, reg, eta=0.1, numItermax=10, + numInnerItermax=200, stopInnerThr=1e-9, verbose=False, + log=False): """ - Solve the entropic regularization optimal transport problem with nonconvex group lasso regularization + Solve the entropic regularization optimal transport problem with nonconvex + group lasso regularization The function solves the following optimization problem: .. math:: - \gamma = arg\min_\gamma <\gamma,M>_F + reg\cdot\Omega_e(\gamma)+ \eta \Omega_g(\gamma) + \gamma = arg\min_\gamma <\gamma,M>_F + reg\cdot\Omega_e(\gamma) + + \eta \Omega_g(\gamma) s.t. \gamma 1 = a @@ -28,11 +40,16 @@ def sinkhorn_lpl1_mm(a,labels_a, b, M, reg, eta=0.1,numItermax = 10,numInnerIter where : - M is the (ns,nt) metric cost matrix - - :math:`\Omega_e` is the entropic regularization term :math:`\Omega_e(\gamma)=\sum_{i,j} \gamma_{i,j}\log(\gamma_{i,j})` - - :math:`\Omega_g` is the group lasso regulaization term :math:`\Omega_g(\gamma)=\sum_{i,c} \|\gamma_{i,\mathcal{I}_c}\|^{1/2}_1` where :math:`\mathcal{I}_c` are the index of samples from class c in the source domain. + - :math:`\Omega_e` is the entropic regularization term + :math:`\Omega_e(\gamma)=\sum_{i,j} \gamma_{i,j}\log(\gamma_{i,j})` + - :math:`\Omega_g` is the group lasso regulaization term + :math:`\Omega_g(\gamma)=\sum_{i,c} \|\gamma_{i,\mathcal{I}_c}\|^{1/2}_1` + where :math:`\mathcal{I}_c` are the index of samples from class c + in the source domain. - a and b are source and target weights (sum to 1) - The algorithm used for solving the problem is the generalised conditional gradient as proposed in [5]_ [7]_ + The algorithm used for solving the problem is the generalised conditional + gradient as proposed in [5]_ [7]_ Parameters @@ -72,8 +89,13 @@ def sinkhorn_lpl1_mm(a,labels_a, b, M, reg, eta=0.1,numItermax = 10,numInnerIter References ---------- - .. [5] N. Courty; R. Flamary; D. Tuia; A. Rakotomamonjy, "Optimal Transport for Domain Adaptation," in IEEE Transactions on Pattern Analysis and Machine Intelligence , vol.PP, no.99, pp.1-1 - .. [7] Rakotomamonjy, A., Flamary, R., & Courty, N. (2015). Generalized conditional gradient: analysis of convergence and applications. arXiv preprint arXiv:1510.06567. + .. [5] N. Courty; R. Flamary; D. Tuia; A. Rakotomamonjy, + "Optimal Transport for Domain Adaptation," in IEEE + Transactions on Pattern Analysis and Machine Intelligence , + vol.PP, no.99, pp.1-1 + .. [7] Rakotomamonjy, A., Flamary, R., & Courty, N. (2015). + Generalized conditional gradient: analysis of convergence + and applications. arXiv preprint arXiv:1510.06567. See Also -------- @@ -94,7 +116,7 @@ def sinkhorn_lpl1_mm(a,labels_a, b, M, reg, eta=0.1,numItermax = 10,numInnerIter W = np.zeros(M.shape) for cpt in range(numItermax): - Mreg = M + eta*W + Mreg = M + eta * W transp = sinkhorn(a, b, Mreg, reg, numItermax=numInnerItermax, stopThr=stopInnerThr) # the transport has been computed. Check if classes are really @@ -102,19 +124,24 @@ def sinkhorn_lpl1_mm(a,labels_a, b, M, reg, eta=0.1,numItermax = 10,numInnerIter W = np.ones(M.shape) for (i, c) in enumerate(classes): majs = np.sum(transp[indices_labels[i]], axis=0) - majs = p*((majs+epsilon)**(p-1)) + majs = p * ((majs + epsilon)**(p - 1)) W[indices_labels[i]] = majs return transp -def sinkhorn_l1l2_gl(a,labels_a, b, M, reg, eta=0.1,numItermax = 10,numInnerItermax = 200,stopInnerThr=1e-9,verbose=False,log=False): + +def sinkhorn_l1l2_gl(a, labels_a, b, M, reg, eta=0.1, numItermax=10, + numInnerItermax=200, stopInnerThr=1e-9, verbose=False, + log=False): """ - Solve the entropic regularization optimal transport problem with group lasso regularization + Solve the entropic regularization optimal transport problem with group + lasso regularization The function solves the following optimization problem: .. math:: - \gamma = arg\min_\gamma <\gamma,M>_F + reg\cdot\Omega_e(\gamma)+ \eta \Omega_g(\gamma) + \gamma = arg\min_\gamma <\gamma,M>_F + reg\cdot\Omega_e(\gamma)+ + \eta \Omega_g(\gamma) s.t. \gamma 1 = a @@ -124,11 +151,16 @@ def sinkhorn_l1l2_gl(a,labels_a, b, M, reg, eta=0.1,numItermax = 10,numInnerIter where : - M is the (ns,nt) metric cost matrix - - :math:`\Omega_e` is the entropic regularization term :math:`\Omega_e(\gamma)=\sum_{i,j} \gamma_{i,j}\log(\gamma_{i,j})` - - :math:`\Omega_g` is the group lasso regulaization term :math:`\Omega_g(\gamma)=\sum_{i,c} \|\gamma_{i,\mathcal{I}_c}\|^2` where :math:`\mathcal{I}_c` are the index of samples from class c in the source domain. + - :math:`\Omega_e` is the entropic regularization term + :math:`\Omega_e(\gamma)=\sum_{i,j} \gamma_{i,j}\log(\gamma_{i,j})` + - :math:`\Omega_g` is the group lasso regulaization term + :math:`\Omega_g(\gamma)=\sum_{i,c} \|\gamma_{i,\mathcal{I}_c}\|^2` + where :math:`\mathcal{I}_c` are the index of samples from class + c in the source domain. - a and b are source and target weights (sum to 1) - The algorithm used for solving the problem is the generalised conditional gradient as proposed in [5]_ [7]_ + The algorithm used for solving the problem is the generalised conditional + gradient as proposed in [5]_ [7]_ Parameters @@ -168,46 +200,54 @@ def sinkhorn_l1l2_gl(a,labels_a, b, M, reg, eta=0.1,numItermax = 10,numInnerIter References ---------- - .. [5] N. Courty; R. Flamary; D. Tuia; A. Rakotomamonjy, "Optimal Transport for Domain Adaptation," in IEEE Transactions on Pattern Analysis and Machine Intelligence , vol.PP, no.99, pp.1-1 - .. [7] Rakotomamonjy, A., Flamary, R., & Courty, N. (2015). Generalized conditional gradient: analysis of convergence and applications. arXiv preprint arXiv:1510.06567. + .. [5] N. Courty; R. Flamary; D. Tuia; A. Rakotomamonjy, + "Optimal Transport for Domain Adaptation," in IEEE Transactions + on Pattern Analysis and Machine Intelligence , vol.PP, no.99, pp.1-1 + .. [7] Rakotomamonjy, A., Flamary, R., & Courty, N. (2015). + Generalized conditional gradient: analysis of convergence and + applications. arXiv preprint arXiv:1510.06567. See Also -------- ot.optim.gcg : Generalized conditional gradient for OT problems """ - lstlab=np.unique(labels_a) + lstlab = np.unique(labels_a) def f(G): - res=0 + res = 0 for i in range(G.shape[1]): for lab in lstlab: - temp=G[labels_a==lab,i] - res+=np.linalg.norm(temp) + temp = G[labels_a == lab, i] + res += np.linalg.norm(temp) return res def df(G): - W=np.zeros(G.shape) + W = np.zeros(G.shape) for i in range(G.shape[1]): for lab in lstlab: - temp=G[labels_a==lab,i] - n=np.linalg.norm(temp) + temp = G[labels_a == lab, i] + n = np.linalg.norm(temp) if n: - W[labels_a==lab,i]=temp/n + W[labels_a == lab, i] = temp / n return W - - return gcg(a,b,M,reg,eta,f,df,G0=None,numItermax = numItermax,numInnerItermax=numInnerItermax, stopThr=stopInnerThr,verbose=verbose,log=log) + return gcg(a, b, M, reg, eta, f, df, G0=None, numItermax=numItermax, + numInnerItermax=numInnerItermax, stopThr=stopInnerThr, + verbose=verbose, log=log) - -def joint_OT_mapping_linear(xs,xt,mu=1,eta=0.001,bias=False,verbose=False,verbose2=False,numItermax = 100,numInnerItermax = 10,stopInnerThr=1e-6,stopThr=1e-5,log=False,**kwargs): +def joint_OT_mapping_linear(xs, xt, mu=1, eta=0.001, bias=False, verbose=False, + verbose2=False, numItermax=100, numInnerItermax=10, + stopInnerThr=1e-6, stopThr=1e-5, log=False, + **kwargs): """Joint OT and linear mapping estimation as proposed in [8] The function solves the following optimization problem: .. math:: - \min_{\gamma,L}\quad \|L(X_s) -n_s\gamma X_t\|^2_F + \mu<\gamma,M>_F + \eta \|L -I\|^2_F + \min_{\gamma,L}\quad \|L(X_s) -n_s\gamma X_t\|^2_F + + \mu<\gamma,M>_F + \eta \|L -I\|^2_F s.t. \gamma 1 = a @@ -216,8 +256,10 @@ def joint_OT_mapping_linear(xs,xt,mu=1,eta=0.001,bias=False,verbose=False,verbos \gamma\geq 0 where : - - M is the (ns,nt) squared euclidean cost matrix between samples in Xs and Xt (scaled by ns) - - :math:`L` is a dxd linear operator that approximates the barycentric mapping + - M is the (ns,nt) squared euclidean cost matrix between samples in + Xs and Xt (scaled by ns) + - :math:`L` is a dxd linear operator that approximates the barycentric + mapping - :math:`I` is the identity matrix (neutral linear mapping) - a and b are uniform source and target weights @@ -272,7 +314,9 @@ def joint_OT_mapping_linear(xs,xt,mu=1,eta=0.001,bias=False,verbose=False,verbos References ---------- - .. [8] M. Perrot, N. Courty, R. Flamary, A. Habrard, "Mapping estimation for discrete optimal transport", Neural Information Processing Systems (NIPS), 2016. + .. [8] M. Perrot, N. Courty, R. Flamary, A. Habrard, + "Mapping estimation for discrete optimal transport", + Neural Information Processing Systems (NIPS), 2016. See Also -------- @@ -281,103 +325,116 @@ def joint_OT_mapping_linear(xs,xt,mu=1,eta=0.001,bias=False,verbose=False,verbos """ - ns,nt,d=xs.shape[0],xt.shape[0],xt.shape[1] + ns, nt, d = xs.shape[0], xt.shape[0], xt.shape[1] if bias: - xs1=np.hstack((xs,np.ones((ns,1)))) - xstxs=xs1.T.dot(xs1) - I=np.eye(d+1) - I[-1]=0 - I0=I[:,:-1] - sel=lambda x : x[:-1,:] + xs1 = np.hstack((xs, np.ones((ns, 1)))) + xstxs = xs1.T.dot(xs1) + I = np.eye(d + 1) + I[-1] = 0 + I0 = I[:, :-1] + + def sel(x): + return x[:-1, :] else: - xs1=xs - xstxs=xs1.T.dot(xs1) - I=np.eye(d) - I0=I - sel=lambda x : x + xs1 = xs + xstxs = xs1.T.dot(xs1) + I = np.eye(d) + I0 = I + + def sel(x): + return x if log: - log={'err':[]} + log = {'err': []} - a,b=unif(ns),unif(nt) - M=dist(xs,xt)*ns - G=emd(a,b,M) + a, b = unif(ns), unif(nt) + M = dist(xs, xt) * ns + G = emd(a, b, M) - vloss=[] + vloss = [] - def loss(L,G): + def loss(L, G): """Compute full loss""" - return np.sum((xs1.dot(L)-ns*G.dot(xt))**2)+mu*np.sum(G*M)+eta*np.sum(sel(L-I0)**2) + return np.sum((xs1.dot(L) - ns * G.dot(xt))**2) + mu * np.sum(G * M) + eta * np.sum(sel(L - I0)**2) def solve_L(G): """ solve L problem with fixed G (least square)""" - xst=ns*G.dot(xt) - return np.linalg.solve(xstxs+eta*I,xs1.T.dot(xst)+eta*I0) + xst = ns * G.dot(xt) + return np.linalg.solve(xstxs + eta * I, xs1.T.dot(xst) + eta * I0) - def solve_G(L,G0): + def solve_G(L, G0): """Update G with CG algorithm""" - xsi=xs1.dot(L) + xsi = xs1.dot(L) + def f(G): - return np.sum((xsi-ns*G.dot(xt))**2) + return np.sum((xsi - ns * G.dot(xt))**2) + def df(G): - return -2*ns*(xsi-ns*G.dot(xt)).dot(xt.T) - G=cg(a,b,M,1.0/mu,f,df,G0=G0,numItermax=numInnerItermax,stopThr=stopInnerThr) + return -2 * ns * (xsi - ns * G.dot(xt)).dot(xt.T) + G = cg(a, b, M, 1.0 / mu, f, df, G0=G0, + numItermax=numInnerItermax, stopThr=stopInnerThr) return G + L = solve_L(G) - L=solve_L(G) - - vloss.append(loss(L,G)) + vloss.append(loss(L, G)) if verbose: - print('{:5s}|{:12s}|{:8s}'.format('It.','Loss','Delta loss')+'\n'+'-'*32) - print('{:5d}|{:8e}|{:8e}'.format(0,vloss[-1],0)) - + print('{:5s}|{:12s}|{:8s}'.format( + 'It.', 'Loss', 'Delta loss') + '\n' + '-' * 32) + print('{:5d}|{:8e}|{:8e}'.format(0, vloss[-1], 0)) # init loop - if numItermax>0: - loop=1 + if numItermax > 0: + loop = 1 else: - loop=0 - it=0 + loop = 0 + it = 0 while loop: - it+=1 + it += 1 # update G - G=solve_G(L,G) + G = solve_G(L, G) - #update L - L=solve_L(G) + # update L + L = solve_L(G) - vloss.append(loss(L,G)) + vloss.append(loss(L, G)) - if it>=numItermax: - loop=0 + if it >= numItermax: + loop = 0 - if abs(vloss[-1]-vloss[-2])/abs(vloss[-2])<stopThr: - loop=0 + if abs(vloss[-1] - vloss[-2]) / abs(vloss[-2]) < stopThr: + loop = 0 if verbose: - if it%20==0: - print('{:5s}|{:12s}|{:8s}'.format('It.','Loss','Delta loss')+'\n'+'-'*32) - print('{:5d}|{:8e}|{:8e}'.format(it,vloss[-1],(vloss[-1]-vloss[-2])/abs(vloss[-2]))) + if it % 20 == 0: + print('{:5s}|{:12s}|{:8s}'.format( + 'It.', 'Loss', 'Delta loss') + '\n' + '-' * 32) + print('{:5d}|{:8e}|{:8e}'.format( + it, vloss[-1], (vloss[-1] - vloss[-2]) / abs(vloss[-2]))) if log: - log['loss']=vloss - return G,L,log + log['loss'] = vloss + return G, L, log else: - return G,L + return G, L -def joint_OT_mapping_kernel(xs,xt,mu=1,eta=0.001,kerneltype='gaussian',sigma=1,bias=False,verbose=False,verbose2=False,numItermax = 100,numInnerItermax = 10,stopInnerThr=1e-6,stopThr=1e-5,log=False,**kwargs): +def joint_OT_mapping_kernel(xs, xt, mu=1, eta=0.001, kerneltype='gaussian', + sigma=1, bias=False, verbose=False, verbose2=False, + numItermax=100, numInnerItermax=10, + stopInnerThr=1e-6, stopThr=1e-5, log=False, + **kwargs): """Joint OT and nonlinear mapping estimation with kernels as proposed in [8] The function solves the following optimization problem: .. math:: - \min_{\gamma,L\in\mathcal{H}}\quad \|L(X_s) -n_s\gamma X_t\|^2_F + \mu<\gamma,M>_F + \eta \|L\|^2_\mathcal{H} + \min_{\gamma,L\in\mathcal{H}}\quad \|L(X_s) - + n_s\gamma X_t\|^2_F + \mu<\gamma,M>_F + \eta \|L\|^2_\mathcal{H} s.t. \gamma 1 = a @@ -386,8 +443,10 @@ def joint_OT_mapping_kernel(xs,xt,mu=1,eta=0.001,kerneltype='gaussian',sigma=1,b \gamma\geq 0 where : - - M is the (ns,nt) squared euclidean cost matrix between samples in Xs and Xt (scaled by ns) - - :math:`L` is a ns x d linear operator on a kernel matrix that approximates the barycentric mapping + - M is the (ns,nt) squared euclidean cost matrix between samples in + Xs and Xt (scaled by ns) + - :math:`L` is a ns x d linear operator on a kernel matrix that + approximates the barycentric mapping - a and b are uniform source and target weights The problem consist in solving jointly an optimal transport matrix @@ -445,7 +504,9 @@ def joint_OT_mapping_kernel(xs,xt,mu=1,eta=0.001,kerneltype='gaussian',sigma=1,b References ---------- - .. [8] M. Perrot, N. Courty, R. Flamary, A. Habrard, "Mapping estimation for discrete optimal transport", Neural Information Processing Systems (NIPS), 2016. + .. [8] M. Perrot, N. Courty, R. Flamary, A. Habrard, + "Mapping estimation for discrete optimal transport", + Neural Information Processing Systems (NIPS), 2016. See Also -------- @@ -454,161 +515,170 @@ def joint_OT_mapping_kernel(xs,xt,mu=1,eta=0.001,kerneltype='gaussian',sigma=1,b """ - ns,nt=xs.shape[0],xt.shape[0] + ns, nt = xs.shape[0], xt.shape[0] - K=kernel(xs,xs,method=kerneltype,sigma=sigma) + K = kernel(xs, xs, method=kerneltype, sigma=sigma) if bias: - K1=np.hstack((K,np.ones((ns,1)))) - I=np.eye(ns+1) - I[-1]=0 - Kp=np.eye(ns+1) - Kp[:ns,:ns]=K + K1 = np.hstack((K, np.ones((ns, 1)))) + I = np.eye(ns + 1) + I[-1] = 0 + Kp = np.eye(ns + 1) + Kp[:ns, :ns] = K # ls regu - #K0 = K1.T.dot(K1)+eta*I - #Kreg=I + # K0 = K1.T.dot(K1)+eta*I + # Kreg=I # RKHS regul - K0 = K1.T.dot(K1)+eta*Kp - Kreg=Kp + K0 = K1.T.dot(K1) + eta * Kp + Kreg = Kp else: - K1=K - I=np.eye(ns) + K1 = K + I = np.eye(ns) # ls regul - #K0 = K1.T.dot(K1)+eta*I - #Kreg=I + # K0 = K1.T.dot(K1)+eta*I + # Kreg=I # proper kernel ridge - K0=K+eta*I - Kreg=K - - - + K0 = K + eta * I + Kreg = K if log: - log={'err':[]} + log = {'err': []} - a,b=unif(ns),unif(nt) - M=dist(xs,xt)*ns - G=emd(a,b,M) + a, b = unif(ns), unif(nt) + M = dist(xs, xt) * ns + G = emd(a, b, M) - vloss=[] + vloss = [] - def loss(L,G): + def loss(L, G): """Compute full loss""" - return np.sum((K1.dot(L)-ns*G.dot(xt))**2)+mu*np.sum(G*M)+eta*np.trace(L.T.dot(Kreg).dot(L)) + return np.sum((K1.dot(L) - ns * G.dot(xt))**2) + mu * np.sum(G * M) + eta * np.trace(L.T.dot(Kreg).dot(L)) def solve_L_nobias(G): """ solve L problem with fixed G (least square)""" - xst=ns*G.dot(xt) - return np.linalg.solve(K0,xst) + xst = ns * G.dot(xt) + return np.linalg.solve(K0, xst) def solve_L_bias(G): """ solve L problem with fixed G (least square)""" - xst=ns*G.dot(xt) - return np.linalg.solve(K0,K1.T.dot(xst)) + xst = ns * G.dot(xt) + return np.linalg.solve(K0, K1.T.dot(xst)) - def solve_G(L,G0): + def solve_G(L, G0): """Update G with CG algorithm""" - xsi=K1.dot(L) + xsi = K1.dot(L) + def f(G): - return np.sum((xsi-ns*G.dot(xt))**2) + return np.sum((xsi - ns * G.dot(xt))**2) + def df(G): - return -2*ns*(xsi-ns*G.dot(xt)).dot(xt.T) - G=cg(a,b,M,1.0/mu,f,df,G0=G0,numItermax=numInnerItermax,stopThr=stopInnerThr) + return -2 * ns * (xsi - ns * G.dot(xt)).dot(xt.T) + G = cg(a, b, M, 1.0 / mu, f, df, G0=G0, + numItermax=numInnerItermax, stopThr=stopInnerThr) return G if bias: - solve_L=solve_L_bias + solve_L = solve_L_bias else: - solve_L=solve_L_nobias + solve_L = solve_L_nobias - L=solve_L(G) + L = solve_L(G) - vloss.append(loss(L,G)) + vloss.append(loss(L, G)) if verbose: - print('{:5s}|{:12s}|{:8s}'.format('It.','Loss','Delta loss')+'\n'+'-'*32) - print('{:5d}|{:8e}|{:8e}'.format(0,vloss[-1],0)) - + print('{:5s}|{:12s}|{:8s}'.format( + 'It.', 'Loss', 'Delta loss') + '\n' + '-' * 32) + print('{:5d}|{:8e}|{:8e}'.format(0, vloss[-1], 0)) # init loop - if numItermax>0: - loop=1 + if numItermax > 0: + loop = 1 else: - loop=0 - it=0 + loop = 0 + it = 0 while loop: - it+=1 + it += 1 # update G - G=solve_G(L,G) + G = solve_G(L, G) - #update L - L=solve_L(G) + # update L + L = solve_L(G) - vloss.append(loss(L,G)) + vloss.append(loss(L, G)) - if it>=numItermax: - loop=0 + if it >= numItermax: + loop = 0 - if abs(vloss[-1]-vloss[-2])/abs(vloss[-2])<stopThr: - loop=0 + if abs(vloss[-1] - vloss[-2]) / abs(vloss[-2]) < stopThr: + loop = 0 if verbose: - if it%20==0: - print('{:5s}|{:12s}|{:8s}'.format('It.','Loss','Delta loss')+'\n'+'-'*32) - print('{:5d}|{:8e}|{:8e}'.format(it,vloss[-1],(vloss[-1]-vloss[-2])/abs(vloss[-2]))) + if it % 20 == 0: + print('{:5s}|{:12s}|{:8s}'.format( + 'It.', 'Loss', 'Delta loss') + '\n' + '-' * 32) + print('{:5d}|{:8e}|{:8e}'.format( + it, vloss[-1], (vloss[-1] - vloss[-2]) / abs(vloss[-2]))) if log: - log['loss']=vloss - return G,L,log + log['loss'] = vloss + return G, L, log else: - return G,L + return G, L +@deprecated("The class OTDA is deprecated in 0.3.1 and will be " + "removed in 0.5" + "\n\tfor standard transport use class EMDTransport instead.") class OTDA(object): + """Class for domain adaptation with optimal transport as proposed in [5] References ---------- - .. [5] N. Courty; R. Flamary; D. Tuia; A. Rakotomamonjy, "Optimal Transport for Domain Adaptation," in IEEE Transactions on Pattern Analysis and Machine Intelligence , vol.PP, no.99, pp.1-1 + .. [5] N. Courty; R. Flamary; D. Tuia; A. Rakotomamonjy, + "Optimal Transport for Domain Adaptation," in IEEE Transactions on + Pattern Analysis and Machine Intelligence , vol.PP, no.99, pp.1-1 """ - def __init__(self,metric='sqeuclidean'): + def __init__(self, metric='sqeuclidean', norm=None): """ Class initialization""" - self.xs=0 - self.xt=0 - self.G=0 - self.metric=metric - self.computed=False - - - def fit(self,xs,xt,ws=None,wt=None,norm=None): - """ Fit domain adaptation between samples is xs and xt (with optional weights)""" - self.xs=xs - self.xt=xt + self.xs = 0 + self.xt = 0 + self.G = 0 + self.metric = metric + self.norm = norm + self.computed = False + + def fit(self, xs, xt, ws=None, wt=None, max_iter=100000): + """Fit domain adaptation between samples is xs and xt + (with optional weights)""" + self.xs = xs + self.xt = xt if wt is None: - wt=unif(xt.shape[0]) + wt = unif(xt.shape[0]) if ws is None: - ws=unif(xs.shape[0]) + ws = unif(xs.shape[0]) - self.ws=ws - self.wt=wt + self.ws = ws + self.wt = wt - self.M=dist(xs,xt,metric=self.metric) - self.normalizeM(norm) - self.G=emd(ws,wt,self.M) - self.computed=True + self.M = dist(xs, xt, metric=self.metric) + self.M = cost_normalization(self.M, self.norm) + self.G = emd(ws, wt, self.M, max_iter) + self.computed = True - def interp(self,direction=1): + def interp(self, direction=1): """Barycentric interpolation for the source (1) or target (-1) samples This Barycentric interpolation solves for each source (resp target) @@ -623,28 +693,28 @@ class OTDA(object): metric could be used in the future. """ - if direction>0: # >0 then source to target - G=self.G - w=self.ws.reshape((self.xs.shape[0],1)) - x=self.xt + if direction > 0: # >0 then source to target + G = self.G + w = self.ws.reshape((self.xs.shape[0], 1)) + x = self.xt else: - G=self.G.T - w=self.wt.reshape((self.xt.shape[0],1)) - x=self.xs + G = self.G.T + w = self.wt.reshape((self.xt.shape[0], 1)) + x = self.xs if self.computed: - if self.metric=='sqeuclidean': - return np.dot(G/w,x) # weighted mean + if self.metric == 'sqeuclidean': + return np.dot(G / w, x) # weighted mean else: - print("Warning, metric not handled yet, using weighted average") - return np.dot(G/w,x) # weighted mean + print( + "Warning, metric not handled yet, using weighted average") + return np.dot(G / w, x) # weighted mean return None else: print("Warning, model not fitted yet, returning None") return None - - def predict(self,x,direction=1): + def predict(self, x, direction=1): """ Out of sample mapping using the formulation from [6] For each sample x to map, it finds the nearest source sample xs and @@ -654,183 +724,1019 @@ class OTDA(object): References ---------- - .. [6] Ferradans, S., Papadakis, N., Peyré, G., & Aujol, J. F. (2014). Regularized discrete optimal transport. SIAM Journal on Imaging Sciences, 7(3), 1853-1882. + .. [6] Ferradans, S., Papadakis, N., Peyré, G., & Aujol, J. F. (2014). + Regularized discrete optimal transport. SIAM Journal on Imaging + Sciences, 7(3), 1853-1882. """ - if direction>0: # >0 then source to target - xf=self.xt - x0=self.xs + if direction > 0: # >0 then source to target + xf = self.xt + x0 = self.xs else: - xf=self.xs - x0=self.xt - - D0=dist(x,x0) # dist netween new samples an source - idx=np.argmin(D0,1) # closest one - xf=self.interp(direction)# interp the source samples - return xf[idx,:]+x-x0[idx,:] # aply the delta to the interpolation - - def normalizeM(self, norm): - """ Apply normalization to the loss matrix - - - Parameters - ---------- - norm : str - type of normalization from 'median','max','log','loglog' - - """ - - if norm == "median": - self.M /= float(np.median(self.M)) - elif norm == "max": - self.M /= float(np.max(self.M)) - elif norm == "log": - self.M = np.log(1 + self.M) - elif norm == "loglog": - self.M = np.log(1 + np.log(1 + self.M)) - + xf = self.xs + x0 = self.xt + D0 = dist(x, x0) # dist netween new samples an source + idx = np.argmin(D0, 1) # closest one + xf = self.interp(direction) # interp the source samples + # aply the delta to the interpolation + return xf[idx, :] + x - x0[idx, :] + +@deprecated("The class OTDA_sinkhorn is deprecated in 0.3.1 and will be" + " removed in 0.5 \nUse class SinkhornTransport instead.") class OTDA_sinkhorn(OTDA): - """Class for domain adaptation with optimal transport with entropic regularization""" - def fit(self,xs,xt,reg=1,ws=None,wt=None,norm=None,**kwargs): - """ Fit regularized domain adaptation between samples is xs and xt (with optional weights)""" - self.xs=xs - self.xt=xt + """Class for domain adaptation with optimal transport with entropic + regularization + + + """ + + def fit(self, xs, xt, reg=1, ws=None, wt=None, **kwargs): + """Fit regularized domain adaptation between samples is xs and xt + (with optional weights)""" + self.xs = xs + self.xt = xt if wt is None: - wt=unif(xt.shape[0]) + wt = unif(xt.shape[0]) if ws is None: - ws=unif(xs.shape[0]) + ws = unif(xs.shape[0]) - self.ws=ws - self.wt=wt + self.ws = ws + self.wt = wt - self.M=dist(xs,xt,metric=self.metric) - self.normalizeM(norm) - self.G=sinkhorn(ws,wt,self.M,reg,**kwargs) - self.computed=True + self.M = dist(xs, xt, metric=self.metric) + self.M = cost_normalization(self.M, self.norm) + self.G = sinkhorn(ws, wt, self.M, reg, **kwargs) + self.computed = True +@deprecated("The class OTDA_lpl1 is deprecated in 0.3.1 and will be" + " removed in 0.5 \nUse class SinkhornLpl1Transport instead.") class OTDA_lpl1(OTDA): - """Class for domain adaptation with optimal transport with entropic and group regularization""" + """Class for domain adaptation with optimal transport with entropic and + group regularization""" - def fit(self,xs,ys,xt,reg=1,eta=1,ws=None,wt=None,norm=None,**kwargs): - """ Fit regularized domain adaptation between samples is xs and xt (with optional weights), See ot.da.sinkhorn_lpl1_mm for fit parameters""" - self.xs=xs - self.xt=xt + def fit(self, xs, ys, xt, reg=1, eta=1, ws=None, wt=None, **kwargs): + """Fit regularized domain adaptation between samples is xs and xt + (with optional weights), See ot.da.sinkhorn_lpl1_mm for fit + parameters""" + self.xs = xs + self.xt = xt if wt is None: - wt=unif(xt.shape[0]) + wt = unif(xt.shape[0]) if ws is None: - ws=unif(xs.shape[0]) + ws = unif(xs.shape[0]) + + self.ws = ws + self.wt = wt - self.ws=ws - self.wt=wt + self.M = dist(xs, xt, metric=self.metric) + self.M = cost_normalization(self.M, self.norm) + self.G = sinkhorn_lpl1_mm(ws, ys, wt, self.M, reg, eta, **kwargs) + self.computed = True - self.M=dist(xs,xt,metric=self.metric) - self.normalizeM(norm) - self.G=sinkhorn_lpl1_mm(ws,ys,wt,self.M,reg,eta,**kwargs) - self.computed=True +@deprecated("The class OTDA_l1L2 is deprecated in 0.3.1 and will be" + " removed in 0.5 \nUse class SinkhornL1l2Transport instead.") class OTDA_l1l2(OTDA): - """Class for domain adaptation with optimal transport with entropic and group lasso regularization""" + """Class for domain adaptation with optimal transport with entropic + and group lasso regularization""" - def fit(self,xs,ys,xt,reg=1,eta=1,ws=None,wt=None,norm=None,**kwargs): - """ Fit regularized domain adaptation between samples is xs and xt (with optional weights), See ot.da.sinkhorn_lpl1_gl for fit parameters""" - self.xs=xs - self.xt=xt + def fit(self, xs, ys, xt, reg=1, eta=1, ws=None, wt=None, **kwargs): + """Fit regularized domain adaptation between samples is xs and xt + (with optional weights), See ot.da.sinkhorn_lpl1_gl for fit + parameters""" + self.xs = xs + self.xt = xt if wt is None: - wt=unif(xt.shape[0]) + wt = unif(xt.shape[0]) if ws is None: - ws=unif(xs.shape[0]) + ws = unif(xs.shape[0]) + + self.ws = ws + self.wt = wt - self.ws=ws - self.wt=wt + self.M = dist(xs, xt, metric=self.metric) + self.M = cost_normalization(self.M, self.norm) + self.G = sinkhorn_l1l2_gl(ws, ys, wt, self.M, reg, eta, **kwargs) + self.computed = True - self.M=dist(xs,xt,metric=self.metric) - self.normalizeM(norm) - self.G=sinkhorn_l1l2_gl(ws,ys,wt,self.M,reg,eta,**kwargs) - self.computed=True +@deprecated("The class OTDA_mapping_linear is deprecated in 0.3.1 and will be" + " removed in 0.5 \nUse class MappingTransport instead.") class OTDA_mapping_linear(OTDA): - """Class for optimal transport with joint linear mapping estimation as in [8]""" + """Class for optimal transport with joint linear mapping estimation as in + [8] + """ def __init__(self): """ Class initialization""" + self.xs = 0 + self.xt = 0 + self.G = 0 + self.L = 0 + self.bias = False + self.computed = False + self.metric = 'sqeuclidean' - self.xs=0 - self.xt=0 - self.G=0 - self.L=0 - self.bias=False - self.computed=False - self.metric='sqeuclidean' - - def fit(self,xs,xt,mu=1,eta=1,bias=False,**kwargs): + def fit(self, xs, xt, mu=1, eta=1, bias=False, **kwargs): """ Fit domain adaptation between samples is xs and xt (with optional weights)""" - self.xs=xs - self.xt=xt - self.bias=bias + self.xs = xs + self.xt = xt + self.bias = bias + self.ws = unif(xs.shape[0]) + self.wt = unif(xt.shape[0]) - self.ws=unif(xs.shape[0]) - self.wt=unif(xt.shape[0]) - - self.G,self.L=joint_OT_mapping_linear(xs,xt,mu=mu,eta=eta,bias=bias,**kwargs) - self.computed=True + self.G, self.L = joint_OT_mapping_linear( + xs, xt, mu=mu, eta=eta, bias=bias, **kwargs) + self.computed = True def mapping(self): return lambda x: self.predict(x) - - def predict(self,x): + def predict(self, x): """ Out of sample mapping estimated during the call to fit""" if self.computed: if self.bias: - x=np.hstack((x,np.ones((x.shape[0],1)))) - return x.dot(self.L) # aply the delta to the interpolation + x = np.hstack((x, np.ones((x.shape[0], 1)))) + return x.dot(self.L) # aply the delta to the interpolation else: print("Warning, model not fitted yet, returning None") return None -class OTDA_mapping_kernel(OTDA_mapping_linear): - """Class for optimal transport with joint nonlinear mapping estimation as in [8]""" +@deprecated("The class OTDA_mapping_kernel is deprecated in 0.3.1 and will be" + " removed in 0.5 \nUse class MappingTransport instead.") +class OTDA_mapping_kernel(OTDA_mapping_linear): + """Class for optimal transport with joint nonlinear mapping + estimation as in [8]""" - def fit(self,xs,xt,mu=1,eta=1,bias=False,kerneltype='gaussian',sigma=1,**kwargs): + def fit(self, xs, xt, mu=1, eta=1, bias=False, kerneltype='gaussian', + sigma=1, **kwargs): """ Fit domain adaptation between samples is xs and xt """ - self.xs=xs - self.xt=xt - self.bias=bias - - self.ws=unif(xs.shape[0]) - self.wt=unif(xt.shape[0]) - self.kernel=kerneltype - self.sigma=sigma - self.kwargs=kwargs - + self.xs = xs + self.xt = xt + self.bias = bias - self.G,self.L=joint_OT_mapping_kernel(xs,xt,mu=mu,eta=eta,bias=bias,**kwargs) - self.computed=True + self.ws = unif(xs.shape[0]) + self.wt = unif(xt.shape[0]) + self.kernel = kerneltype + self.sigma = sigma + self.kwargs = kwargs + self.G, self.L = joint_OT_mapping_kernel( + xs, xt, mu=mu, eta=eta, bias=bias, **kwargs) + self.computed = True - def predict(self,x): + def predict(self, x): """ Out of sample mapping estimated during the call to fit""" if self.computed: - K=kernel(x,self.xs,method=self.kernel,sigma=self.sigma,**self.kwargs) + K = kernel( + x, self.xs, method=self.kernel, sigma=self.sigma, + **self.kwargs) if self.bias: - K=np.hstack((K,np.ones((x.shape[0],1)))) + K = np.hstack((K, np.ones((x.shape[0], 1)))) return K.dot(self.L) else: print("Warning, model not fitted yet, returning None") - return None
\ No newline at end of file + return None + + +def distribution_estimation_uniform(X): + """estimates a uniform distribution from an array of samples X + + Parameters + ---------- + X : array-like, shape (n_samples, n_features) + The array of samples + + Returns + ------- + mu : array-like, shape (n_samples,) + The uniform distribution estimated from X + """ + + return unif(X.shape[0]) + + +class BaseTransport(BaseEstimator): + """Base class for OTDA objects + + Notes + ----- + All estimators should specify all the parameters that can be set + at the class level in their ``__init__`` as explicit keyword + arguments (no ``*args`` or ``**kwargs``). + + fit method should: + - estimate a cost matrix and store it in a `cost_` attribute + - estimate a coupling matrix and store it in a `coupling_` + attribute + - estimate distributions from source and target data and store them in + mu_s and mu_t attributes + - store Xs and Xt in attributes to be used later on in transform and + inverse_transform methods + + transform method should always get as input a Xs parameter + inverse_transform method should always get as input a Xt parameter + """ + + def fit(self, Xs=None, ys=None, Xt=None, yt=None): + """Build a coupling matrix from source and target sets of samples + (Xs, ys) and (Xt, yt) + + Parameters + ---------- + Xs : array-like, shape (n_source_samples, n_features) + The training input samples. + ys : array-like, shape (n_source_samples,) + The class labels + Xt : array-like, shape (n_target_samples, n_features) + The training input samples. + yt : array-like, shape (n_labeled_target_samples,) + The class labels + + Returns + ------- + self : object + Returns self. + """ + + # check the necessary inputs parameters are here + if check_params(Xs=Xs, Xt=Xt): + + # pairwise distance + self.cost_ = dist(Xs, Xt, metric=self.metric) + self.cost_ = cost_normalization(self.cost_, self.norm) + + 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_) + + # 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 + + # 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 + + return self + + def fit_transform(self, Xs=None, ys=None, Xt=None, yt=None): + """Build a coupling matrix from source and target sets of samples + (Xs, ys) and (Xt, yt) and transports source samples Xs onto target + ones Xt + + Parameters + ---------- + Xs : array-like, shape (n_source_samples, n_features) + The training input samples. + ys : array-like, shape (n_source_samples,) + The class labels + Xt : array-like, shape (n_target_samples, n_features) + The training input samples. + yt : array-like, shape (n_labeled_target_samples,) + The class labels + + Returns + ------- + transp_Xs : array-like, shape (n_source_samples, n_features) + The source samples samples. + """ + + return self.fit(Xs, ys, Xt, yt).transform(Xs, ys, Xt, yt) + + def transform(self, Xs=None, ys=None, Xt=None, yt=None, batch_size=128): + """Transports source samples Xs onto target ones Xt + + Parameters + ---------- + Xs : array-like, shape (n_source_samples, n_features) + The training input samples. + ys : array-like, shape (n_source_samples,) + The class labels + Xt : array-like, shape (n_target_samples, n_features) + The training input samples. + yt : array-like, shape (n_labeled_target_samples,) + The class labels + batch_size : int, optional (default=128) + The batch size for out of sample inverse transform + + Returns + ------- + transp_Xs : array-like, shape (n_source_samples, n_features) + The transport source samples. + """ + + # check the necessary inputs parameters are here + if check_params(Xs=Xs): + + 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 + + # 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: + + # 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_) + + # define the transported points + transp_Xs_ = transp_Xs_[idx, :] + Xs[bi] - self.xs_[idx, :] + + transp_Xs.append(transp_Xs_) + + transp_Xs = np.concatenate(transp_Xs, axis=0) + + return transp_Xs + + def inverse_transform(self, Xs=None, ys=None, Xt=None, yt=None, + batch_size=128): + """Transports target samples Xt onto target samples Xs + + Parameters + ---------- + Xs : array-like, shape (n_source_samples, n_features) + The training input samples. + ys : array-like, shape (n_source_samples,) + The class labels + Xt : array-like, shape (n_target_samples, n_features) + The training input samples. + yt : array-like, shape (n_labeled_target_samples,) + The class labels + batch_size : int, optional (default=128) + The batch size for out of sample inverse transform + + Returns + ------- + transp_Xt : array-like, shape (n_source_samples, n_features) + The transported target samples. + """ + + # check the necessary inputs parameters are here + if check_params(Xt=Xt): + + 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 + + # 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: + + 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_) + + # define the transported points + transp_Xt_ = transp_Xt_[idx, :] + Xt[bi] - self.xt_[idx, :] + + transp_Xt.append(transp_Xt_) + + transp_Xt = np.concatenate(transp_Xt, axis=0) + + return transp_Xt + + +class SinkhornTransport(BaseTransport): + """Domain Adapatation OT method based on Sinkhorn Algorithm + + Parameters + ---------- + reg_e : float, optional (default=1) + Entropic regularization parameter + max_iter : int, float, optional (default=1000) + The minimum number of iteration before stopping the optimization + algorithm if no it has not converged + tol : float, optional (default=10e-9) + The precision required to stop the optimization algorithm. + mapping : string, optional (default="barycentric") + The kind of mapping to apply to transport samples from a domain into + another one. + if "barycentric" only the samples used to estimate the coupling can + be transported from a domain to another one. + metric : string, optional (default="sqeuclidean") + The ground metric for the Wasserstein problem + norm : string, optional (default=None) + If given, normalize the ground metric to avoid numerical errors that + can occur with large metric values. + distribution : string, optional (default="uniform") + The kind of distribution estimation to employ + verbose : int, optional (default=0) + Controls the verbosity of the optimization algorithm + log : int, optional (default=0) + Controls the logs of the optimization algorithm + limit_max: float, optional (defaul=np.infty) + Controls the semi supervised mode. Transport between labeled source + and target samples of different classes will exhibit an infinite cost + + Attributes + ---------- + coupling_ : array-like, shape (n_source_samples, n_target_samples) + The optimal coupling + log_ : dictionary + The dictionary of log, empty dic if parameter log is not True + + References + ---------- + .. [1] N. Courty; R. Flamary; D. Tuia; A. Rakotomamonjy, + "Optimal Transport for Domain Adaptation," in IEEE Transactions + on Pattern Analysis and Machine Intelligence , vol.PP, no.99, pp.1-1 + .. [2] M. Cuturi, Sinkhorn Distances : Lightspeed Computation of Optimal + Transport, Advances in Neural Information Processing Systems (NIPS) + 26, 2013 + """ + + def __init__(self, reg_e=1., max_iter=1000, + tol=10e-9, verbose=False, log=False, + metric="sqeuclidean", norm=None, + distribution_estimation=distribution_estimation_uniform, + out_of_sample_map='ferradans', limit_max=np.infty): + + self.reg_e = reg_e + self.max_iter = max_iter + self.tol = tol + self.verbose = verbose + self.log = log + self.metric = metric + self.norm = norm + self.limit_max = limit_max + self.distribution_estimation = distribution_estimation + self.out_of_sample_map = out_of_sample_map + + def fit(self, Xs=None, ys=None, Xt=None, yt=None): + """Build a coupling matrix from source and target sets of samples + (Xs, ys) and (Xt, yt) + + Parameters + ---------- + Xs : array-like, shape (n_source_samples, n_features) + The training input samples. + ys : array-like, shape (n_source_samples,) + The class labels + Xt : array-like, shape (n_target_samples, n_features) + The training input samples. + yt : array-like, shape (n_labeled_target_samples,) + The class labels + + Returns + ------- + self : object + Returns self. + """ + + super(SinkhornTransport, self).fit(Xs, ys, Xt, yt) + + # coupling estimation + returned_ = sinkhorn( + a=self.mu_s, b=self.mu_t, M=self.cost_, reg=self.reg_e, + numItermax=self.max_iter, stopThr=self.tol, + verbose=self.verbose, log=self.log) + + # deal with the value of log + if self.log: + self.coupling_, self.log_ = returned_ + else: + self.coupling_ = returned_ + self.log_ = dict() + + return self + + +class EMDTransport(BaseTransport): + """Domain Adapatation OT method based on Earth Mover's Distance + + Parameters + ---------- + mapping : string, optional (default="barycentric") + The kind of mapping to apply to transport samples from a domain into + another one. + if "barycentric" only the samples used to estimate the coupling can + be transported from a domain to another one. + metric : string, optional (default="sqeuclidean") + The ground metric for the Wasserstein problem + norm : string, optional (default=None) + If given, normalize the ground metric to avoid numerical errors that + can occur with large metric values. + distribution : string, optional (default="uniform") + The kind of distribution estimation to employ + verbose : int, optional (default=0) + Controls the verbosity of the optimization algorithm + log : int, optional (default=0) + Controls the logs of the optimization algorithm + limit_max: float, optional (default=10) + Controls the semi supervised mode. Transport between labeled source + and target samples of different classes will exhibit an infinite cost + (10 times the maximum value of the cost matrix) + max_iter : int, optional (default=100000) + The maximum number of iterations before stopping the optimization + algorithm if it has not converged. + + Attributes + ---------- + coupling_ : array-like, shape (n_source_samples, n_target_samples) + The optimal coupling + + References + ---------- + .. [1] N. Courty; R. Flamary; D. Tuia; A. Rakotomamonjy, + "Optimal Transport for Domain Adaptation," in IEEE Transactions + on Pattern Analysis and Machine Intelligence , vol.PP, no.99, pp.1-1 + """ + + def __init__(self, metric="sqeuclidean", norm=None, + distribution_estimation=distribution_estimation_uniform, + out_of_sample_map='ferradans', limit_max=10, + max_iter=100000): + + self.metric = metric + self.norm = norm + self.limit_max = limit_max + self.distribution_estimation = distribution_estimation + self.out_of_sample_map = out_of_sample_map + self.max_iter = max_iter + + def fit(self, Xs, ys=None, Xt=None, yt=None): + """Build a coupling matrix from source and target sets of samples + (Xs, ys) and (Xt, yt) + + Parameters + ---------- + Xs : array-like, shape (n_source_samples, n_features) + The training input samples. + ys : array-like, shape (n_source_samples,) + The class labels + Xt : array-like, shape (n_target_samples, n_features) + The training input samples. + yt : array-like, shape (n_labeled_target_samples,) + The class labels + + Returns + ------- + self : object + Returns self. + """ + + super(EMDTransport, self).fit(Xs, ys, Xt, yt) + + # coupling estimation + self.coupling_ = emd( + a=self.mu_s, b=self.mu_t, M=self.cost_, numItermax=self.max_iter + ) + + return self + + +class SinkhornLpl1Transport(BaseTransport): + """Domain Adapatation OT method based on sinkhorn algorithm + + LpL1 class regularization. + + Parameters + ---------- + reg_e : float, optional (default=1) + Entropic regularization parameter + reg_cl : float, optional (default=0.1) + Class regularization parameter + mapping : string, optional (default="barycentric") + The kind of mapping to apply to transport samples from a domain into + another one. + if "barycentric" only the samples used to estimate the coupling can + be transported from a domain to another one. + metric : string, optional (default="sqeuclidean") + The ground metric for the Wasserstein problem + norm : string, optional (default=None) + If given, normalize the ground metric to avoid numerical errors that + can occur with large metric values. + distribution : string, optional (default="uniform") + The kind of distribution estimation to employ + max_iter : int, float, optional (default=10) + The minimum number of iteration before stopping the optimization + algorithm if no it has not converged + max_inner_iter : int, float, optional (default=200) + The number of iteration in the inner loop + verbose : int, optional (default=0) + Controls the verbosity of the optimization algorithm + limit_max: float, optional (defaul=np.infty) + Controls the semi supervised mode. Transport between labeled source + and target samples of different classes will exhibit an infinite cost + + Attributes + ---------- + coupling_ : array-like, shape (n_source_samples, n_target_samples) + The optimal coupling + + References + ---------- + + .. [1] N. Courty; R. Flamary; D. Tuia; A. Rakotomamonjy, + "Optimal Transport for Domain Adaptation," in IEEE + Transactions on Pattern Analysis and Machine Intelligence , + vol.PP, no.99, pp.1-1 + .. [2] Rakotomamonjy, A., Flamary, R., & Courty, N. (2015). + Generalized conditional gradient: analysis of convergence + and applications. arXiv preprint arXiv:1510.06567. + + """ + + def __init__(self, reg_e=1., reg_cl=0.1, + max_iter=10, max_inner_iter=200, + tol=10e-9, verbose=False, + metric="sqeuclidean", norm=None, + distribution_estimation=distribution_estimation_uniform, + out_of_sample_map='ferradans', limit_max=np.infty): + + self.reg_e = reg_e + self.reg_cl = reg_cl + self.max_iter = max_iter + self.max_inner_iter = max_inner_iter + self.tol = tol + self.verbose = verbose + self.metric = metric + self.norm = norm + self.distribution_estimation = distribution_estimation + self.out_of_sample_map = out_of_sample_map + self.limit_max = limit_max + + def fit(self, Xs, ys=None, Xt=None, yt=None): + """Build a coupling matrix from source and target sets of samples + (Xs, ys) and (Xt, yt) + + Parameters + ---------- + Xs : array-like, shape (n_source_samples, n_features) + The training input samples. + ys : array-like, shape (n_source_samples,) + The class labels + Xt : array-like, shape (n_target_samples, n_features) + The training input samples. + yt : array-like, shape (n_labeled_target_samples,) + The class labels + + Returns + ------- + self : object + Returns self. + """ + + # check the necessary inputs parameters are here + if check_params(Xs=Xs, Xt=Xt, ys=ys): + + super(SinkhornLpl1Transport, self).fit(Xs, ys, Xt, yt) + + 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 + + +class SinkhornL1l2Transport(BaseTransport): + """Domain Adapatation OT method based on sinkhorn algorithm + + l1l2 class regularization. + + Parameters + ---------- + reg_e : float, optional (default=1) + Entropic regularization parameter + reg_cl : float, optional (default=0.1) + Class regularization parameter + mapping : string, optional (default="barycentric") + The kind of mapping to apply to transport samples from a domain into + another one. + if "barycentric" only the samples used to estimate the coupling can + be transported from a domain to another one. + metric : string, optional (default="sqeuclidean") + The ground metric for the Wasserstein problem + norm : string, optional (default=None) + If given, normalize the ground metric to avoid numerical errors that + can occur with large metric values. + distribution : string, optional (default="uniform") + The kind of distribution estimation to employ + max_iter : int, float, optional (default=10) + The minimum number of iteration before stopping the optimization + algorithm if no it has not converged + max_inner_iter : int, float, optional (default=200) + The number of iteration in the inner loop + verbose : int, optional (default=0) + Controls the verbosity of the optimization algorithm + log : int, optional (default=0) + Controls the logs of the optimization algorithm + limit_max: float, optional (default=10) + Controls the semi supervised mode. Transport between labeled source + and target samples of different classes will exhibit an infinite cost + (10 times the maximum value of the cost matrix) + + Attributes + ---------- + coupling_ : array-like, shape (n_source_samples, n_target_samples) + The optimal coupling + log_ : dictionary + The dictionary of log, empty dic if parameter log is not True + + References + ---------- + + .. [1] N. Courty; R. Flamary; D. Tuia; A. Rakotomamonjy, + "Optimal Transport for Domain Adaptation," in IEEE + Transactions on Pattern Analysis and Machine Intelligence , + vol.PP, no.99, pp.1-1 + .. [2] Rakotomamonjy, A., Flamary, R., & Courty, N. (2015). + Generalized conditional gradient: analysis of convergence + and applications. arXiv preprint arXiv:1510.06567. + + """ + + def __init__(self, reg_e=1., reg_cl=0.1, + max_iter=10, max_inner_iter=200, + tol=10e-9, verbose=False, log=False, + metric="sqeuclidean", norm=None, + distribution_estimation=distribution_estimation_uniform, + out_of_sample_map='ferradans', limit_max=10): + + self.reg_e = reg_e + self.reg_cl = reg_cl + self.max_iter = max_iter + self.max_inner_iter = max_inner_iter + self.tol = tol + self.verbose = verbose + self.log = log + self.metric = metric + self.norm = norm + self.distribution_estimation = distribution_estimation + self.out_of_sample_map = out_of_sample_map + self.limit_max = limit_max + + def fit(self, Xs, ys=None, Xt=None, yt=None): + """Build a coupling matrix from source and target sets of samples + (Xs, ys) and (Xt, yt) + + Parameters + ---------- + Xs : array-like, shape (n_source_samples, n_features) + The training input samples. + ys : array-like, shape (n_source_samples,) + The class labels + Xt : array-like, shape (n_target_samples, n_features) + The training input samples. + yt : array-like, shape (n_labeled_target_samples,) + The class labels + + Returns + ------- + self : object + Returns self. + """ + + # check the necessary inputs parameters are here + if check_params(Xs=Xs, Xt=Xt, ys=ys): + + super(SinkhornL1l2Transport, self).fit(Xs, ys, Xt, yt) + + 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) + + # deal with the value of log + if self.log: + self.coupling_, self.log_ = returned_ + else: + self.coupling_ = returned_ + self.log_ = dict() + + return self + + +class MappingTransport(BaseEstimator): + """MappingTransport: DA methods that aims at jointly estimating a optimal + transport coupling and the associated mapping + + Parameters + ---------- + mu : float, optional (default=1) + Weight for the linear OT loss (>0) + eta : float, optional (default=0.001) + Regularization term for the linear mapping L (>0) + bias : bool, optional (default=False) + Estimate linear mapping with constant bias + metric : string, optional (default="sqeuclidean") + The ground metric for the Wasserstein problem + norm : string, optional (default=None) + If given, normalize the ground metric to avoid numerical errors that + can occur with large metric values. + kernel : string, optional (default="linear") + The kernel to use either linear or gaussian + sigma : float, optional (default=1) + The gaussian kernel parameter + max_iter : int, optional (default=100) + Max number of BCD iterations + tol : float, optional (default=1e-5) + Stop threshold on relative loss decrease (>0) + max_inner_iter : int, optional (default=10) + Max number of iterations (inner CG solver) + inner_tol : float, optional (default=1e-6) + Stop threshold on error (inner CG solver) (>0) + verbose : bool, optional (default=False) + Print information along iterations + log : bool, optional (default=False) + record log if True + + Attributes + ---------- + coupling_ : array-like, shape (n_source_samples, n_target_samples) + The optimal coupling + mapping_ : array-like, shape (n_features (+ 1), n_features) + (if bias) for kernel == linear + The associated mapping + array-like, shape (n_source_samples (+ 1), n_features) + (if bias) for kernel == gaussian + log_ : dictionary + The dictionary of log, empty dic if parameter log is not True + + References + ---------- + + .. [8] M. Perrot, N. Courty, R. Flamary, A. Habrard, + "Mapping estimation for discrete optimal transport", + Neural Information Processing Systems (NIPS), 2016. + + """ + + def __init__(self, mu=1, eta=0.001, bias=False, metric="sqeuclidean", + norm=None, kernel="linear", sigma=1, max_iter=100, tol=1e-5, + max_inner_iter=10, inner_tol=1e-6, log=False, verbose=False, + verbose2=False): + + self.metric = metric + self.norm = norm + self.mu = mu + self.eta = eta + self.bias = bias + self.kernel = kernel + self.sigma = sigma + self.max_iter = max_iter + self.tol = tol + self.max_inner_iter = max_inner_iter + self.inner_tol = inner_tol + self.log = log + self.verbose = verbose + self.verbose2 = verbose2 + + def fit(self, Xs=None, ys=None, Xt=None, yt=None): + """Builds an optimal coupling and estimates the associated mapping + from source and target sets of samples (Xs, ys) and (Xt, yt) + + Parameters + ---------- + Xs : array-like, shape (n_source_samples, n_features) + The training input samples. + ys : array-like, shape (n_source_samples,) + The class labels + Xt : array-like, shape (n_target_samples, n_features) + The training input samples. + yt : array-like, shape (n_labeled_target_samples,) + The class labels + + Returns + ------- + self : object + Returns self + """ + + # check the necessary inputs parameters are here + if check_params(Xs=Xs, Xt=Xt): + + self.xs_ = Xs + self.xt_ = Xt + + if self.kernel == "linear": + returned_ = joint_OT_mapping_linear( + Xs, Xt, mu=self.mu, eta=self.eta, bias=self.bias, + verbose=self.verbose, verbose2=self.verbose2, + numItermax=self.max_iter, + numInnerItermax=self.max_inner_iter, stopThr=self.tol, + stopInnerThr=self.inner_tol, log=self.log) + + elif self.kernel == "gaussian": + returned_ = joint_OT_mapping_kernel( + Xs, Xt, mu=self.mu, eta=self.eta, bias=self.bias, + sigma=self.sigma, verbose=self.verbose, + verbose2=self.verbose, numItermax=self.max_iter, + numInnerItermax=self.max_inner_iter, + stopInnerThr=self.inner_tol, stopThr=self.tol, + log=self.log) + + # deal with the value of log + if self.log: + self.coupling_, self.mapping_, self.log_ = returned_ + else: + self.coupling_, self.mapping_ = returned_ + self.log_ = dict() + + return self + + def transform(self, Xs): + """Transports source samples Xs onto target ones Xt + + Parameters + ---------- + Xs : array-like, shape (n_source_samples, n_features) + The training input samples. + + Returns + ------- + transp_Xs : array-like, shape (n_source_samples, n_features) + The transport source samples. + """ + + # check the necessary inputs parameters are here + if check_params(Xs=Xs): + + 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 + + # compute transported samples + transp_Xs = np.dot(transp, self.xt_) + else: + if self.kernel == "gaussian": + K = kernel(Xs, self.xs_, method=self.kernel, + sigma=self.sigma) + elif self.kernel == "linear": + K = Xs + if self.bias: + K = np.hstack((K, np.ones((Xs.shape[0], 1)))) + transp_Xs = K.dot(self.mapping_) + + return transp_Xs diff --git a/ot/datasets.py b/ot/datasets.py index 7816833..e4fe118 100644 --- a/ot/datasets.py +++ b/ot/datasets.py @@ -2,12 +2,16 @@ Simple example datasets for OT """ +# Author: Remi Flamary <remi.flamary@unice.fr> +# +# License: MIT License + import numpy as np import scipy as sp -def get_1D_gauss(n,m,s): +def get_1D_gauss(n, m, s): """return a 1D histogram for a gaussian distribution (n bins, mean m and std s) Parameters @@ -27,12 +31,12 @@ def get_1D_gauss(n,m,s): 1D histogram for a gaussian distribution """ - x=np.arange(n,dtype=np.float64) - h=np.exp(-(x-m)**2/(2*s**2)) - return h/h.sum() + x = np.arange(n, dtype=np.float64) + h = np.exp(-(x - m)**2 / (2 * s**2)) + return h / h.sum() -def get_2D_samples_gauss(n,m,sigma): +def get_2D_samples_gauss(n, m, sigma): """return n samples drawn from 2D gaussian N(m,sigma) Parameters @@ -52,17 +56,17 @@ def get_2D_samples_gauss(n,m,sigma): n samples drawn from N(m,sigma) """ - if np.isscalar(sigma): - sigma=np.array([sigma,]) - if len(sigma)>1: - P=sp.linalg.sqrtm(sigma) - res= np.random.randn(n,2).dot(P)+m + if np.isscalar(sigma): + sigma = np.array([sigma, ]) + if len(sigma) > 1: + P = sp.linalg.sqrtm(sigma) + res = np.random.randn(n, 2).dot(P) + m else: - res= np.random.randn(n,2)*np.sqrt(sigma)+m + res = np.random.randn(n, 2) * np.sqrt(sigma) + m return res -def get_data_classif(dataset,n,nz=.5,theta=0,**kwargs): +def get_data_classif(dataset, n, nz=.5, theta=0, **kwargs): """ dataset generation for classification problems Parameters @@ -84,48 +88,53 @@ def get_data_classif(dataset,n,nz=.5,theta=0,**kwargs): labels of the samples """ - if dataset.lower()=='3gauss': - y=np.floor((np.arange(n)*1.0/n*3))+1 - x=np.zeros((n,2)) + if dataset.lower() == '3gauss': + y = np.floor((np.arange(n) * 1.0 / n * 3)) + 1 + x = np.zeros((n, 2)) # class 1 - x[y==1,0]=-1.; x[y==1,1]=-1. - x[y==2,0]=-1.; x[y==2,1]=1. - x[y==3,0]=1. ; x[y==3,1]=0 - - x[y!=3,:]+=1.5*nz*np.random.randn(sum(y!=3),2) - x[y==3,:]+=2*nz*np.random.randn(sum(y==3),2) - - elif dataset.lower()=='3gauss2': - y=np.floor((np.arange(n)*1.0/n*3))+1 - x=np.zeros((n,2)) - y[y==4]=3 + x[y == 1, 0] = -1. + x[y == 1, 1] = -1. + x[y == 2, 0] = -1. + x[y == 2, 1] = 1. + x[y == 3, 0] = 1. + x[y == 3, 1] = 0 + + x[y != 3, :] += 1.5 * nz * np.random.randn(sum(y != 3), 2) + x[y == 3, :] += 2 * nz * np.random.randn(sum(y == 3), 2) + + elif dataset.lower() == '3gauss2': + y = np.floor((np.arange(n) * 1.0 / n * 3)) + 1 + x = np.zeros((n, 2)) + y[y == 4] = 3 # class 1 - x[y==1,0]=-2.; x[y==1,1]=-2. - x[y==2,0]=-2.; x[y==2,1]=2. - x[y==3,0]=2. ; x[y==3,1]=0 - - x[y!=3,:]+=nz*np.random.randn(sum(y!=3),2) - x[y==3,:]+=2*nz*np.random.randn(sum(y==3),2) - - elif dataset.lower()=='gaussrot' : - rot=np.array([[np.cos(theta),np.sin(theta)],[-np.sin(theta),np.cos(theta)]]) - m1=np.array([-1,1]) - m2=np.array([1,-1]) - y=np.floor((np.arange(n)*1.0/n*2))+1 - n1=np.sum(y==1) - n2=np.sum(y==2) - x=np.zeros((n,2)) - - x[y==1,:]=get_2D_samples_gauss(n1,m1,nz) - x[y==2,:]=get_2D_samples_gauss(n2,m2,nz) - - x=x.dot(rot) - - + x[y == 1, 0] = -2. + x[y == 1, 1] = -2. + x[y == 2, 0] = -2. + x[y == 2, 1] = 2. + x[y == 3, 0] = 2. + x[y == 3, 1] = 0 + + x[y != 3, :] += nz * np.random.randn(sum(y != 3), 2) + x[y == 3, :] += 2 * nz * np.random.randn(sum(y == 3), 2) + + elif dataset.lower() == 'gaussrot': + rot = np.array( + [[np.cos(theta), np.sin(theta)], [-np.sin(theta), np.cos(theta)]]) + m1 = np.array([-1, 1]) + m2 = np.array([1, -1]) + y = np.floor((np.arange(n) * 1.0 / n * 2)) + 1 + n1 = np.sum(y == 1) + n2 = np.sum(y == 2) + x = np.zeros((n, 2)) + + x[y == 1, :] = get_2D_samples_gauss(n1, m1, nz) + x[y == 2, :] = get_2D_samples_gauss(n2, m2, nz) + + x = x.dot(rot) else: - x=np.array(0) - y=np.array(0) + x = np.array(0) + y = np.array(0) print("unknown dataset") - return x,y.astype(int)
\ No newline at end of file + return x, y.astype(int) @@ -3,43 +3,50 @@ Dimension reduction with optimal transport """ +# Author: Remi Flamary <remi.flamary@unice.fr> +# +# License: MIT License + +from scipy import linalg import autograd.numpy as np from pymanopt.manifolds import Stiefel from pymanopt import Problem from pymanopt.solvers import SteepestDescent, TrustRegions -import scipy.linalg as la -def dist(x1,x2): + +def dist(x1, x2): """ Compute squared euclidean distance between samples (autograd) """ - x1p2=np.sum(np.square(x1),1) - x2p2=np.sum(np.square(x2),1) - return x1p2.reshape((-1,1))+x2p2.reshape((1,-1))-2*np.dot(x1,x2.T) + x1p2 = np.sum(np.square(x1), 1) + x2p2 = np.sum(np.square(x2), 1) + return x1p2.reshape((-1, 1)) + x2p2.reshape((1, -1)) - 2 * np.dot(x1, x2.T) + -def sinkhorn(w1,w2,M,reg,k): +def sinkhorn(w1, w2, M, reg, k): """Sinkhorn algorithm with fixed number of iteration (autograd) """ - K=np.exp(-M/reg) - ui=np.ones((M.shape[0],)) - vi=np.ones((M.shape[1],)) + K = np.exp(-M / reg) + ui = np.ones((M.shape[0],)) + vi = np.ones((M.shape[1],)) for i in range(k): - vi=w2/(np.dot(K.T,ui)) - ui=w1/(np.dot(K,vi)) - G=ui.reshape((M.shape[0],1))*K*vi.reshape((1,M.shape[1])) + vi = w2 / (np.dot(K.T, ui)) + ui = w1 / (np.dot(K, vi)) + G = ui.reshape((M.shape[0], 1)) * K * vi.reshape((1, M.shape[1])) return G -def split_classes(X,y): + +def split_classes(X, y): """split samples in X by classes in y """ - lstsclass=np.unique(y) - return [X[y==i,:].astype(np.float32) for i in lstsclass] + lstsclass = np.unique(y) + return [X[y == i, :].astype(np.float32) for i in lstsclass] + + +def fda(X, y, p=2, reg=1e-16): + """ + Fisher Discriminant Analysis -def fda(X,y,p=2,reg=1e-16): - """ - Fisher Discriminant Analysis - - Parameters ---------- X : numpy.ndarray (n,d) @@ -59,62 +66,62 @@ def fda(X,y,p=2,reg=1e-16): proj : fun projection function including mean centering - - """ - - mx=np.mean(X) - X-=mx.reshape((1,-1)) - + + """ + + mx = np.mean(X) + X -= mx.reshape((1, -1)) + # data split between classes - d=X.shape[1] - xc=split_classes(X,y) - nc=len(xc) - - p=min(nc-1,p) - - Cw=0 + d = X.shape[1] + xc = split_classes(X, y) + nc = len(xc) + + p = min(nc - 1, p) + + Cw = 0 for x in xc: - Cw+=np.cov(x,rowvar=False) - Cw/=nc - - mxc=np.zeros((d,nc)) - + Cw += np.cov(x, rowvar=False) + Cw /= nc + + mxc = np.zeros((d, nc)) + for i in range(nc): - mxc[:,i]=np.mean(xc[i]) - - mx0=np.mean(mxc,1) - Cb=0 + mxc[:, i] = np.mean(xc[i]) + + mx0 = np.mean(mxc, 1) + Cb = 0 for i in range(nc): - Cb+=(mxc[:,i]-mx0).reshape((-1,1))*(mxc[:,i]-mx0).reshape((1,-1)) - - w,V=la.eig(Cb,Cw+reg*np.eye(d)) - - idx=np.argsort(w.real) - - Popt=V[:,idx[-p:]] - - - + Cb += (mxc[:, i] - mx0).reshape((-1, 1)) * \ + (mxc[:, i] - mx0).reshape((1, -1)) + + w, V = linalg.eig(Cb, Cw + reg * np.eye(d)) + + idx = np.argsort(w.real) + + Popt = V[:, idx[-p:]] + def proj(X): - return (X-mx.reshape((1,-1))).dot(Popt) - + return (X - mx.reshape((1, -1))).dot(Popt) + return Popt, proj -def wda(X,y,p=2,reg=1,k=10,solver = None,maxiter=100,verbose=0,P0=None): - """ + +def wda(X, y, p=2, reg=1, k=10, solver=None, maxiter=100, verbose=0, P0=None): + """ Wasserstein Discriminant Analysis [11]_ - + The function solves the following optimization problem: .. math:: P = \\text{arg}\min_P \\frac{\\sum_i W(PX^i,PX^i)}{\\sum_{i,j\\neq i} W(PX^i,PX^j)} where : - + - :math:`P` is a linear projection operator in the Stiefel(p,d) manifold - :math:`W` is entropic regularized Wasserstein distances - - :math:`X^i` are samples in the dataset corresponding to class i - + - :math:`X^i` are samples in the dataset corresponding to class i + Parameters ---------- X : numpy.ndarray (n,d) @@ -147,54 +154,50 @@ def wda(X,y,p=2,reg=1,k=10,solver = None,maxiter=100,verbose=0,P0=None): ---------- .. [11] Flamary, R., Cuturi, M., Courty, N., & Rakotomamonjy, A. (2016). Wasserstein Discriminant Analysis. arXiv preprint arXiv:1608.08063. - - """ - - mx=np.mean(X) - X-=mx.reshape((1,-1)) - + + """ # noqa + + mx = np.mean(X) + X -= mx.reshape((1, -1)) + # data split between classes - d=X.shape[1] - xc=split_classes(X,y) + d = X.shape[1] + xc = split_classes(X, y) # compute uniform weighs - wc=[np.ones((x.shape[0]),dtype=np.float32)/x.shape[0] for x in xc] - + wc = [np.ones((x.shape[0]), dtype=np.float32) / x.shape[0] for x in xc] + def cost(P): # wda loss - loss_b=0 - loss_w=0 - - for i,xi in enumerate(xc): - xi=np.dot(xi,P) - for j,xj in enumerate(xc[i:]): - xj=np.dot(xj,P) - M=dist(xi,xj) - G=sinkhorn(wc[i],wc[j+i],M,reg,k) - if j==0: - loss_w+=np.sum(G*M) + loss_b = 0 + loss_w = 0 + + for i, xi in enumerate(xc): + xi = np.dot(xi, P) + for j, xj in enumerate(xc[i:]): + xj = np.dot(xj, P) + M = dist(xi, xj) + G = sinkhorn(wc[i], wc[j + i], M, reg, k) + if j == 0: + loss_w += np.sum(G * M) else: - loss_b+=np.sum(G*M) - - # loss inversed because minimization - return loss_w/loss_b - - + loss_b += np.sum(G * M) + + # loss inversed because minimization + return loss_w / loss_b + # declare manifold and problem - manifold = Stiefel(d, p) + manifold = Stiefel(d, p) problem = Problem(manifold=manifold, cost=cost) - + # declare solver and solve if solver is None: - solver= SteepestDescent(maxiter=maxiter,logverbosity=verbose) - elif solver in ['tr','TrustRegions']: - solver= TrustRegions(maxiter=maxiter,logverbosity=verbose) - - Popt = solver.solve(problem,x=P0) - - def proj(X): - return (X-mx.reshape((1,-1))).dot(Popt) - - return Popt, proj + solver = SteepestDescent(maxiter=maxiter, logverbosity=verbose) + elif solver in ['tr', 'TrustRegions']: + solver = TrustRegions(maxiter=maxiter, logverbosity=verbose) + Popt = solver.solve(problem, x=P0) + def proj(X): + return (X - mx.reshape((1, -1))).dot(Popt) + return Popt, proj diff --git a/ot/gpu/__init__.py b/ot/gpu/__init__.py index 40b11c0..c8f9433 100644 --- a/ot/gpu/__init__.py +++ b/ot/gpu/__init__.py @@ -4,4 +4,9 @@ from . import bregman from . import da from .bregman import sinkhorn +# Author: Remi Flamary <remi.flamary@unice.fr> +# Leo Gautheron <https://github.com/aje> +# +# License: MIT License + __all__ = ["bregman", "da", "sinkhorn"] diff --git a/ot/gpu/bregman.py b/ot/gpu/bregman.py index 2c3e317..47939c4 100644 --- a/ot/gpu/bregman.py +++ b/ot/gpu/bregman.py @@ -3,13 +3,18 @@ Bregman projections for regularized OT with GPU """ +# Author: Remi Flamary <remi.flamary@unice.fr> +# Leo Gautheron <https://github.com/aje> +# +# License: MIT License + import numpy as np import cudamat def sinkhorn(a, b, M_GPU, reg, numItermax=1000, stopThr=1e-9, verbose=False, log=False, returnAsGPU=False): - """ + r""" Solve the entropic regularization optimal transport problem on GPU The function solves the following optimization problem: @@ -82,7 +87,7 @@ def sinkhorn(a, b, M_GPU, reg, numItermax=1000, stopThr=1e-9, verbose=False, ot.lp.emd : Unregularized OT ot.optim.cg : General regularized OT - """ + """ # init data Nini = len(a) Nfin = len(b) @@ -92,11 +97,11 @@ def sinkhorn(a, b, M_GPU, reg, numItermax=1000, stopThr=1e-9, verbose=False, # we assume that no distances are null except those of the diagonal of # distances - u = (np.ones(Nini)/Nini).reshape((Nini, 1)) + u = (np.ones(Nini) / Nini).reshape((Nini, 1)) u_GPU = cudamat.CUDAMatrix(u) a_GPU = cudamat.CUDAMatrix(a.reshape((Nini, 1))) ones_GPU = cudamat.empty(u_GPU.shape).assign(1) - v = (np.ones(Nfin)/Nfin).reshape((Nfin, 1)) + v = (np.ones(Nfin) / Nfin).reshape((Nfin, 1)) v_GPU = cudamat.CUDAMatrix(v) b_GPU = cudamat.CUDAMatrix(b.reshape((Nfin, 1))) @@ -121,7 +126,7 @@ def sinkhorn(a, b, M_GPU, reg, numItermax=1000, stopThr=1e-9, verbose=False, ones_GPU.divide(Kp_GPU.dot(v_GPU), target=u_GPU) if (np.any(KtransposeU_GPU.asarray() == 0) or - not u_GPU.allfinite() or not v_GPU.allfinite()): + not u_GPU.allfinite() or not v_GPU.allfinite()): # we have reached the machine precision # come back to previous solution and quit loop print('Warning: numerical errors at iteration', cpt) @@ -142,7 +147,8 @@ def sinkhorn(a, b, M_GPU, reg, numItermax=1000, stopThr=1e-9, verbose=False, if verbose: if cpt % 200 == 0: - print('{:5s}|{:12s}'.format('It.', 'Err')+'\n'+'-'*19) + print( + '{:5s}|{:12s}'.format('It.', 'Err') + '\n' + '-' * 19) print('{:5d}|{:8e}|'.format(cpt, err)) cpt += 1 if log: diff --git a/ot/gpu/da.py b/ot/gpu/da.py index b05ff70..05c580f 100644 --- a/ot/gpu/da.py +++ b/ot/gpu/da.py @@ -3,6 +3,14 @@ Domain adaptation with optimal transport with GPU implementation """ +# Author: Remi Flamary <remi.flamary@unice.fr> +# Nicolas Courty <ncourty@irisa.fr> +# Michael Perrot <michael.perrot@univ-st-etienne.fr> +# Leo Gautheron <https://github.com/aje> +# +# License: MIT License + + import numpy as np from ..utils import unif from ..da import OTDA @@ -167,7 +175,7 @@ def sinkhorn_lpl1_mm(a, labels_a, b, M_GPU, reg, eta=0.1, numItermax=10, tmpC_GPU = cudamat.empty((Nfin, nbRow)).assign(0) transp_GPU.transpose().select_columns(indices_labels[i], tmpC_GPU) majs_GPU = tmpC_GPU.sum(axis=1).add(epsilon) - cudamat.pow(majs_GPU, (p-1)) + cudamat.pow(majs_GPU, (p - 1)) majs_GPU.mult(p) tmpC_GPU.assign(0) diff --git a/ot/lp/EMD.h b/ot/lp/EMD.h index fb7feca..15e9115 100644 --- a/ot/lp/EMD.h +++ b/ot/lp/EMD.h @@ -23,8 +23,12 @@ using namespace lemon; typedef unsigned int node_id_type; +enum ProblemType { + INFEASIBLE, + OPTIMAL, + UNBOUNDED +}; -void EMD_wrap(int n1,int n2, double *X, double *Y,double *D, double *G, - double* alpha, double* beta, double *cost, int max_iter); +int EMD_wrap(int n1,int n2, double *X, double *Y,double *D, double *G, double* alpha, double* beta, double *cost, int max_iter); #endif diff --git a/ot/lp/EMD_wrapper.cpp b/ot/lp/EMD_wrapper.cpp index 0977e75..8ac43c7 100644 --- a/ot/lp/EMD_wrapper.cpp +++ b/ot/lp/EMD_wrapper.cpp @@ -15,10 +15,11 @@ #include "EMD.h" -void EMD_wrap(int n1, int n2, double *X, double *Y, double *D, double *G, +int EMD_wrap(int n1, int n2, double *X, double *Y, double *D, double *G, double* alpha, double* beta, double *cost, int max_iter) { // beware M and C anre strored in row major C style!!! - int n, m, i,cur; + int n, m, i, cur; + double max; typedef FullBipartiteDigraph Digraph; DIGRAPH_TYPEDEFS(FullBipartiteDigraph); @@ -44,7 +45,7 @@ void EMD_wrap(int n1, int n2, double *X, double *Y, double *D, double *G, std::vector<int> indI(n), indJ(m); std::vector<double> weights1(n), weights2(m); Digraph di(n, m); - NetworkSimplexSimple<Digraph,double,double, node_id_type> net(di, true, n+m, n*m,max_iter); + NetworkSimplexSimple<Digraph,double,double, node_id_type> net(di, true, n+m, n*m, max_iter); // Set supply and demand, don't account for 0 values (faster) @@ -108,5 +109,5 @@ void EMD_wrap(int n1, int n2, double *X, double *Y, double *D, double *G, }; - + return ret; } diff --git a/ot/lp/__init__.py b/ot/lp/__init__.py index 915a18c..a14d4e4 100644 --- a/ot/lp/__init__.py +++ b/ot/lp/__init__.py @@ -3,6 +3,10 @@ Solvers for the original linear program OT problem """ +# Author: Remi Flamary <remi.flamary@unice.fr> +# +# License: MIT License + import numpy as np # import compiled emd from .emd_wrap import emd_c, emd2_c @@ -10,8 +14,7 @@ from ..utils import parmap import multiprocessing - -def emd(a, b, M, dual_variables=False, max_iter=-1): +def emd(a, b, M, numItermax=100000, dual_variables=False): """Solves the Earth Movers distance problem and returns the OT matrix @@ -36,6 +39,9 @@ def emd(a, b, M, dual_variables=False, max_iter=-1): Target histogram (uniform weigth if empty list) M : (ns,nt) ndarray, float64 loss matrix + numItermax : int, optional (default=100000) + The maximum number of iterations before stopping the optimization + algorithm if it has not converged. Returns ------- @@ -48,7 +54,7 @@ def emd(a, b, M, dual_variables=False, max_iter=-1): Simple example with obvious solution. The function emd accepts lists and perform automatic conversion to numpy arrays - + >>> import ot >>> a=[.5,.5] >>> b=[.5,.5] @@ -80,13 +86,13 @@ def emd(a, b, M, dual_variables=False, max_iter=-1): if len(b) == 0: b = np.ones((M.shape[1], ), dtype=np.float64)/M.shape[1] - G, alpha, beta = emd_c(a, b, M, max_iter) + G, alpha, beta = emd_c(a, b, M, numItermax) if dual_variables: return G, alpha, beta return G -def emd2(a, b, M,processes=multiprocessing.cpu_count(), max_iter=-1): - """Solves the Earth Movers distance problem and returns the loss +def emd2(a, b, M, processes=multiprocessing.cpu_count(), numItermax=100000): + """Solves the Earth Movers distance problem and returns the loss .. math:: \gamma = arg\min_\gamma <\gamma,M>_F @@ -109,6 +115,9 @@ def emd2(a, b, M,processes=multiprocessing.cpu_count(), max_iter=-1): Target histogram (uniform weigth if empty list) M : (ns,nt) ndarray, float64 loss matrix + numItermax : int, optional (default=100000) + The maximum number of iterations before stopping the optimization + algorithm if it has not converged. Returns ------- @@ -121,15 +130,15 @@ def emd2(a, b, M,processes=multiprocessing.cpu_count(), max_iter=-1): Simple example with obvious solution. The function emd accepts lists and perform automatic conversion to numpy arrays - - + + >>> import ot >>> a=[.5,.5] >>> b=[.5,.5] >>> M=[[0.,1.],[1.,0.]] >>> ot.emd2(a,b,M) 0.0 - + References ---------- @@ -152,16 +161,13 @@ def emd2(a, b, M,processes=multiprocessing.cpu_count(), max_iter=-1): a = np.ones((M.shape[0], ), dtype=np.float64)/M.shape[0] if len(b) == 0: b = np.ones((M.shape[1], ), dtype=np.float64)/M.shape[1] - + if len(b.shape)==1: - return emd2_c(a, b, M, max_iter)[0] - else: - nb=b.shape[1] - #res=[emd2_c(a,b[:,i].copy(),M) for i in range(nb)] - def f(b): - return emd2_c(a,b,M, max_iter)[0] - res= parmap(f, [b[:,i] for i in range(nb)],processes) - return np.array(res) - - -
\ No newline at end of file + return emd2_c(a, b, M, numItermax)[0] + nb = b.shape[1] + # res = [emd2_c(a, b[:, i].copy(), M, numItermax) for i in range(nb)] + + def f(b): + return emd2_c(a,b,M, max_iter)[0] + res= parmap(f, [b[:,i] for i in range(nb)],processes) + return np.array(res) diff --git a/ot/lp/emd_wrap.pyx b/ot/lp/emd_wrap.pyx index 4febb32..7056e0e 100644 --- a/ot/lp/emd_wrap.pyx +++ b/ot/lp/emd_wrap.pyx @@ -1,9 +1,12 @@ # -*- coding: utf-8 -*- """ -Created on Thu Sep 11 08:42:08 2014 - -@author: rflamary +Cython linker with C solver """ + +# Author: Remi Flamary <remi.flamary@unice.fr> +# +# License: MIT License + import numpy as np cimport numpy as np @@ -12,47 +15,50 @@ cimport cython cdef extern from "EMD.h": - void EMD_wrap(int n1,int n2, double *X, double *Y,double *D, double *G, double *cost, - double* alpha, double* beta, int max_iter) + int EMD_wrap(int n1,int n2, double *X, double *Y,double *D, double *G, double* alpha, double* beta, double *cost, int max_iter) + cdef enum ProblemType: INFEASIBLE, OPTIMAL, UNBOUNDED @cython.boundscheck(False) @cython.wraparound(False) -def emd_c( np.ndarray[double, ndim=1, mode="c"] a,np.ndarray[double, ndim=1, mode="c"] b,np.ndarray[double, ndim=2, mode="c"] M, int maxiter): +def emd_c( np.ndarray[double, ndim=1, mode="c"] a,np.ndarray[double, ndim=1, mode="c"] b,np.ndarray[double, ndim=2, mode="c"] M, int max_iter): """ Solves the Earth Movers distance problem and returns the optimal transport matrix - + gamm=emd(a,b,M) - + .. math:: - \gamma = arg\min_\gamma <\gamma,M>_F - + \gamma = arg\min_\gamma <\gamma,M>_F + s.t. \gamma 1 = a - - \gamma^T 1= b - + + \gamma^T 1= b + \gamma\geq 0 where : - + - M is the metric cost matrix - a and b are the sample weights - + Parameters ---------- a : (ns,) ndarray, float64 - source histogram + source histogram b : (nt,) ndarray, float64 target histogram M : (ns,nt) ndarray, float64 - loss matrix - - + loss matrix + max_iter : int + The maximum number of iterations before stopping the optimization + algorithm if it has not converged. + + Returns ------- gamma: (ns x nt) ndarray Optimal transportation matrix for the given parameters - + """ cdef int n1= M.shape[0] cdef int n2= M.shape[1] @@ -61,7 +67,8 @@ def emd_c( np.ndarray[double, ndim=1, mode="c"] a,np.ndarray[double, ndim=1, mod cdef np.ndarray[double, ndim=2, mode="c"] G=np.zeros([n1, n2]) cdef np.ndarray[double, ndim=1, mode="c"] alpha=np.zeros(n1) cdef np.ndarray[double, ndim=1, mode="c"] beta=np.zeros(n2) - + + if not len(a): a=np.ones((n1,))/n1 @@ -69,47 +76,54 @@ def emd_c( np.ndarray[double, ndim=1, mode="c"] a,np.ndarray[double, ndim=1, mod b=np.ones((n2,))/n2 # calling the function - EMD_wrap(n1,n2,<double*> a.data,<double*> b.data,<double*> M.data,<double*> G.data, - <double*> alpha.data, <double*> beta.data, <double*> &cost, maxiter) + cdef int resultSolver = EMD_wrap(n1,n2,<double*> a.data,<double*> b.data,<double*> M.data,<double*> G.data, <double*> alpha.data, <double*> beta.data, <double*> &cost, max_iter) + if resultSolver != OPTIMAL: + if resultSolver == INFEASIBLE: + print("Problem infeasible. Try to increase numItermax.") + elif resultSolver == UNBOUNDED: + print("Problem unbounded") return G, alpha, beta @cython.boundscheck(False) @cython.wraparound(False) -def emd2_c( np.ndarray[double, ndim=1, mode="c"] a,np.ndarray[double, ndim=1, mode="c"] b,np.ndarray[double, ndim=2, mode="c"] M, int maxiter): +def emd2_c( np.ndarray[double, ndim=1, mode="c"] a,np.ndarray[double, ndim=1, mode="c"] b,np.ndarray[double, ndim=2, mode="c"] M, int max_iter): """ Solves the Earth Movers distance problem and returns the optimal transport loss - + gamm=emd(a,b,M) - + .. math:: - \gamma = arg\min_\gamma <\gamma,M>_F - + \gamma = arg\min_\gamma <\gamma,M>_F + s.t. \gamma 1 = a - - \gamma^T 1= b - + + \gamma^T 1= b + \gamma\geq 0 where : - + - M is the metric cost matrix - a and b are the sample weights - + Parameters ---------- a : (ns,) ndarray, float64 - source histogram + source histogram b : (nt,) ndarray, float64 target histogram M : (ns,nt) ndarray, float64 - loss matrix - - + loss matrix + max_iter : int + The maximum number of iterations before stopping the optimization + algorithm if it has not converged. + + Returns ------- gamma: (ns x nt) ndarray Optimal transportation matrix for the given parameters - + """ cdef int n1= M.shape[0] cdef int n2= M.shape[1] @@ -119,16 +133,19 @@ def emd2_c( np.ndarray[double, ndim=1, mode="c"] a,np.ndarray[double, ndim=1, mo cdef np.ndarray[double, ndim = 1, mode = "c"] alpha = np.zeros([n1]) cdef np.ndarray[double, ndim = 1, mode = "c"] beta = np.zeros([n2]) - + if not len(a): a=np.ones((n1,))/n1 if not len(b): b=np.ones((n2,))/n2 # calling the function - EMD_wrap(n1,n2,<double*> a.data,<double*> b.data,<double*> M.data,<double*> G.data, - <double*> alpha.data, <double*> beta.data, <double*> &cost, maxiter) - + cdef int resultSolver = EMD_wrap(n1,n2,<double*> a.data,<double*> b.data,<double*> M.data,<double*> G.data, <double*> alpha.data, <double*> beta.data, <double*> &cost, max_iter) + if resultSolver != OPTIMAL: + if resultSolver == INFEASIBLE: + print("Problem infeasible. Try to inscrease numItermax.") + elif resultSolver == UNBOUNDED: + print("Problem unbounded") return cost, alpha, beta diff --git a/ot/optim.py b/ot/optim.py index 79f4f66..1d09adc 100644 --- a/ot/optim.py +++ b/ot/optim.py @@ -3,13 +3,19 @@ Optimization algorithms for OT """ +# Author: Remi Flamary <remi.flamary@unice.fr> +# +# License: MIT License + import numpy as np from scipy.optimize.linesearch import scalar_search_armijo from .lp import emd from .bregman import sinkhorn # The corresponding scipy function does not work for matrices -def line_search_armijo(f,xk,pk,gfk,old_fval,args=(),c1=1e-4,alpha0=0.99): + + +def line_search_armijo(f, xk, pk, gfk, old_fval, args=(), c1=1e-4, alpha0=0.99): """ Armijo linesearch function that works with matrices @@ -51,20 +57,21 @@ def line_search_armijo(f,xk,pk,gfk,old_fval,args=(),c1=1e-4,alpha0=0.99): def phi(alpha1): fc[0] += 1 - return f(xk + alpha1*pk, *args) + return f(xk + alpha1 * pk, *args) if old_fval is None: phi0 = phi(0.) else: phi0 = old_fval - derphi0 = np.sum(pk*gfk) # Quickfix for matrices - alpha,phi1 = scalar_search_armijo(phi,phi0,derphi0,c1=c1,alpha0=alpha0) + derphi0 = np.sum(pk * gfk) # Quickfix for matrices + alpha, phi1 = scalar_search_armijo( + phi, phi0, derphi0, c1=c1, alpha0=alpha0) - return alpha,fc[0],phi1 + return alpha, fc[0], phi1 -def cg(a,b,M,reg,f,df,G0=None,numItermax = 200,stopThr=1e-9,verbose=False,log=False): +def cg(a, b, M, reg, f, df, G0=None, numItermax=200, stopThr=1e-9, verbose=False, log=False): """ Solve the general regularized OT problem with conditional gradient @@ -128,74 +135,74 @@ def cg(a,b,M,reg,f,df,G0=None,numItermax = 200,stopThr=1e-9,verbose=False,log=Fa """ - loop=1 + loop = 1 if log: - log={'loss':[]} + log = {'loss': []} if G0 is None: - G=np.outer(a,b) + G = np.outer(a, b) else: - G=G0 + G = G0 def cost(G): - return np.sum(M*G)+reg*f(G) + return np.sum(M * G) + reg * f(G) - f_val=cost(G) + f_val = cost(G) if log: log['loss'].append(f_val) - it=0 + it = 0 if verbose: - print('{:5s}|{:12s}|{:8s}'.format('It.','Loss','Delta loss')+'\n'+'-'*32) - print('{:5d}|{:8e}|{:8e}'.format(it,f_val,0)) + print('{:5s}|{:12s}|{:8s}'.format( + 'It.', 'Loss', 'Delta loss') + '\n' + '-' * 32) + print('{:5d}|{:8e}|{:8e}'.format(it, f_val, 0)) while loop: - it+=1 - old_fval=f_val - + it += 1 + old_fval = f_val # problem linearization - Mi=M+reg*df(G) + Mi = M + reg * df(G) # set M positive - Mi+=Mi.min() + Mi += Mi.min() # solve linear program - Gc=emd(a,b,Mi) + Gc = emd(a, b, Mi) - deltaG=Gc-G + deltaG = Gc - G # line search - alpha,fc,f_val = line_search_armijo(cost,G,deltaG,Mi,f_val) + alpha, fc, f_val = line_search_armijo(cost, G, deltaG, Mi, f_val) - G=G+alpha*deltaG + G = G + alpha * deltaG # test convergence - if it>=numItermax: - loop=0 - - delta_fval=(f_val-old_fval)/abs(f_val) - if abs(delta_fval)<stopThr: - loop=0 + if it >= numItermax: + loop = 0 + delta_fval = (f_val - old_fval) / abs(f_val) + if abs(delta_fval) < stopThr: + loop = 0 if log: log['loss'].append(f_val) if verbose: - if it%20 ==0: - print('{:5s}|{:12s}|{:8s}'.format('It.','Loss','Delta loss')+'\n'+'-'*32) - print('{:5d}|{:8e}|{:8e}'.format(it,f_val,delta_fval)) - + if it % 20 == 0: + print('{:5s}|{:12s}|{:8s}'.format( + 'It.', 'Loss', 'Delta loss') + '\n' + '-' * 32) + print('{:5d}|{:8e}|{:8e}'.format(it, f_val, delta_fval)) if log: - return G,log + return G, log else: return G -def gcg(a,b,M,reg1,reg2,f,df,G0=None,numItermax = 10,numInnerItermax = 200,stopThr=1e-9,verbose=False,log=False): + +def gcg(a, b, M, reg1, reg2, f, df, G0=None, numItermax=10, numInnerItermax=200, stopThr=1e-9, verbose=False, log=False): """ Solve the general regularized OT problem with the generalized conditional gradient @@ -264,70 +271,68 @@ def gcg(a,b,M,reg1,reg2,f,df,G0=None,numItermax = 10,numInnerItermax = 200,stopT """ - loop=1 + loop = 1 if log: - log={'loss':[]} + log = {'loss': []} if G0 is None: - G=np.outer(a,b) + G = np.outer(a, b) else: - G=G0 + G = G0 def cost(G): - return np.sum(M*G)+ reg1*np.sum(G*np.log(G)) + reg2*f(G) + return np.sum(M * G) + reg1 * np.sum(G * np.log(G)) + reg2 * f(G) - f_val=cost(G) + f_val = cost(G) if log: log['loss'].append(f_val) - it=0 + it = 0 if verbose: - print('{:5s}|{:12s}|{:8s}'.format('It.','Loss','Delta loss')+'\n'+'-'*32) - print('{:5d}|{:8e}|{:8e}'.format(it,f_val,0)) + print('{:5s}|{:12s}|{:8s}'.format( + 'It.', 'Loss', 'Delta loss') + '\n' + '-' * 32) + print('{:5d}|{:8e}|{:8e}'.format(it, f_val, 0)) while loop: - it+=1 - old_fval=f_val - + it += 1 + old_fval = f_val # problem linearization - Mi=M+reg2*df(G) + Mi = M + reg2 * df(G) # solve linear program with Sinkhorn - #Gc = sinkhorn_stabilized(a,b, Mi, reg1, numItermax = numInnerItermax) - Gc = sinkhorn(a,b, Mi, reg1, numItermax = numInnerItermax) + # Gc = sinkhorn_stabilized(a,b, Mi, reg1, numItermax = numInnerItermax) + Gc = sinkhorn(a, b, Mi, reg1, numItermax=numInnerItermax) - deltaG=Gc-G + deltaG = Gc - G # line search - dcost=Mi+reg1*(1+np.log(G)) #?? - alpha,fc,f_val = line_search_armijo(cost,G,deltaG,dcost,f_val) + dcost = Mi + reg1 * (1 + np.log(G)) # ?? + alpha, fc, f_val = line_search_armijo(cost, G, deltaG, dcost, f_val) - G=G+alpha*deltaG + G = G + alpha * deltaG # test convergence - if it>=numItermax: - loop=0 - - delta_fval=(f_val-old_fval)/abs(f_val) - if abs(delta_fval)<stopThr: - loop=0 + if it >= numItermax: + loop = 0 + delta_fval = (f_val - old_fval) / abs(f_val) + if abs(delta_fval) < stopThr: + loop = 0 if log: log['loss'].append(f_val) if verbose: - if it%20 ==0: - print('{:5s}|{:12s}|{:8s}'.format('It.','Loss','Delta loss')+'\n'+'-'*32) - print('{:5d}|{:8e}|{:8e}'.format(it,f_val,delta_fval)) - + if it % 20 == 0: + print('{:5s}|{:12s}|{:8s}'.format( + 'It.', 'Loss', 'Delta loss') + '\n' + '-' * 32) + print('{:5d}|{:8e}|{:8e}'.format(it, f_val, delta_fval)) if log: - return G,log + return G, log else: return G - @@ -2,91 +2,84 @@ Functions for plotting OT matrices """ +# Author: Remi Flamary <remi.flamary@unice.fr> +# +# License: MIT License import numpy as np import matplotlib.pylab as pl from matplotlib import gridspec -def plot1D_mat(a,b,M,title=''): - """ Plot matrix M with the source and target 1D distribution - - Creates a subplot with the source distribution a on the left and +def plot1D_mat(a, b, M, title=''): + """ Plot matrix M with the source and target 1D distribution + + Creates a subplot with the source distribution a on the left and target distribution b on the tot. The matrix M is shown in between. - - + + Parameters ---------- - - a : np.array (na,) + a : np.array, shape (na,) Source distribution - b : np.array (nb,) - Target distribution - M : np.array (na,nb) + b : np.array, shape (nb,) + Target distribution + M : np.array, shape (na,nb) Matrix to plot - - - """ - - na=M.shape[0] - nb=M.shape[1] - + na, nb = M.shape + gs = gridspec.GridSpec(3, 3) - - - xa=np.arange(na) - xb=np.arange(nb) - - - ax1=pl.subplot(gs[0,1:]) - pl.plot(xb,b,'r',label='Target distribution') + + xa = np.arange(na) + xb = np.arange(nb) + + ax1 = pl.subplot(gs[0, 1:]) + pl.plot(xb, b, 'r', label='Target distribution') pl.yticks(()) pl.title(title) - - #pl.axis('off') - - ax2=pl.subplot(gs[1:,0]) - pl.plot(a,xa,'b',label='Source distribution') + + ax2 = pl.subplot(gs[1:, 0]) + pl.plot(a, xa, 'b', label='Source distribution') pl.gca().invert_xaxis() pl.gca().invert_yaxis() pl.xticks(()) - #pl.ylim((0,n)) - #pl.axis('off') - - pl.subplot(gs[1:,1:],sharex=ax1,sharey=ax2) - pl.imshow(M,interpolation='nearest') - - pl.xlim((0,nb)) + pl.subplot(gs[1:, 1:], sharex=ax1, sharey=ax2) + pl.imshow(M, interpolation='nearest') + pl.axis('off') + + pl.xlim((0, nb)) + pl.tight_layout() + pl.subplots_adjust(wspace=0., hspace=0.2) -def plot2D_samples_mat(xs,xt,G,thr=1e-8,**kwargs): + +def plot2D_samples_mat(xs, xt, G, thr=1e-8, **kwargs): """ Plot matrix M in 2D with lines using alpha values - - Plot lines between source and target 2D samples with a color + + Plot lines between source and target 2D samples with a color proportional to the value of the matrix G between samples. - - + + Parameters ---------- - - xs : np.array (ns,2) + xs : ndarray, shape (ns,2) Source samples positions - b : np.array (nt,2) + b : ndarray, shape (nt,2) Target samples positions - G : np.array (na,nb) + G : ndarray, shape (na,nb) OT matrix thr : float, optional threshold above which the line is drawn **kwargs : dict - paameters given to the plot functions (default color is black if nothing given) - + paameters given to the plot functions (default color is black if + nothing given) """ - if ('color' not in kwargs) and ('c' not in kwargs): - kwargs['color']='k' - mx=G.max() + if ('color' not in kwargs) and ('c' not in kwargs): + kwargs['color'] = 'k' + mx = G.max() for i in range(xs.shape[0]): for j in range(xt.shape[0]): - if G[i,j]/mx>thr: - pl.plot([xs[i,0],xt[j,0]],[xs[i,1],xt[j,1]],alpha=G[i,j]/mx,**kwargs) -
\ No newline at end of file + if G[i, j] / mx > thr: + pl.plot([xs[i, 0], xt[j, 0]], [xs[i, 1], xt[j, 1]], + alpha=G[i, j] / mx, **kwargs) diff --git a/ot/utils.py b/ot/utils.py index 7ad7637..31a002b 100644 --- a/ot/utils.py +++ b/ot/utils.py @@ -2,36 +2,50 @@ """ Various function that can be usefull """ + +# Author: Remi Flamary <remi.flamary@unice.fr> +# +# License: MIT License + +import multiprocessing +from functools import reduce +import time + import numpy as np from scipy.spatial.distance import cdist -import multiprocessing +import sys +import warnings + + +__time_tic_toc = time.time() -import time -__time_tic_toc=time.time() def tic(): """ Python implementation of Matlab tic() function """ global __time_tic_toc - __time_tic_toc=time.time() + __time_tic_toc = time.time() + def toc(message='Elapsed time : {} s'): """ Python implementation of Matlab toc() function """ - t=time.time() - print(message.format(t-__time_tic_toc)) - return t-__time_tic_toc + t = time.time() + print(message.format(t - __time_tic_toc)) + return t - __time_tic_toc + def toq(): """ Python implementation of Julia toc() function """ - t=time.time() - return t-__time_tic_toc + t = time.time() + return t - __time_tic_toc -def kernel(x1,x2,method='gaussian',sigma=1,**kwargs): +def kernel(x1, x2, method='gaussian', sigma=1, **kwargs): """Compute kernel matrix""" - if method.lower() in ['gaussian','gauss','rbf']: - K=np.exp(-dist(x1,x2)/(2*sigma**2)) + if method.lower() in ['gaussian', 'gauss', 'rbf']: + K = np.exp(-dist(x1, x2) / (2 * sigma**2)) return K + def unif(n): """ return a uniform histogram of length n (simplex) @@ -48,17 +62,19 @@ def unif(n): """ - return np.ones((n,))/n + return np.ones((n,)) / n + -def clean_zeros(a,b,M): - """ Remove all components with zeros weights in a and b +def clean_zeros(a, b, M): + """ Remove all components with zeros weights in a and b """ - M2=M[a>0,:][:,b>0].copy() # copy force c style matrix (froemd) - a2=a[a>0] - b2=b[b>0] - return a2,b2,M2 + M2 = M[a > 0, :][:, b > 0].copy() # copy force c style matrix (froemd) + a2 = a[a > 0] + b2 = b[b > 0] + return a2, b2, M2 -def dist(x1,x2=None,metric='sqeuclidean'): + +def dist(x1, x2=None, metric='sqeuclidean'): """Compute distance between samples in x1 and x2 using function scipy.spatial.distance.cdist Parameters @@ -84,12 +100,12 @@ def dist(x1,x2=None,metric='sqeuclidean'): """ if x2 is None: - x2=x1 + x2 = x1 - return cdist(x1,x2,metric=metric) + return cdist(x1, x2, metric=metric) -def dist0(n,method='lin_square'): +def dist0(n, method='lin_square'): """Compute standard cost matrices of size (n,n) for OT problems Parameters @@ -111,16 +127,50 @@ def dist0(n,method='lin_square'): """ - res=0 - if method=='lin_square': - x=np.arange(n,dtype=np.float64).reshape((n,1)) - res=dist(x,x) + res = 0 + if method == 'lin_square': + x = np.arange(n, dtype=np.float64).reshape((n, 1)) + res = dist(x, x) return res +def cost_normalization(C, norm=None): + """ Apply normalization to the loss matrix + + + Parameters + ---------- + C : np.array (n1, n2) + The cost matrix to normalize. + norm : str + type of normalization from 'median','max','log','loglog'. Any other + value do not normalize. + + + Returns + ------- + + C : np.array (n1, n2) + The input cost matrix normalized according to given norm. + + """ + + if norm == "median": + C /= float(np.median(C)) + elif norm == "max": + C /= float(np.max(C)) + elif norm == "log": + C = np.log(1 + C) + elif norm == "loglog": + C = np.log(1 + np.log(1 + C)) + + return C + + def dots(*args): """ dots function for multiple matrix multiply """ - return reduce(np.dot,args) + return reduce(np.dot, args) + def fun(f, q_in, q_out): """ Utility function for parmap with no serializing problems """ @@ -130,6 +180,7 @@ def fun(f, q_in, q_out): break q_out.put((i, f(x))) + def parmap(f, X, nprocs=multiprocessing.cpu_count()): """ paralell map for multiprocessing """ q_in = multiprocessing.Queue(1) @@ -147,4 +198,241 @@ def parmap(f, X, nprocs=multiprocessing.cpu_count()): [p.join() for p in proc] - return [x for i, x in sorted(res)]
\ No newline at end of file + return [x for i, x in sorted(res)] + + +def check_params(**kwargs): + """check_params: check whether some parameters are missing + """ + + missing_params = [] + check = True + + for param in kwargs: + if kwargs[param] is None: + missing_params.append(param) + + if len(missing_params) > 0: + print("POT - Warning: following necessary parameters are missing") + for p in missing_params: + print("\n", p) + + check = False + + return check + + +class deprecated(object): + """Decorator to mark a function or class as deprecated. + + deprecated class from scikit-learn package + https://github.com/scikit-learn/scikit-learn/blob/master/sklearn/utils/deprecation.py + Issue a warning when the function is called/the class is instantiated and + adds a warning to the docstring. + The optional extra argument will be appended to the deprecation message + and the docstring. Note: to use this with the default value for extra, put + in an empty of parentheses: + >>> from ot.deprecation import deprecated + >>> @deprecated() + ... def some_function(): pass + + Parameters + ---------- + extra : string + to be added to the deprecation messages + """ + + # Adapted from http://wiki.python.org/moin/PythonDecoratorLibrary, + # but with many changes. + + def __init__(self, extra=''): + self.extra = extra + + def __call__(self, obj): + """Call method + Parameters + ---------- + obj : object + """ + if isinstance(obj, type): + return self._decorate_class(obj) + else: + return self._decorate_fun(obj) + + def _decorate_class(self, cls): + msg = "Class %s is deprecated" % cls.__name__ + if self.extra: + msg += "; %s" % self.extra + + # FIXME: we should probably reset __new__ for full generality + init = cls.__init__ + + def wrapped(*args, **kwargs): + warnings.warn(msg, category=DeprecationWarning) + return init(*args, **kwargs) + + cls.__init__ = wrapped + + wrapped.__name__ = '__init__' + wrapped.__doc__ = self._update_doc(init.__doc__) + wrapped.deprecated_original = init + + return cls + + def _decorate_fun(self, fun): + """Decorate function fun""" + + msg = "Function %s is deprecated" % fun.__name__ + if self.extra: + msg += "; %s" % self.extra + + def wrapped(*args, **kwargs): + warnings.warn(msg, category=DeprecationWarning) + return fun(*args, **kwargs) + + wrapped.__name__ = fun.__name__ + wrapped.__dict__ = fun.__dict__ + wrapped.__doc__ = self._update_doc(fun.__doc__) + + return wrapped + + def _update_doc(self, olddoc): + newdoc = "DEPRECATED" + if self.extra: + newdoc = "%s: %s" % (newdoc, self.extra) + if olddoc: + newdoc = "%s\n\n%s" % (newdoc, olddoc) + return newdoc + + +def _is_deprecated(func): + """Helper to check if func is wraped by our deprecated decorator""" + if sys.version_info < (3, 5): + raise NotImplementedError("This is only available for python3.5 " + "or above") + closures = getattr(func, '__closure__', []) + if closures is None: + closures = [] + is_deprecated = ('deprecated' in ''.join([c.cell_contents + for c in closures + if isinstance(c.cell_contents, str)])) + return is_deprecated + + +class BaseEstimator(object): + """Base class for most objects in POT + adapted from sklearn BaseEstimator class + + Notes + ----- + All estimators should specify all the parameters that can be set + at the class level in their ``__init__`` as explicit keyword + arguments (no ``*args`` or ``**kwargs``). + """ + + @classmethod + def _get_param_names(cls): + """Get parameter names for the estimator""" + try: + from inspect import signature + except ImportError: + from .externals.funcsigs import signature + # fetch the constructor or the original constructor before + # deprecation wrapping if any + init = getattr(cls.__init__, 'deprecated_original', cls.__init__) + if init is object.__init__: + # No explicit constructor to introspect + return [] + + # introspect the constructor arguments to find the model parameters + # to represent + init_signature = signature(init) + # Consider the constructor parameters excluding 'self' + parameters = [p for p in init_signature.parameters.values() + if p.name != 'self' and p.kind != p.VAR_KEYWORD] + for p in parameters: + if p.kind == p.VAR_POSITIONAL: + raise RuntimeError("POT estimators should always " + "specify their parameters in the signature" + " of their __init__ (no varargs)." + " %s with constructor %s doesn't " + " follow this convention." + % (cls, init_signature)) + # Extract and sort argument names excluding 'self' + return sorted([p.name for p in parameters]) + + def get_params(self, deep=True): + """Get parameters for this estimator. + + Parameters + ---------- + deep : boolean, optional + If True, will return the parameters for this estimator and + contained subobjects that are estimators. + + Returns + ------- + params : mapping of string to any + Parameter names mapped to their values. + """ + out = dict() + for key in self._get_param_names(): + # We need deprecation warnings to always be on in order to + # catch deprecated param values. + # This is set in utils/__init__.py but it gets overwritten + # when running under python3 somehow. + warnings.simplefilter("always", DeprecationWarning) + try: + with warnings.catch_warnings(record=True) as w: + value = getattr(self, key, None) + if len(w) and w[0].category == DeprecationWarning: + # if the parameter is deprecated, don't show it + continue + finally: + warnings.filters.pop(0) + + # XXX: should we rather test if instance of estimator? + if deep and hasattr(value, 'get_params'): + deep_items = value.get_params().items() + out.update((key + '__' + k, val) for k, val in deep_items) + out[key] = value + return out + + def set_params(self, **params): + """Set the parameters of this estimator. + + The method works on simple estimators as well as on nested objects + (such as pipelines). The latter have parameters of the form + ``<component>__<parameter>`` so that it's possible to update each + component of a nested object. + + Returns + ------- + self + """ + if not params: + # Simple optimisation to gain speed (inspect is slow) + return self + valid_params = self.get_params(deep=True) + # for key, value in iteritems(params): + for key, value in params.items(): + split = key.split('__', 1) + if len(split) > 1: + # nested objects case + name, sub_name = split + if name not in valid_params: + raise ValueError('Invalid parameter %s for estimator %s. ' + 'Check the list of available parameters ' + 'with `estimator.get_params().keys()`.' % + (name, self)) + sub_object = valid_params[name] + sub_object.set_params(**{sub_name: value}) + else: + # simple objects case + if key not in valid_params: + raise ValueError('Invalid parameter %s for estimator %s. ' + 'Check the list of available parameters ' + 'with `estimator.get_params().keys()`.' % + (key, self.__class__.__name__)) + setattr(self, key, value) + return self diff --git a/requirements.txt b/requirements.txt index 319f5e1..9bfca43 100644 --- a/requirements.txt +++ b/requirements.txt @@ -3,3 +3,5 @@ scipy cython matplotlib sphinx-gallery +autograd +pymanopt @@ -1,2 +1,6 @@ [metadata] description-file = README.md + +[flake8] +exclude = __init__.py +ignore = E265,E501 diff --git a/test/test_bregman.py b/test/test_bregman.py new file mode 100644 index 0000000..4a800fd --- /dev/null +++ b/test/test_bregman.py @@ -0,0 +1,137 @@ +"""Tests for module bregman on OT with bregman projections """ + +# Author: Remi Flamary <remi.flamary@unice.fr> +# +# License: MIT License + +import numpy as np +import ot + + +def test_sinkhorn(): + # test sinkhorn + n = 100 + rng = np.random.RandomState(0) + + x = rng.randn(n, 2) + u = ot.utils.unif(n) + + M = ot.dist(x, x) + + G = ot.sinkhorn(u, u, M, 1, stopThr=1e-10) + + # check constratints + np.testing.assert_allclose( + u, G.sum(1), atol=1e-05) # cf convergence sinkhorn + np.testing.assert_allclose( + u, G.sum(0), atol=1e-05) # cf convergence sinkhorn + + +def test_sinkhorn_empty(): + # test sinkhorn + n = 100 + rng = np.random.RandomState(0) + + x = rng.randn(n, 2) + u = ot.utils.unif(n) + + M = ot.dist(x, x) + + G, log = ot.sinkhorn([], [], M, 1, stopThr=1e-10, verbose=True, log=True) + # check constratints + np.testing.assert_allclose(u, G.sum(1), atol=1e-05) + np.testing.assert_allclose(u, G.sum(0), atol=1e-05) + + G, log = ot.sinkhorn([], [], M, 1, stopThr=1e-10, + method='sinkhorn_stabilized', verbose=True, log=True) + # check constratints + np.testing.assert_allclose(u, G.sum(1), atol=1e-05) + np.testing.assert_allclose(u, G.sum(0), atol=1e-05) + + G, log = ot.sinkhorn( + [], [], M, 1, stopThr=1e-10, method='sinkhorn_epsilon_scaling', + verbose=True, log=True) + # check constratints + np.testing.assert_allclose(u, G.sum(1), atol=1e-05) + np.testing.assert_allclose(u, G.sum(0), atol=1e-05) + + +def test_sinkhorn_variants(): + # test sinkhorn + n = 100 + rng = np.random.RandomState(0) + + x = rng.randn(n, 2) + u = ot.utils.unif(n) + + M = ot.dist(x, x) + + G0 = ot.sinkhorn(u, u, M, 1, method='sinkhorn', stopThr=1e-10) + Gs = ot.sinkhorn(u, u, M, 1, method='sinkhorn_stabilized', stopThr=1e-10) + Ges = ot.sinkhorn( + u, u, M, 1, method='sinkhorn_epsilon_scaling', stopThr=1e-10) + Gerr = ot.sinkhorn(u, u, M, 1, method='do_not_exists', stopThr=1e-10) + + # check values + np.testing.assert_allclose(G0, Gs, atol=1e-05) + np.testing.assert_allclose(G0, Ges, atol=1e-05) + np.testing.assert_allclose(G0, Gerr) + + +def test_bary(): + + n_bins = 100 # nb bins + + # Gaussian distributions + a1 = ot.datasets.get_1D_gauss(n_bins, m=30, s=10) # m= mean, s= std + a2 = ot.datasets.get_1D_gauss(n_bins, m=40, s=10) + + # creating matrix A containing all distributions + A = np.vstack((a1, a2)).T + + # loss matrix + normalization + M = ot.utils.dist0(n_bins) + M /= M.max() + + alpha = 0.5 # 0<=alpha<=1 + weights = np.array([1 - alpha, alpha]) + + # wasserstein + reg = 1e-3 + bary_wass = ot.bregman.barycenter(A, M, reg, weights) + + np.testing.assert_allclose(1, np.sum(bary_wass)) + + ot.bregman.barycenter(A, M, reg, log=True, verbose=True) + + +def test_unmix(): + + n_bins = 50 # nb bins + + # Gaussian distributions + a1 = ot.datasets.get_1D_gauss(n_bins, m=20, s=10) # m= mean, s= std + a2 = ot.datasets.get_1D_gauss(n_bins, m=40, s=10) + + a = ot.datasets.get_1D_gauss(n_bins, m=30, s=10) + + # creating matrix A containing all distributions + D = np.vstack((a1, a2)).T + + # loss matrix + normalization + M = ot.utils.dist0(n_bins) + M /= M.max() + + M0 = ot.utils.dist0(2) + M0 /= M0.max() + h0 = ot.unif(2) + + # wasserstein + reg = 1e-3 + um = ot.bregman.unmix(a, D, M, M0, h0, reg, 1, alpha=0.01,) + + np.testing.assert_allclose(1, np.sum(um), rtol=1e-03, atol=1e-03) + np.testing.assert_allclose([0.5, 0.5], um, rtol=1e-03, atol=1e-03) + + ot.bregman.unmix(a, D, M, M0, h0, reg, + 1, alpha=0.01, log=True, verbose=True) diff --git a/test/test_da.py b/test/test_da.py new file mode 100644 index 0000000..104a798 --- /dev/null +++ b/test/test_da.py @@ -0,0 +1,469 @@ +"""Tests for module da on Domain Adaptation """ + +# Author: Remi Flamary <remi.flamary@unice.fr> +# +# License: MIT License + +import numpy as np +from numpy.testing.utils import assert_allclose, assert_equal + +import ot +from ot.datasets import get_data_classif +from ot.utils import unif + + +def test_sinkhorn_lpl1_transport_class(): + """test_sinkhorn_transport + """ + + ns = 150 + nt = 200 + + Xs, ys = get_data_classif('3gauss', ns) + Xt, yt = get_data_classif('3gauss2', nt) + + clf = ot.da.SinkhornLpl1Transport() + + # test its computed + clf.fit(Xs=Xs, ys=ys, Xt=Xt) + assert hasattr(clf, "cost_") + assert hasattr(clf, "coupling_") + + # test dimensions of coupling + assert_equal(clf.cost_.shape, ((Xs.shape[0], Xt.shape[0]))) + assert_equal(clf.coupling_.shape, ((Xs.shape[0], Xt.shape[0]))) + + # test margin constraints + mu_s = unif(ns) + mu_t = unif(nt) + assert_allclose(np.sum(clf.coupling_, axis=0), mu_t, rtol=1e-3, atol=1e-3) + assert_allclose(np.sum(clf.coupling_, axis=1), mu_s, rtol=1e-3, atol=1e-3) + + # test transform + transp_Xs = clf.transform(Xs=Xs) + assert_equal(transp_Xs.shape, Xs.shape) + + Xs_new, _ = get_data_classif('3gauss', ns + 1) + transp_Xs_new = clf.transform(Xs_new) + + # check that the oos method is working + assert_equal(transp_Xs_new.shape, Xs_new.shape) + + # test inverse transform + transp_Xt = clf.inverse_transform(Xt=Xt) + assert_equal(transp_Xt.shape, Xt.shape) + + Xt_new, _ = get_data_classif('3gauss2', nt + 1) + transp_Xt_new = clf.inverse_transform(Xt=Xt_new) + + # check that the oos method is working + assert_equal(transp_Xt_new.shape, Xt_new.shape) + + # test fit_transform + transp_Xs = clf.fit_transform(Xs=Xs, ys=ys, Xt=Xt) + assert_equal(transp_Xs.shape, Xs.shape) + + # test semi supervised mode + clf = ot.da.SinkhornLpl1Transport() + clf.fit(Xs=Xs, ys=ys, Xt=Xt) + n_unsup = np.sum(clf.cost_) + + # test semi supervised mode + clf = ot.da.SinkhornLpl1Transport() + clf.fit(Xs=Xs, ys=ys, Xt=Xt, yt=yt) + assert_equal(clf.cost_.shape, ((Xs.shape[0], Xt.shape[0]))) + n_semisup = np.sum(clf.cost_) + + assert n_unsup != n_semisup, "semisupervised mode not working" + + +def test_sinkhorn_l1l2_transport_class(): + """test_sinkhorn_transport + """ + + ns = 150 + nt = 200 + + Xs, ys = get_data_classif('3gauss', ns) + Xt, yt = get_data_classif('3gauss2', nt) + + clf = ot.da.SinkhornL1l2Transport() + + # test its computed + clf.fit(Xs=Xs, ys=ys, Xt=Xt) + assert hasattr(clf, "cost_") + assert hasattr(clf, "coupling_") + assert hasattr(clf, "log_") + + # test dimensions of coupling + assert_equal(clf.cost_.shape, ((Xs.shape[0], Xt.shape[0]))) + assert_equal(clf.coupling_.shape, ((Xs.shape[0], Xt.shape[0]))) + + # test margin constraints + mu_s = unif(ns) + mu_t = unif(nt) + assert_allclose(np.sum(clf.coupling_, axis=0), mu_t, rtol=1e-3, atol=1e-3) + assert_allclose(np.sum(clf.coupling_, axis=1), mu_s, rtol=1e-3, atol=1e-3) + + # test transform + transp_Xs = clf.transform(Xs=Xs) + assert_equal(transp_Xs.shape, Xs.shape) + + Xs_new, _ = get_data_classif('3gauss', ns + 1) + transp_Xs_new = clf.transform(Xs_new) + + # check that the oos method is working + assert_equal(transp_Xs_new.shape, Xs_new.shape) + + # test inverse transform + transp_Xt = clf.inverse_transform(Xt=Xt) + assert_equal(transp_Xt.shape, Xt.shape) + + Xt_new, _ = get_data_classif('3gauss2', nt + 1) + transp_Xt_new = clf.inverse_transform(Xt=Xt_new) + + # check that the oos method is working + assert_equal(transp_Xt_new.shape, Xt_new.shape) + + # test fit_transform + transp_Xs = clf.fit_transform(Xs=Xs, ys=ys, Xt=Xt) + assert_equal(transp_Xs.shape, Xs.shape) + + # test semi supervised mode + clf = ot.da.SinkhornL1l2Transport() + clf.fit(Xs=Xs, ys=ys, Xt=Xt) + n_unsup = np.sum(clf.cost_) + + # test semi supervised mode + clf = ot.da.SinkhornL1l2Transport() + clf.fit(Xs=Xs, ys=ys, Xt=Xt, yt=yt) + assert_equal(clf.cost_.shape, ((Xs.shape[0], Xt.shape[0]))) + n_semisup = np.sum(clf.cost_) + + assert n_unsup != n_semisup, "semisupervised mode not working" + + # check everything runs well with log=True + clf = ot.da.SinkhornL1l2Transport(log=True) + clf.fit(Xs=Xs, ys=ys, Xt=Xt) + assert len(clf.log_.keys()) != 0 + + +def test_sinkhorn_transport_class(): + """test_sinkhorn_transport + """ + + ns = 150 + nt = 200 + + Xs, ys = get_data_classif('3gauss', ns) + Xt, yt = get_data_classif('3gauss2', nt) + + clf = ot.da.SinkhornTransport() + + # test its computed + clf.fit(Xs=Xs, Xt=Xt) + assert hasattr(clf, "cost_") + assert hasattr(clf, "coupling_") + assert hasattr(clf, "log_") + + # test dimensions of coupling + assert_equal(clf.cost_.shape, ((Xs.shape[0], Xt.shape[0]))) + assert_equal(clf.coupling_.shape, ((Xs.shape[0], Xt.shape[0]))) + + # test margin constraints + mu_s = unif(ns) + mu_t = unif(nt) + assert_allclose(np.sum(clf.coupling_, axis=0), mu_t, rtol=1e-3, atol=1e-3) + assert_allclose(np.sum(clf.coupling_, axis=1), mu_s, rtol=1e-3, atol=1e-3) + + # test transform + transp_Xs = clf.transform(Xs=Xs) + assert_equal(transp_Xs.shape, Xs.shape) + + Xs_new, _ = get_data_classif('3gauss', ns + 1) + transp_Xs_new = clf.transform(Xs_new) + + # check that the oos method is working + assert_equal(transp_Xs_new.shape, Xs_new.shape) + + # test inverse transform + transp_Xt = clf.inverse_transform(Xt=Xt) + assert_equal(transp_Xt.shape, Xt.shape) + + Xt_new, _ = get_data_classif('3gauss2', nt + 1) + transp_Xt_new = clf.inverse_transform(Xt=Xt_new) + + # check that the oos method is working + assert_equal(transp_Xt_new.shape, Xt_new.shape) + + # test fit_transform + transp_Xs = clf.fit_transform(Xs=Xs, Xt=Xt) + assert_equal(transp_Xs.shape, Xs.shape) + + # test semi supervised mode + clf = ot.da.SinkhornTransport() + clf.fit(Xs=Xs, Xt=Xt) + n_unsup = np.sum(clf.cost_) + + # test semi supervised mode + clf = ot.da.SinkhornTransport() + clf.fit(Xs=Xs, ys=ys, Xt=Xt, yt=yt) + assert_equal(clf.cost_.shape, ((Xs.shape[0], Xt.shape[0]))) + n_semisup = np.sum(clf.cost_) + + assert n_unsup != n_semisup, "semisupervised mode not working" + + # check everything runs well with log=True + clf = ot.da.SinkhornTransport(log=True) + clf.fit(Xs=Xs, ys=ys, Xt=Xt) + assert len(clf.log_.keys()) != 0 + + +def test_emd_transport_class(): + """test_sinkhorn_transport + """ + + ns = 150 + nt = 200 + + Xs, ys = get_data_classif('3gauss', ns) + Xt, yt = get_data_classif('3gauss2', nt) + + clf = ot.da.EMDTransport() + + # test its computed + clf.fit(Xs=Xs, Xt=Xt) + assert hasattr(clf, "cost_") + assert hasattr(clf, "coupling_") + + # test dimensions of coupling + assert_equal(clf.cost_.shape, ((Xs.shape[0], Xt.shape[0]))) + assert_equal(clf.coupling_.shape, ((Xs.shape[0], Xt.shape[0]))) + + # test margin constraints + mu_s = unif(ns) + mu_t = unif(nt) + assert_allclose(np.sum(clf.coupling_, axis=0), mu_t, rtol=1e-3, atol=1e-3) + assert_allclose(np.sum(clf.coupling_, axis=1), mu_s, rtol=1e-3, atol=1e-3) + + # test transform + transp_Xs = clf.transform(Xs=Xs) + assert_equal(transp_Xs.shape, Xs.shape) + + Xs_new, _ = get_data_classif('3gauss', ns + 1) + transp_Xs_new = clf.transform(Xs_new) + + # check that the oos method is working + assert_equal(transp_Xs_new.shape, Xs_new.shape) + + # test inverse transform + transp_Xt = clf.inverse_transform(Xt=Xt) + assert_equal(transp_Xt.shape, Xt.shape) + + Xt_new, _ = get_data_classif('3gauss2', nt + 1) + transp_Xt_new = clf.inverse_transform(Xt=Xt_new) + + # check that the oos method is working + assert_equal(transp_Xt_new.shape, Xt_new.shape) + + # test fit_transform + transp_Xs = clf.fit_transform(Xs=Xs, Xt=Xt) + assert_equal(transp_Xs.shape, Xs.shape) + + # test semi supervised mode + clf = ot.da.EMDTransport() + clf.fit(Xs=Xs, Xt=Xt) + n_unsup = np.sum(clf.cost_) + + # test semi supervised mode + clf = ot.da.EMDTransport() + clf.fit(Xs=Xs, ys=ys, Xt=Xt, yt=yt) + assert_equal(clf.cost_.shape, ((Xs.shape[0], Xt.shape[0]))) + n_semisup = np.sum(clf.cost_) + + assert n_unsup != n_semisup, "semisupervised mode not working" + + +def test_mapping_transport_class(): + """test_mapping_transport + """ + + ns = 150 + nt = 200 + + Xs, ys = get_data_classif('3gauss', ns) + Xt, yt = get_data_classif('3gauss2', nt) + Xs_new, _ = get_data_classif('3gauss', ns + 1) + + ########################################################################## + # kernel == linear mapping tests + ########################################################################## + + # check computation and dimensions if bias == False + clf = ot.da.MappingTransport(kernel="linear", bias=False) + clf.fit(Xs=Xs, Xt=Xt) + assert hasattr(clf, "coupling_") + assert hasattr(clf, "mapping_") + assert hasattr(clf, "log_") + + assert_equal(clf.coupling_.shape, ((Xs.shape[0], Xt.shape[0]))) + assert_equal(clf.mapping_.shape, ((Xs.shape[1], Xt.shape[1]))) + + # test margin constraints + mu_s = unif(ns) + mu_t = unif(nt) + assert_allclose(np.sum(clf.coupling_, axis=0), mu_t, rtol=1e-3, atol=1e-3) + assert_allclose(np.sum(clf.coupling_, axis=1), mu_s, rtol=1e-3, atol=1e-3) + + # test transform + transp_Xs = clf.transform(Xs=Xs) + assert_equal(transp_Xs.shape, Xs.shape) + + transp_Xs_new = clf.transform(Xs_new) + + # check that the oos method is working + assert_equal(transp_Xs_new.shape, Xs_new.shape) + + # check computation and dimensions if bias == True + clf = ot.da.MappingTransport(kernel="linear", bias=True) + clf.fit(Xs=Xs, Xt=Xt) + assert_equal(clf.coupling_.shape, ((Xs.shape[0], Xt.shape[0]))) + assert_equal(clf.mapping_.shape, ((Xs.shape[1] + 1, Xt.shape[1]))) + + # test margin constraints + mu_s = unif(ns) + mu_t = unif(nt) + assert_allclose(np.sum(clf.coupling_, axis=0), mu_t, rtol=1e-3, atol=1e-3) + assert_allclose(np.sum(clf.coupling_, axis=1), mu_s, rtol=1e-3, atol=1e-3) + + # test transform + transp_Xs = clf.transform(Xs=Xs) + assert_equal(transp_Xs.shape, Xs.shape) + + transp_Xs_new = clf.transform(Xs_new) + + # check that the oos method is working + assert_equal(transp_Xs_new.shape, Xs_new.shape) + + ########################################################################## + # kernel == gaussian mapping tests + ########################################################################## + + # check computation and dimensions if bias == False + clf = ot.da.MappingTransport(kernel="gaussian", bias=False) + clf.fit(Xs=Xs, Xt=Xt) + + assert_equal(clf.coupling_.shape, ((Xs.shape[0], Xt.shape[0]))) + assert_equal(clf.mapping_.shape, ((Xs.shape[0], Xt.shape[1]))) + + # test margin constraints + mu_s = unif(ns) + mu_t = unif(nt) + assert_allclose(np.sum(clf.coupling_, axis=0), mu_t, rtol=1e-3, atol=1e-3) + assert_allclose(np.sum(clf.coupling_, axis=1), mu_s, rtol=1e-3, atol=1e-3) + + # test transform + transp_Xs = clf.transform(Xs=Xs) + assert_equal(transp_Xs.shape, Xs.shape) + + transp_Xs_new = clf.transform(Xs_new) + + # check that the oos method is working + assert_equal(transp_Xs_new.shape, Xs_new.shape) + + # check computation and dimensions if bias == True + clf = ot.da.MappingTransport(kernel="gaussian", bias=True) + clf.fit(Xs=Xs, Xt=Xt) + assert_equal(clf.coupling_.shape, ((Xs.shape[0], Xt.shape[0]))) + assert_equal(clf.mapping_.shape, ((Xs.shape[0] + 1, Xt.shape[1]))) + + # test margin constraints + mu_s = unif(ns) + mu_t = unif(nt) + assert_allclose(np.sum(clf.coupling_, axis=0), mu_t, rtol=1e-3, atol=1e-3) + assert_allclose(np.sum(clf.coupling_, axis=1), mu_s, rtol=1e-3, atol=1e-3) + + # test transform + transp_Xs = clf.transform(Xs=Xs) + assert_equal(transp_Xs.shape, Xs.shape) + + transp_Xs_new = clf.transform(Xs_new) + + # check that the oos method is working + assert_equal(transp_Xs_new.shape, Xs_new.shape) + + # check everything runs well with log=True + clf = ot.da.MappingTransport(kernel="gaussian", log=True) + clf.fit(Xs=Xs, Xt=Xt) + assert len(clf.log_.keys()) != 0 + + +def test_otda(): + + n_samples = 150 # nb samples + np.random.seed(0) + + xs, ys = ot.datasets.get_data_classif('3gauss', n_samples) + xt, yt = ot.datasets.get_data_classif('3gauss2', n_samples) + + a, b = ot.unif(n_samples), ot.unif(n_samples) + + # LP problem + da_emd = ot.da.OTDA() # init class + da_emd.fit(xs, xt) # fit distributions + da_emd.interp() # interpolation of source samples + da_emd.predict(xs) # interpolation of source samples + + np.testing.assert_allclose(a, np.sum(da_emd.G, 1)) + np.testing.assert_allclose(b, np.sum(da_emd.G, 0)) + + # sinkhorn regularization + lambd = 1e-1 + da_entrop = ot.da.OTDA_sinkhorn() + da_entrop.fit(xs, xt, reg=lambd) + da_entrop.interp() + da_entrop.predict(xs) + + np.testing.assert_allclose(a, np.sum(da_entrop.G, 1), rtol=1e-3, atol=1e-3) + np.testing.assert_allclose(b, np.sum(da_entrop.G, 0), rtol=1e-3, atol=1e-3) + + # 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) + da_lpl1.interp() + da_lpl1.predict(xs) + + np.testing.assert_allclose(a, np.sum(da_lpl1.G, 1), rtol=1e-3, atol=1e-3) + np.testing.assert_allclose(b, np.sum(da_lpl1.G, 0), rtol=1e-3, atol=1e-3) + + # 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) + da_l1l2.interp() + da_l1l2.predict(xs) + + np.testing.assert_allclose(a, np.sum(da_l1l2.G, 1), rtol=1e-3, atol=1e-3) + np.testing.assert_allclose(b, np.sum(da_l1l2.G, 0), rtol=1e-3, atol=1e-3) + + # linear mapping + da_emd = ot.da.OTDA_mapping_linear() # init class + da_emd.fit(xs, xt, numItermax=10) # fit distributions + da_emd.predict(xs) # interpolation of source samples + + # nonlinear mapping + da_emd = ot.da.OTDA_mapping_kernel() # init class + da_emd.fit(xs, xt, numItermax=10) # fit distributions + da_emd.predict(xs) # interpolation of source samples + + +# if __name__ == "__main__": + +# test_sinkhorn_transport_class() +# test_emd_transport_class() +# test_sinkhorn_l1l2_transport_class() +# test_sinkhorn_lpl1_transport_class() +# test_mapping_transport_class() diff --git a/test/test_dr.py b/test/test_dr.py new file mode 100644 index 0000000..915012d --- /dev/null +++ b/test/test_dr.py @@ -0,0 +1,59 @@ +"""Tests for module dr on Dimensionality Reduction """ + +# Author: Remi Flamary <remi.flamary@unice.fr> +# +# License: MIT License + +import numpy as np +import ot +import pytest + +try: # test if autograd and pymanopt are installed + import ot.dr + nogo = False +except ImportError: + nogo = True + + +@pytest.mark.skipif(nogo, reason="Missing modules (autograd or pymanopt)") +def test_fda(): + + n_samples = 90 # nb samples in source and target datasets + np.random.seed(0) + + # generate gaussian dataset + xs, ys = ot.datasets.get_data_classif('gaussrot', n_samples) + + n_features_noise = 8 + + xs = np.hstack((xs, np.random.randn(n_samples, n_features_noise))) + + p = 1 + + Pfda, projfda = ot.dr.fda(xs, ys, p) + + projfda(xs) + + np.testing.assert_allclose(np.sum(Pfda**2, 0), np.ones(p)) + + +@pytest.mark.skipif(nogo, reason="Missing modules (autograd or pymanopt)") +def test_wda(): + + n_samples = 100 # nb samples in source and target datasets + np.random.seed(0) + + # generate gaussian dataset + xs, ys = ot.datasets.get_data_classif('gaussrot', n_samples) + + n_features_noise = 8 + + xs = np.hstack((xs, np.random.randn(n_samples, n_features_noise))) + + p = 2 + + Pwda, projwda = ot.dr.wda(xs, ys, p, maxiter=10) + + projwda(xs) + + np.testing.assert_allclose(np.sum(Pwda**2, 0), np.ones(p)) diff --git a/test/test_emd_multi.py b/test/test_emd_multi.py deleted file mode 100644 index ee0a20e..0000000 --- a/test/test_emd_multi.py +++ /dev/null @@ -1,48 +0,0 @@ -#!/usr/bin/env python2 -# -*- coding: utf-8 -*- -""" -Created on Fri Mar 10 09:56:06 2017 - -@author: rflamary -""" - -import numpy as np -import pylab as pl -import ot - -from ot.datasets import get_1D_gauss as gauss -reload(ot.lp) - -#%% parameters - -n=5000 # nb bins - -# bin positions -x=np.arange(n,dtype=np.float64) - -# Gaussian distributions -a=gauss(n,m=20,s=5) # m= mean, s= std - -ls= range(20,1000,10) -nb=len(ls) -b=np.zeros((n,nb)) -for i in range(nb): - b[:,i]=gauss(n,m=ls[i],s=10) - -# loss matrix -M=ot.dist(x.reshape((n,1)),x.reshape((n,1))) -#M/=M.max() - -#%% - -print('Computing {} EMD '.format(nb)) - -# emd loss 1 proc -ot.tic() -emd_loss4=ot.emd2(a,b,M,1) -ot.toc('1 proc : {} s') - -# emd loss multipro proc -ot.tic() -emd_loss4=ot.emd2(a,b,M) -ot.toc('multi proc : {} s') diff --git a/test/test_gpu.py b/test/test_gpu.py new file mode 100644 index 0000000..615c2a7 --- /dev/null +++ b/test/test_gpu.py @@ -0,0 +1,79 @@ +"""Tests for module gpu for gpu acceleration """ + +# Author: Remi Flamary <remi.flamary@unice.fr> +# +# License: MIT License + +import numpy as np +import ot +import time +import pytest + +try: # test if cudamat installed + import ot.gpu + nogpu = False +except ImportError: + nogpu = True + + +@pytest.mark.skipif(nogpu, reason="No GPU available") +def test_gpu_sinkhorn(): + + rng = np.random.RandomState(0) + + def describe_res(r): + print("min:{:.3E}, max::{:.3E}, mean::{:.3E}, std::{:.3E}".format( + np.min(r), np.max(r), np.mean(r), np.std(r))) + + for n_samples in [50, 100, 500, 1000]: + print(n_samples) + a = rng.rand(n_samples // 4, 100) + b = rng.rand(n_samples, 100) + time1 = time.time() + transport = ot.da.OTDA_sinkhorn() + transport.fit(a, b) + G1 = transport.G + time2 = time.time() + transport = ot.gpu.da.OTDA_sinkhorn() + transport.fit(a, b) + G2 = transport.G + time3 = time.time() + print("Normal sinkhorn, time: {:6.2f} sec ".format(time2 - time1)) + describe_res(G1) + print(" GPU sinkhorn, time: {:6.2f} sec ".format(time3 - time2)) + describe_res(G2) + + np.testing.assert_allclose(G1, G2, rtol=1e-5, atol=1e-5) + + +@pytest.mark.skipif(nogpu, reason="No GPU available") +def test_gpu_sinkhorn_lpl1(): + + rng = np.random.RandomState(0) + + def describe_res(r): + print("min:{:.3E}, max:{:.3E}, mean:{:.3E}, std:{:.3E}" + .format(np.min(r), np.max(r), np.mean(r), np.std(r))) + + for n_samples in [50, 100, 500]: + print(n_samples) + a = rng.rand(n_samples // 4, 100) + labels_a = np.random.randint(10, size=(n_samples // 4)) + b = rng.rand(n_samples, 100) + time1 = time.time() + transport = ot.da.OTDA_lpl1() + transport.fit(a, labels_a, b) + G1 = transport.G + time2 = time.time() + transport = ot.gpu.da.OTDA_lpl1() + transport.fit(a, labels_a, b) + G2 = transport.G + time3 = time.time() + print("Normal sinkhorn lpl1, time: {:6.2f} sec ".format( + time2 - time1)) + describe_res(G1) + print(" GPU sinkhorn lpl1, time: {:6.2f} sec ".format( + time3 - time2)) + describe_res(G2) + + np.testing.assert_allclose(G1, G2, rtol=1e-5, atol=1e-5) diff --git a/test/test_gpu_sinkhorn.py b/test/test_gpu_sinkhorn.py deleted file mode 100644 index bfa2cd2..0000000 --- a/test/test_gpu_sinkhorn.py +++ /dev/null @@ -1,26 +0,0 @@ -import ot -import numpy as np -import time -import ot.gpu - -def describeRes(r): - print("min:{:.3E}, max::{:.3E}, mean::{:.3E}, std::{:.3E}".format(np.min(r),np.max(r),np.mean(r),np.std(r))) - - -for n in [5000, 10000, 15000, 20000]: - print(n) - a = np.random.rand(n // 4, 100) - b = np.random.rand(n, 100) - time1 = time.time() - transport = ot.da.OTDA_sinkhorn() - transport.fit(a, b) - G1 = transport.G - time2 = time.time() - transport = ot.gpu.da.OTDA_sinkhorn() - transport.fit(a, b) - G2 = transport.G - time3 = time.time() - print("Normal sinkhorn, time: {:6.2f} sec ".format(time2 - time1)) - describeRes(G1) - print(" GPU sinkhorn, time: {:6.2f} sec ".format(time3 - time2)) - describeRes(G2)
\ No newline at end of file diff --git a/test/test_gpu_sinkhorn_lpl1.py b/test/test_gpu_sinkhorn_lpl1.py deleted file mode 100644 index e6cdd31..0000000 --- a/test/test_gpu_sinkhorn_lpl1.py +++ /dev/null @@ -1,28 +0,0 @@ -import ot -import numpy as np -import time -import ot.gpu - - -def describeRes(r): - print("min:{:.3E}, max:{:.3E}, mean:{:.3E}, std:{:.3E}" - .format(np.min(r), np.max(r), np.mean(r), np.std(r))) - -for n in [5000, 10000, 15000, 20000]: - print(n) - a = np.random.rand(n // 4, 100) - labels_a = np.random.randint(10, size=(n // 4)) - b = np.random.rand(n, 100) - time1 = time.time() - transport = ot.da.OTDA_lpl1() - transport.fit(a, labels_a, b) - G1 = transport.G - time2 = time.time() - transport = ot.gpu.da.OTDA_lpl1() - transport.fit(a, labels_a, b) - G2 = transport.G - time3 = time.time() - print("Normal sinkhorn lpl1, time: {:6.2f} sec ".format(time2 - time1)) - describeRes(G1) - print(" GPU sinkhorn lpl1, time: {:6.2f} sec ".format(time3 - time2)) - describeRes(G2) diff --git a/test/test_load_module.py b/test/test_load_module.py deleted file mode 100644 index a04c5df..0000000 --- a/test/test_load_module.py +++ /dev/null @@ -1,10 +0,0 @@ - - -import ot -import doctest - -# test lp solver -doctest.testmod(ot.lp,verbose=True) - -# test bregman solver -doctest.testmod(ot.bregman,verbose=True) diff --git a/test/test_optim.py b/test/test_optim.py new file mode 100644 index 0000000..69496a5 --- /dev/null +++ b/test/test_optim.py @@ -0,0 +1,67 @@ +"""Tests for module optim fro OT optimization """ + +# Author: Remi Flamary <remi.flamary@unice.fr> +# +# License: MIT License + +import numpy as np +import ot + + +def test_conditional_gradient(): + + n_bins = 100 # nb bins + np.random.seed(0) + # bin positions + x = np.arange(n_bins, dtype=np.float64) + + # Gaussian distributions + a = ot.datasets.get_1D_gauss(n_bins, m=20, s=5) # m= mean, s= std + b = ot.datasets.get_1D_gauss(n_bins, m=60, s=10) + + # loss matrix + M = ot.dist(x.reshape((n_bins, 1)), x.reshape((n_bins, 1))) + M /= M.max() + + def f(G): + return 0.5 * np.sum(G**2) + + def df(G): + return G + + reg = 1e-1 + + G, log = ot.optim.cg(a, b, M, reg, f, df, verbose=True, log=True) + + np.testing.assert_allclose(a, G.sum(1)) + np.testing.assert_allclose(b, G.sum(0)) + + +def test_generalized_conditional_gradient(): + + n_bins = 100 # nb bins + np.random.seed(0) + # bin positions + x = np.arange(n_bins, dtype=np.float64) + + # Gaussian distributions + a = ot.datasets.get_1D_gauss(n_bins, m=20, s=5) # m= mean, s= std + b = ot.datasets.get_1D_gauss(n_bins, m=60, s=10) + + # loss matrix + M = ot.dist(x.reshape((n_bins, 1)), x.reshape((n_bins, 1))) + M /= M.max() + + def f(G): + return 0.5 * np.sum(G**2) + + def df(G): + return G + + reg1 = 1e-3 + reg2 = 1e-1 + + G, log = ot.optim.gcg(a, b, M, reg1, reg2, f, df, verbose=True, log=True) + + np.testing.assert_allclose(a, G.sum(1), atol=1e-05) + np.testing.assert_allclose(b, G.sum(0), atol=1e-05) diff --git a/test/test_ot.py b/test/test_ot.py new file mode 100644 index 0000000..acd8718 --- /dev/null +++ b/test/test_ot.py @@ -0,0 +1,102 @@ +"""Tests for main module ot """ + +# Author: Remi Flamary <remi.flamary@unice.fr> +# +# License: MIT License + +import numpy as np +import ot + + +def test_doctest(): + + import doctest + + # test lp solver + doctest.testmod(ot.lp, verbose=True) + + # test bregman solver + doctest.testmod(ot.bregman, verbose=True) + + +def test_emd_emd2(): + # test emd and emd2 for simple identity + n = 100 + rng = np.random.RandomState(0) + + x = rng.randn(n, 2) + u = ot.utils.unif(n) + + M = ot.dist(x, x) + + G = ot.emd(u, u, M) + + # check G is identity + np.testing.assert_allclose(G, np.eye(n) / n) + # check constratints + np.testing.assert_allclose(u, G.sum(1)) # cf convergence sinkhorn + np.testing.assert_allclose(u, G.sum(0)) # cf convergence sinkhorn + + w = ot.emd2(u, u, M) + # check loss=0 + np.testing.assert_allclose(w, 0) + + +def test_emd_empty(): + # test emd and emd2 for simple identity + n = 100 + rng = np.random.RandomState(0) + + x = rng.randn(n, 2) + u = ot.utils.unif(n) + + M = ot.dist(x, x) + + G = ot.emd([], [], M) + + # check G is identity + np.testing.assert_allclose(G, np.eye(n) / n) + # check constratints + np.testing.assert_allclose(u, G.sum(1)) # cf convergence sinkhorn + np.testing.assert_allclose(u, G.sum(0)) # cf convergence sinkhorn + + w = ot.emd2([], [], M) + # check loss=0 + np.testing.assert_allclose(w, 0) + + +def test_emd2_multi(): + + from ot.datasets import get_1D_gauss as gauss + + n = 1000 # nb bins + + # bin positions + x = np.arange(n, dtype=np.float64) + + # Gaussian distributions + a = gauss(n, m=20, s=5) # m= mean, s= std + + ls = np.arange(20, 1000, 20) + nb = len(ls) + b = np.zeros((n, nb)) + for i in range(nb): + b[:, i] = gauss(n, m=ls[i], s=10) + + # loss matrix + M = ot.dist(x.reshape((n, 1)), x.reshape((n, 1))) + # M/=M.max() + + print('Computing {} EMD '.format(nb)) + + # emd loss 1 proc + ot.tic() + emd1 = ot.emd2(a, b, M, 1) + ot.toc('1 proc : {} s') + + # emd loss multipro proc + ot.tic() + emdn = ot.emd2(a, b, M) + ot.toc('multi proc : {} s') + + np.testing.assert_allclose(emd1, emdn) diff --git a/test/test_plot.py b/test/test_plot.py new file mode 100644 index 0000000..f7debee --- /dev/null +++ b/test/test_plot.py @@ -0,0 +1,49 @@ +"""Tests for module plot for visualization """ + +# Author: Remi Flamary <remi.flamary@unice.fr> +# +# License: MIT License + +import numpy as np +import matplotlib +matplotlib.use('Agg') + + +def test_plot1D_mat(): + + import ot + + n_bins = 100 # nb bins + + # bin positions + x = np.arange(n_bins, dtype=np.float64) + + # Gaussian distributions + a = ot.datasets.get_1D_gauss(n_bins, m=20, s=5) # m= mean, s= std + b = ot.datasets.get_1D_gauss(n_bins, m=60, s=10) + + # loss matrix + M = ot.dist(x.reshape((n_bins, 1)), x.reshape((n_bins, 1))) + M /= M.max() + + ot.plot.plot1D_mat(a, b, M, 'Cost matrix M') + + +def test_plot2D_samples_mat(): + + import ot + + n_bins = 50 # nb samples + + mu_s = np.array([0, 0]) + cov_s = np.array([[1, 0], [0, 1]]) + + mu_t = np.array([4, 4]) + cov_t = np.array([[1, -.8], [-.8, 1]]) + + xs = ot.datasets.get_2D_samples_gauss(n_bins, mu_s, cov_s) + xt = ot.datasets.get_2D_samples_gauss(n_bins, mu_t, cov_t) + + G = 1.0 * (np.random.rand(n_bins, n_bins) < 0.01) + + ot.plot.plot2D_samples_mat(xs, xt, G, thr=1e-5) diff --git a/test/test_utils.py b/test/test_utils.py new file mode 100644 index 0000000..1bd37cd --- /dev/null +++ b/test/test_utils.py @@ -0,0 +1,125 @@ +"""Tests for module utils for timing and parallel computation """ + +# Author: Remi Flamary <remi.flamary@unice.fr> +# +# License: MIT License + + +import ot +import numpy as np + + +def test_parmap(): + + n = 100 + + def f(i): + return 1.0 * i * i + + a = np.arange(n) + + l1 = list(map(f, a)) + + l2 = list(ot.utils.parmap(f, a)) + + np.testing.assert_allclose(l1, l2) + + +def test_tic_toc(): + + import time + + ot.tic() + time.sleep(0.5) + t = ot.toc() + t2 = ot.toq() + + # test timing + np.testing.assert_allclose(0.5, t, rtol=1e-2, atol=1e-2) + + # test toc vs toq + np.testing.assert_allclose(t, t2, rtol=1e-2, atol=1e-2) + + +def test_kernel(): + + n = 100 + + x = np.random.randn(n, 2) + + K = ot.utils.kernel(x, x) + + # gaussian kernel has ones on the diagonal + np.testing.assert_allclose(np.diag(K), np.ones(n)) + + +def test_unif(): + + n = 100 + + u = ot.unif(n) + + np.testing.assert_allclose(1, np.sum(u)) + + +def test_dist(): + + n = 100 + + x = np.random.randn(n, 2) + + D = np.zeros((n, n)) + for i in range(n): + for j in range(n): + D[i, j] = np.sum(np.square(x[i, :] - x[j, :])) + + D2 = ot.dist(x, x) + D3 = ot.dist(x) + + # dist shoul return squared euclidean + np.testing.assert_allclose(D, D2) + np.testing.assert_allclose(D, D3) + + +def test_dist0(): + + n = 100 + M = ot.utils.dist0(n, method='lin_square') + + # dist0 default to linear sampling with quadratic loss + np.testing.assert_allclose(M[0, -1], (n - 1) * (n - 1)) + + +def test_dots(): + + n1, n2, n3, n4 = 100, 50, 200, 100 + + A = np.random.randn(n1, n2) + B = np.random.randn(n2, n3) + C = np.random.randn(n3, n4) + + X1 = ot.utils.dots(A, B, C) + + X2 = A.dot(B.dot(C)) + + np.testing.assert_allclose(X1, X2) + + +def test_clean_zeros(): + + n = 100 + nz = 50 + nz2 = 20 + u1 = ot.unif(n) + u1[:nz] = 0 + u1 = u1 / u1.sum() + u2 = ot.unif(n) + u2[:nz2] = 0 + u2 = u2 / u2.sum() + + M = ot.utils.dist0(n) + + a, b, M2 = ot.utils.clean_zeros(u1, u2, M) + + assert len(a) == n - nz + assert len(b) == n - nz2 |