Repository: alan-turing-institute/mogp-emulator (language: python)

Issue:
Remove unused MCMC code
Remove the old code used for MCMC sampling as it is unused at the moment. Also need to remove some of the benchmarks and update those files to match the current version of the code.
Remove MCMC benchmark from docs
The MCMC benchmark docs pages remain in place despite the benchmark being removed (due to temporary removal of the MCMC capabilities while the GP class was refactored). These should be removed in an update of the documentation.
Feature/pivotcholesky
Implements pivoted Cholesky decomposition as a method for handling the emulator nugget. If this option is chosen, the Cholesky factorization uses pivoting to avoid problems with singular matrices and automatically ignores points that are numerically similar enough to others such that the covariance matrix is nearly singular. Addresses issue #9.

Changes in this PR:

* Implements pivoted Cholesky decomposition in the `linalg/cholesky.py` file, plus a function for inverting pivot matrices.
* Modifications to the `GaussianProcess` class to allow pivoting as an option through setting the `nugget` attribute.
* Additional changes were required to re-order the inputs/targets to use the cached factorized covariance matrix (as the rows/columns have been re-arranged due to pivoting). This is done using the pivot matrix and its transpose as indices where appropriate in all calls to the `cho_solve` function -- once to rearrange the RHS before solving using the pivot, and again after solving with the transpose to return to the ordering in the original data.
* Includes a benchmark to illustrate how pivoting improves performance on problems with duplicate input/target values.
* Updates to documentation to reflect the above changes.


Discussion / hints:


**TODO:** `cho_solve` function combined with array reordering using the pivot should be wrapped into a function to improve legibility and ensure it hasn't been forgotten anywhere. (**Update:** Fixed)

Fix patch:
diff --git a/docs/benchmarks/benchmarks.rst b/docs/benchmarks/benchmarks.rst
index 14b881a2..9a28190b 100644
--- a/docs/benchmarks/benchmarks.rst
+++ b/docs/benchmarks/benchmarks.rst
@@ -8,65 +8,70 @@ Benchmarks
    rosenbrock
    branin
    tsunami
-   mcmc_benchmark
    mice_benchmark
    pivot_benchmark
    gkdr_benchmark
    histmatch_benchmark
 
-The code includes a series of benchmarks that illustrate various pieces of the implementation. Benchmarks
-can be run from the ``mogp_emulator/tests`` directory by entering ``make all`` or ``make benchmarks`` to
-run all benchmarks, or ``make rosenbrock``, ``make branin``, ``make tsunami``, ``make mcmc``, ``make pivot``, ``make mice``, ``make gKDR``, or ``make histmatch`` to run the individual benchmarks.
+The code includes a series of benchmarks that illustrate various
+pieces of the implementation. Benchmarks can be run from the
+``mogp_emulator/benchmarks`` directory by entering ``make all`` or ``make
+benchmarks`` to run all benchmarks, or ``make rosenbrock``, ``make
+branin``, ``make tsunami``, ``make mice``, ``make pivot``, ``make
+gKDR``, or ``make histmatch`` to run the individual benchmarks.
 
 Single Emulator Convergence Tests
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-The first benchmark examines the convergence of a single emulator applied to the Rosenbrock function in
-several different dimensions (more details can be found `here <https://www.sfu.ca/~ssurjano/rosen.html>`__).
-This illustrates how the emulator predictions improve as the number of training points is increased
-for different numbers of input parameters. The benchmark evaluates the Rosenbrock function in 4, 6, and
-8 dimensions and shows that the mean squared prediction error and the mean variance improve with the
-number of training points used. Matplotlib can optionally be used to visualize the results.
+The first benchmark examines the convergence of a single emulator
+applied to the Rosenbrock function in several different dimensions
+(more details can be found `here
+<https://www.sfu.ca/~ssurjano/rosen.html>`__).  This illustrates how
+the emulator predictions improve as the number of training points is
+increased for different numbers of input parameters. The benchmark
+evaluates the Rosenbrock function in 4, 6, and 8 dimensions and shows
+that the mean squared prediction error and the mean variance improve
+with the number of training points used. Matplotlib can optionally be
+used to visualize the results.
 
 Multi-Output Convergence Tests
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-The second benchmark examines the convergence of multiple emulators derived from the same input values.
-This benchmark is based on the 2D Branin function (more details on this function can be found `here <https://www.sfu.ca/~ssurjano/branin.html>`__)
-. The code uses 8 different realizations of the Branin
-function using different parameter values, and then examines the convergence of the 8 different
-emulators fit using different number of parameter values based on the prediction errors and
-variance values. The results can optionally be visualized using Matplotlib.
+The second benchmark examines the convergence of multiple emulators
+derived from the same input values.  This benchmark is based on the 2D
+Branin function (more details on this function can be found `here
+<https://www.sfu.ca/~ssurjano/branin.html>`__) . The code uses 8
+different realizations of the Branin function using different
+parameter values, and then examines the convergence of the 8 different
+emulators fit using different number of parameter values based on the
+prediction errors and variance values. The results can optionally be
+visualized using Matplotlib.
 
 Multi-Output Performance Benchmark
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-A performance benchmark is included that uses a set of Tsunami simulation results to examine the
-speed at which the code fits multiple emulators in parallel. The code fits 8, 16, 32, and 64 emulators
-using 1, 2, 4, and 8 processess and notes the time required to perform the fitting. Note that the results
-will depend on the number of cores on the computer -- once you exceed the number of cores, the performance
-will degrade. As with the other benchmarks, Matplotlib can optionally be used to plot the results.
-
-MCMC Benchmark
-~~~~~~~~~~~~~~
-
-A benchmark applying the software to fitting an emulator with MCMC sampling is included. The code
-draws hyperparameter samples and compares the resulting posterior distributions with the values
-found via maximum likelihood estimation. If Matplotlib is installed, a histogram of the parameter
-samples is shown.
+A performance benchmark is included that uses a set of Tsunami
+simulation results to examine the speed at which the code fits
+multiple emulators in parallel. The code fits 8, 16, 32, and 64
+emulators using 1, 2, 4, and 8 processess and notes the time required
+to perform the fitting. Note that the results will depend on the
+number of cores on the computer -- once you exceed the number of
+cores, the performance will degrade. As with the other benchmarks,
+Matplotlib can optionally be used to plot the results.
 
 MICE Benchmark
 ~~~~~~~~~~~~~~
 
-A benchmark comparing the MICE Sequential design method to Latin Hypercube sampling is also available.
-This creates designs of a variety of sizes and computes the error on unseen data for the 2D Branin
-function. It compares the accuracy of the sequential design to the Latin Hypercube for both the
-predictions and uncertainties.
+A benchmark comparing the MICE Sequential design method to Latin
+Hypercube sampling is also available.  This creates designs of a
+variety of sizes and computes the error on unseen data for the 2D
+Branin function. It compares the accuracy of the sequential design to
+the Latin Hypercube for both the predictions and uncertainties.
 
 Pivoted Cholesky Benchmark
 ~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-A benchmark shows how a duplicated input can degrade emulator performance
+This benchmark shows how a duplicated input can degrade emulator performance
 and how using pivoting can fix this more effectively than adding a nugget.
 It compares the accuracy of the emulator and its uncertainty for pivoting
 and adaptive nugget handling to the 2D Branin function.
@@ -74,12 +79,15 @@ and adaptive nugget handling to the 2D Branin function.
 Dimension Reduction Benchmark
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-A benchmark illustrating gradient-based kernel dimension reduction is available in ``tests/benchmark_kdr_GP.py``
-This problem contains a 100 dimension function with a single active dimension and shows how the loss depends
-on the number of dimensions included in the reduced dimension space.
+A benchmark illustrating gradient-based kernel dimension reduction is
+available in ``tests/benchmark_kdr_GP.py`` This problem contains a 100
+dimension function with a single active dimension and shows how the
+loss depends on the number of dimensions included in the reduced
+dimension space.
 
 History Matching Benchmark
 ~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-A benchmark illustrating use of the ``HistoryMatching`` class to rule out parts of the input space
-using a GP emulator to make predictions. The benchmark contains a 1D and a 2D example.
+A benchmark illustrating use of the ``HistoryMatching`` class to rule
+out parts of the input space using a GP emulator to make
+predictions. The benchmark contains a 1D and a 2D example.
diff --git a/docs/benchmarks/mcmc_benchmark.rst b/docs/benchmarks/mcmc_benchmark.rst
deleted file mode 100644
index 5e77b938..00000000
--- a/docs/benchmarks/mcmc_benchmark.rst
+++ /dev/null
@@ -1,7 +0,0 @@
-.. _mcmc_benchmark:
-
-**********************************
-MCMC Benchmark
-**********************************
-
-.. automodule:: mogp_emulator.benchmarks.benchmark_MCMC
\ No newline at end of file
diff --git a/docs/implementation/MCMC.rst b/docs/implementation/MCMC.rst
deleted file mode 100644
index d7e9df84..00000000
--- a/docs/implementation/MCMC.rst
+++ /dev/null
@@ -1,8 +0,0 @@
-.. _MCMC:
-
-**********************************
-The ``MCMC`` Module
-**********************************
-
-.. automodule:: mogp_emulator.MCMC
-    :members:
diff --git a/docs/implementation/implementation.rst b/docs/implementation/implementation.rst
index 523dcf3d..2ac42695 100644
--- a/docs/implementation/implementation.rst
+++ b/docs/implementation/implementation.rst
@@ -18,5 +18,4 @@ mogp_emulator Implementation Details
    ExperimentalDesign
    SequentialDesign
    HistoryMatching
-   MCMC
    linalg
diff --git a/mogp_emulator/MCMC.py b/mogp_emulator/MCMC.py
deleted file mode 100644
index 79225adf..00000000
--- a/mogp_emulator/MCMC.py
+++ /dev/null
@@ -1,261 +0,0 @@
-import numpy as np
-from numpy.linalg import LinAlgError
-from inspect import signature
-import warnings
-
-def MH_proposal(current_params, step_sizes):
-    """
-    Propose an MCMC step using a Metropolis-Hastings method
-    
-    Proposes the next point in an MCMC sampler using the Metropolis-Hastings method.
-    Inputs are the current point and a covariance matrix. The next point is drawn from
-    a multivariate normal distribution centered around the current point. The covariance
-    matrix must be a 2D ``n`` by ``n`` array, where ``n`` is the number of parameters, and
-    must be positive definite. Returns a 1D array of length ``n`` holding the new proposed
-    parameter values.
-    
-    :param current_params: Current value of the parameters. Must be a 1D array.
-    :type current_params: ndarray
-    :param step_sizes: Covariance matrix from which steps are drawn. Must be a 2D array
-                       with both dimensions the same length as ``current_params``, and
-                       must be positive definite.
-    :type step_sizes: ndarray
-    :returns: New value of parameters, a 1D array with the same length as ``current_params``
-    :rtype: ndarray
-    """
-    
-    current_params = np.array(current_params)
-    step_sizes = np.array(step_sizes)
-    
-    assert len(current_params.shape) == 1, "current parameters must be a 1D array"
-    assert len(step_sizes.shape) == 2, "step sizes must be a 2D array"
-    assert step_sizes.shape[0] == step_sizes.shape[1], "step sizes must be a square array"
-    assert len(current_params) == step_sizes.shape[0], "length of current parameters must match length of step sizes"
-    assert np.all(np.diag(step_sizes) > 0.), "step sizes must be a positive definite matrix"
-    
-    return np.random.multivariate_normal(mean=current_params, cov = step_sizes)
-
-
-def MCMC_step(loglikelihood, current_params, step_sizes, loglike_sign = 1.):
-    """
-    Method to take a weak prior Metropolis-Hastings MCMC step
-    
-    Take a Metropolis-Hastings MCMC step with given log-likelihood function and step sizes.
-    Method uses a multivariate normal distribution centered around the current parameter
-    values from which the next point is drawn, and evaluates the log-likelihood for both
-    points. If the next point has a larger log-likelihood, the step is always accepted,
-    and if the log-likelihood is less for the proposed step it is accepted with a
-    probability based on the difference between the two. If the method encounters an
-    error when evalutaing the log-likelihood, the step is rejected. Returns the next point
-    and a boolean indicating whether or not the step was accepted.
-    
-    An optional parameter ``loglike_sign`` can be passed that must be a float with the
-    value +/- 1. This is multiplied by the log-likelihood and thus allows methods that
-    compute the negative log-likelihood to be used in this routine. If values other than
-    +/- 1 are passed, the method will raise an error.
-    
-    :param loglikelihood: Log-likelihood function to be used in the MCMC step. Must be
-                          callable and must accept a single argument, which is the array
-                          holding the parameters. If this function computes the negative
-                          log-likelihood, pass ``loglike_sign = -1.`` to the function as well.
-    :type loglikelihood: function or other callable
-    :param current_params: Current value of the parameters. Must be a 1D array.
-    :type current_params: ndarray
-    :param step_sizes: Covariance matrix from which steps are drawn. Must be a 2D array
-                       with both dimensions the same length as ``current_params``, and
-                       must be positive definite.
-    :type step_sizes: ndarray
-    :param loglike_sign: Sign for the log-likelihood function. If the provided
-                         ``loglikelihood`` function computes the negative log-likelihood,
-                         pass ``-1.`` for this parameter. Optional, default value is ``1.``
-    :type loglike_sign: float
-    :returns: Proposed next point and whether or not the point is accepted as a tuple.
-              The first return item is the next point (as a 1D array with the same length\
-              as ``current_params``) and the second item is a boolean indicating whether
-              or not the step was accepted.
-    :rtype: tuple containing a ndarray and a bool
-    """
-    
-    assert callable(loglikelihood), "loglikelihood must be a callable function"
-    assert len(signature(loglikelihood).parameters) == 1
-    assert loglike_sign == 1. or loglike_sign == -1., "loglikelihood sign must be +/- 1"
-    
-    next_point = MH_proposal(current_params, step_sizes)
-
-    try:
-        H = loglike_sign*(-loglikelihood(current_params) + loglikelihood(next_point))
-    except (FloatingPointError, AssertionError, LinAlgError):
-        H = np.nan
-
-    if H >= np.log(np.random.random()) and np.isfinite(H):
-        accept = True
-    else:
-        accept = False
-
-    return next_point, accept
-
-def sample_MCMC(loglikelihood, start, step_sizes, n_samples = 1000, thin = 0, loglike_sign = 1.):
-    """
-    Draw MCMC samples for a given log-likelihood function with weak priors
-    
-    Compute an MCMC chain for a given log-likelihood function with weak priors. Function
-    requires the log-likelihood function, the starting point for the MCMC chain, and an
-    array describing the step sizes. Optional parameters are the number of steps to take,
-    how the thin the samples, and a sign for the log-likelihood function.
-    
-    The log-likelihood function must be a function or other callable that accepts a single
-    argument, which is a 1D array holding the current parameter values. The starting point
-    must be a 1D array holding the starting parameter values, which must match the length of
-    the input to the log-likelihood function. The step sizes must be a 2D array with each
-    dimension having the same length as the number of parameters, and must be positive
-    definite. The step size array is used as the covariance matrix for a multivariate
-    normal distribution, from which steps are drawn.
-    
-    Optional parameters are the number of steps to be taken (must be a positive integer).
-    Note that if the chain is thinned, the number of points in the final MCMC chain will
-    differ from the number of steps taken.
-    
-    Thinning may be specified with a non-negative integer. If a positive integer is
-    given, the chain will be thinned by only keeping every ``thin`` steps. Note that
-    ``thin = 1`` means that the chain will not be thinned. If ``thin = 0`` is given
-    (the default value), the chain will automatically be thinned by computing the
-    autocorrelation of the chain for each parameter separately and estimating the value
-    needed to eliminate correlations in the chain. If the autothinning method fails
-    (usually occurrs if the posterior is multimodal), the chain will not be thinned
-    and a warning will be given. More details on the autothinning procedure are
-    described in the corresponding function.
-    
-    An optional parameter ``loglike_sign`` can be passed that must be a float with the
-    value +/- 1. This is multiplied by the log-likelihood and thus allows methods that
-    compute the negative log-likelihood to be used in this routine.
-    
-    Returns the final thinned MCMC chain (a 2D array, where the first dimension indicates
-    the different samples and the second dimension indicates the different parameters),
-    an array holding all rejected steps (also a 2D array like the MCMC chain, useful for
-    diagnosing problems with convergence), the fraction of steps that are accepted (also
-    useful for diagnosing problems with convergence), and the first lag autocorrelation of
-    the thinned MCMC chain.
-    
-    :param loglikelihood: Log-likelihood function to be used in the MCMC step. Must be
-                          callable and must accept a single argument, which is the array
-                          holding the parameters. If this function computes the negative
-                          log-likelihood, pass ``loglike_sign = -1.`` to the function as well.
-    :type loglikelihood: function or other callable
-    :param start: Starting value of the parameters. Must be a 1D array.
-    :type start: ndarray
-    :param step_sizes: Covariance matrix from which steps are drawn. Must be a 2D array
-                       with both dimensions the same length as ``current_params``, and
-                       must be positive definite.
-    :type step_sizes: ndarray
-    :param n_samples: Number of steps to be taken. Must be a positive integer. Optional,
-                      default value is 1000. Note that if the chain is thinned, the
-                      final MCMC chain will be shorter than this.
-    :type n_samples: int
-    :param thin: Integer describing how to thin the MCMC chain. A positive integer
-                 indicates manual thinning by keeping every ``thin`` steps (note
-                 that ``thin = 1`` means the chain will not be thinned). A value
-                 of ``0`` will attempt to autothin the chain.
-    :type thin: int
-    :param loglike_sign: Sign for the log-likelihood function. If the provided
-                         ``loglikelihood`` function computes the negative log-likelihood,
-                         pass ``-1.`` for this parameter. Optional, default value is ``1.``
-    :type loglike_sign: float
-    :returns: MCMC chain (2D array), array of rejected points (2D array), acceptance rate
-              (float), and first lag autocorrelation of the thinned MCMC chain (float)
-    :rtype: tuple containing (ndarray, ndarray, float, float)
-    """
-
-    n_samples = int(n_samples)
-    thin = int(thin)
-
-    assert n_samples > 0, "number of samples must be a positive integer"
-    assert thin >= 0, "thin must be a non-negative integer"
-    
-    start = np.array(start)
-    
-    assert start.ndim == 1, "starting point must be a 1d array"
-    n_params = len(start)
-
-    samples = np.zeros((n_samples, n_params))
-    samples[0] = start
-    rejected = []
-
-    for i in range(n_samples - 1):
-        next_point, accept = MCMC_step(loglikelihood, samples[i], step_sizes, loglike_sign)
-        if accept:
-            samples[i+1] = np.copy(next_point)
-        else:
-            samples[i+1] = np.copy(samples[i])
-            rejected.append(np.array(next_point))
-        
-    acceptance = float(n_samples - len(rejected))/float(n_samples)
-
-    if thin == 0:
-        thin_freq = autothin_samples(samples)
-    else:
-        thin_freq = thin
-    
-    thinned = samples[::thin_freq]
-
-    first_lag = np.zeros(n_params)
-
-    if n_samples > 1:
-        for i in range(n_params):
-            autocorr = np.correlate(thinned[:,i]-np.mean(thinned[:,i]), thinned[:,i]-np.mean(thinned[:,i]), mode="full")
-            if np.max(autocorr) > 0.:
-                first_lag[i] = autocorr[np.argmax(autocorr)+1]/np.max(autocorr)
-
-    return thinned, np.array(rejected), acceptance, first_lag
-
-def autothin_samples(signal):
-    """
-    Automatically estimate thinning needed to obtain uncorrelated samples
-    
-    This function attempts to estimate the thinning needed to obtain uncorrelated samples in an
-    MCMC chain. For each separate parameter, the function computes the autocorrelation and
-    estimates the lag needed to obtain uncorrelated samples. This is done by recognizing that
-    the standard deviation of the autocorrelation of a random signal scales inversely with the
-    square root of the number of samples, so when the autocorrelation drops below three times
-    this, we use this as a guess to when the signal is uncorrelated. The maximum lag across all
-    parameters is returned. If the chain contains fewer than 10 points, or the autocorrelation
-    never drops below the target value, the method gives a warning and
-    returns 1.
-    
-    :param signal: MCMC chain to be thinned. Must be a 1D or 2D array. If 2D, the first dimension
-                   indicates the MCMC samples, and the second dimension indicates the different
-                   parameter values. If 1D, the array is assumed to hold samples of a single
-                   parameter, and the array is reshaped to have a singleton second dimension.
-    :type signal: ndarray
-    :returns: Estimated stride needed to thin the samples to obtain uncorrelated samples. If
-              the method does not succeed, gives a warning and returns ``1``.
-    :rtype: int
-    """
-    
-    signal = np.array(signal)
-    if signal.ndim == 1:
-        signal = np.reshape(signal, (len(signal), 1))
-
-    assert signal.ndim == 2, "input signal to be thinned must be a 1d or 2d array"
-    
-    n_samples, n_params = signal.shape
-    
-    maxthin = 0
-    
-    if n_samples >= 10:
-        for i in range(n_params):
-            autocorr = np.correlate(signal[:,i]-np.mean(signal[:,i]), signal[:,i]-np.mean(signal[:,i]), mode="full")
-            start = np.argmax(autocorr)
-            if np.max(autocorr) > 0.:
-                autocorr = autocorr/np.max(autocorr)
-            for j in range(1, len(autocorr) - start):
-                if np.abs(autocorr[start + j]) < 3./np.sqrt(float(n_samples)):
-                    if j > maxthin:
-                        maxthin = j
-                    break
-    
-    if maxthin == 0:
-        warnings.warn("automatic thinning failed, posterior distribution may be multimodal")
-        maxthin = 1
-        
-    return maxthin
-    
\ No newline at end of file
diff --git a/setup.py b/setup.py
index f3039e86..9b77e53a 100644
--- a/setup.py
+++ b/setup.py
@@ -4,7 +4,7 @@
 MAJOR = 0
 MINOR = 5
 MICRO = 0
-PRERELEASE = 2
+PRERELEASE = 3
 ISRELEASED = False
 version = "{}.{}.{}".format(MAJOR, MINOR, MICRO)
 


Test patch:
diff --git a/mogp_emulator/tests/test_MCMC.py b/mogp_emulator/tests/test_MCMC.py
deleted file mode 100644
index 9b282883..00000000
--- a/mogp_emulator/tests/test_MCMC.py
+++ /dev/null
@@ -1,162 +0,0 @@
-import numpy as np
-import pytest
-from numpy.testing import assert_allclose
-from ..MCMC import MH_proposal, MCMC_step, sample_MCMC, autothin_samples
-
-def test_MH_proposal():
-    "test the Metropolis-Hastings proposal distribution"
-    
-    np.random.seed(438)
-    
-    current_param = np.zeros(2)
-    cov = np.eye(2)
-    
-    assert_allclose(MH_proposal(current_param, cov), np.array([-0.392811129690724 ,  0.4478392215027822]))
-    
-    with pytest.raises(AssertionError):
-        MH_proposal(np.zeros((2,2)), cov)
-        
-    with pytest.raises(AssertionError):
-        MH_proposal(np.zeros(2), np.ones((2, 2, 2)))
-        
-    with pytest.raises(AssertionError):
-        MH_proposal(np.zeros(2), np.ones((2, 3)))
-        
-    with pytest.raises(AssertionError):
-        MH_proposal(np.zeros(3), np.ones((2, 2)))
-        
-    with pytest.raises(AssertionError):
-        MH_proposal(np.zeros(2), -1.*np.eye(2))
-
-def test_MCMC_step():
-    "test the MCMC step routine"
-    
-    def loglikelihood(x):
-        return 1.
-        
-    current_param = np.zeros(2)
-    cov = np.eye(2)
-    
-    next_point_expected = np.array([-0.392811129690724 ,  0.4478392215027822])
-    
-    np.random.seed(438)
-    
-    next_point, accept = MCMC_step(loglikelihood, current_param, cov)
-    
-    assert_allclose(next_point, next_point_expected)
-    assert accept
-    
-    def loglikelihood(x):
-        return 100.*np.sum(x)
-
-    np.random.seed(438)
-    
-    next_point, accept = MCMC_step(loglikelihood, current_param, cov, 1.)
-    
-    assert_allclose(next_point, next_point_expected)
-    assert accept
-    
-    np.random.seed(438)
-    
-    next_point, accept = MCMC_step(loglikelihood, current_param, cov, -1.)
-    
-    assert_allclose(next_point, next_point_expected)
-    assert not accept
-    
-    with pytest.raises(AssertionError):
-        MCMC_step(loglikelihood, current_param, cov, 0.)
-        
-    def loglikelihood(x, y):
-        return -1.
-    
-    with pytest.raises(AssertionError):
-        MCMC_step(loglikelihood, current_param, cov)
-        
-    with pytest.raises(AssertionError):
-        MCMC_step(2., current_param, cov)
-        
-def test_sample_MCMC():
-    "test MCMC sampling"
-    
-    def loglikelihood(x):
-        return 100.*np.sum(x)
-    
-    current_param = np.zeros(2)
-    cov = np.eye(2)
-    
-    samples_expected = np.array([[0., 0.],
-                                [-0.392811129690724 ,  0.4478392215027822],
-                                [-1.6681438506139665,  2.697940030289671 ],
-                                [-1.6681438506139665,  2.697940030289671 ]])
-    rejected_expected = np.array([[-3.4166078360234584,  2.5824123019889043]])
-    acceptance_expected = 0.75
-    first_lag_expected = np.array([0.2886356417500504, 0.2824351401835904])
-    
-    np.random.seed(438)
-    
-    with pytest.warns(Warning):
-        samples, rejected, acceptance, first_lag = sample_MCMC(loglikelihood, current_param, cov, n_samples = 4, thin = 0)
-    
-    assert_allclose(samples, samples_expected)
-    assert_allclose(rejected, rejected_expected)
-    assert_allclose(acceptance, acceptance_expected)
-    assert_allclose(first_lag, first_lag_expected)
-    
-    np.random.seed(438)
-    
-    samples, rejected, acceptance, first_lag = sample_MCMC(loglikelihood, current_param, cov, n_samples = 4, thin = 1)
-    
-    assert_allclose(samples, samples_expected)
-    assert_allclose(rejected, rejected_expected)
-    assert_allclose(acceptance, acceptance_expected)
-    assert_allclose(first_lag, first_lag_expected)
-    
-    samples_expected = np.array([[0., 0.],
-                                [-1.6681438506139665,  2.697940030289671 ]])
-    rejected_expected = np.array([[-3.4166078360234584,  2.5824123019889043]])
-    acceptance_expected = 0.75
-    first_lag_expected = np.array([-0.5, -0.5])
-    
-    np.random.seed(438)
-    
-    samples, rejected, acceptance, first_lag = sample_MCMC(loglikelihood, current_param, cov, n_samples = 4, thin = 2)
-    
-    assert_allclose(samples, samples_expected)
-    assert_allclose(rejected, rejected_expected)
-    assert_allclose(acceptance, acceptance_expected)
-    assert_allclose(first_lag, first_lag_expected)
-    
-    with pytest.raises(AssertionError):
-        sample_MCMC(loglikelihood, current_param, cov, n_samples = -1)
-        
-    with pytest.raises(AssertionError):
-        sample_MCMC(loglikelihood, current_param, cov, 1000, -1)
-        
-    with pytest.raises(AssertionError):
-        sample_MCMC(loglikelihood, np.zeros((2,3)), cov, 1000, -1)
-
-def test_autothin_samples():
-    "test the autothinning routine"
-    
-    a = np.zeros(100)
-    a[45:55] = 1.
-    
-    assert autothin_samples(a) == 7
-    
-    a = np.zeros((100, 2))
-    a[45:55, 0] = 1.
-    a[40:60, 1] = 1.
-    
-    assert autothin_samples(a) == 11
-    
-    a = np.zeros(9)
-    
-    with pytest.warns(Warning):
-        assert autothin_samples(a) == 1
-    
-    a = np.zeros((2, 2, 2))
-    
-    with pytest.raises(AssertionError):
-        autothin_samples(a)
-    
-