Initial commit of new Nomad Analytics Toolkit Tutorial

.gitignore

+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+wheels/
+pip-wheel-metadata/
+share/python-wheels/
+*.egg-info/
+.installed.cfg
+*.egg
+MANIFEST
+
+# PyInstaller
+#  Usually these files are written by a python script from a template
+#  before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.nox/
+.coverage
+.coverage.*
+nosetests.xml
+coverage.xml
+*.cover
+*.py,cover
+.hypothesis/
+.pytest_cache/
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+local_settings.py
+db.sqlite3
+db.sqlite3-journal
+
+# Flask stuff:
+instance/
+.webassets-cache
+
+# Scrapy stuff:
+.scrapy
+
+# Sphinx documentation
+docs/_build/
+
+# PyBuilder
+target/
+
+# Jupyter Notebook
+.ipynb_checkpoints
+
+# IPython
+profile_default/
+ipython_config.py
+
+# pyenv
+.python-version
+
+# pipenv
+#   According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control.
+#   However, in case of collaboration, if having platform-specific dependencies or dependencies
+#   having no cross-platform support, pipenv may install dependencies that don't work, or not
+#   install all needed dependencies.
+#Pipfile.lock
+
+# PEP 582; used by e.g. github.com/David-OConnor/pyflow
+__pypackages__/
+
+# Celery stuff
+celerybeat-schedule
+celerybeat.pid
+
+# SageMath parsed files
+*.sage.py
+
+# Environments
+.env
+.venv
+env/
+venv/
+ENV/
+env.bak/
+venv.bak/
+
+# Spyder project settings
+.spyderproject
+.spyproject
+
+# Rope project settings
+.ropeproject
+
+# mkdocs documentation
+/site
+
+# mypy
+.mypy_cache/
+.dmypy.json
+dmypy.json
+
+# Pyre type checker
+.pyre/
+
+# ctags
+tags
+.vscode
+
+# vim swap files
+*.swp
+
+# osx
+.DS_Store
+
+# Extra folder
+.idea

CHANGELOG.md

+# v1.0.0
+## 14/01/2020
+
+1. [](#new)
+    * ChangeLog started...

LICENSE

+
+                                 Apache License
+                           Version 2.0, January 2004
+                        http://www.apache.org/licenses/
+
+   TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION
+
+   1. Definitions.
+
+      "License" shall mean the terms and conditions for use, reproduction,
+      and distribution as defined by Sections 1 through 9 of this document.
+
+      "Licensor" shall mean the copyright owner or entity authorized by
+      the copyright owner that is granting the License.
+
+      "Legal Entity" shall mean the union of the acting entity and all
+      other entities that control, are controlled by, or are under common
+      control with that entity. For the purposes of this definition,
+      "control" means (i) the power, direct or indirect, to cause the
+      direction or management of such entity, whether by contract or
+      otherwise, or (ii) ownership of fifty percent (50%) or more of the
+      outstanding shares, or (iii) beneficial ownership of such entity.
+
+      "You" (or "Your") shall mean an individual or Legal Entity
+      exercising permissions granted by this License.
+
+      "Source" form shall mean the preferred form for making modifications,
+      including but not limited to software source code, documentation
+      source, and configuration files.
+
+      "Object" form shall mean any form resulting from mechanical
+      transformation or translation of a Source form, including but
+      not limited to compiled object code, generated documentation,
+      and conversions to other media types.
+
+      "Work" shall mean the work of authorship, whether in Source or
+      Object form, made available under the License, as indicated by a
+      copyright notice that is included in or attached to the work
+      (an example is provided in the Appendix below).
+
+      "Derivative Works" shall mean any work, whether in Source or Object
+      form, that is based on (or derived from) the Work and for which the
+      editorial revisions, annotations, elaborations, or other modifications
+      represent, as a whole, an original work of authorship. For the purposes
+      of this License, Derivative Works shall not include works that remain
+      separable from, or merely link (or bind by name) to the interfaces of,
+      the Work and Derivative Works thereof.
+
+      "Contribution" shall mean any work of authorship, including
+      the original version of the Work and any modifications or additions
+      to that Work or Derivative Works thereof, that is intentionally
+      submitted to Licensor for inclusion in the Work by the copyright owner
+      or by an individual or Legal Entity authorized to submit on behalf of
+      the copyright owner. For the purposes of this definition, "submitted"
+      means any form of electronic, verbal, or written communication sent
+      to the Licensor or its representatives, including but not limited to
+      communication on electronic mailing lists, source code control systems,
+      and issue tracking systems that are managed by, or on behalf of, the
+      Licensor for the purpose of discussing and improving the Work, but
+      excluding communication that is conspicuously marked or otherwise
+      designated in writing by the copyright owner as "Not a Contribution."
+
+      "Contributor" shall mean Licensor and any individual or Legal Entity
+      on behalf of whom a Contribution has been received by Licensor and
+      subsequently incorporated within the Work.
+
+   2. Grant of Copyright License. Subject to the terms and conditions of
+      this License, each Contributor hereby grants to You a perpetual,
+      worldwide, non-exclusive, no-charge, royalty-free, irrevocable
+      copyright license to reproduce, prepare Derivative Works of,
+      publicly display, publicly perform, sublicense, and distribute the
+      Work and such Derivative Works in Source or Object form.
+
+   3. Grant of Patent License. Subject to the terms and conditions of
+      this License, each Contributor hereby grants to You a perpetual,
+      worldwide, non-exclusive, no-charge, royalty-free, irrevocable
+      (except as stated in this section) patent license to make, have made,
+      use, offer to sell, sell, import, and otherwise transfer the Work,
+      where such license applies only to those patent claims licensable
+      by such Contributor that are necessarily infringed by their
+      Contribution(s) alone or by combination of their Contribution(s)
+      with the Work to which such Contribution(s) was submitted. If You
+      institute patent litigation against any entity (including a
+      cross-claim or counterclaim in a lawsuit) alleging that the Work
+      or a Contribution incorporated within the Work constitutes direct
+      or contributory patent infringement, then any patent licenses
+      granted to You under this License for that Work shall terminate
+      as of the date such litigation is filed.
+
+   4. Redistribution. You may reproduce and distribute copies of the
+      Work or Derivative Works thereof in any medium, with or without
+      modifications, and in Source or Object form, provided that You
+      meet the following conditions:
+
+      (a) You must give any other recipients of the Work or
+          Derivative Works a copy of this License; and
+
+      (b) You must cause any modified files to carry prominent notices
+          stating that You changed the files; and
+
+      (c) You must retain, in the Source form of any Derivative Works
+          that You distribute, all copyright, patent, trademark, and
+          attribution notices from the Source form of the Work,
+          excluding those notices that do not pertain to any part of
+          the Derivative Works; and
+
+      (d) If the Work includes a "NOTICE" text file as part of its
+          distribution, then any Derivative Works that You distribute must
+          include a readable copy of the attribution notices contained
+          within such NOTICE file, excluding those notices that do not
+          pertain to any part of the Derivative Works, in at least one
+          of the following places: within a NOTICE text file distributed
+          as part of the Derivative Works; within the Source form or
+          documentation, if provided along with the Derivative Works; or,
+          within a display generated by the Derivative Works, if and
+          wherever such third-party notices normally appear. The contents
+          of the NOTICE file are for informational purposes only and
+          do not modify the License. You may add Your own attribution
+          notices within Derivative Works that You distribute, alongside
+          or as an addendum to the NOTICE text from the Work, provided
+          that such additional attribution notices cannot be construed
+          as modifying the License.
+
+      You may add Your own copyright statement to Your modifications and
+      may provide additional or different license terms and conditions
+      for use, reproduction, or distribution of Your modifications, or
+      for any such Derivative Works as a whole, provided Your use,
+      reproduction, and distribution of the Work otherwise complies with
+      the conditions stated in this License.
+
+   5. Submission of Contributions. Unless You explicitly state otherwise,
+      any Contribution intentionally submitted for inclusion in the Work
+      by You to the Licensor shall be under the terms and conditions of
+      this License, without any additional terms or conditions.
+      Notwithstanding the above, nothing herein shall supersede or modify
+      the terms of any separate license agreement you may have executed
+      with Licensor regarding such Contributions.
+
+   6. Trademarks. This License does not grant permission to use the trade
+      names, trademarks, service marks, or product names of the Licensor,
+      except as required for reasonable and customary use in describing the
+      origin of the Work and reproducing the content of the NOTICE file.
+
+   7. Disclaimer of Warranty. Unless required by applicable law or
+      agreed to in writing, Licensor provides the Work (and each
+      Contributor provides its Contributions) on an "AS IS" BASIS,
+      WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
+      implied, including, without limitation, any warranties or conditions
+      of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A
+      PARTICULAR PURPOSE. You are solely responsible for determining the
+      appropriateness of using or redistributing the Work and assume any
+      risks associated with Your exercise of permissions under this License.
+
+   8. Limitation of Liability. In no event and under no legal theory,
+      whether in tort (including negligence), contract, or otherwise,
+      unless required by applicable law (such as deliberate and grossly
+      negligent acts) or agreed to in writing, shall any Contributor be
+      liable to You for damages, including any direct, indirect, special,
+      incidental, or consequential damages of any character arising as a
+      result of this License or out of the use or inability to use the
+      Work (including but not limited to damages for loss of goodwill,
+      work stoppage, computer failure or malfunction, or any and all
+      other commercial damages or losses), even if such Contributor
+      has been advised of the possibility of such damages.
+
+   9. Accepting Warranty or Additional Liability. While redistributing
+      the Work or Derivative Works thereof, You may choose to offer,
+      and charge a fee for, acceptance of support, warranty, indemnity,
+      or other liability obligations and/or rights consistent with this
+      License. However, in accepting such obligations, You may act only
+      on Your own behalf and on Your sole responsibility, not on behalf
+      of any other Contributor, and only if You agree to indemnify,
+      defend, and hold each Contributor harmless for any liability
+      incurred by, or claims asserted against, such Contributor by reason
+      of your accepting any such warranty or additional liability.
+
+   END OF TERMS AND CONDITIONS
+
+   APPENDIX: How to apply the Apache License to your work.
+
+      To apply the Apache License to your work, attach the following
+      boilerplate notice, with the fields enclosed by brackets "[]"
+      replaced with your own identifying information. (Don't include
+      the brackets!)  The text should be enclosed in the appropriate
+      comment syntax for the file format. We also recommend that a
+      file or class name and description of purpose be included on the
+      same "printed page" as the copyright notice for easier
+      identification within third-party archives.
+
+   Copyright [yyyy] [name of copyright owner]
+
+   Licensed under the Apache License, Version 2.0 (the "License");
+   you may not use this file except in compliance with the License.
+   You may obtain a copy of the License at
+
+       http://www.apache.org/licenses/LICENSE-2.0
+
+   Unless required by applicable law or agreed to in writing, software
+   distributed under the License is distributed on an "AS IS" BASIS,
+   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+   See the License for the specific language governing permissions and
+   limitations under the License.

NOTICE

+Max Planck Society - Fritz Haber Institute
+
+Copyright (c) 2018+ Fritz Haber Institute of the Max Planck Society
+Contributors: Benjamin Regler <regler@fhi-berlin.mpg.de> - https://github.com/Sommerregen
+              (Big-Data Analytics for Materials Science, Luca M. Ghiringhelli)
+Affiliation:  Theory Department, Faradayweg 4-6, 14195 Berlin, Germany
+URL:          http://th.fhi-berlin.mpg.de/site/ 
+Citation:     <>
+              
+This software is part of the publication <>. An interactive online tutorial can be found in the NOMAD Analytics Toolkit <>.
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+   http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.

README.md

+# Total cumulative mutual information (NOMAD Analytics Toolkit)
+
+The identification of relevant features, i.e., the driving variables that determine a process or the property of a system, is an essential part of the analysis of data sets whose entries are described by a large number of variables. The preferred measure for quantifying the relevance of nonlinear statistical dependencies is mutual information, which requires as input probability distributions. Probability distributions cannot be reliably sampled and estimated from limited data, especially for real-valued data samples such as lengths or energies.
+
+This interactive notebook includes the original implementation of total cumulative mutual information (TCMI) to reproduce the main results presented in the publication:
+
+<div style="padding: 1ex; margin-top: 1ex; margin-bottom: 1ex; border-style: dotted; border-width: 1pt; border-color: blue; border-radius: 3px;">
+B. Regler, M. Scheffler, and L. M. Ghiringhelli: "TCMI: a non-parametric mutual-dependence estimator for multivariate continuous distributions"
+</div>
+
+TCMI is a measure of the relevance of mutual dependencies based on cumulative probability distributions. TCMI can be estimated directly from sample data and is a non-parametric, robust and deterministic measure that facilitates comparisons and rankings between feature sets with different cardinality. The ranking induced by TCMI allows for feature selection, i.e. the identification of the set of relevant features that are statistical related to the process or the property of a system, while taking into account the number of data samples as well as the cardinality of the feature subsets.
+
+It is compared to [Cumulative mutual information (CMI)](https://dx.doi.org/10.1137/1.9781611972832.22), [Multivariate maximal correlation analysis (MAC)](http://proceedings.mlr.press/v32/nguyenc14.html), [Universal dependency analysis (UDS)](https://dx.doi.org/10.1137/1.9781611974348.89), and [Monte Carlo dependency estimation (MCDE)](https://dx.doi.org/10.1145/3335783.3335795).
+
+This repository (notebook and code) is released under the [Apache License, Version 2.0](http://www.apache.org/licenses/). Please see the [LICENSE](LICENSE) file.
+
+**Important notes:**
+<ul style="color: #8b0000; font-style: italic;">
+<li>All comparisons have been computed with the Java package <code>MCDE</code> written in Scala, which is not part of the repository. To use the most recent and maintained implementation, please visit <a href="https://github.com/edouardfouche/MCDE">https://github.com/edouardfouche/MCDE</a> and run all examples with 50,000 iterations.</li>
+<li>For the sake of simplicity, all results have been cached. However, results can be recalculated after adjusting the respective test sections. Depending on the test, the calculation time ranges from minutes to days.</li>
+</ul>
+
+#### Project Information:
+
+**Maintainer:** <a href="https://gitlab.rzg.mpg.de/rbenj" style="color: #808080;" title="Maintainer">&#x1F464; Benjamin Regler</a>
+
+**Status:** <span style="color: #008000;">&#10004; Under active development</span>
+
+## License
+
+Copyright (c) 2018+ Fritz Haber Institute of the Max Planck Society ([Benjamin Regler][github]).
+
+[github]: https://github.com/Sommerregen "GitLab account of Benjamin Regler"

data/octet-binary-compound-semiconductors.csv

+Combination,Delta E,EA(A),EN(A),H(A),IP(A),L(A),rd(A),rp(A),rs(A),EA(B),EN(B),H(B),IP(B),L(B),rd(B),rp(B),rs(B)
+Li+F,-0.0594925,-0.698,3.0135,-2.874,-5.329,-0.978,6.93,1.995,1.652,-4.273,11.8385,-11.294,-19.404,1.251,1.428,0.371,0.406
+Li+Cl,-0.038388,-0.698,3.0135,-2.874,-5.329,-0.978,6.93,1.995,1.652,-3.971,8.9365,-8.7,-13.902,0.574,1.666,0.756,0.679
+Li+Br,-0.032847,-0.698,3.0135,-2.874,-5.329,-0.978,6.93,1.995,1.652,-3.739,8.1945,-8.001,-12.65,0.708,1.869,0.882,0.749
+Li+I,-0.022172,-0.698,3.0135,-2.874,-5.329,-0.978,6.93,1.995,1.652,-3.513,7.385,-7.236,-11.257,0.213,1.722,1.071,0.896
+Be+O,0.430056,0.631,4.414,-5.6,-9.459,-2.098,2.877,1.211,1.078,-3.006,9.7195,-9.197,-16.433,2.541,2.219,0.427,0.462
+Be+S,0.5064495,0.631,4.414,-5.6,-9.459,-2.098,2.877,1.211,1.078,-2.845,7.32,-7.106,-11.795,0.642,2.366,0.847,0.742
+Be+Se,0.49488,0.631,4.414,-5.6,-9.459,-2.098,2.877,1.211,1.078,-2.751,6.8485,-6.654,-10.946,1.316,2.177,0.952,0.798
+Be+Te,0.466361,0.631,4.414,-5.6,-9.459,-2.098,2.877,1.211,1.078,-2.666,6.2665,-6.109,-9.867,0.099,1.827,1.141,0.945
+B+N,1.712778,-0.107,4.1485,-3.715,-8.19,2.248,1.946,0.826,0.805,-1.867,7.726,-7.239,-13.585,3.057,1.54,0.511,0.539
+B+P,1.0196875,-0.107,4.1485,-3.715,-8.19,2.248,1.946,0.826,0.805,-1.92,5.8355,-5.596,-9.751,0.183,1.771,0.966,0.826
+B+As,0.8789035,-0.107,4.1485,-3.715,-8.19,2.248,1.946,0.826,0.805,-1.839,5.5505,-5.341,-9.262,0.064,2.023,1.043,0.847
+C+C,2.6381375,-0.872,5.862,-5.416,-10.852,1.992,1.631,0.63,0.644,-0.872,5.862,-5.416,-10.852,1.992,1.631,0.63,0.644
+Na+F,-0.1458045,-0.716,2.9695,-2.819,-5.223,-0.718,6.566,2.597,1.715,-4.273,11.8385,-11.294,-19.404,1.251,1.428,0.371,0.406
+Na+Cl,-0.1330525,-0.716,2.9695,-2.819,-5.223,-0.718,6.566,2.597,1.715,-3.971,8.9365,-8.7,-13.902,0.574,1.666,0.756,0.679
+Na+Br,-0.1266025,-0.716,2.9695,-2.819,-5.223,-0.718,6.566,2.597,1.715,-3.739,8.1945,-8.001,-12.65,0.708,1.869,0.882,0.749
+Na+I,-0.115454,-0.716,2.9695,-2.819,-5.223,-0.718,6.566,2.597,1.715,-3.513,7.385,-7.236,-11.257,0.213,1.722,1.071,0.896
+Mg+O,-0.177796,0.693,3.672,-4.782,-8.037,-1.358,3.171,1.897,1.33,-3.006,9.7195,-9.197,-16.433,2.541,2.219,0.427,0.462
+Mg+S,-0.086812,0.693,3.672,-4.782,-8.037,-1.358,3.171,1.897,1.33,-2.845,7.32,-7.106,-11.795,0.642,2.366,0.847,0.742
+Mg+Se,-0.055127,0.693,3.672,-4.782,-8.037,-1.358,3.171,1.897,1.33,-2.751,6.8485,-6.654,-10.946,1.316,2.177,0.952,0.798
+Mg+Te,-0.0049405,0.693,3.672,-4.782,-8.037,-1.358,3.171,1.897,1.33,-2.666,6.2665,-6.109,-9.867,0.099,1.827,1.141,0.945
+Al+N,0.0729255,-0.313,3.0465,-2.784,-5.78,0.695,1.939,1.393,1.092,-1.867,7.726,-7.239,-13.585,3.057,1.54,0.511,0.539
+Al+P,0.2188715,-0.313,3.0465,-2.784,-5.78,0.695,1.939,1.393,1.092,-1.92,5.8355,-5.596,-9.751,0.183,1.771,0.966,0.826
+Al+As,0.212353,-0.313,3.0465,-2.784,-5.78,0.695,1.939,1.393,1.092,-1.839,5.5505,-5.341,-9.262,0.064,2.023,1.043,0.847
+Al+Sb,0.1503865,-0.313,3.0465,-2.784,-5.78,0.695,1.939,1.393,1.092,-1.847,5.1575,-4.991,-8.468,0.105,2.065,1.232,1.001
+Si+C,0.6686675,-0.993,4.3755,-4.163,-7.758,0.44,1.89,1.134,0.938,-0.872,5.862,-5.416,-10.852,1.992,1.631,0.63,0.644
+Si+Si,0.275475,-0.993,4.3755,-4.163,-7.758,0.44,1.89,1.134,0.938,-0.993,4.3755,-4.163,-7.758,0.44,1.89,1.134,0.938
+K+F,-0.146464,-0.621,2.527,-2.426,-4.433,-0.697,1.785,2.443,2.128,-4.273,11.8385,-11.294,-19.404,1.251,1.428,0.371,0.406
+K+Cl,-0.164521,-0.621,2.527,-2.426,-4.433,-0.697,1.785,2.443,2.128,-3.971,8.9365,-8.7,-13.902,0.574,1.666,0.756,0.679
+K+Br,-0.1664765,-0.621,2.527,-2.426,-4.433,-0.697,1.785,2.443,2.128,-3.739,8.1945,-8.001,-12.65,0.708,1.869,0.882,0.749
+K+I,-0.168284,-0.621,2.527,-2.426,-4.433,-0.697,1.785,2.443,2.128,-3.513,7.385,-7.236,-11.257,0.213,1.722,1.071,0.896
+Ca+O,-0.2660115,0.304,3.062,-3.864,-6.428,-2.133,0.679,2.324,1.757,-3.006,9.7195,-9.197,-16.433,2.541,2.219,0.427,0.462
+Ca+S,-0.3692695,0.304,3.062,-3.864,-6.428,-2.133,0.679,2.324,1.757,-2.845,7.32,-7.106,-11.795,0.642,2.366,0.847,0.742
+Ca+Se,-0.3613305,0.304,3.062,-3.864,-6.428,-2.133,0.679,2.324,1.757,-2.751,6.8485,-6.654,-10.946,1.316,2.177,0.952,0.798
+Ca+Te,-0.350466,0.304,3.062,-3.864,-6.428,-2.133,0.679,2.324,1.757,-2.666,6.2665,-6.109,-9.867,0.099,1.827,1.141,0.945
+Cu+F,-0.0187695,-1.638,5.0135,-4.856,-8.389,-0.641,2.576,1.68,1.197,-4.273,11.8385,-11.294,-19.404,1.251,1.428,0.371,0.406
+Cu+Cl,0.156381,-1.638,5.0135,-4.856,-8.389,-0.641,2.576,1.68,1.197,-3.971,8.9365,-8.7,-13.902,0.574,1.666,0.756,0.679
+Cu+Br,0.1520705,-1.638,5.0135,-4.856,-8.389,-0.641,2.576,1.68,1.197,-3.739,8.1945,-8.001,-12.65,0.708,1.869,0.882,0.749
+Cu+I,0.203405,-1.638,5.0135,-4.856,-8.389,-0.641,2.576,1.68,1.197,-3.513,7.385,-7.236,-11.257,0.213,1.722,1.071,0.896
+Zn+O,0.1018155,1.081,4.5275,-6.217,-10.136,-1.194,2.254,1.547,1.099,-3.006,9.7195,-9.197,-16.433,2.541,2.219,0.427,0.462
+Zn+S,0.275334,1.081,4.5275,-6.217,-10.136,-1.194,2.254,1.547,1.099,-2.845,7.32,-7.106,-11.795,0.642,2.366,0.847,0.742
+Zn+Se,0.259251,1.081,4.5275,-6.217,-10.136,-1.194,2.254,1.547,1.099,-2.751,6.8485,-6.654,-10.946,1.316,2.177,0.952,0.798
+Zn+Te,0.2411885,1.081,4.5275,-6.217,-10.136,-1.194,2.254,1.547,1.099,-2.666,6.2665,-6.109,-9.867,0.099,1.827,1.141,0.945
+Ga+N,0.4332885,-0.108,2.963,-2.732,-5.818,0.13,2.163,1.33,0.994,-1.867,7.726,-7.239,-13.585,3.057,1.54,0.511,0.539
+Ga+P,0.340923,-0.108,2.963,-2.732,-5.818,0.13,2.163,1.33,0.994,-1.92,5.8355,-5.596,-9.751,0.183,1.771,0.966,0.826
+Ga+As,0.2714115,-0.108,2.963,-2.732,-5.818,0.13,2.163,1.33,0.994,-1.839,5.5505,-5.341,-9.262,0.064,2.023,1.043,0.847
+Ga+Sb,0.158226,-0.108,2.963,-2.732,-5.818,0.13,2.163,1.33,0.994,-1.847,5.1575,-4.991,-8.468,0.105,2.065,1.232,1.001
+Ge+Ge,0.202076,-0.949,4.258,-4.046,-7.567,2.175,2.373,1.162,0.917,-0.949,4.258,-4.046,-7.567,2.175,2.373,1.162,0.917
+Rb+F,-0.136152,-0.59,2.4395,-2.36,-4.289,-0.705,1.96,3.199,2.24,-4.273,11.8385,-11.294,-19.404,1.251,1.428,0.371,0.406
+Rb+Cl,-0.160783,-0.59,2.4395,-2.36,-4.289,-0.705,1.96,3.199,2.24,-3.971,8.9365,-8.7,-13.902,0.574,1.666,0.756,0.679
+Rb+Br,-0.164282,-0.59,2.4395,-2.36,-4.289,-0.705,1.96,3.199,2.24,-3.739,8.1945,-8.001,-12.65,0.708,1.869,0.882,0.749
+Rb+I,-0.168527,-0.59,2.4395,-2.36,-4.289,-0.705,1.96,3.199,2.24,-3.513,7.385,-7.236,-11.257,0.213,1.722,1.071,0.896
+Sr+O,-0.220613,0.343,2.8445,-3.641,-6.032,-1.379,1.204,2.548,1.911,-3.006,9.7195,-9.197,-16.433,2.541,2.219,0.427,0.462
+Sr+S,-0.368853,0.343,2.8445,-3.641,-6.032,-1.379,1.204,2.548,1.911,-2.845,7.32,-7.106,-11.795,0.642,2.366,0.847,0.742
+Sr+Se,-0.3751365,0.343,2.8445,-3.641,-6.032,-1.379,1.204,2.548,1.911,-2.751,6.8485,-6.654,-10.946,1.316,2.177,0.952,0.798
+Sr+Te,-0.380649,0.343,2.8445,-3.641,-6.032,-1.379,1.204,2.548,1.911,-2.666,6.2665,-6.109,-9.867,0.099,1.827,1.141,0.945
+Ag+F,-0.1557595,-1.667,4.8625,-4.71,-8.058,-0.479,2.968,1.883,1.316,-4.273,11.8385,-11.294,-19.404,1.251,1.428,0.371,0.406
+Ag+Cl,-0.043813,-1.667,4.8625,-4.71,-8.058,-0.479,2.968,1.883,1.316,-3.971,8.9365,-8.7,-13.902,0.574,1.666,0.756,0.679
+Ag+Br,-0.0303675,-1.667,4.8625,-4.71,-8.058,-0.479,2.968,1.883,1.316,-3.739,8.1945,-8.001,-12.65,0.708,1.869,0.882,0.749
+Ag+I,0.036523,-1.667,4.8625,-4.71,-8.058,-0.479,2.968,1.883,1.316,-3.513,7.385,-7.236,-11.257,0.213,1.722,1.071,0.896
+Cd+O,-0.0871575,0.839,4.371,-5.952,-9.581,-1.309,2.604,1.736,1.232,-3.006,9.7195,-9.197,-16.433,2.541,2.219,0.427,0.462
+Cd+S,0.0709195,0.839,4.371,-5.952,-9.581,-1.309,2.604,1.736,1.232,-2.845,7.32,-7.106,-11.795,0.642,2.366,0.847,0.742
+Cd+Se,0.0831555,0.839,4.371,-5.952,-9.581,-1.309,2.604,1.736,1.232,-2.751,6.8485,-6.654,-10.946,1.316,2.177,0.952,0.798
+Cd+Te,0.1127115,0.839,4.371,-5.952,-9.581,-1.309,2.604,1.736,1.232,-2.666,6.2665,-6.109,-9.867,0.099,1.827,1.141,0.945
+In+N,0.150714,-0.256,2.8965,-2.697,-5.537,0.368,3.108,1.498,1.134,-1.867,7.726,-7.239,-13.585,3.057,1.54,0.511,0.539
+In+P,0.170062,-0.256,2.8965,-2.697,-5.537,0.368,3.108,1.498,1.134,-1.92,5.8355,-5.596,-9.751,0.183,1.771,0.966,0.826
+In+As,0.122301,-0.256,2.8965,-2.697,-5.537,0.368,3.108,1.498,1.134,-1.839,5.5505,-5.341,-9.262,0.064,2.023,1.043,0.847
+In+Sb,0.0795525,-0.256,2.8965,-2.697,-5.537,0.368,3.108,1.498,1.134,-1.847,5.1575,-4.991,-8.468,0.105,2.065,1.232,1.001
+Sn+Sn,0.0163115,-1.039,4.041,-3.866,-7.043,0.008,2.03,1.344,1.057,-1.039,4.041,-3.866,-7.043,0.008,2.03,1.344,1.057
+B+Sb,0.5808505,-0.107,4.1485,-3.715,-8.19,2.248,1.946,0.826,0.805,-1.847,5.1575,-4.991,-8.468,0.105,2.065,1.232,1.001
+Cs+F,-0.11158,-0.57,2.288,-2.22,-4.006,-0.548,1.974,3.164,2.464,-4.273,11.8385,-11.294,-19.404,1.251,1.428,0.371,0.406
+Cs+Cl,-0.1524795,-0.57,2.288,-2.22,-4.006,-0.548,1.974,3.164,2.464,-3.971,8.9365,-8.7,-13.902,0.574,1.666,0.756,0.679
+Cs+Br,-0.1579615,-0.57,2.288,-2.22,-4.006,-0.548,1.974,3.164,2.464,-3.739,8.1945,-8.001,-12.65,0.708,1.869,0.882,0.749
+Cs+I,-0.164926,-0.57,2.288,-2.22,-4.006,-0.548,1.974,3.164,2.464,-3.513,7.385,-7.236,-11.257,0.213,1.722,1.071,0.896
+Ba+O,-0.095162,0.278,2.619,-3.346,-5.516,-2.129,1.351,2.632,2.149,-3.006,9.7195,-9.197,-16.433,2.541,2.219,0.427,0.462
+Ba+S,-0.3264895,0.278,2.619,-3.346,-5.516,-2.129,1.351,2.632,2.149,-2.845,7.32,-7.106,-11.795,0.642,2.366,0.847,0.742
+Ba+Se,-0.3501065,0.278,2.619,-3.346,-5.516,-2.129,1.351,2.632,2.149,-2.751,6.8485,-6.654,-10.946,1.316,2.177,0.952,0.798
+Ba+Te,-0.3810935,0.278,2.619,-3.346,-5.516,-2.129,1.351,2.632,2.149,-2.666,6.2665,-6.109,-9.867,0.099,1.827,1.141,0.945
+Ge+C,0.807408,-0.949,4.258,-4.046,-7.567,2.175,2.373,1.162,0.917,-0.872,5.862,-5.416,-10.852,1.992,1.631,0.63,0.644
+Sn+C,0.45098,-1.039,4.041,-3.866,-7.043,0.008,2.03,1.344,1.057,-0.872,5.862,-5.416,-10.852,1.992,1.631,0.63,0.644
+Ge+Si,0.263938,-0.949,4.258,-4.046,-7.567,2.175,2.373,1.162,0.917,-0.993,4.3755,-4.163,-7.758,0.44,1.89,1.134,0.938
+Sn+Si,0.136481,-1.039,4.041,-3.866,-7.043,0.008,2.03,1.344,1.057,-0.993,4.3755,-4.163,-7.758,0.44,1.89,1.134,0.938
+Sn+Ge,0.086873,-1.039,4.041,-3.866,-7.043,0.008,2.03,1.344,1.057,-0.949,4.258,-4.046,-7.567,2.175,2.373,1.162,0.917

requirements.txt

+numpy==1.16.2
+scipy==1.2.1
+pandas==0.24.2
+sklearn==0.0
+joblib==0.13.2
+diskcache==4.1.0

tcmi/__init__.py

+# -*- coding: utf-8 -*-
+"""
+@package    tcmi
+
+@copyright  Copyright (c) 2018+ Fritz Haber Institute of the Max Planck Society,
+            Benjamin Regler <regler@fhi-berlin.mpg.de>
+@license    See LICENSE file for details.
+
+Licensed under the Apache License, Version 2.0 (the 'License').
+You may not use this file except in compliance with the License.
+"""
+
+# Metadata
+__name__ = 'tcmi'
+__version__ = '1.0.0'
+__description__ = 'A Python package for estimating mutual dependencies of multivariate continuous distributions'
+__author__ = 'Benjamin Regler <regler@fhi-berlin.mpg.de>'
+__url__ = 'https://github.com/sommerregen/tcmi'
+__email__ = 'regler@fhi-berlin.mpg.de'
+__copyright__ = 'Copyright 2018+, Fritz Haber Institute of the Max Planck ' \
+                'Society'
+__license__ = 'Apache License, Version 2.0'

tcmi/entropy.py

+# -*- coding: utf-8 -*-
+"""
+@package    tcmi.entropy
+
+@copyright  Copyright (c) 2018+ Fritz Haber Institute of the Max Planck Society,
+            Benjamin Regler <regler@fhi-berlin.mpg.de>
+@license    See LICENSE file for details.
+
+Licensed under the Apache License, Version 2.0 (the "License").
+You may not use this file except in compliance with the License.
+"""
+
+import joblib
+import functools
+import itertools
+import collections
+import numpy as np
+
+from .utils import ndindex, chunk_iterator
+from scipy.stats import hypergeom as _hypergeom
+
+
+def cumulative_mutual_information(y, x, adjust=0, cache=False, n_jobs=None,
+                                  verbose=0, pre_dispatch='2*n_jobs',
+                                  return_scores=False):
+    """Compute the cumulative mutual information shared by X and Y.
+    """
+    eps = np.finfo(np.float_).eps
+    kwargs = dict(n_jobs=n_jobs, verbose=verbose, pre_dispatch=pre_dispatch,
+                  return_vector=True)
+
+    x = clean_data(x, allow_duplicates=True)
+
+    hce, hce_corr, ce, counts = conditional_cumulative_entropy(y, *x,
+                                                               return_counts=True,
+                                                               **kwargs)
+    hce_corr = hce_corr.sum(axis=1)
+    hce = hce.sum(axis=1)
+
+    hce0 = cumulative_baseline_correction(y, *x, correction=counts,
+                                          cache=cache, **kwargs)
+    hce0 = hce0.sum(axis=1)
+    ce0 = cumulative_entropy(y)
+
+    score = 1 - hce / ce.sum(axis=1)
+    score0 = 1 - hce0 / ce0
+
+    score_corr = np.zeros_like(score)
+    if np.any(hce_corr > eps):
+        score_corr = np.clip(1 - hce / np.maximum(hce_corr, eps), 0, 1)
+
+    total_score = np.clip(score - score0 - score_corr, 0, 1)
+
+    result = total_score.mean()
+    if return_scores:
+        result = (result, (total_score, score, score_corr, score0))
+    return result
+
+
+def cumulative_entropy(y, return_vector=False, return_inverse=False,
+                       direction=0):
+    """Computes the cumulative entropy of y.
+    """
+    index = _direction_to_index(direction)
+    
+    inverse = False
+    if return_inverse:
+        y, inverse, probability = np.unique(y, return_inverse=True,
+                                            return_counts=True)
+    else:
+        y, probability = np.unique(y, return_counts=True)
+
+    dy = np.diff(y)
+    size = y.size
+
+    entropies = []
+    directions = ((index, ) if isinstance(index, int) else (-1, 1))
+    for direction in directions:
+        ce = np.zeros(size, dtype=np.float_)
+
+        if direction == -1:
+            p = np.add.accumulate(probability, dtype=np.float_)
+            p /= p[-1]
+
+            start = (1 if return_vector else size - 1)
+
+            for i in range(start, size):
+                vector = p[:i]
+                ce[i] -= np.inner(dy[:i], vector * np.log2(vector))
+            entropies.append(ce if return_vector else ce[-1])
+
+        elif direction == 1:
+            counts = probability.sum()
+            p = np.add.accumulate(probability[::-1], dtype=np.float_)[::-1]
+            p /= counts
+
+            stop = (size if return_vector else 2)
+
+            for i in range(1, stop):
+                vector = p[i:]
+                ce[i - 1] -= np.inner(dy[i - 1:], vector * np.log2(vector))
+            entropies.append(ce if return_vector else ce[0])
+
+        else:
+            raise RuntimeError('Unknown direction "{:d}".'.format(direction))
+
+    result = (np.array(entropies) if len(entropies) > 1 else entropies[0])
+    if return_inverse:
+        result = (result, inverse)
+    return result
+
+
+def _compute_entropy(y, x, x_sorted, dimensions, masks, offset=0):
+    """Compute entropy.
+    """
+    cache = (None, True)
+    shape = (2, y.size)
+    counts = np.zeros(shape, dtype=np.int_)
+    entropy = np.zeros_like(counts, dtype=np.float)
+    entropy0 = np.zeros_like(counts, dtype=np.float)
+    entropy_corr = np.zeros_like(counts, dtype=np.float)
+
+    local_storage = {}
+    for i, dimension in enumerate(itertools.product(*dimensions)):
+        key, mask = cache
+
+        if key != dimension[:-1]:
+            mask = np.logical_and.reduce(
+                [masks[i][j] for i, j in enumerate(dimension[:-1])])
+            cache = (dimension[:-1], mask)
+
+        mask = np.logical_and(mask, masks[-1][dimension[-1]])
+
+        key = np.packbits(mask).tobytes()
+
+        bucket = local_storage.get(key, [])
+        if bucket:
+            for slot, index, value, value0, value_corr in bucket:
+                counts[slot, index] += 1
+                entropy[slot, index] -= value
+                entropy0[slot, index] -= value0
+                entropy_corr[slot, index] -= value_corr
+
+            continue
+
+        bucket = []
+        for slot, submask in enumerate(mask, offset):
+            size = submask.sum()
+            index = max(size - 1, 0)
+
+            counts[slot, index] += 1
+            if size < 2:
+                bucket.append((slot, index, np.float(0), np.float(0), np.float(0)))
+                continue
+
+            yc = np.compress(submask, y)
+
+            umask = np.ones(size + 1, dtype=np.bool_)
+            umask[1:-1] = (yc[1:] != yc[:-1])
+            yu = np.compress(umask[:-1], yc)
+
+            px = []
+            for variable in x[slot]:
+                value = np.compress(submask, variable)
+
+                if np.any(value != value[0]):
+                    px.append(
+                        np.searchsorted(np.sort(value), value, side='right'))
+
+            if not px:
+                px.append(np.full(size, size, dtype=np.uintp))
+            px = np.minimum.reduce(px)
+
+            if slot == 1:
+                py = np.searchsorted(yc, yc, side='left')
+                np.subtract(size, py, out=py)
+
+                pxy = np.minimum(px, py)
+                py = py / py[0]
+
+                dy = np.ediff1d(yu, to_begin=0)
+            else:
+                py = np.searchsorted(yc, yc, side='right')
+
+                pxy = np.minimum(px, py)
+                py = py / py[-1]
+
+                dy = np.ediff1d(yu, to_end=0)
+
+            if yu.size < yc.size:
+                delta_y = np.zeros(size, dtype=np.float_)
+                idx = np.nonzero(umask[:-1])[0]
+                delta_y[idx] = dy
+                dy = delta_y
+                del delta_y
+
+            value0 = np.einsum('i,i,i', dy, py, np.log2(py))
+            value = np.einsum('i,i,i', dy, pxy / size, np.log2(pxy / px))
+
+            px_s = np.sort(px)
+            pxy_s = np.sort(pxy)
+            if slot:
+                px_s = px_s[::-1]
+                pxy_s = pxy_s[::-1]
+            value_corr = np.einsum('i,i,i', dy, pxy_s / size, np.log2(pxy_s / px_s))
+
+            value_corr = np.minimum(value_corr, value)
+
+            entropy[slot, index] -= value
+            entropy0[slot, index] -= value0
+            entropy_corr[slot, index] -= value_corr
+
+            bucket.append((slot, index, value, value0, value_corr))
+
+        local_storage[key] = bucket
+
+    return entropy, entropy0, entropy_corr, counts
+
+
+def conditional_cumulative_entropy(y, *x, direction=0, cache=False,
+                                   early_stopping_rounds=10, threshold=1e-3,
+                                   n_jobs=None, verbose=0, pre_dispatch='2*n_jobs',
+                                   return_counts=False, return_vector=False):
+    """Computes the conditional cumulative entropy of y given x.
+    """
+    y = np.asarray(y)
+
+    order = y.argsort(kind='mergesort')
+
+    y = np.take(y, order)
+    x = tuple(np.take(value, order) for value in x)
+    x = clean_data(x, allow_duplicates=True, return_original=False,
+                   return_reverse=True)
+
+    x_sorted = []
+    for variables in x:
+        variables = tuple(np.sort(v) for v in variables)
+        x_sorted.append(variables)
+
+    masks, dimensions = _get_selection_masks(x[0], direction=direction,
+                                             return_dimension=True)
+
+    offset = _direction_to_index(direction)
+    if not isinstance(offset, int):
+        offset = 0
+
+    shape = (2, y.size)
+    counts = np.zeros(shape, dtype=np.int_)
+    entropy = np.zeros_like(counts, dtype=np.float)
+    entropy0 = np.zeros_like(counts, dtype=np.float)
+    entropy_corr = np.zeros_like(counts, dtype=np.float)
+
+    processor = joblib.Parallel(n_jobs=n_jobs, backend='multiprocessing',
+                                verbose=verbose, pre_dispatch=pre_dispatch)
+    callback = joblib.delayed(_compute_entropy)
+
+    chunks = processor._effective_n_jobs()
+    iterator = chunk_iterator(ndindex(*dimensions, raw=True, grouped=True,
+                                      multi_index=True), chunks)
+
+    buffer = collections.deque(maxlen=early_stopping_rounds)
+    buffer.extend((i, 0, -1) for i in range(2))
+
+    counter = 0
+    total = np.prod(dimensions)
+    eps = np.finfo(np.float).eps
+
+    for i, indices in enumerate(iterator):
+        levels, indices = tuple(zip(*indices))
+
+        results = processor(callback(y, x, x_sorted, index, masks, offset=offset)
+                            for index in indices)
+
+        for level, values in zip(levels, results):
+            entropy += values[0]
+            entropy0 += values[1]
+            entropy_corr += values[2]
+            counts += values[3]
+
+            total_counts = np.count_nonzero(counts, axis=1)
+            ce = entropy0.sum(axis=1) / (total_counts + eps)
+            hce = entropy.sum(axis=1) / (total_counts + eps)
+
+            xx, yy, scores = tuple(zip(*buffer))
+
+            progress = counts.sum() / total
+            score = np.mean(hce / (ce + eps))
+            buffer.append((progress, score, np.mean(yy)))
+
+            model = np.poly1d(np.polyfit(xx, yy, 1))
+            value = model(1)
+
+            model = np.poly1d(np.polyfit(xx, scores, 1))
+            value0 = model(1)
+
+            if 0 <= value <= 1 and 0 <= value0 <= 1 and \
+                    abs(value - value0) < threshold:
+                counter += 1
+
+                if counter > early_stopping_rounds > 0:
+                    break
+            else:
+                counter = 0
+
+        if counter > early_stopping_rounds > 0:
+            break
+
+    results = (entropy, entropy_corr, entropy0)
+    if not return_vector:
+        results = tuple(vector.sum(axis=1) for vector in results)
+    if return_counts:
+        results += (counts, )
+
+    index = _direction_to_index(direction)
+    return tuple(result[index] for result in results)
+
+
+def cumulative_baseline_correction(y, *x, direction=0, correction=None, cache=False,
+                                   n_jobs=None, verbose=0, pre_dispatch='2*n_jobs',
+                                   return_vector=False, return_baseline=False):
+    """Computes the correction of chance baseline correction of y given x.
+    """
+    index = _direction_to_index(direction)
+
+    y = np.asarray(y)
+    weights = []
+
+    order = y.argsort(kind='mergesort')
+    x = tuple(np.take(value, order) for value in x)
+    variables = clean_data(x, allow_duplicates=True, return_original=False,
+                           return_reverse=True)
+
+    for variable in variables:
+        weight = _get_permutation_weights(*variable)
+        weights.append(weight)
+
+    if correction is not None:
+        counts = correction / correction.T[-1][:, np.newaxis]
+
+        weight = np.arange(y.size + 1, dtype=np.int_)
+        np.power(weight, len(x), out=weight)
+        weight = np.diff(weight)[::-1]
+
+        mask = np.logical_and(counts == 0, weights)
+
+        weights = counts / weight
+        weights[mask] = 1
+
+    processor = joblib.Parallel(n_jobs=n_jobs, backend='multiprocessing',
+                                verbose=verbose, pre_dispatch=pre_dispatch)
+    callback = joblib.delayed(compute_cumulative_baseline)
+
+    baseline = processor(callback(y, i + 1, return_vector=False, direction=0)
+                         for i in range(y.size))
+    baseline = np.column_stack(baseline)
+
+    corrections = []
+    for i, weight in enumerate(weights):
+        correction = weight[1:] * baseline[i][1:]
+        corrections.append(correction / weight.sum())
+
+    if return_vector:
+        corrections = np.atleast_2d(corrections)
+    else:
+        corrections = [np.sum(correction) for correction in corrections]
+
+    results = corrections[index]
+    if return_baseline:
+        results = (results, baseline)
+    return results
+
+
+def compute_cumulative_baseline(y, component, direction=0, return_entropy=False,
+                                return_vector=False):
+    """Computes the correction of chance baseline correction of response.
+    """
+    index = _direction_to_index(direction)
+
+    y = np.sort(y, kind='mergesort')
+    if component < 2:
+        score = np.zeros(2, dtype=np.float_)
+
+        result = (score[index], )
+        if return_entropy:
+            entropy = np.zeros((2, component), dtype=np.float_)
+            result += (entropy[index], )
+        if return_vector:
+            entropy = np.zeros((2, component), dtype=np.float_)
+            result += (entropy, )
+
+        return result[0] if len(result) == 1 else result
+
+    n = y.size
+    b = component
+
+    numbers = np.arange(b + 1, dtype=np.int_)
+    ce = np.zeros((2, numbers.size), dtype=np.float_)
+
+    ce[0, 1:] = (numbers[1:] / b) * np.log2(numbers[1:] / b)
+    ce[1, 1:] = ce[0, 1:][::-1]
+
+    entropy = np.zeros((2, n), dtype=np.float_, order='f')
+    fraction = b / n
+
+    size = n - b + 1
+
+    ij = np.arange(n, b - 2, -1, dtype=np.int_)
+    ij[0] = 0
+
+    np.add.accumulate(ij, out=ij)
+
+    matrix = np.zeros((ij[-1], 2), dtype=np.float_)
+    y = np.column_stack((y, -y[::-1]))
+
+    partial_entropy = np.zeros((n, 2), dtype=np.float_)
+    vector = np.zeros((2, b), dtype=np.float_)
+    scores = np.zeros(2, dtype=np.float_)
+
+    for i in range(1, n):
+        left = max(1, i + b - n)
+        right = 1 + min(i, b)
+        nij = numbers[left:right]
+
+        probability = _hypergeom.pmf(b - nij, n - 1, n - i, b - 1)
+
+
+        cutoff = (-1 if b < right else None)
+        jx = max(0, i - (n - b + 1))
+
+        for j, mij in enumerate(b - nij[:cutoff]):
+            number = (n - i) - (mij - 1)
+
+            m = np.int_(mij - 1)
+            weights = 1 - m / np.arange(n - i, m, -1)
+            weights[0] = mij / (n - i)
+
+            np.multiply.accumulate(weights, out=weights)
+
+            delta = y[i:number + i] - y[i - 1]
+
+            value = delta * weights[:, None] * ce[0, left + j] * probability[j]
+            offset = min(number, n - i)
+            s = slice(i, i + offset)
+
+            partial_entropy[s] -= value[:offset]
+            scores -= value.sum(axis=0)
+
+            k = j + jx
+            cursor = i - k - 1
+            l, r = ij[cursor:cursor + 2]
+
+            if k > 0:
+                factor = matrix[l + i - cursor - 2]
+                partial_entropy[s] += np.outer(weights[:offset], factor)
+
+                if i - k < size:
+                    matrix[r:ij[cursor + 2]] += matrix[l + 1:r]
+
+            matrix[l:r] = partial_entropy[cursor + 1:]
+
+            value = partial_entropy[s].T
+            vector[:, k + 1] += value.sum(axis=-1)
+
+            entropy[:, s] += value
+            partial_entropy.fill(0)
+
+    entropy *= fraction
+    vector *= fraction
+    scores *= fraction
+
+    result = (scores[index], )
+    if return_entropy:
+        result += (entropy[index], )
+    if return_vector:
+        result += (vector[index], )
+    return result if len(result) > 1 else result[0]
+
+
+def clean_data(x, allow_duplicates=False, return_original=True,
+               return_reverse=False):
+    """Clean data.
+    """
+    order = []
+    values = []
+    reverse = []
+    duplicates = set()
+
+    length = len(x)
+    for i, vector in enumerate(x, 1):
+        unique_vector = np.unique(vector)
+        size = unique_vector.size
+
+        if size > 1 or length == i:
+            vs = np.sort(vector)
+            unique_vector = np.searchsorted(vs, vector, side='right')
+
+            fingerprint = unique_vector.tobytes()
+            reversed_fingerprint = unique_vector[::-1].tobytes()
+
+            if fingerprint not in duplicates:
+                if not allow_duplicates:
+                    duplicates.add(fingerprint)
+                    duplicates.add(reversed_fingerprint)
+
+                values.append(vector if return_original else unique_vector)
+                order.append(size)
+
+                if return_reverse:
+                    unique_vector = np.searchsorted(vs, vector, side='left')
+                    reverse.append(vector[::-1] if return_original
+                                   else vs.size - unique_vector)
+
+    order = np.argsort(order, kind='mergesort')[::-1]
+
+    result = (tuple(values[i] for i in order), )
+    if return_reverse:
+        result += (tuple(reverse[i] for i in order), )
+    return result if len(result) > 1 else result[0]
+
+
+def _direction_to_index(direction):
+    """Map direction identifier to index.
+    """
+    directions = {-1: 0, 0: slice(None), 1: 1, '<=': 0, '<=>': slice(None), '>=': 1}
+
+    if direction not in directions:
+        raise RuntimeError('Unknown direction "{:d}".'.format(direction))
+
+    return directions[direction]
+
+
+def _get_permutation_weights(*x):
+    """Get permutation weights.
+    """
+    x = [np.sort(v, kind='mergesort') for v in x]
+    size = x[0].size
+
+    #
+    #
+    dtype = 'S' + str(len(x) * x[0].itemsize)
+
+    data = np.column_stack(x)
+    data = data.view(dtype=dtype).flatten()
+
+    block = np.ones(size + 1, dtype=np.bool_)
+    block[1:-1] = data[1:] != data[:-1]
+    block = np.nonzero(block)[0]
+
+    weights = np.zeros(size, dtype=np.bool_)
+    for i in range(1, block.size):
+        left, right = block[i - 1:i + 1]
+        width = right - left
+
+        if width > 1 and min(width, size - right) > 0:
+            for j in range(left + 1, right):
+                limits = []
+                for variable in x:
+                    limit = np.searchsorted(variable, variable[j], side='right')
+                    limits.append(size - limit)
+
+                if len(limits) > 1 and min(limits) > 0:
+                    weights[j - 1] = 1
+
+        weights[right - 1] = 1
+
+    return weights
+
+
+def _get_selection_masks(x, direction=0, return_dimension=False):
+    """Get selection masks of x.
+    """
+    index = _direction_to_index(direction)
+    operators = [np.less_equal, np.greater_equal]
+    if isinstance(index, int):
+        operators = (operators[index], )
+
+    masks = []
+    dimensions = []
+
+    for variable in x:
+        values = np.unique(variable)
+
+        mask = []
+        for operator in operators:
+            mask.append(tuple(operator(variable, value) for value in values))
+
+        mask = np.stack(mask, axis=1)
+
+        dimensions.append(values.size)
+        masks.append(mask)
+
+    return (masks, dimensions) if return_dimension else masks

tcmi/estimators.py

+# -*- coding: utf-8 -*-
+"""
+@package    tcmi.estimators
+
+@copyright  Copyright (c) 2018+ Fritz Haber Institute of the Max Planck Society,
+            Benjamin Regler <regler@fhi-berlin.mpg.de>
+@license    See LICENSE file for details.
+
+Licensed under the Apache License, Version 2.0 (the "License").
+You may not use this file except in compliance with the License.
+"""
+
+import os
+import re
+import tempfile
+import subprocess
+
+import numpy as np
+import sklearn as sk
+
+from . import entropy
+
+
+class DependenceEstimator(sk.base.BaseEstimator, sk.base.RegressorMixin):
+    """A general-purpose dependence estimator for measuring mutual information.
+    """
+    
+    def __init__(self, method='tcmi', n_iter=50000, cache=False, n_jobs=None,
+                 verbose=0, fit_params=None, pre_dispatch='2*n_jobs'):
+        """Initialize dependence estimator.
+        """
+        self.method = method
+        self.n_iter = n_iter
+        self.cache = cache
+        
+        self.n_jobs = n_jobs
+        self.verbose = verbose
+        self.fit_params = fit_params
+        self.pre_dispatch = pre_dispatch
+
+    def fit(self, x, y):
+        """Fit dependence estimator.
+        """
+        x, y = sk.utils.check_X_y(x, y, multi_output=False, y_numeric=True)
+        size = y.shape[-1]
+
+        self._variables = (x, y)
+
+    def predict(self, x):
+        """Predict dependence using mutual information.
+        """
+        if not self._variables:
+            raise RuntimeError('Estimator not fitted yet.')
+        return self.score(*self._variables)
+
+    def score(self, x, y=None):
+        """Score mutual dependence.
+        """
+        method = self.method.lower()
+
+        score = 0
+        if method == 'tcmi':
+            if not isinstance(x, (tuple, list, np.ndarray)):
+                x = x.to_numpy()
+            if isinstance(x, np.ndarray):
+                x = x.T
+                
+            x = tuple(np.asarray(variable) for variable in x)
+            kwargs = dict(n_jobs=self.n_jobs, verbose=self.verbose, adjust=0.1,
+                          pre_dispatch=self.pre_dispatch, cache=self.cache)
+            score = entropy.cumulative_mutual_information(y, x, **kwargs)
+
+        return np.clip(score, 0, 1)

tcmi/subspace_search.py

+# -*- coding: utf-8 -*-
+"""
+@package    python.subspace_search
+
+@copyright  Copyright (c) 2018+ Fritz Haber Institute of the Max Planck Society,
+            Benjamin Regler <regler@fhi-berlin.mpg.de>
+@license    See LICENSE file for details.
+
+Licensed under the Apache License, Version 2.0 (the "License").
+You may not use this file except in compliance with the License.
+"""
+
+import time
+import heapq
+import bisect
+import joblib
+import pickle
+import hashlib
+import warnings
+import itertools
+import collections
+
+import numpy as np
+import scipy.stats
+
+from . import utils
+
+from sklearn import metrics
+from sklearn.metrics.scorer import _check_multimetric_scoring
+from sklearn.model_selection._validation import _fit_and_score
+from sklearn.model_selection import cross_validate as _cross_validate
+
+
+def cross_validate(estimator, x, y, scoring=None, cv='warn', **kwargs):
+    """Cross-validate model.
+    """
+    if cv is None or cv == 1:
+        scorers, _ = _check_multimetric_scoring(estimator, scoring=scoring)
+        indices = np.arange(len(y))
+
+        scores, = _fit_and_score(estimator, x, y, scorers, indices, indices,
+                                kwargs.get('verbose', False),
+                                None, kwargs.get('fit_params', None),
+                                return_train_score=False, return_times=False,
+                                return_estimator=False, error_score=False)
+
+        stats = {}
+        for name in scorers:
+            stats['test_%s' % name] = np.array([scores[name], ])
+    else:
+        stats = _cross_validate(estimator, x, y, cv=cv, scoring=scoring,
+                                **kwargs)
+
+    return stats
+
+
+def evaluate(model, x, y, cv='warn', scoring='', fit_params=None, n_jobs=None,
+             key=''):
+    """Evaluate model and compute statistics.
+    """
+    if not scoring:
+        if np.issubdtype(y.dtype, np.integer):
+            scoring =_SCORINGS['classification']
+        else:
+            scoring = _SCORINGS['regression']
+    elif not isinstance(scoring, (list, tuple)):
+        scoring = (scoring, )
+
+    scorers = {}
+    for scorer in scoring:
+        name = _get_name(scorer)
+        scorers[name] = metrics.make_scorer(_mutual_information)
+
+    with warnings.catch_warnings():
+        warnings.simplefilter(action='ignore', category=FutureWarning)
+        stats = cross_validate(model, x, y, cv=cv, scoring=scorers,
+                               fit_params=fit_params, n_jobs=n_jobs)
+
+    for key in list(stats):
+        if not key.startswith('test_'):
+            del stats[key]
+            continue
+
+        raw_value = stats.pop(key)
+
+        key = key[5:]
+        scorer = scorers[key]
+        raw_value *= scorer._sign
+
+        stat_value = raw_value
+        if stat_value.size == 1:
+            stat_value = np.concatenate((stat_value, stat_value))
+
+        stats.update({
+            key: raw_value,
+            '{:s}_mean'.format(key): stat_value.mean(),
+            '{:s}_std'.format(key): stat_value.std(),
+        })
+
+    return stats
+
+
+def check_subspace(x, subspace, threshold=0.95):
+    """Check subspace before computation.
+    """
+    filtered_subspace = [feature.strip('+-|') for feature in subspace]
+    alias = filtered_subspace[-1]
+    key = subspace[-1]
+
+    if alias in filtered_subspace[:-1]:
+        return False
+
+    rejection = 1 - threshold
+    data = utils.prepare_data(x[[key]], [], copy=True)
+    for column, reference in data.items():
+        for feature in subspace[:-1]:
+            if np.corrcoef(x[feature], reference)[0, 1]**2 > threshold:
+                return False
+
+            value = scipy.stats.ks_2samp(x[feature], reference).pvalue
+            if value > rejection:
+                False
+
+    return True
+
+
+def _compute_subspace(task, settings):
+    """Simple wrapper to compute a subspace task.
+    """
+    kwargs = settings.copy()
+    i, key, features, params = task
+    callback, model, x, y = [kwargs.pop(k) for k in ('callback', 'model', 'x', 'y')]
+
+    stats = callback(model, x[features], y, key=key, **kwargs)
+    return i, features, stats
+
+
+def _mutual_information(y_true, y_pred):
+    """Mutual information score.
+    """
+
+    return y_pred
+
+
+def _get_name(scoring):
+    """Get name of score.
+    """
+    return 'mutual_information'
+
+
+
+def get_subspaces(data, target, model, cv='warn', alpha=0.95, beta=0, tol=5e-3, 
+                  subspace=(), depth=-1, scoring='mutual_information_score', fit_params=None,
+                  return_score=False, verbose=0, n_jobs=1, pre_dispatch='2*n_jobs',
+                  cache=False):
+    """Get subspaces.
+    """
+    x, y = data.drop(target, axis=1), data[target]
+    scorer = metrics.make_scorer(_mutual_information)
+    sign = scorer._sign
+
+    if depth < 0:
+        depth = x.shape[-1]
+
+    if not isinstance(fit_params, dict):
+        fit_params = {}
+
+    keys = [key for key in sorted(x) if key not in subspace]
+    queue = collections.deque([{
+        'stats': None,
+        'subspace': subspace,
+        'keys': keys,
+        'path': (np.float(0), )
+    }])
+
+    nodes = []
+    bounds = {}
+    safe_features = dict()
+    subspaces = dict()
+    optimal = (np.finfo(np.float).min, np.float_(0), len(subspace), queue[0])
+
+    counter = 0
+    eps = np.finfo(np.float_).eps
+
+    evaluate_model = evaluate
+
+    if subspace:
+        features = list(subspace)
+        key = _get_key(subspace)
+        size = len(subspace)
+        
+        params = fit_params.copy()
+        if 'eval_set' in params:
+            params['eval_set'] = [(x_test[features], y_test)
+                                  for x_test, y_test in params['eval_set']]
+
+        stats = evaluate_model(model, x[features], y, cv=cv, key=key,
+                               fit_params=params, scoring=scoring,
+                               n_jobs=n_jobs)
+        score_mean, score_std = _get_score(scoring, scorer, stats, alpha=alpha)
+
+        queue[0].update({
+            'stats': stats,
+            'path': queue[0]['path'] + (sign * score_mean, )
+        })
+        optimal = (score_mean, score_std, len(subspace), queue[0])
+        bounds[size] = score_mean - score_std - eps
+
+    n_cores = 1
+    if cv is None or cv == 1:
+        n_cores = n_jobs
+
+    processor = joblib.Parallel(n_jobs=n_cores, pre_dispatch=pre_dispatch)
+    callback = joblib.delayed(_compute_subspace)
+
+    kwargs = dict(callback=evaluate_model, model=model, x=x, y=y, cv=cv,
+                  scoring=scoring, n_jobs=n_jobs)
+
+    while queue:
+        node = queue.popleft()
+        node_mean, node_std = _get_score(scoring, scorer, node['stats'],
+                                         alpha=alpha)
+        
+        subspace = list(node['subspace'])
+        key = _get_key(subspace)
+        size = len(subspace)
+
+        delta = (node_mean - node_std) - (optimal[0] - optimal[1])
+        if key in subspaces or (delta < 0 and size > optimal[2]):
+            continue
+
+        if verbose:
+            print('Subspace: {{{:s}}}  ->  [{:.2f} +- {:.2f}]'
+                  .format(','.join(subspace), node_mean, node_std))
+
+        subspaces[key] = safe_features.pop(key, 0)
+        counter += 1
+        
+        keys = node['keys']
+        scores = np.zeros((len(keys), 2), dtype=np.float_)
+        statistics = [None for i in range(scores.shape[0])]
+
+        tasks = []
+        for i, key in enumerate(keys):
+            features = subspace + [key, ]
+
+            key = _get_key(features)
+            if key in subspaces or key in safe_features or not check_subspace(x, features):
+                scores[i] = (-np.inf, np.float_(0))
+                continue
+
+            safe_features[key] = -1
+            params = fit_params.copy()
+            if 'eval_set' in params:
+                params['eval_set'] = [(x_test[features], y_test)
+                                      for x_test, y_test in params['eval_set']]
+
+            task = (i, key, features, params)
+            tasks.append(task)
+            safe_features[key] = -1
+
+        for i, features, stats in processor(callback(task, kwargs)
+                                            for task in tasks):
+            if verbose:
+                print('{:<85s}'.format(','.join(features)), end=' ')
+
+            scores[i] = _get_score(scoring, scorer, stats, alpha=alpha)
+            statistics[i] = stats
+
+            if verbose:
+                print('[Score: {:5.2f} +- {:5.2f}]'
+                      .format(sign * scores[i][0], scores[i][1]))
+    
+        indices = sorted(_get_node(score[1] - score[0], keys[i], i)
+                         for i, score in enumerate(scores))
+        keys, indices = tuple(zip(*[(keys[i], i) for *_, i in indices]))
+
+        next_bound = bounds.get(size, bounds.get(size - 1, -np.inf))
+        lower_bound = max(bounds.get(size - 1, -np.inf),
+                          node_mean - node_std + tol)
+        
+        if verbose:
+            print('Bound: {:.2f} -> {:.2f}'
+                  .format(sign * lower_bound, sign * next_bound))
+
+        upper_bound = np.subtract(*scores[indices[0]])
+        if np.isfinite(lower_bound) and np.isfinite(upper_bound):
+            bound = lower_bound + beta * (upper_bound - lower_bound)
+        else:
+            bound = lower_bound
+        
+        for i, j in enumerate(indices):
+            score_mean, score_std = scores[j]
+            
+            if score_mean < bound:
+                break
+
+            delta = (score_mean - score_std) - (optimal[0] - optimal[1])
+            difference = (score_mean - score_std) - (node_mean - node_std)
+            if difference < 0 or np.abs(difference) < tol:
+                continue
+
+            child_node = {
+                'stats': statistics[j],
+                'keys': keys[:i] + keys[i + 1:],
+                'subspace': tuple(subspace + [keys[i], ]),
+                'path': node['path'] + (sign * score_mean, )
+            }
+
+            key = _get_key(child_node['subspace'])
+            if key in subspaces:
+                print('Something went wrong.\n', key, keys, scores)
+                continue
+            
+            nodes.append(_get_node(score_std - score_mean, counter, child_node))
+            
+            length = len(child_node['subspace'])
+            if delta > tol or (np.abs(delta) < tol and length < optimal[2]):
+                optimal = (score_mean, score_std, length, child_node.copy())
+            
+            index = len(nodes) - 1
+            if child_node['keys'] and length < depth and difference > tol and \
+                   not np.isclose(score_mean - score_std, 1):
+                queue.append(child_node)
+                safe_features[key] = index
+            else:
+                subspaces[key] = index
+            counter += 1
+        
+        score_mean, score_std = scores[indices[0]]
+        bounds[size] = max(next_bound, score_mean - score_std - eps)
+
+        if verbose:
+            print('[{:s} - {:d} subspaces remaining]\n'
+                  .format(time.strftime("%Y-%m-%d %H:%M:%S"), len(queue)))
+
+    if verbose:
+        print('\nBalance search tree')
+
+    paths = {}
+    size = len(nodes)
+    safe_features = set()
+
+    cursor = 0
+    while cursor < len(nodes):
+        *_, index, node = nodes[cursor]
+        subspace = node['subspace']
+        cursor += 1
+
+        key = _get_key(subspace)
+        if key in safe_features:
+            continue
+
+        weight = None
+        depth = len(subspace)
+        if depth > 1:
+            node_mean, node_std = _get_score(scoring, scorer, node['stats'],
+                                             alpha=alpha)
+            
+            for keys in itertools.combinations(subspace, depth - 1):
+                features = sorted(keys)
+                key = _get_key(features)
+
+                if key not in subspaces:
+                    params = fit_params.copy()                    
+                    if 'eval_set' in params:
+                        params['eval_set'] = [(x_test[features], y_test)
+                                              for x_test, y_test in params['eval_set']]
+                    
+                    if verbose:
+                        print('{:<90s}'.format(','.join(features)), end=' ')
+
+                    stats = evaluate_model(model, x[features], y, cv=cv, key=key,
+                                           fit_params=params, scoring=scoring,
+                                           n_jobs=n_jobs)
+
+                    score_mean, score_std = _get_score(scoring, scorer, stats,
+                                                       alpha=alpha)
+                    if verbose:
+                        print('[Score: {:5.2f} +- {:5.2f}]'
+                              .format(sign * score_mean, score_std))
+                    
+                    path = paths.get(key, (np.float(0), ))
+                    
+                    child_node = {
+                        'stats': stats,
+                        'subspace': keys,
+                        'path': path + (sign * score_mean, )
+                    }
+
+                    delta = (score_mean - score_std) - (node_mean - node_std)
+                    if np.abs(delta) < tol:
+                        nodes.append(_get_node(score_std - score_mean, counter,
+                                               child_node))
+                        paths[key] = child_node['path']
+                        subspaces[key] = len(nodes) - 1
+                        weight = np.inf
+                        counter += 1
+                else:
+                    subset = nodes[subspaces.get(key)]
+                    stats = subset[-1]['stats']
+
+                    score_mean, score_std = _get_score(scoring, scorer, stats,
+                                                       alpha=alpha)
+
+                    delta = (score_mean - score_std) - (node_mean - node_std)
+                    if np.abs(delta) < tol:
+                        weight = np.inf
+
+        if weight is not None:
+            nodes[cursor - 1] = (weight, index, node)
+        elif depth > 1:
+            for keys in itertools.combinations(subspace, depth - 1):
+                key = _get_key(keys)
+                safe_features.add(key)
+
+    heapq.heapify(nodes)
+    safe_features.clear()
+    subspaces = []
+    subsets = []
+
+    cursor = np.inf
+    while nodes:
+        *score, node = heapq.heappop(nodes)
+        node.pop('path', None)
+        node.pop('keys', None)
+
+        subspace = sorted(node['subspace'])
+        key = _get_key(subspace)
+        depth = len(subspace)
+
+        if key in safe_features or score[0] > 0 or depth > cursor:
+            continue
+
+        bucket = None
+        for size, subset in subsets:
+            if depth <= size:
+                break
+
+            if subset.issubset(subspace):
+                bucket = subset
+                break
+
+        if bucket is None:
+            bisect.insort_left(subsets, (depth, set(subspace)))            
+            subspaces.append(node)
+            safe_features.add(key)
+            cursor = depth
+            
+    del subsets
+    del nodes
+    
+    result = (subspaces, )    
+    if return_score:
+        del optimal[-1]['keys']
+        result += (optimal[-1], )
+    return result if len(result) > 1 else result[0]
+
+
+def _get_score(key, scorer, stats=None, alpha=0.95):
+    """Get score.
+    """
+    if stats is None:
+        stats = {}
+
+    key = _get_name(key)
+    mean = scorer._sign * stats.get(key + '_mean', -scorer._sign * np.inf)
+
+    values = stats.get(key, None)
+    std = np.float_(0)
+    
+    return mean, std
+
+
+def _get_node(value, counter, node, start=1, stop=4):
+    """Get node.
+    """
+    return tuple(np.around(value, decimals=num)
+                 for num in range(start, stop + 1)) + (counter, node)
+
+
+def _get_key(subspace):
+    """Get key.
+    """
+    return ','.join(sorted(subspace, key=lambda k: (k.strip('-+|'), len(k), k)))
+
+
+_SCORINGS = {
+    'classification': (
+        'accuracy',
+        'f1_micro',
+        'precision_micro',
+        'recall_micro'
+    ),
+    'regression': (
+        'r2',
+        'neg_root_mean_squared_error',
+        'neg_mean_absolute_error',
+        'neg_maximum_absolute_error',
+        'neg_relative_root_mean_squared_error',
+        'neg_relative_mean_absolute_error',
+        'neg_relative_maximum_absolute_error'
+    )
+}

tcmi/utils.py

+# -*- coding: utf-8 -*-
+"""
+@package    tcmi.utils
+
+@copyright  Copyright (c) 2018+ Fritz Haber Institute of the Max Planck Society,
+            Benjamin Regler <regler@fhi-berlin.mpg.de>
+@license    See LICENSE file for details.
+
+Licensed under the Apache License, Version 2.0 (the "License").
+You may not use this file except in compliance with the License.
+"""
+
+import itertools
+import numpy as np
+import pandas as pd
+
+from joblib import hashing
+
+# Mappings
+_MAPPINGS = {
+    '|{}|': [np.abs],
+    '-{}': [np.negative],
+    '-|{}|': [np.abs, np.negative]
+}
+
+
+def get_fingerprint(x):
+    """Computes the fingerprint of a Numpy vector.
+    """
+    fingerprint = np.searchsorted(np.sort(x), x, side='left')
+    return compute_hash(fingerprint)
+
+
+def compute_hash(x):
+    """Hash object.
+    """
+    return hashing.hash(x, hash_name='md5')
+
+
+def is_numeric(obj):
+    """Check if object is numeric.
+    """
+    
+    flag = bool(isinstance(obj, np.ndarray) and obj.dtype.kind in 'OSU')
+    if isinstance(obj, pd.Series):
+        flag |= pd.api.types.is_categorical_dtype(obj)
+    elif isinstance(obj, pd.DataFrame):
+        for key in obj.columns:
+            flag |= is_numeric(obj[key])
+    
+    attrs = ['__add__', '__sub__', '__mul__', '__truediv__', '__pow__']
+    return all(hasattr(obj, attr) for attr in attrs) and not flag
+
+
+def prepare_data(data, target, copy=False):
+    """Prepare data by agumenting feature space.
+    """   
+    if copy:
+        data = data.copy()
+
+    if isinstance(target, str):
+        target = [target]
+    
+    keys = sorted(data)
+    hashes = set()
+    
+    for key in keys:
+        value = data[key]
+        if key in target:
+            continue
+
+        fingerprint = get_fingerprint(value)
+        hashes.add(fingerprint)
+
+        for label, chain in _MAPPINGS.items():
+            alias = label.replace('{}', '')
+            label = label.format(key)
+
+            result = value
+            if is_numeric(value):
+                for func in chain:
+                    result = func(result)
+                
+            fingerprint = get_fingerprint(result)
+            if fingerprint not in hashes:
+                data[label] = result
+                hashes.add(fingerprint)
+    
+    return data
+
+
+def filter_subsets(subsets, remove_duplicates=False):
+    """Filter subsets.
+    """
+    mappings = sorted(_MAPPINGS, key=lambda x: (-len(x), x))
+    
+    results = []
+    duplicates = set()
+    for subset in subsets:        
+        subspace_original = subset['subspace']
+        size = len(subspace_original)
+
+        normalized_subspace = subspace_original
+        for mapping in mappings:
+            prefix, suffix = mapping.split('{}', 1)
+            normalized_subspace = tuple(strip(k, prefix, suffix)
+                                        for k in normalized_subspace)
+
+        if remove_duplicates:
+            subspace = []
+            for k in normalized_subspace:
+                if k not in subspace:
+                    subspace.append(k)
+            normalized_subspace = tuple(subspace)
+
+        key = ','.join(sorted(normalized_subspace))
+        if key in duplicates:            continue
+
+        duplicates.add(key)
+        subset = subset.copy()
+        subset.update({
+            'subspace': normalized_subspace,
+            'subspace_original': subspace_original
+        })
+        results.append(subset)
+
+    return results
+
+
+def strip(text, prefix='', suffix=''):
+    """Remove substring from the left and right side of the text.
+    """
+    return strip_right(strip_left(text, prefix), suffix)
+
+
+def strip_left(text, prefix):
+    """Remove substring from the left side of the text.
+    """
+    if prefix and text.startswith(prefix):
+        text = text[len(prefix):]
+    return text
+
+
+def strip_right(text, suffix):
+    """Remove substring from the right side of the text.
+    """
+    if suffix and text.endswith(suffix):
+        text = text[:-len(suffix)]
+    return text
+
+
+def wrap_iterator(iterator, wrap=False, index=None):
+    """Return wrapped iterator for `yield from` syntax.
+    """
+    if wrap:
+        iterator = iter([iterator]) 
+    return (iterator if index is None
+            else itertools.zip_longest((), iterator, fillvalue=index))
+
+
+def chunk_iterator(iterable, n):
+    """Group an iterator in chunks of n without padding.
+    """
+    iterator = iter(iterable)
+    return iter(lambda: list(itertools.islice(iterable, n)), [])
+
+
+def index_split(index, dimension=1, method='symmetric'):
+    """Splits an index into parts.
+    """
+    splits = []
+    size = len(index)
+
+    if method == 'adaptive':
+        stack = [index]
+        split = []
+
+        while stack:
+            indices = stack
+            stack = []
+
+            for index in indices:
+                size = len(index)
+                if size < 3:
+                    split.extend(index)
+                    continue
+
+                divider = size // 2
+                split.append(index[divider])
+
+                stack.append(index[:divider])
+                stack.append(index[divider + 1:])
+
+            splits.append(np.array(split))
+            split = []
+
+    elif method == 'symmetric':
+        for i in range(np.math.ceil(size / 2)):
+            a, b = index[i], index[-1-i]
+            split = np.array((a, b))
+
+            splits.append(split if a < b else split[0:1])
+
+    elif method == 'interleave':
+        divider = max(2, np.sqrt(size).astype(np.int) // dimension)
+        step = size // divider + 1
+
+        splits = []
+        for i in range(step):
+            splits.append(index[i::step])
+
+    else:
+        raise KeyError('Unknown split method "{:s}".'.format(method))
+
+    return splits
+
+
+def ndindex(*indices, method='symmetric', raw=False, grouped=False,
+            multi_index=False):
+    """An N-dimensional iterator object to index arrays.
+    """
+    dimension = len(indices)
+    splits = [index_split(np.arange(index), dimension=dimension,
+                          method=method) for index in indices]
+
+    pool = [split.pop(0) for split in splits]
+    empty = np.array([], dtype=np.int_)
+
+    iteration = None
+    if multi_index:
+        iteration = 0
+
+    iterator = (pool if raw else itertools.product(*pool))
+    yield from wrap_iterator(iterator, wrap=grouped, index=iteration)
+
+    while True:
+        loop = False
+        if multi_index:
+            iteration += 1
+
+        staging = []
+        for split in splits:
+            flag = len(split) > 0
+            loop |= flag
+
+            staging.append(split.pop(0) if flag else empty)
+
+        if not loop:
+            break
+
+        for i in range(dimension):
+            if staging[i].size == 0:
+                continue
+
+            indices = pool.copy()
+            for j in range(i + 1, dimension):
+                indices[j] = np.concatenate((pool[j], staging[j]))
+            indices[i] = staging[i]
+
+            iterator = (indices if raw else itertools.product(*indices))
+            yield from wrap_iterator(iterator, wrap=grouped, index=iteration)
+
+        pool = [np.concatenate(v) for v in zip(pool, staging)]