Added kde tension code, general code restructuring

.gitignore

@@ -1,3 +1,4 @@
+# standard:
 *.DS_Store
 __pycache__/
 *.py[cod]

.material.py

@@ -2,6 +2,33 @@
 This is random material, do not read it :)
 """
 
+
+def _vec_to_log_pdm(vec, d):
+    """
+    """
+    # get indexes:
+    ind = np.tril_indices(d, 0)
+    # initialize:
+    mat = np.zeros((d, d))
+    mat[ind] = vec
+    # take exponential of the diagonal to ensure positivity:
+    mat[np.diag_indices(d)] = np.exp(np.diagonal(mat))
+    #
+    return np.dot(mat, mat.T)
+
+
+def _log_pdm_to_vec(pdm, d):
+    """
+    """
+    # decompose:
+    mat = np.linalg.cholesky(pdm)
+    # take log of diagonal:
+    mat[np.diag_indices(d)] = np.log(np.diagonal(mat))
+    #
+    return mat[np.tril_indices(d, 0)]
+
+
+
 def _temp_vec_kde_pdf(x, samples, weights):
     """
     Utility function to compute the KDE

.travis.yml

@@ -6,16 +6,16 @@ git:
 
 jobs:
   include:
-    - env:
+    - name: "Anaconda Python 3.7"
+      env:
         - CHANNEL="defaults"
         - PYDIST="ANACONDA"
-      python: "3.6"
-    - env:
+      python: "3.7"
+    - name: "Anaconda Python 3.8"
+      env:
         - CHANNEL="conda-forge"
         - PYDIST="ANACONDA"
       python: "3.8"
-    - name: "3.7 with bionic"
-      python: "3.7"
 
 install:
   # Setup anaconda following http://conda.pydata.org/docs/travis.html
@@ -29,6 +29,7 @@ install:
       conda create -q -n test-environment -c $CHANNEL python=$TRAVIS_PYTHON_VERSION scipy matplotlib;
       source activate test-environment;
     fi
+  - pip install --upgrade pip
   - python --version
   - pip install .
   - pip install coveralls

Makefile

@@ -1,18 +1,16 @@
-#
-# Simple make file to store commands that are used often
-#
+# Simple make file to store commands that are used often:
+
+install:
+	@python setup.py install --user
 
 test:
 	@python -m unittest discover tensiometer/tests
 
-install:
-	@python setup.py install
-
 coverage_report:
 	@coverage report
 
 test_with_coverage:
-	@coverage run -m unittest discover tensiometer/tests
+	@coverage run -m unittest discover tensiometer/tests -vvv -f
 	@coverage report
 
 run_examples:
@@ -21,7 +19,13 @@ run_examples:
 		jupyter nbconvert --execute --to html $$i --ExecutePreprocessor.timeout=-1; \
 	done;
 
-doc:
+prepare_examples:
+	@cd docs/example_notebooks && \
+	for i in *.ipynb ; do \
+		jupyter nbconvert --to html $$i; \
+	done;
+
+documentation:
 	@sphinx-build -b html docs/source build
 
 release:

README.rst

@@ -24,8 +24,8 @@ Tensiometer
 Description
 ============
 
-The tensiometer package is a collection of tools that extends `GetDist <https://pypi.org/project/GetDist/>`_ capabilities
-to test the level of agreement/disagreement between different posterior distributions.
+The tensiometer package is a collection of tools to test the level of
+agreement/disagreement between different posterior distributions.
 
 The best way to get up to speed is to read through the worked example
 `full worked example <https://tensiometer.readthedocs.org/en/latest/tension_example.html>`_
@@ -42,17 +42,18 @@ Alternatively one can download the source code from github::
 
   git clone https://github.com/mraveri/tensiometer.git
 
-and install with::
+and install it locally with the shortcut::
 
-  python setup.py install
+  make install
 
 You can test that the code is working properly by using::
 
   make test
 
+
 Dependencies
 =============
-* Python 3.6+
-* matplotlib 2.2+ (3.1+ recommended)
-* scipy
-* GetDist
+
+Uses mostly standard python packages and GetDist.
+
+For the full list of requirements see the `requirements.txt` file.

docs/source/conf.py

@@ -18,6 +18,7 @@
 autoclass_content = 'both'
 
 nitpicky = True
+nitpick_ignore = [('py:class', 'optional')]
 
 
 # Prevent spurious errors for every field ivar (not sure why..)
@@ -57,7 +58,9 @@ def setup(app):
                        'numpy': ('https://docs.scipy.org/doc/numpy/', None),
                        'scipy': ('https://docs.scipy.org/doc/scipy/reference/', None),
                        'matplotlib': ('https://matplotlib.org/', None),
-                       'getdist': ('https://getdist.readthedocs.io/en/latest/', None)}
+                       'getdist': ('https://getdist.readthedocs.io/en/latest/', None),
+                       "tensorflow": ("https://www.tensorflow.org/api_docs/python", "https://github.com/GPflow/tensorflow-intersphinx/raw/master/tf2_py_objects.inv"),
+                       "tensorflow_probability": ("https://www.tensorflow.org/probability/api_docs/python", "https://github.com/GPflow/tensorflow-intersphinx/raw/master/tfp_py_objects.inv")}
 
 # plot_formats = [('png', 80)]
 plot_html_show_formats = False
@@ -118,7 +121,7 @@ def setup(app):
 # further.  For a list of options available for each theme, see the
 # documentation.
 html_theme_options = {
-    'collapse_navigation': True,    
+    'collapse_navigation': True,
     }
 # Add any paths that contain custom themes here, relative to this directory.
 # html_theme_path = []

docs/source/index.rst

@@ -19,6 +19,7 @@ The best place to start is by looking at the worked examples below.
 
    Quantify measured parameters <https://tensiometer.readthedocs.org/en/latest/measured_parameters.html>
    Quantify convergence of chains <https://tensiometer.readthedocs.org/en/latest/chains_convergence_test.html>
+   Compute non-Gaussian tension estimates <https://tensiometer.readthedocs.org/en/latest/non_gaussian_tension.html>
 
 Other examples are provided for specific tasks.
 

docs/source/mcmc_tension.rst

@@ -3,5 +3,13 @@ tensiometer.mcmc_tension
 
 .. currentmodule:: tensiometer.mcmc_tension
 
-.. automodule:: tensiometer.mcmc_tension
+.. automodule:: tensiometer.mcmc_tension.__init__
+
+.. automodule:: tensiometer.mcmc_tension.param_diff
+   :members:
+
+.. automodule:: tensiometer.mcmc_tension.kde
+   :members:
+
+.. automodule:: tensiometer.mcmc_tension.flow
    :members:

requirements.txt

@@ -1,9 +1,12 @@
 GetDist>=1.1.2
-numpy
+numpy>=1.19.1
 matplotlib>=2.2.0
-scipy>=1.0.0
-joblib
-coverage
+scipy>=1.6.2
+joblib>=1.0.1
+coverage>=5.5
 tqdm
 autograd
-pymanopt
+pymanopt>=0.2.5
+tensorflow>=2.0.0
+tensorflow_probability>=0.11.1
+numba>=0.53.0

setup.py

@@ -35,12 +35,17 @@ def get_long_description():
 
 # get requirements:
 def get_requirements():
-    with open('README.rst',  encoding='utf-8-sig') as f:
+    requirements = []
+    with open('requirements.txt',  encoding='utf-8-sig') as f:
         lines = f.readlines()
-        i = -1
-        while '=====' not in lines[i]:
-            i -= 1
-        return ''.join(lines[:i])
+        for line in lines:
+            line = line.rstrip()
+            if '>=' in line:
+                pc, ver = line.split('>=')
+                requirements.append(pc+' (>='+ver+')')
+            else:
+                requirements.append(line)
+    return requirements
 
 
 # setup:
@@ -59,27 +64,16 @@ def get_requirements():
                     },
                  packages=setuptools.find_packages(),
                  platforms='any',
-                 install_requires=[
-                    'GetDist (>=1.1.2)',
-                    'numpy',
-                    'matplotlib (>=2.2.0)',
-                    'scipy (>=1.0.0)',
-                    'joblib',
-                    'coverage',
-                    'tqdm',
-                    'autograd',
-                    'pymanopt',
-                    ],
+                 install_requires=get_requirements(),
                  classifiers=[
-                    'Development Status :: 2 - Pre-Alpha',
+                    'Development Status :: 3 - Alpha',
                     'Operating System :: OS Independent',
                     'Intended Audience :: Science/Research',
                     'Programming Language :: Python :: 3',
-                    'Programming Language :: Python :: 3.6',
                     'Programming Language :: Python :: 3.7',
                     'Programming Language :: Python :: 3.8',
                     ],
-                 python_requires='>=3.6',
+                 python_requires='>=3.7',
                  zip_safe=False,
                  keywords=['MCMC']
                  )

tensiometer/__init__.py

@@ -1,3 +1,5 @@
 __author__ = 'Marco Raveri'
 __version__ = "0.0.8"
 __url__ = "https://tensiometer.readthedocs.io"
+
+from . import gaussian_tension, mcmc_tension, cosmosis_interface

tensiometer/cosmosis_interface.py

@@ -43,7 +43,7 @@ def MCSamplesFromCosmosis(chain_root, chain_min_root=None,
     :param chain_min_root: (optional) name of the file containing the
         explicit best fit.
     :param param_name_dict: (optional) a dictionary with the mapping between
-        cosmosis names and reasonable names.
+        cosmosis names and reasonable parameter names.
     :param param_label_dict: (optional) dictionary with the mapping between
         parameter names and parameter labels, since Cosmosis does not save
         the labels in the chain.
@@ -80,7 +80,7 @@ def MCSamplesFromCosmosis(chain_root, chain_min_root=None,
     if sampler == 'nested':
         # get number of samples to use:
         nsamples = int(list(filter(lambda x: 'nsample=' in x, info))
-                       [0].split('=')[1])
+                       [0].replace(' ', '').split('=')[1])
         # get the chain:
         chain = chain[-nsamples:]
         # get all quantities:
@@ -96,6 +96,19 @@ def MCSamplesFromCosmosis(chain_root, chain_min_root=None,
             param_labels.pop(param_names.index('post'))
         param_names.pop(param_names.index('weight'))
         param_names.pop(param_names.index('post'))
+    elif sampler == 'mcmc':
+        # get all quantities:
+        indexes = [i for i in range(len(param_names))
+                   if i != param_names.index('post')]
+        samples = chain[:, indexes]
+        loglike = chain[:, param_names.index('post')]
+        # Cosmosis does not weight samples:
+        samples, idx, weights = np.unique(samples, return_index=True, return_counts=True, axis=0)
+        loglike = loglike[idx]
+        # delete the weights and likelihood from names:
+        if param_labels is not None:
+            param_labels.pop(param_names.index('post'))
+        param_names.pop(param_names.index('post'))
     elif sampler == 'uncorrelated':
         # get all quantities:
         indexes = [i for i in range(len(param_names))
@@ -138,9 +151,6 @@ def MCSamplesFromCosmosis(chain_root, chain_min_root=None,
     mc_samples = polish_samples(mc_samples)
     # get the best fit:
     if chain_min_root is not None:
-        #mc_samples.bestfit = get_maximum_likelihood(chain_min_root,
-        #                                            param_name_dict,
-        #                                            param_label_dict)
         # since getdist does not cache the best fit we have to override the
         # method in this brute way:
         funcType = types.MethodType
@@ -221,13 +231,14 @@ def get_sampler_type(info):
                     'polychord': 'nested',
                     'multinest': 'nested',
                     'apriori': 'uncorrelated',
+                    'emcee': 'mcmc',
                     'pmaxlike': 'max_like',
                     'maxlike': 'max_like'
                     }
     # find the sampler in the parameters:
     temp = list(filter(lambda x: 'sampler' in x, info))
     if len(temp) > 0:
-        sampler = temp[0].split('=')[1].lower()
+        sampler = temp[0].replace(' ', '').split('=')[1].lower()
         if sampler in sampler_dict.keys():
             sampler = sampler_dict[sampler]
         else:
@@ -248,7 +259,7 @@ def get_name_tag(info):
     """
     temp = list(filter(lambda x: 'run_name' in x, info))
     if len(temp) > 0:
-        name_tag = temp[0].split('=')[1].lower()
+        name_tag = temp[0].replace(' ', '').split('=')[1].lower()
         name_tag = name_tag.rstrip().lstrip()
     else:
         name_tag = None