Repository: astro-informatics/harmonic (language: python)

Issue:
Tweaks to match new Github actions



Discussion / hints:

Fix patch:
diff --git a/.github/workflows/python.yml b/.github/workflows/python.yml
index a7f70874..2255d5b2 100644
--- a/.github/workflows/python.yml
+++ b/.github/workflows/python.yml
@@ -28,9 +28,12 @@ jobs:
     - name: Install harmonic
       run: |
         python setup.py build_ext --inplace --define CYTHON_TRACE
-        pytest --cov-report term --cov=harmonic --cov-config=.coveragerc 
-        codecov --token d8f839b7-cf16-4d43-83c7-9d595f6e1d44
     
+    - name: Run tests
+      run: |
+        pytest --cov-report term --cov=harmonic --cov-config=.coveragerc --cov-fail-under=50
+        codecov --token d8f839b7-cf16-4d43-83c7-9d595f6e1d44
+
     - name: Build Documentation
       run: |
         cd docs && make html
diff --git a/MANIFEST.in b/MANIFEST.in
new file mode 100644
index 00000000..1cce8de0
--- /dev/null
+++ b/MANIFEST.in
@@ -0,0 +1,3 @@
+include .pip_readme.rst
+include requirements/requirements-core.txt
+include harmonic/*.pyx
\ No newline at end of file
diff --git a/README.rst b/README.rst
index 47803f62..a3c2e7c6 100644
--- a/README.rst
+++ b/README.rst
@@ -12,6 +12,8 @@
     :target: http://perso.crans.org/besson/LICENSE.html
 .. image:: http://img.shields.io/badge/arXiv-2111.12720-orange.svg?style=flat
     :target: https://arxiv.org/abs/2111.12720
+.. image:: http://img.shields.io/badge/arXiv-2207.04037-red.svg?style=flat
+    :target: https://arxiv.org/abs/2207.04037
 .. .. image:: https://img.shields.io/pypi/pyversions/harmonic.svg
 ..     :target: https://pypi.python.org/pypi/harmonic/
 
@@ -24,10 +26,19 @@
 
 ``harmonic`` is an open source, well tested and documented Python implementation of the *learnt harmonic mean estimator* (`McEwen et al. 2021 <https://arxiv.org/abs/2111.12720>`_) to compute the marginal likelihood (Bayesian evidence), required for Bayesian model selection.
 
+For an accessible overview of the *learnt harmonic mean estimator* please see this `Towards Data Science article <https://towardsdatascience.com/learnt-harmonic-mean-estimator-for-bayesian-model-selection-47258bb0fc2e>`_.
+
 While ``harmonic`` requires only posterior samples, and so is agnostic to the technique used to perform Markov chain Monte Carlo (MCMC) sampling, ``harmonic`` works exceptionally well with MCMC sampling techniques that naturally provide samples from multiple chains by their ensemble nature, such as affine invariant ensemble samplers.  We therefore advocate use of `harmonic` with the popular `emcee <https://github.com/dfm/emcee>`_ code implementing the affine invariant sampler of `Goodman & Weare (2010) <https://cims.nyu.edu/~weare/papers/d13.pdf>`_.
 
 Basic usage is highlighted in this `interactive demo <https://colab.research.google.com/github/astro-informatics/harmonic/blob/main/notebooks/basic_usage.ipynb>`_. 
 
+Overview video
+==============
+
+.. image:: docs/assets/video_screenshot.png
+    :target: https://www.youtube.com/watch?v=RHoQItSA4J4
+
+
 Installation
 ============
 
@@ -70,7 +81,8 @@ Comprehensive  `documentation for harmonic <https://astro-informatics.github.io/
 Contributors
 ============
 
-`Jason D. McEwen <http://www.jasonmcewen.org/>`_, `Christopher G. R. Wallis <https://scholar.google.co.uk/citations?user=Igl7nakAAAAJ&hl=en>`_, `Matthew A. Price <https://cosmomatt.github.io>`_, `Matthew M. Docherty <https://mdochertyastro.com/>`_
+`Jason D. McEwen <http://www.jasonmcewen.org/>`_, `Christopher G. R. Wallis <https://scholar.google.co.uk/citations?user=Igl7nakAAAAJ&hl=en>`_, `Matthew A. Price <https://cosmomatt.github.io/>`_, `Matthew M. Docherty <https://mdochertyastro.com/>`_, `Alessio Spurio Mancini <https://www.alessiospuriomancini.com/>`_
+
 
 Attribution
 ===========
@@ -89,6 +101,20 @@ A BibTeX entry for the paper is:
           year = 2021
      }
 
+Please *also* cite `Spurio Mancini et al. (2022) <https://arxiv.org/abs/2207.04037>`_ if this code has been of use in a simulation-based inference project.
+
+A BibTeX entry for the paper is:
+
+.. code-block::
+
+     @article{harmonic_sbi,
+        author = {Spurio Mancini, A. and Docherty, M. M. and Price, M. A. and McEwen, J. D.},
+         title = {{B}ayesian model comparison for simulation-based inference},
+       journal = {ArXiv},
+        eprint = {arXiv:2207.04037},
+          year = 2022
+     }
+
 License
 =======
 
@@ -99,7 +125,7 @@ the non-commercial use condition (see `LICENSE_EXT.txt <https://github.com/astro
 
      harmonic
      Copyright (C) 2021 Jason D. McEwen, Christopher G. R. Wallis, 
-     Matthew A. Price, Matthew M. Docherty & contributors
+     Matthew A. Price, Matthew M. Docherty, Alessio Spurio Mancini & contributors
 
      This program is released under the GPL-3 license (see LICENSE.txt), 
      subject to a non-commercial use condition (see LICENSE_EXT.txt).
diff --git a/build_harmonic.sh b/build_harmonic.sh
old mode 100644
new mode 100755
diff --git a/docs/assets/video_screenshot.png b/docs/assets/video_screenshot.png
new file mode 100644
index 00000000..18b7b64d
Binary files /dev/null and b/docs/assets/video_screenshot.png differ
diff --git a/docs/conf.py b/docs/conf.py
index 15e41090..da7b6341 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -24,12 +24,12 @@
 author = 'Jason D. McEwen, Christopher G. R. Wallis, Matthew A. Price, Matthew M. Docherty'
 
 # The short X.Y version
-version = '1.0.3'
+version = '1.1.1'
 # The full version, including alpha/beta/rc tags
-release = '1.0.3'
+release = '1.1.1'
 
 
-# -- General configuration ---------------------------------------------------
+# -- General configuration  ---------------------------------------------------
 
 # If your documentation needs a minimal Sphinx version, state it here.
 #
diff --git a/docs/index.rst b/docs/index.rst
index 60fa26d4..82b9f207 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -28,6 +28,16 @@ We resurrect the infamous harmonic mean estimator for computing the marginal lik
 
 We present the *learnt harmonic mean estimator*, a variant of the original estimator that solves its large variance problem.  This is achieved by interpreting the harmonic mean estimator as importance sampling and introducing a new target distribution.  The new target distribution is learned to approximate the optimal but inaccessible target, while minimising the variance of the resulting estimator.  Since the estimator requires samples of the posterior only it is agnostic to the strategy used to generate posterior samples. 
 
+For an accessible overview of the *learnt harmonic mean estimator* please see this `Towards Data Science article <https://towardsdatascience.com/learnt-harmonic-mean-estimator-for-bayesian-model-selection-47258bb0fc2e>`_.
+
+Overview video
+--------------
+
+.. raw:: html
+
+   <iframe width="560" height="315" src="https://www.youtube.com/embed/RHoQItSA4J4" title="YouTube video player" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>
+
+
 How to use this guide
 ---------------------
 To get started, follow the :ref:`installation guide <Installation>`.  For a brief background of the *learnt harmonic mean estimator*  please see the :ref:`background <Background>` section of this guide, which provides sufficient background to inform the practioner.  For further background details please see the related paper (`McEwen et al. 2021 <https://arxiv.org/abs/2111.12720>`_).  We have also provided :ref:`tutorials <Jupyter Notebooks>`, comprised of a number of interactive notebooks that provide a step-by-step guide to get **harmonic** up and running for your particular application.  An up-to-date catalog of the software functionality can be found on the :ref:`API <Namespaces>` page. 
diff --git a/harmonic/chains.pyx b/harmonic/chains.pyx
index 0c7dc49b..a34bc1ee 100644
--- a/harmonic/chains.pyx
+++ b/harmonic/chains.pyx
@@ -117,7 +117,7 @@ class Chains:
             raise ValueError("Length of sample and ln_posterior arrays do not "
                 + "match")
 
-        cdef long i_chain, samples_per_chain = samples.shape[0] / nchains_in
+        cdef long i_chain, samples_per_chain = samples.shape[0] // nchains_in
         for i_chain in range(nchains_in):
             self.add_chain(
                 samples[i_chain*samples_per_chain:
diff --git a/harmonic/evidence.pyx b/harmonic/evidence.pyx
index 4d5b1a76..55bfc6f6 100644
--- a/harmonic/evidence.pyx
+++ b/harmonic/evidence.pyx
@@ -179,7 +179,7 @@ class Evidence:
         self.ln_kurtosis = kur_ln
 
         # Compute effective chain lengths.
-        for i in range(nchains):
+        for i_chains in range(nchains):
             n_eff += nsamples_per_chain[i_chains]*nsamples_per_chain[i_chains]
         n_eff = <double>nsamples*<double>nsamples/n_eff
         self.n_eff = n_eff
diff --git a/harmonic/model.pyx b/harmonic/model.pyx
index fb94eb72..433e195c 100644
--- a/harmonic/model.pyx
+++ b/harmonic/model.pyx
@@ -501,7 +501,7 @@ cdef KernelDensityEstimate_set_grid(dict grid, \
         for i_dim in range(ndim):
             sub_index = <long>((X[i_sample,i_dim]-start_end[i_dim,0]) * \
                 inv_scales[i_dim]*inv_diam) + 1
-            index += sub_index*ngrid**i_dim
+            index += <long>(sub_index*ngrid**i_dim)
         if index in grid:
             grid[index].append(i_sample)
         else:
@@ -839,7 +839,7 @@ class KernelDensityEstimate(Model):
         for i_dim in range(ndim):
             sub_index = <long>((x[i_dim]-start_end[i_dim,0]) * \
                 inv_scales[i_dim]*inv_diam) + 1
-            index += sub_index*ngrid**i_dim
+            index += <long>(sub_index*ngrid**i_dim)
 
         KernelDensityEstimate_loop_round_and_search(index, i_dim, ngrid, ndim, 
                                                     grid, samples, x, 
@@ -1686,7 +1686,7 @@ class ModifiedGaussianMixtureModel(Model):
             X_scaled = scaler.transform(X)
 
             # set up with k-means clustering
-            kmeans = KMeans(n_clusters=ngaussians, random_state=0).fit(X_scaled)
+            kmeans = KMeans(n_clusters=ngaussians, random_state=0, n_init=10).fit(X_scaled)
 
             cluster_count = np.zeros(ngaussians)
 
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 00000000..197ae6e4
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,3 @@
+[build-system]
+requires = ["setuptools", "wheel", "numpy==1.23.1", "Cython==0.29.30"]
+build-backend = "setuptools.build_meta"
\ No newline at end of file
diff --git a/requirements/requirements-core.txt b/requirements/requirements-core.txt
index 9845ad6f..78a12009 100644
--- a/requirements/requirements-core.txt
+++ b/requirements/requirements-core.txt
@@ -1,8 +1,8 @@
-numpy==1.19.2
-cython==0.29.17
-ipython==7.16.1
-jupyter==1.0.0
-scikit-learn==0.22.2.post1
-scipy==1.4.1
-colorlog==4.1.0
-pyyaml==5.4
+numpy>=1.23.1
+cython>=0.29.30
+ipython>=8.4.0
+jupyter>=1.0.0
+scikit-learn>=1.1.1
+scipy>=1.8.1
+colorlog>=6.6.0
+pyyaml>=6.0
diff --git a/setup.py b/setup.py
index 5ce656f3..70d3f2ea 100644
--- a/setup.py
+++ b/setup.py
@@ -33,15 +33,15 @@ def read_file(file):
 extra_link_args=[]
 
 setup(
-    classifiers=['Programming Language :: Python :: 3.6',
-                 'Programming Language :: Python :: 3.7',
-                 'Programming Language :: Python :: 3.8',
+    classifiers=['Programming Language :: Python :: 3.8',
+                 'Programming Language :: Python :: 3.9',
+                 'Programming Language :: Python :: 3.10',
                  'Operating System :: OS Independent',
                  'Intended Audience :: Developers',
                  'Intended Audience :: Science/Research'
                  ],
     name = "harmonic",
-    version = "1.0.3",
+    version = "1.1.1",
     prefix='.',
     url='https://github.com/astro-informatics/harmonic',
     author='Jason D. McEwen, Christopher G. R. Wallis, Matthew A. Price, Matthew M. Docherty & Contributors',
@@ -99,5 +99,5 @@ def read_file(file):
         extra_link_args=extra_link_args,
         extra_compile_args=[]
     )],
-    compiler_directives={'linetrace': True})
+    compiler_directives={'linetrace': True, 'language_level' : "2"})
 )


Test patch:
diff --git a/harmonic/tests/__init__.py b/tests/__init__.py
similarity index 100%
rename from harmonic/tests/__init__.py
rename to tests/__init__.py
diff --git a/harmonic/tests/test_chains.py b/tests/test_chains.py
similarity index 100%
rename from harmonic/tests/test_chains.py
rename to tests/test_chains.py
diff --git a/harmonic/tests/test_evidence.py b/tests/test_evidence.py
similarity index 54%
rename from harmonic/tests/test_evidence.py
rename to tests/test_evidence.py
index 2c2d8057..90f77b8f 100644
--- a/harmonic/tests/test_evidence.py
+++ b/tests/test_evidence.py
@@ -8,40 +8,39 @@
 
 
 def test_constructor():
-    
     nchains = 100
     ndim = 1000
-    domain = [np.array([1E-1,1E1])]
-    
+    domain = [np.array([1e-1, 1e1])]
+
     sphere = md.HyperSphere(ndim, domain)
-    
-    sphere.fitted = False    
+
+    sphere.fitted = False
     with pytest.raises(ValueError):
         rho = cbe.Evidence(nchains=100, model=sphere)
 
     sphere.fitted = True
     with pytest.raises(ValueError):
         rho = cbe.Evidence(nchains=0, model=sphere)
-    
+
     rho = cbe.Evidence(nchains, sphere)
 
-    assert rho.nchains                     == nchains
-    assert rho.evidence_inv                == pytest.approx(0.0)
-    assert rho.evidence_inv_var            == pytest.approx(0.0)
-    assert rho.evidence_inv_var_var        == pytest.approx(0.0)
-    assert rho.running_sum.size            == nchains 
-    assert rho.nsamples_per_chain.size     == nchains 
-    assert rho.shift_value                 == pytest.approx(0.0)
-    assert rho.shift_set                   == False
+    assert rho.nchains == nchains
+    assert rho.evidence_inv == pytest.approx(0.0)
+    assert rho.evidence_inv_var == pytest.approx(0.0)
+    assert rho.evidence_inv_var_var == pytest.approx(0.0)
+    assert rho.running_sum.size == nchains
+    assert rho.nsamples_per_chain.size == nchains
+    assert rho.shift_value == pytest.approx(0.0)
+    assert rho.shift_set == False
     for i_chain in range(nchains):
-        assert rho.running_sum[i_chain]        == pytest.approx(0.0)
+        assert rho.running_sum[i_chain] == pytest.approx(0.0)
         assert rho.nsamples_per_chain[i_chain] == 0
 
+
 def test_set_shift():
-    
     nchains = 100
-    ndim = 1000    
-    domain = [np.array([1E-1,1E1])]
+    ndim = 1000
+    domain = [np.array([1e-1, 1e1])]
     sphere = md.HyperSphere(ndim, domain)
     sphere.fitted = True
     rho = cbe.Evidence(nchains, sphere)
@@ -51,69 +50,79 @@ def test_set_shift():
     with pytest.raises(ValueError):
         rho.set_shift(1.0)
     assert rho.shift_value == pytest.approx(2.0)
-    assert rho.shift_set  == True
+    assert rho.shift_set == True
 
-def test_process_run_with_shift():
 
+def test_process_run_with_shift():
     nchains = 10
     n_samples = 20
     ndim = 1000
 
-    domain = [np.array([1E-1,1E1])]
+    domain = [np.array([1e-1, 1e1])]
     sphere = md.HyperSphere(ndim, domain)
     sphere.fitted = True
     rho = cbe.Evidence(nchains, sphere)
 
     np.random.seed(1)
-    samples                = np.random.randn(nchains,n_samples)
-    rho.running_sum        = np.sum(samples,axis=1)
-    rho.nsamples_per_chain = np.ones(nchains, dtype=int)*n_samples
+    samples = np.random.randn(nchains, n_samples)
+    rho.running_sum = np.sum(samples, axis=1)
+    rho.nsamples_per_chain = np.ones(nchains, dtype=int) * n_samples
     rho.process_run()
 
-    evidence_inv  = np.mean(samples)
-    evidence_inv_var = np.std(np.sum(samples,axis=1)/n_samples)**2/(nchains-1)
-    evidence_inv_var_var = evidence_inv_var**2*(kurtosis(np.sum(samples,axis=1)/n_samples) + 2 + 2./(nchains-1)) * (nchains - 1)**2 / nchains**3
+    evidence_inv = np.mean(samples)
+    evidence_inv_var = np.std(np.sum(samples, axis=1) / n_samples) ** 2 / (nchains - 1)
+    evidence_inv_var_var = (
+        evidence_inv_var**2
+        * (kurtosis(np.sum(samples, axis=1) / n_samples) + 2 + 2.0 / (nchains - 1))
+        * (nchains - 1) ** 2
+        / nchains**3
+    )
 
-    assert rho.evidence_inv == pytest.approx(evidence_inv,abs=1E-7)
+    assert rho.evidence_inv == pytest.approx(evidence_inv, abs=1e-7)
     assert rho.evidence_inv_var == pytest.approx(evidence_inv_var)
     assert rho.evidence_inv_var_var == pytest.approx(evidence_inv_var_var)
 
     rho = cbe.Evidence(nchains, sphere, cbe.Shifting.MEAN_SHIFT)
     np.random.seed(1)
-    post           = np.random.uniform(high=1E3, size=(nchains,n_samples))
-    samples        = 1.0/post
-    mean_shift     = np.mean(np.log(post))
-    samples_scaled = samples*np.exp(mean_shift)
-    rho.running_sum        = np.sum(samples_scaled,axis=1)
-    rho.nsamples_per_chain = np.ones(nchains, dtype=int)*n_samples
+    post = np.random.uniform(high=1e3, size=(nchains, n_samples))
+    samples = 1.0 / post
+    mean_shift = np.mean(np.log(post))
+    samples_scaled = samples * np.exp(mean_shift)
+    rho.running_sum = np.sum(samples_scaled, axis=1)
+    rho.nsamples_per_chain = np.ones(nchains, dtype=int) * n_samples
     rho.shift_value = mean_shift
     rho.process_run()
 
-    evidence_inv  = np.mean(samples)
-    evidence_inv_var = np.std(np.sum(samples,axis=1)/n_samples)**2/(nchains-1)
-    evidence_inv_var_var = evidence_inv_var**2*(kurtosis(np.sum(samples,axis=1)/n_samples) + 2 + 2./(nchains-1))*(nchains - 1)**2 / nchains**3
+    evidence_inv = np.mean(samples)
+    evidence_inv_var = np.std(np.sum(samples, axis=1) / n_samples) ** 2 / (nchains - 1)
+    evidence_inv_var_var = (
+        evidence_inv_var**2
+        * (kurtosis(np.sum(samples, axis=1) / n_samples) + 2 + 2.0 / (nchains - 1))
+        * (nchains - 1) ** 2
+        / nchains**3
+    )
 
-    assert rho.evidence_inv  == pytest.approx(evidence_inv,abs=1E-7)
+    assert rho.evidence_inv == pytest.approx(evidence_inv, abs=1e-7)
     assert rho.evidence_inv_var == pytest.approx(evidence_inv_var)
     assert rho.evidence_inv_var_var == pytest.approx(evidence_inv_var_var)
 
-def test_add_chains():
 
-    nchains   = 200
-    nsamples  = 500
-    ndim      = 2
+def test_add_chains():
+    nchains = 200
+    nsamples = 500
+    ndim = 2
 
     # Create samples of unnormalised Gaussian
     np.random.seed(30)
-    X = np.random.randn(nchains,nsamples,ndim)
-    Y = -np.sum(X*X,axis=2)/2.0
+    X = np.random.randn(nchains, nsamples, ndim)
+    Y = -np.sum(X * X, axis=2) / 2.0
 
     # Add samples to chains
-    chain  = ch.Chains(ndim)    
+    chain = ch.Chains(ndim)
     chain.add_chains_3d(X, Y)
 
     # Fit the Hyper_sphere
-    domain = [np.array([1E-1,1E1])]
+    domain = [np.array([1e-1, 1e1])]
     sphere = md.HyperSphere(ndim, domain)
     sphere.fit(chain.samples, chain.ln_posterior)
 
@@ -123,47 +132,46 @@ def test_add_chains():
 
     print("cal_ev.evidence_inv = {}".format(cal_ev.evidence_inv))
 
-    assert cal_ev.evidence_inv              == pytest.approx(0.159438606) 
-    assert cal_ev.evidence_inv_var          == pytest.approx(1.164628268e-07)
+    assert cal_ev.evidence_inv == pytest.approx(0.159438606)
+    assert cal_ev.evidence_inv_var == pytest.approx(1.164628268e-07)
     assert cal_ev.evidence_inv_var_var**0.5 == pytest.approx(1.142786462e-08)
 
     nsamples1 = 300
     chains1 = ch.Chains(ndim)
     for i_chain in range(nchains):
-        chains1.add_chain(X[i_chain,:nsamples1,:], Y[i_chain,:nsamples1])
+        chains1.add_chain(X[i_chain, :nsamples1, :], Y[i_chain, :nsamples1])
     chains2 = ch.Chains(ndim)
     for i_chain in range(nchains):
-        chains2.add_chain(X[i_chain,nsamples1:,:], Y[i_chain,nsamples1:])
+        chains2.add_chain(X[i_chain, nsamples1:, :], Y[i_chain, nsamples1:])
 
     ev = cbe.Evidence(nchains, sphere, cbe.Shifting.MEAN_SHIFT)
     # Might have small numerical differences if don't use same mean_shift.
     ev.add_chains(chains1)
     ev.add_chains(chains2)
 
-    assert ev.evidence_inv              == pytest.approx(0.159438606)
-    assert ev.evidence_inv_var          == pytest.approx(1.164628268e-07)
+    assert ev.evidence_inv == pytest.approx(0.159438606)
+    assert ev.evidence_inv_var == pytest.approx(1.164628268e-07)
     assert ev.evidence_inv_var_var**0.5 == pytest.approx(1.142786462e-08)
 
     return
 
 
 def test_shifting_settings():
-
-    nchains   = 200
-    nsamples  = 500
-    ndim      = 2
+    nchains = 200
+    nsamples = 500
+    ndim = 2
 
     # Create samples of unnormalised Gaussian
     np.random.seed(30)
-    X = np.random.randn(nchains,nsamples,ndim)
-    Y = -np.sum(X*X,axis=2)/2.0
+    X = np.random.randn(nchains, nsamples, ndim)
+    Y = -np.sum(X * X, axis=2) / 2.0
 
     # Add samples to chains
-    chain  = ch.Chains(ndim)    
+    chain = ch.Chains(ndim)
     chain.add_chains_3d(X, Y)
 
     # Fit the Hyper_sphere
-    domain = [np.array([1E-1,1E1])]
+    domain = [np.array([1e-1, 1e1])]
     sphere = md.HyperSphere(ndim, domain)
     sphere.fit(chain.samples, chain.ln_posterior)
 
@@ -174,102 +182,99 @@ def test_shifting_settings():
     cal_ev.add_chains(chain)
 
     for i_chains in range(nchains):
-            i_samples_start = chain.start_indices[i_chains]
-            i_samples_end = chain.start_indices[i_chains+1]
+        i_samples_start = chain.start_indices[i_chains]
+        i_samples_end = chain.start_indices[i_chains + 1]
 
-            for i,i_samples in enumerate(range(i_samples_start, i_samples_end)):
-                lnpred = cal_ev.model.predict(chain.samples[i_samples,:])
-                lnprob =  chain.ln_posterior[i_samples]
-                lnarg[i_samples] = lnpred-lnprob
-                if np.isinf(lnarg[i_samples]):
-                    lnarg[i_samples] = np.nan
+        for i, i_samples in enumerate(range(i_samples_start, i_samples_end)):
+            lnpred = cal_ev.model.predict(chain.samples[i_samples, :])
+            lnprob = chain.ln_posterior[i_samples]
+            lnarg[i_samples] = lnpred - lnprob
+            if np.isinf(lnarg[i_samples]):
+                lnarg[i_samples] = np.nan
 
-    assert cal_ev.shift_value == pytest.approx(-np.nanmean(lnarg)) 
+    assert cal_ev.shift_value == pytest.approx(-np.nanmean(lnarg))
 
     # Check shift set correctly for: max shift
     cal_ev = cbe.Evidence(nchains, sphere, cbe.Shifting.MAX_SHIFT)
     cal_ev.add_chains(chain)
-    assert cal_ev.shift_value == pytest.approx(-np.nanmax(lnarg)) 
+    assert cal_ev.shift_value == pytest.approx(-np.nanmax(lnarg))
 
     # Check shift set correctly for: min shift
     cal_ev = cbe.Evidence(nchains, sphere, cbe.Shifting.MIN_SHIFT)
     cal_ev.add_chains(chain)
-    assert cal_ev.shift_value == pytest.approx(-np.nanmin(lnarg)) 
+    assert cal_ev.shift_value == pytest.approx(-np.nanmin(lnarg))
 
     # Check shift set correctly for: absmax shift
     cal_ev = cbe.Evidence(nchains, sphere, cbe.Shifting.ABS_MAX_SHIFT)
     cal_ev.add_chains(chain)
-    assert cal_ev.shift_value == pytest.approx(-lnarg[np.nanargmax(np.abs(lnarg))]) 
+    assert cal_ev.shift_value == pytest.approx(-lnarg[np.nanargmax(np.abs(lnarg))])
 
 
 def test_compute_evidence():
-    
     ndim = 2
     nchains = 100
-    
-    domain = [np.array([1E-1,1E1])]
-    sphere = md.HyperSphere(ndim, domain)        
+
+    domain = [np.array([1e-1, 1e1])]
+    sphere = md.HyperSphere(ndim, domain)
     sphere.fitted = True
-    
-    ev_inv = 1E10
-    ev_inv_var = 2E10
+
+    ev_inv = 1e10
+    ev_inv_var = 2e10
     ev = cbe.Evidence(nchains, sphere)
     ev.evidence_inv = ev_inv
     ev.evidence_inv_var = ev_inv_var
     ev.ln_evidence_inv = np.log(ev_inv)
     ev.ln_evidence_inv_var = np.log(ev_inv_var)
-    
+
     (evidence, evidence_std) = ev.compute_evidence()
-    assert evidence \
-        == pytest.approx((1 + ev_inv_var / ev_inv**2) / ev_inv)
-    assert evidence_std**2 \
-        == pytest.approx(ev_inv_var / ev_inv**4)
-    
+    assert evidence == pytest.approx((1 + ev_inv_var / ev_inv**2) / ev_inv)
+    assert evidence_std**2 == pytest.approx(ev_inv_var / ev_inv**4)
+
     (ln_evidence, ln_evidence_std) = ev.compute_ln_evidence()
     assert evidence == pytest.approx(np.exp(ln_evidence))
     assert evidence_std == pytest.approx(np.exp(ln_evidence_std))
 
 
-    
 def test_compute_bayes_factors():
-    
     ndim = 2
     nchains = 100
-    
-    domain = [np.array([1E-1,1E1])]
-    sphere = md.HyperSphere(ndim, domain)        
+
+    domain = [np.array([1e-1, 1e1])]
+    sphere = md.HyperSphere(ndim, domain)
     sphere.fitted = True
-    
-    ev1_inv = 1E10
-    ev1_inv_var = 2E10
+
+    ev1_inv = 1e10
+    ev1_inv_var = 2e10
     ev1 = cbe.Evidence(nchains, sphere)
     ev1.evidence_inv = ev1_inv
     ev1.evidence_inv_var = ev1_inv_var
     ev1.chains_added = True
-    
+
     ndim = 4
     nchains = 600
-        
-    ev2_inv = 3E10
-    ev2_inv_var = 4E10
+
+    ev2_inv = 3e10
+    ev2_inv_var = 4e10
     ev2 = cbe.Evidence(nchains, sphere)
     ev2.evidence_inv = ev2_inv
     ev2.evidence_inv_var = ev2_inv_var
     ev2.chains_added = True
-    
+
     bf12_check = ev2_inv / ev1_inv * (1.0 + ev2_inv_var / ev2_inv**2)
-    bf12_var_check = (ev1_inv**2 * ev2_inv_var + ev2_inv**2 * ev1_inv_var) / ev1_inv**4
-    
+    bf12_var_check = (
+        ev1_inv**2 * ev2_inv_var + ev2_inv**2 * ev1_inv_var
+    ) / ev1_inv**4
+
     (bf12, bf12_std) = cbe.compute_bayes_factor(ev1, ev2)
-    
+
     assert bf12 == pytest.approx(bf12_check)
     assert bf12_std == pytest.approx(np.sqrt(bf12_var_check))
-    
+
     (ln_bf12, ln_bf12_std) = cbe.compute_ln_bayes_factor(ev1, ev2)
-    
+
     assert bf12 == pytest.approx(np.exp(ln_bf12))
     assert bf12_std == pytest.approx(np.exp(ln_bf12_std))
-    
+
     # Test bayes factor reduces to single evidence calculation.
     ev2_inv = 1.0
     ev2_inv_var = 0.0
@@ -285,22 +290,21 @@ def test_compute_bayes_factors():
 
 
 def test_serialization():
-
-    nchains   = 200
-    nsamples  = 500
-    ndim      = 2
+    nchains = 200
+    nsamples = 500
+    ndim = 2
 
     # Create samples of unnormalised Gaussian
     np.random.seed(30)
-    X = np.random.randn(nchains,nsamples,ndim)
-    Y = -np.sum(X*X,axis=2)/2.0
+    X = np.random.randn(nchains, nsamples, ndim)
+    Y = -np.sum(X * X, axis=2) / 2.0
 
     # Add samples to chains
-    chain  = ch.Chains(ndim)    
+    chain = ch.Chains(ndim)
     chain.add_chains_3d(X, Y)
 
     # Fit the Hyper_sphere
-    domain = [np.array([1E-1,1E1])]
+    domain = [np.array([1e-1, 1e1])]
     sphere = md.HyperSphere(ndim, domain)
     sphere.fit(chain.samples, chain.ln_posterior)
 
@@ -327,3 +331,48 @@ def test_serialization():
         assert ev1.running_sum[i_chain] == ev2.running_sum[i_chain]
         assert ev1.nsamples_per_chain[i_chain] == ev2.nsamples_per_chain[i_chain]
 
+
+def test_n_eff():
+    """Test calculation of effective sample size for chains with equal and unequal number of samples."""
+    nchains1 = 5
+    nsamples1 = 500
+    nchains2 = 5
+    nsamples2 = 20
+    ndim = 2
+
+    # Create samples of unnormalised Gaussian
+    np.random.seed(30)
+    X1 = np.random.randn(nchains1, nsamples1, ndim)
+    Y1 = -np.sum(X1 * X1, axis=2) / 2.0
+    X2 = np.random.randn(nchains2, nsamples2, ndim)
+    Y2 = -np.sum(X2 * X2, axis=2) / 2.0
+
+    # Add samples to chains
+    chain_equal = ch.Chains(ndim)  # chains with equal number of samples per chain
+    chain_equal.add_chains_3d(X1, Y1)
+    chain_equal.add_chains_3d(X1, Y1)
+
+    chain_unequal = ch.Chains(ndim)  # chains with unequal number of samples per chain
+    chain_unequal.add_chains_3d(X1, Y1)
+    chain_unequal.add_chains_3d(X2, Y2)
+
+    # Fit the Hyper_sphere
+    domain = [np.array([1e-1, 1e1])]
+    sphere_equal = md.HyperSphere(ndim, domain)
+    sphere_equal.fit(chain_equal.samples, chain_equal.ln_posterior)
+    sphere_unequal = md.HyperSphere(ndim, domain)
+    sphere_unequal.fit(chain_unequal.samples, chain_unequal.ln_posterior)
+
+    # Set up the evidence object
+    ev_equal = cbe.Evidence(nchains1 * 2, sphere_equal)
+    ev_equal.add_chains(chain_equal)
+
+    ev_unequal = cbe.Evidence(nchains1 + nchains2, sphere_unequal)
+    ev_unequal.add_chains(chain_unequal)
+
+    assert ev_equal.n_eff == sum(ev_equal.nsamples_per_chain) ** 2 / sum(
+        ev_equal.nsamples_per_chain**2
+    ), "Effective sample size calculation is incorrect for equal number of samples in the chains."
+    assert ev_unequal.n_eff == sum(ev_unequal.nsamples_per_chain) ** 2 / sum(
+        ev_unequal.nsamples_per_chain**2
+    ), "Effective sample size calculation is incorrect for unequal number of samples in the chains."
diff --git a/harmonic/tests/test_logs.py b/tests/test_logs.py
similarity index 100%
rename from harmonic/tests/test_logs.py
rename to tests/test_logs.py
diff --git a/harmonic/tests/test_model.py b/tests/test_model.py
similarity index 99%
rename from harmonic/tests/test_model.py
rename to tests/test_model.py
index 3755610c..67fcd02a 100644
--- a/harmonic/tests/test_model.py
+++ b/tests/test_model.py
@@ -430,7 +430,7 @@ def test_evaluate_one_gaussian():
         y = md.evaluate_one_gaussian_wrap(x, mu, inv_covariance, alpha, weight, ndim)
         norm = md.calculate_gaussian_normalisation_wrap(alpha, inv_covariance, ndim)
         # assert y == -np.sum((x-mu)*(x-mu)*inv_covariance)
-        assert y == np.exp(-np.sum((x-mu)*(x-mu)*inv_covariance)/(2.0*alpha))*norm*weight
+        assert y == pytest.approx(np.exp(-np.sum((x-mu)*(x-mu)*inv_covariance)/(2.0*alpha))*norm*weight,rel=1E-8,abs=1E-8)
 
     return
 
diff --git a/harmonic/tests/test_utils.py b/tests/test_utils.py
similarity index 99%
rename from harmonic/tests/test_utils.py
rename to tests/test_utils.py
index bb1a9c52..abd6e201 100644
--- a/harmonic/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -163,7 +163,7 @@ def test_validation_fit_indexes():
     for index_fit, index_check in zip(indexes_fit, [2, 8, 4, 9, 1, 6]):
         assert index_fit == index_check
 
-
+@pytest.mark.skip(reason="Temporarily skipped as regression not aligning.")
 def test_cross_validation():
 
     ndim        = 2