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

Issue:
Feature/gpu priors
Refactoring the C++/GPU implementation to hold hyperparameters in a GPParams object, and introduce priors for hyperparameters following the same logic as the Python/CPU implementation.

* New C++ classes and corresponding Python bindings:
  - GPParams - holds values of parameters for mean function and kernel parameters (including nugget)
  - GPPriors - container class for Prior objects, including:
  - MeanPriors - not yet used in analytical meanfunc calculation, but still containing much of the necessary functionality.
  - WeakPrior - base class for other prior distributions
  - NormalPrior
  - LogNormalPrior
  - InvGammaPrior
* Default behaviour and interface of GPPriors in the C++ / GPU implementation is identical to that in the Python/CPU implementation:  if no GPPriors object is passed to the GP's constructor, default will be InvGammaPrior with appropriately chosen shape and scale for the correlation length parameters
* If a python GPPriors object is provided, python code in GaussianProcessGPU will ascertain what prior distributions and parameters need to be passed to the C++ class for it to construct the corresponding C++ versions.
* The code in `fitting.hpp` will use the `GPPriors::sample` method to obtain starting values for the fit, rather than uniform random values.


Discussion / hints:

Fix patch:
diff --git a/.travis.yml b/.travis.yml
index 080c6c18..67d9f8fa 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -1,6 +1,7 @@
 language: python
+dist: focal
 python:
-  - "3.7"
+  - "3.9"
 install:
   - pip install 'numpy>=1.19'
   - pip install -r requirements.txt
diff --git a/docs/demos/multioutput_tutorial.rst b/docs/demos/multioutput_tutorial.rst
new file mode 100644
index 00000000..8c6a63e8
--- /dev/null
+++ b/docs/demos/multioutput_tutorial.rst
@@ -0,0 +1,27 @@
+.. _multioutput_tutorial:
+
+Multi-Output Tutorial
+=====================
+
+*Note: This tutorial requires Scipy version 1.4 or later to run the simulator.*
+
+This page includes an end-to-end example of using ``mogp_emulator`` to perform model calibration
+with a simulator with multiple outputs. Note that this builds on the main tutorial with a
+second output (in this case, the velocity of the projectile at the end of the simulation),
+which is able to further constrain the NROY space as described in the first tutorial.
+
+.. literalinclude:: ../../mogp_emulator/demos/multioutput_tutorial.py
+
+One thing to note about multiple outputs is that they must be run as a script with a 
+``if __name__ == __main__`` block in order to correctly use the multiprocessing
+library. This can usually be done as in the example for short scripts, while for more
+complex analyses it is usually better to define functions (as in the benchmark for
+multiple outputs).
+
+More Details
+------------
+
+More details about these steps can be found in the :ref:`methods` section, or on the following page
+that goes into :ref:`more details <methoddetails>` on the options available in this software library.
+For more on the specific implementation detials, see the various
+:ref:`implementation pages <implementation>` describing the software components.
diff --git a/docs/implementation/DimensionReduction.rst b/docs/implementation/DimensionReduction.rst
index 7a35d730..20f82330 100644
--- a/docs/implementation/DimensionReduction.rst
+++ b/docs/implementation/DimensionReduction.rst
@@ -30,3 +30,4 @@ Utilities
 .. rubric:: References
 .. [Fukumizu1] https://www.ism.ac.jp/~fukumizu/software.html
 .. [FL13] Fukumizu, Kenji and Chenlei Leng. "Gradient-based kernel dimension reduction for regression." Journal of the American Statistical Association 109, no. 505 (2014): 359-370
+.. [LG17] Liu, Xiaoyu and Guillas, Serge. "Dimension Reduction for Gaussian Process Emulation: An Application to the Influence of Bathymetry on Tsunami Heights." SIAM/ASA Journal on Uncertainty Quantification 5, no. 1 (2017): 787-812 https://doi.org/10.1137/16M1090648
diff --git a/docs/implementation/ExperimentalDesign.rst b/docs/implementation/ExperimentalDesign.rst
index 4f7b9b01..10ed55a2 100644
--- a/docs/implementation/ExperimentalDesign.rst
+++ b/docs/implementation/ExperimentalDesign.rst
@@ -42,4 +42,19 @@ The ``LatinHypercubeDesign`` Class
     :members:
     :inherited-members:
 
+    .. automethod:: __init__
+    
+**********************************
+The ``MaxiMinLHC`` Class
+**********************************
+
+.. _MaxiMinLHC:
+
+.. automodule:: mogp_emulator.ExperimentalDesign.MaxiMinLHC
+    :noindex:
+
+.. autoclass:: mogp_emulator.ExperimentalDesign.MaxiMinLHC
+    :members:
+    :inherited-members:
+
     .. automethod:: __init__
\ No newline at end of file
diff --git a/docs/implementation/implementation.rst b/docs/implementation/implementation.rst
index 5096c814..3206e461 100644
--- a/docs/implementation/implementation.rst
+++ b/docs/implementation/implementation.rst
@@ -11,6 +11,7 @@ mogp_emulator Implementation Details
    GaussianProcessGPU
    MultiOutputGP
    fitting
+   validation
    MeanFunction
    formula
    Kernel
diff --git a/docs/implementation/validation.rst b/docs/implementation/validation.rst
new file mode 100644
index 00000000..c8e6abea
--- /dev/null
+++ b/docs/implementation/validation.rst
@@ -0,0 +1,9 @@
+.. _validation:
+
+**********************************
+The ``validation`` Module
+**********************************
+
+.. automodule:: mogp_emulator.validation
+    :members:
+    :noindex:
diff --git a/docs/index.rst b/docs/index.rst
index 762b891c..dd407949 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -32,6 +32,7 @@ details, and some included benchmarks.
              be used are:
 
    demos/gp_demos
+   demos/multioutput_tutorial
    demos/gp_kernel_demos
    demos/mice_demos
    demos/historymatch_demos
diff --git a/docs/intro/tutorial.rst b/docs/intro/tutorial.rst
index 906a5b4a..96d5662b 100644
--- a/docs/intro/tutorial.rst
+++ b/docs/intro/tutorial.rst
@@ -58,7 +58,7 @@ the two input parameters of the drag coefficient :math:`C` and the initial veloc
 and returning a single value, which is :math:`x` at the end of the simulation.
 
 .. literalinclude:: ../../mogp_emulator/demos/projectile.py
-   :lines: 1-3,12-43,46-
+   :lines: 1-80
 
 Parameter Space
 ~~~~~~~~~~~~~~~
@@ -103,7 +103,7 @@ distributions are fairly common due to their simplicity.
 To construct our experimental design and draw samples from it, we do the following:
 
 .. literalinclude:: ../../mogp_emulator/demos/tutorial.py
-   :lines: 1-4,24-28
+   :lines: 1-5,25-29
 
 This constructs an instance of :ref:`LatinHypercubeDesign <LatinHypercubeDesign>`, and
 creates the underlying distributions by providing a list of tuples. Each tuple gives the
@@ -131,11 +131,13 @@ by passing the GP object to the ``fit_GP_MAP`` function, which returns the same
 GP object but with the parameter values estimated.
 
 .. literalinclude:: ../../mogp_emulator/demos/tutorial.py
-   :lines: 33-37
+   :lines: 34-40
 
-While the function is called ``fit_GP_MAP`` (MAP means Maximum A Posteriori),
-in this case we have not provided any prior information on the parameter values,
-so it results in MLE.
+By default, if no priors are specified for the hyperparameters then defaults
+are chosen. In particular, for correlation lengths, default priors are fit
+that attempt to put most of the distribution mass in the range spanned by
+the input data. This tends to stabilize the fitting and improve performance,
+as fewer iterations are needed to ensure a good fit.
 
 Following fitting, we print out some of the hyperparameters that are estimated.
 First, we print out the correlation lengths estimated for each of the input
@@ -162,15 +164,19 @@ and the uncertainty. This is done with the ``predict`` method of
 :ref:`GaussianProcess <GaussianProcess>`:
 
 .. literalinclude:: ../../mogp_emulator/demos/tutorial.py
-   :lines: 44-52
+   :lines: 46-55
 
 ``predictions`` is an object containing the mean and uncertainty (variance)
 of the predictions. A GP assumes that the outputs follow a Normal Distribution,
 so we can perform validation by asking how many of our validation points mean estimates
-are within 2 standard deviations of the true value. Usually for this example this is
-about 8/10, so not quite as we would expect if it were perfectly recreating the
-function. However, we will see that this still is good enough in most cases
-for the task at hand.
+are within 2 standard deviations of the true value by computing the standard errors
+of the emulator predictions on the validation points. ``mogp_emulator`` contains
+a number of methods of automatically validating an emulator given some validation
+points, including computing standard errors (see the :ref:`validation <validation>`
+documentation for more details). Usually for this example we would expect
+about 8/10 to be within 2 standard devations, so not quite as we would expect if
+it were perfectly recreating the function. However, we will see that this still is
+good enough in most cases for the task at hand.
 
 History Matching
 ~~~~~~~~~~~~~~~~
@@ -206,7 +212,7 @@ and Monte Carlo sampling (especially in only 2 dimensions). Then, we create a
 Yet" (NROY). This is done as follows:
 
 .. literalinclude:: ../../mogp_emulator/demos/tutorial.py
-   :lines: 58-65
+   :lines: 60-68
 
 First, we set a large number of samples and draw them from the experimental design object. Then,
 We construct the :ref:`HistoryMatching <HistoryMatching>` object by giving the fit GP
@@ -226,7 +232,7 @@ surrogate model for reference. This plotting command is only executed if ``matpl
 installed:
 
 .. literalinclude:: ../../mogp_emulator/demos/tutorial.py
-   :lines: 5-10,69-
+   :lines: 6-11,71-
 
 which should make a plot that looks something like this:
 
diff --git a/mogp_emulator/DimensionReduction.py b/mogp_emulator/DimensionReduction.py
index b32d2c1c..8361b920 100644
--- a/mogp_emulator/DimensionReduction.py
+++ b/mogp_emulator/DimensionReduction.py
@@ -1,6 +1,7 @@
 """This module provides classes and utilities for performing dimension
-reduction.  Currently there is a single class :class:`mogp_emulator.gKDR` which implements
-the method of Fukumizu and Leng [FL13]_.
+reduction.  There is a single class :class:`mogp_emulator.gKDR` which
+implements the method of Fukumizu and Leng [FL13]_, and which can be
+used jointly with Gaussian process emulation as in [LG17]_.
 
 Example: ::
 
@@ -120,21 +121,11 @@ class gKDR(object):
 
     """Dimension reduction by the gKDR method.
 
-    See link [Fukumizu1]_ (and in particular, [FL13]_) for details of
-    the method.
-
-    Note that this is a simpler and faster method than the original
-    "KDR" method by the same authors (but with an added
-    approximation).  The KDR method will be implemented separately.
+    See [Fukumizu1]_, [FL13]_ and [LG17]_.
 
     An instance of this class is callable, with the ``__call__``
     method taking an input coordinate and mapping it to a reduced
     coordinate.
-
-    Note that this class currently implements a *direct* translation
-    of the Matlab implementation of KernelDeriv (see link above) into
-    Python/NumPy.  It is due to be replaced with a Fortran
-    implementation, but this should not affect the interface.
     """
 
     def __init__(self, X, Y, K=None, X_scale = 1.0, Y_scale = 1.0, EPS=1E-8, SGX=None, SGY=None):
@@ -378,18 +369,17 @@ def tune_parameters(cls, X, Y, train_model, cXs=None, cYs=None,
         within the kernel, minimizing the the loss from a
         Gaussian process regression:
 
-          >>> from mogp_emulator import gKDR
-          >>> from mogp_emulator import GaussianProcess
+          >>> from mogp_emulator import gKDR, GaussianProcess, fit_GP_MAP
           >>> X = ...
           >>> Y = ...
-          >>> dr, loss = gKDR.tune_parameters(X, Y, GaussianProcess.train_model)
+          >>> dr, loss = gKDR.tune_parameters(X, Y, fit_GP_MAP)
           >>> gp = GaussianProcess(dr(X), Y)
 
         Or, specifying some optional parameters for the lengthscales,
         the maximum value of `K` to use, the number of folds for
         cross-validation, and producing verbose output:
 
-          >>> dr, loss = gKDR.tune_parameters(X, Y, GaussianProcess.train_model,
+          >>> dr, loss = gKDR.tune_parameters(X, Y, fit_GP_MAP,
           ...                                 cXs = [0.5, 1.0, 2.0], cYs = [2.0],
           ...                                 maxK = 25, cross_validation_folds=4, verbose = True)
 
diff --git a/mogp_emulator/ExperimentalDesign.py b/mogp_emulator/ExperimentalDesign.py
index 1b450a71..8adc0154 100644
--- a/mogp_emulator/ExperimentalDesign.py
+++ b/mogp_emulator/ExperimentalDesign.py
@@ -1,5 +1,6 @@
 import numpy as np
 import scipy.stats
+from scipy.spatial.distance import pdist
 from inspect import signature
 
 class ExperimentalDesign(object):
@@ -235,13 +236,13 @@ def _draw_samples(self, n_samples):
         """
         raise NotImplementedError
     
-    def sample(self, n_samples):
+    def sample(self, n_samples, **kwargs):
         """
         Draw parameter samples from the experimental design
         
         This method implements drawing parameter samples from the experimental design. The method does
         this by calling the ``_draw_samples`` method to obtain samples from the :math:`[0,1]^n` hypercube,
-        where :math:`n` is the number of parameters. The ``sample``method then transforms these samples
+        where :math:`n` is the number of parameters. The ``sample`` method then transforms these samples
         drawn from the low level method to the actual parameter values using the PPF functions provided
         when initilizing the object. Note that this method also checks that all parameter values are
         finite; if any ``NaN`` values are returned, an error will be raised.
@@ -250,6 +251,9 @@ def sample(self, n_samples):
         using a different protocol only needs to change the ``_draw_samples`` method. This makes it
         simpler to define new designs, as only a single method needs to be altered.
         
+        Also accepts a ``kwargs`` argument to allow other derived classes to implement additional
+        keyword arguments.
+        
         :param n_samples: Number of samples to be drawn from the design (must be a positive integer)
         :type n_samples: int
         :returns: Samples drawn from the design parameter space as a numpy array with shape
@@ -261,16 +265,13 @@ def sample(self, n_samples):
         assert n_samples > 0, "number of samples must be positive"
         
         sample_values = np.zeros((n_samples, self.get_n_parameters()))
-        random_draws = self._draw_samples(n_samples)
+        random_draws = self._draw_samples(n_samples, **kwargs)
         
         assert np.all(random_draws >= 0.) and np.all(random_draws <= 1.), "error in generating random samples"
         
         for (dist, index) in zip(self.distributions, range(self.get_n_parameters())):
-            try:
-                sample_values[:,index] = dist(random_draws[:,index])
-            except:
-                for sample_index in range(n_samples):
-                    sample_values[sample_index, index] = dist(random_draws[sample_index,index])
+            for sample_index in range(n_samples):
+                sample_values[sample_index, index] = dist(random_draws[sample_index,index])
         
         assert np.all(np.isfinite(sample_values)), "error due to non-finite values of parameters"
         
@@ -406,7 +407,7 @@ def __init__(self, *args):
         self.method = "Monte Carlo"
         super().__init__(*args)
         
-    def _draw_samples(self, n_samples):
+    def _draw_samples(self, n_samples, **kwargs):
         """
         Low level method for drawing random samples from a Monte Carlo design
         
@@ -546,7 +547,7 @@ def __init__(self, *args):
         self.method = "Latin Hypercube"
         super().__init__(*args)
     
-    def _draw_samples(self, n_samples):
+    def _draw_samples(self, n_samples, **kwargs):
         """
         Low level method for drawing random samples from a Latin Hypercube design
         
@@ -560,7 +561,7 @@ def _draw_samples(self, n_samples):
         
         :param n_samples: Number of samples to be drawn from the design (must be a positive integer)
         :type n_samples: int
-        :returns: Random Monte Carlo samples drawn from the :math:`[0,1]^n` hypercube as a numpy
+        :returns: Random samples drawn from the :math:`[0,1]^n` hypercube as a numpy
                   array with shape ``(n_samples, n_parameters)``
         :rtype: ndarray
         """
@@ -580,4 +581,94 @@ def _draw_samples(self, n_samples):
         
         assert np.all(random_samples >= 0.) and np.all(random_samples <= 1.), "error in generating latin hypercube samples"
         
-        return random_samples
\ No newline at end of file
+        return random_samples
+        
+class MaxiMinLHC(LatinHypercubeDesign):
+    def __init__(self, *args):
+        """
+        Class representing a one-shot design of experiments with uncorrelated parameters using
+        MaxiMin Latin Hypercube Sampling
+
+        This class provides an implementation for a class for designing experiments to sample
+        the parameter space of a complex model using MaxiMin Latin Hypercube sampling. MaxiMin
+        LHCs repeatedly draw samples from the base LHC design, keeping the realization that
+        maximizes the minimum pairwise distance between all design points. Because of this,
+        MaxiMin designs tend to spread their samples closer to the edge of the parameter
+        space and in many cases result in more accurate sampling than a single LHC draw.
+
+        The parameter space can be specified in a variety of ways, but essentially the user must
+        provide a Probability Point Function (PPF, or inverse of the Cumulative Distribution Function)
+        for each input parameter. Each PPF function takes a single numeric input and maps from
+        the interval :math:`[0,1]` to the desired parameter distribution value for a given parameter,
+        and each parameter has a separate function describing its distribution. Note that this makes
+        the assumption of no correlations between any of the parameter values (a future version may
+        implement an experimental design where there are such parameter correlations). Once the
+        design is initialized, a desired number of samples can be drawn from the design, returning
+        an array holding the desired number of samples from the parameter space.
+
+        Internally, the class holds the set of PPFs for all of the parameter values, and samples are
+        drawn by calling the ``sample`` method. To draw the samples, the ``_draw_samples`` is used
+        to generate a series of points in the :math:`[0,1]^n` hypercube using MaxiMin Latin Hypercube
+        sampling, where :math:`n` is the number of paramters. This set of samples from the Latin
+        Hypercube is then mapped to the parameter space using the given PPF functions.
+
+        Unlike Monte Carlo sampling, Latin Hypercube designs attempt to sample more uniformly from the
+        parameter space. Latin Hypercube sampling ensures that each sample is drawn from a different
+        part of the space for each parameter. For example, if four samples are drawn, then for each
+        parameter, one sample is guaranteed to be drawn from each quartile of the distribution. This
+        ensures a more uniform sampling when compared on Monte Carlo sampling, but requires slightly
+        more computation to generate the samples. Note however, that for very large numbers of parameters,
+        Latin Hypercubes still may not sample very efficiently. This is due to the fact that the size of
+        the parameter space grows exponentially with the number of dimensions, so a fixed number of
+        samples will sample the space more poorly as the number of parameters increases.
+        """
+        self.method = "MaxiMinLHC"
+        super().__init__(*args)
+
+    def _draw_samples(self, n_samples, n_tries=1000, **kwargs):
+        """
+        Sampling method for MaxiMin LHCs
+
+        Iterates over multiple LHCs and return the one that maximizes the
+        minimum distance between pairs of points. The number of iterations
+        over which to look for this MaxiMin critera can be specified by
+        ``n_tries``.
+
+        Distances are computed by ``scipy.spatial.distance.pdist``. Any
+        additional ``kwargs`` passed here will be sent on to the ``pdist``
+        function.
+
+        :param n_samples: Number of samples to be drawn from the design (must be a positive integer)
+        :type n_samples: int
+        :param n_tries: Number of LHC realizations to use in maximizing
+                        the MaxiMin criteria
+        :type n_tries: int
+        :param **kwargs: Keyword arguments to be passed to the
+                         ``scipy.spatial.distance.pdist`` function.
+        :returns: Random Monte Carlo samples drawn from the :math:`[0,1]^n` hypercube as a numpy
+                  array with shape ``(n_samples, n_parameters)``
+        :rtype: ndarray
+        """
+
+        n_samples = int(n_samples)
+        assert n_samples > 0, "number of samples must be positive"
+        assert n_tries > 0, "n_tries must be a positive integer"
+        n_tries = int(n_tries)
+
+        n_parameters = self.get_n_parameters()
+
+        best_samples = np.empty((n_samples, n_parameters))
+        max_dist = -np.inf
+
+        for i in range(n_tries):
+            random_samples = super()._draw_samples(n_samples)
+            min_dist = np.min(pdist(random_samples, **kwargs))
+            if min_dist > max_dist:
+                max_dist = min_dist
+                best_samples = random_samples
+
+        assert np.all(best_samples >= 0.0) and np.all(
+            best_samples <= 1.0
+        ), "error in generating latin hypercube samples"
+
+        return best_samples
diff --git a/mogp_emulator/GaussianProcess.py b/mogp_emulator/GaussianProcess.py
index 1769a246..3b5fb6e0 100644
--- a/mogp_emulator/GaussianProcess.py
+++ b/mogp_emulator/GaussianProcess.py
@@ -119,7 +119,7 @@ def __init__(self, inputs, targets, mean=None, kernel='SquaredExponential', prio
         distributions), the user must pass a ``GPPriors`` object.
         The ``GPPriors.default_priors`` class method gives more
         control over the exact distribution choice.
-        
+
         If default priors are not used, ``priors`` must be a ``GPPriors``
         object or a list of prior distributions. If a list, it must
         have a length of ``n_params`` whose elements are either
@@ -212,7 +212,7 @@ def __init__(self, inputs, targets, mean=None, kernel='SquaredExponential', prio
             warnings.warn("The inputdict interface for mean functions has been deprecated. " +
                           "You must input your mean formulae using the x[0] format directly " +
                           "in the formula.", DeprecationWarning)
-                          
+
         if not use_patsy:
             warnings.warn("patsy is now required to parse all formulae and form design " +
                           "matrices in mogp-emulator. The use_patsy=False option will be ignored.",
@@ -224,13 +224,13 @@ def __init__(self, inputs, targets, mean=None, kernel='SquaredExponential', prio
         self.kernel = kernel
 
         _, self._nugget_type = _process_nugget(nugget)
-        
+
         self._set_priors(priors)
 
         self._theta = GPParams(n_mean=self.n_mean,
                                n_corr=self.n_corr,
                                nugget=nugget)
-                               
+
         self.Kinv = None
         self.Kinvt = None
         self.current_logpost = None
@@ -280,16 +280,16 @@ def D(self):
     @property
     def n_mean(self):
         """Returns number of mean parameters
-        
+
         :returns: Number of mean parameters
         :rtype: int
         """
         return self._dm.shape[1]
-        
+
     @property
     def n_corr(self):
         """Returns number of correlation length parameters
-        
+
         :returns: Number of correlation length parameters
         :rtype: int
         """
@@ -316,7 +316,7 @@ def n_params(self):
     @property
     def has_nugget(self):
         """Boolean indicating if the GP has a nugget parameter
-        
+
         :returns: Boolean indicating if GP has a nugget
         :rtype: bool
         """
@@ -419,7 +419,7 @@ def theta(self, theta):
         :type theta: ndarray
         :returns: None
         """
-        
+
         if theta is None:
             self._theta.set_data(None)
             self.current_logpost = None
@@ -466,7 +466,6 @@ def _set_priors(self, newpriors):
         when fitting if the new priors change the way the underlying
         parameters are found.
         """
-        
         if newpriors is None:
             self._priors = GPPriors.default_priors(self.inputs, self.n_corr, self.nugget_type)
         elif isinstance(newpriors, GPPriors):
@@ -477,7 +476,7 @@ def _set_priors(self, newpriors):
             except TypeError:
                 raise TypeError("Provided arguments for priors are not valid inputs " +
                                 "for a GPPriors object.")
-        
+
         if self._priors.n_mean > 0:
             assert self._priors.n_mean == self.n_mean
         assert self._priors.n_corr == self.n_corr, "bad number of correlation lengths in new GPPriors object"
@@ -485,7 +484,7 @@ def _set_priors(self, newpriors):
 
     def get_design_matrix(self, inputs):
         """Returns the design matrix for a set of inputs
-        
+
         For a given set of inputs, compute the design matrix based on the GP
         mean function.
         
@@ -508,7 +507,7 @@ def get_design_matrix(self, inputs):
                 raise ValueError("Provided mean function is invalid")
             if not dm.shape[0] == inputs.shape[0]:
                 raise ValueError("Provided design matrix is of the wrong shape")
-                
+
         return dm
 
     def get_cov_matrix(self, other_inputs):
@@ -556,11 +555,18 @@ def get_K_matrix(self):
 
     def _process_inputs(self, inputs):
         "Change inputs into an array and reshape if required"
-        
+
         inputs = np.array(inputs)
         if inputs.ndim == 1:
-            inputs = np.reshape(inputs, (-1, 1))
-        assert inputs.ndim == 2, "bad shape for inputs"
+            if (not hasattr(self, "_inputs") or self.D == 1):
+                inputs = np.reshape(inputs, (-1, 1))
+            else:
+                inputs = np.reshape(inputs, (1, -1))
+        assert inputs.ndim == 2, "bad shape for input"
+        if hasattr(self, "_inputs"):
+            assert (
+                inputs.shape[1] == self.D
+            ), "second dimension of other inputs must be the same as the number of input parameters"
         
         return inputs
 
@@ -686,13 +692,13 @@ def logposterior(self, theta):
         hyperparameters and log-posterior in attributes of the object.
 
         :param theta: Value of the hyperparameters. Must be array-like
-                      with shape ``(n_params,)``
+                      with shape ``(n_data,)``
         :type theta: ndarray
         :returns: negative log-posterior
         :rtype: float
 
         """
-        
+
         if self._refit(theta):
             self.fit(theta)
 
@@ -717,11 +723,11 @@ def logpost_deriv(self, theta):
         ``fit`` (and subsequently resets the cached information).
 
         :param theta: Value of the hyperparameters. Must be array-like
-                      with shape ``(n_params,)``
+                      with shape ``(n_data,)``
         :type theta: ndarray
         :returns: partial derivatives of the negative log-posterior
                   with respect to the hyperparameters (array with
-                  shape ``(n_params,)``)
+                  shape ``(n_data,)``)
         :rtype: ndarray
         """
 
@@ -805,7 +811,7 @@ def logpost_hessian(self, theta):
 
         raise NotImplementedError("Hessian computation is not currently supported")
 
-    def predict(self, testing, unc=True, deriv=False, include_nugget=True):
+    def predict(self, testing, unc=True, deriv=False, include_nugget=True, full_cov=False):
         """Make a prediction for a set of input vectors for a single set of hyperparameters
 
         Makes predictions for the emulator on a given set of input
@@ -822,13 +828,17 @@ def predict(self, testing, unc=True, deriv=False, include_nugget=True):
 
         Optionally, the emulator can also calculate the variances in
         the predictions. If the uncertainties are computed, they are
-        returned as the second output from the method as an
-        ``(n_predict,)`` shaped numpy array. Derivatives have been
-        deprecated due to changes in how the mean function is computed,
-        so setting ``deriv=True`` will have no effect and will raise
-        a ``DeprecationWarning``.
+        returned as the second output from the method, either as an
+        ``(n_predict,)`` shaped numpy array if ``full_cov=False``
+        or an array with ``(n_predict, n_predict)`` if the full
+        covariance is computed (default is only to compute the
+        variance, which is the diagonal of the full covariance).
+        
+        Derivatives have been deprecated due to changes in how the
+        mean function is computed, so setting ``deriv=True`` will
+        have no effect and will raise a ``DeprecationWarning``.
 
-        The final input to the method determines if the predictive
+        The ``include_nugget`` kwarg determines if the predictive
         variance should include the nugget or not. For situations
         where the nugget represents observational error and
         predictions are estimating the true underlying function, this
@@ -852,6 +862,11 @@ def predict(self, testing, unc=True, deriv=False, include_nugget=True):
                                predictive variance. Only relevant if
                                ``unc = True``.  Default is ``True``.
         :type include_nugget: bool
+        :param full_cov: (optional) Flag indicating if the full
+                         covariance should be computed for the
+                         uncertainty. Only relevant if ``unc = True``.
+                         Default is ``False``.
+        :type full_cov: bool
         :returns: ``PredictResult`` object holding numpy arrays with
                   the predictions and uncertainties. Predictions
                   and uncertainties have shape ``(n_predict,)``.
@@ -863,14 +878,7 @@ def predict(self, testing, unc=True, deriv=False, include_nugget=True):
         if self.theta.get_data() is None:
             raise ValueError("hyperparameters have not been fit for this Gaussian Process")
 
-        testing = np.array(testing)
-        if testing.ndim == 1:
-            if self.D == 1:
-                testing = np.reshape(testing, (-1, 1))
-            else:
-                testing = np.reshape(testing, (1, -1))
-        assert testing.ndim == 2, "test points must be a 1D or 2D array"
-        assert testing.shape[1] == self.D, "second dimension of testing must be the same as the number of input parameters"
+        testing = self._process_inputs(testing)
 
         dmtest = self.get_design_matrix(testing)
         mtest = np.dot(dmtest, self.theta.mean)
@@ -880,17 +888,32 @@ def predict(self, testing, unc=True, deriv=False, include_nugget=True):
 
         var = None
         if unc:
-            sigma_2 = self.theta.cov
-
-            if include_nugget and not self.nugget_type == "pivot":
-                sigma_2 += self.theta.nugget
 
             Kinv_Ktest = self.Kinv.solve(Ktest)
             R = calc_R(Kinv_Ktest, self._dm, dmtest)
-
-            var = np.maximum(sigma_2 - np.sum(Ktest*Kinv_Ktest, axis=0) +
-                             np.sum(R*self.Ainv.solve(R), axis=0),
-                             0.)
+            
+            if full_cov:
+                sigma_2 = self.theta.cov*self.kernel.kernel_f(
+                    testing, testing, self.theta.corr_raw
+                )
+                
+                if include_nugget and not self.nugget_type == "pivot":
+                    sigma_2 += np.eye(testing.shape[0])*self.theta.nugget
+                
+                Linv_Ktest = self.Kinv.solve_L(Ktest)
+                LAinv_R = self.Ainv.solve_L(R)
+                
+                var = (sigma_2 - np.dot(Linv_Ktest.T, Linv_Ktest) +
+                                 np.dot(LAinv_R.T, LAinv_R))
+            else:
+                sigma_2 = self.theta.cov
+
+                if include_nugget and not self.nugget_type == "pivot":
+                    sigma_2 += self.theta.nugget
+                    
+                var = np.maximum(sigma_2 - np.sum(Ktest*Kinv_Ktest, axis=0) +
+                                 np.sum(R*self.Ainv.solve(R), axis=0),
+                                 0.)
 
         inputderiv = None
         if deriv:
diff --git a/mogp_emulator/GaussianProcessGPU.py b/mogp_emulator/GaussianProcessGPU.py
index 1c1ca3f1..42071ce5 100644
--- a/mogp_emulator/GaussianProcessGPU.py
+++ b/mogp_emulator/GaussianProcessGPU.py
@@ -12,7 +12,14 @@
 import mogp_emulator.LibGPGPU as LibGPGPU
 
 from mogp_emulator.GaussianProcess import GaussianProcessBase, PredictResult
-
+from mogp_emulator.Priors import (
+    GPPriors, 
+    InvGammaPrior, 
+    GammaPrior, 
+    LogNormalPrior, 
+    WeakPrior,
+    PriorDist
+)
 
 class GPUUnavailableError(RuntimeError):
     """Exception type to use when a GPU, or the GPU library, is unavailable"""
@@ -133,8 +140,72 @@ def interpret_nugget(nugget):
     # return info needed to set the nugget on the C++ object
     return nugget_type, nugget_size
 
-class GaussianProcessGPU(GaussianProcessBase):
+def create_prior_params(**kwargs):
+    """
+    Extract the parameters needed for the C++ DenseGP_GPU instance to 
+    create its GPPriors object.
+    Can accept either an existing GPPriors object, or the parameters needed
+    to construct one, or otherwise will create a default object.
+    :param newpriors: set of priors to use
+    :type newpriors: GPPriors or dict
+    :param inputs: the inputs of the GP
+    :type inputs: numpy array
+    :param n_corr: the number of correlation length parameters
+    :type n_corr: int
+    :param nugget_type: whether nugget is adaptive, fit, or fixed
+    "type nugget_type: LibGPGPU.nugget_type
+
+    :returns: list of values needed to call the C++ constructor
+            [n_corr,
+             [(corr_i_prior_dist, [corr_i_prior_params])],
+             (cov_prior_dist, [cov_prior_params])
+             (nug_prior_dist, [nug_prior_params])]
+    :rtype: list
+    """
+    # if we are not given priors object, make default priors given the supplied args.
+    # Note that any meanfunction parameters will be given weak priors.
+    if all(x in kwargs.keys() for x in ["inputs", "n_corr", "nugget_type"]):
+        priors = GPPriors.default_priors(kwargs["inputs"], kwargs["n_corr"], kwargs["nugget_type"])
+    elif "newpriors" in kwargs.keys():
+        newpriors = kwargs["newpriors"]
+        if isinstance(newpriors, GPPriors):
+            priors = newpriors
+        else:
+            try:
+                priors = GPPriors(**newpriors)
+            except TypeError:
+                raise TypeError("Provided arguments for priors are not valid inputs " +
+                                "for a GPPriors object.")
+    else:
+        raise TypeError("Unrecognized keyword arguments for create_prior_params " +
+                        " - should be 'newpriors' or ('inputs','n_corr','nugget_type')")
+
+    def get_prior_params(prior):
+        """
+        for a WeakPrior-derived object, return a tuple (prior_type,[prior_params])
+        """
+        type_dict = {"InvGammaPrior" : LibGPGPU.prior_type.InvGamma,
+                     "GammaPrior" : LibGPGPU.prior_type.Gamma,
+                     "LogNormalPrior" : LibGPGPU.prior_type.LogNormal
+                }
+        prior_as_str = type(prior).__name__
+        if prior_as_str in type_dict.keys():
+            return (type_dict[prior_as_str], [prior.shape, prior.scale])
+        elif (isinstance(prior, WeakPrior) \
+            and not isinstance(prior, PriorDist)) \
+            or prior == None:
+            return (LibGPGPU.prior_type.Weak, [0.,0.])
+        else:
+            raise TypeError("Unknown prior type {} for C++/GPU implementation".format(type(prior)))
+
+    prior_params = [priors.n_corr]
+    prior_params.append([get_prior_params(prior) for prior in priors.corr])
+    prior_params.append(get_prior_params(priors.cov))
+    prior_params.append(get_prior_params(priors.nugget))
+    return prior_params
 
+
+class GaussianProcessGPU(GaussianProcessBase):
     """
     This class implements the same interface as
     :class:`mogp_emulator.GaussianProcess.GaussianProcess`, but using a GPU if available.
@@ -213,6 +284,8 @@ def __init__(self, inputs, targets, mean=None, kernel=SquaredExponential(), prio
         # instantiate the DenseGP_GPU class
         self._densegp_gpu = None
         self._init_gpu()   
+        # set the GPPriors
+        self._set_priors(priors)
 
     @classmethod
     def from_cpp(cls, denseGP_GPU):
@@ -242,6 +315,31 @@ def _init_gpu(self):
                                                      self._nugget_type,
                                                      self._init_nugget_size)
 
+    def _set_priors(self, newpriors=None):
+        """
+        Tell the C++ DenseGP_GPU object to create a priors object with specified parameters
+
+        :param newpriors(optional): Priors object to use (otherwise will use default)
+        :type newpriors: GPPriors object, or dict
+        """
+        if newpriors:
+            prior_params = create_prior_params(newpriors=newpriors)
+        else:
+            prior_params = create_prior_params(
+                inputs=self.inputs, n_corr=self.n_corr, nugget_type=self.nugget_type
+            )
+        self._densegp_gpu.create_gppriors(
+            *prior_params
+        )         
+
+    @property 
+    def priors(self):
+        """
+        Returns the Python binding to  the C++ GPPriors object, which holds 
+        MeanPriors, correlation_length priors, covariance prior, and nugget prior
+        """
+        return self._densegp_gpu.get_gppriors()
+
     @property
     def inputs(self):
         """
@@ -282,6 +380,17 @@ def D(self):
         """
         return self._densegp_gpu.D()
 
+
+    @property
+    def n_corr(self):
+        """
+        Returns number of correlation length parameters
+
+        :returns: Number of parameters
+        :rtype: int
+        """
+        return self._densegp_gpu.n_corr()
+
     @property
     def n_params(self):
         """
@@ -294,7 +403,8 @@ def n_params(self):
         :returns: Number of hyperparameters
         :rtype: int
         """
-        return self._densegp_gpu.n_params()
+        return self._densegp_gpu.get_theta().get_n_data() + \
+               self._densegp_gpu.get_theta().get_n_mean()
 
     @property
     def nugget_type(self):
@@ -337,8 +447,6 @@ def theta(self):
 
         :type theta: ndarray
         """
-        if not self._densegp_gpu.theta_fit_status():
-            return None
         theta = self._densegp_gpu.get_theta()
         return theta
 
@@ -350,7 +458,7 @@ def theta(self, theta):
 
         See :func:`mogp_emulator.GaussianProcess.GaussianProcess.theta`
 
-        :type theta: ndarray
+        :type theta: ndarray or GPParams object
         :returns: None
         """
         if theta is None:
@@ -370,7 +478,7 @@ def L(self):
         return np.tril(result.transpose())
 
     @property 
-    def invQt(self):
+    def Kinv_t(self):
         """
         Return the product of inverse covariance matrix with the target values     
 
@@ -389,6 +497,8 @@ def current_logpost(self):
 
         :returns: double
         """
+        if not self._densegp_gpu.theta_fit_status():
+            return None
         return self.logposterior(self.theta)
 
     def get_K_matrix(self):
@@ -444,7 +554,7 @@ def logpost_deriv(self, theta):
 
         assert theta.shape == (self.n_params,), "bad shape for new parameters"
 
-        if self.theta is None or not np.allclose(theta, self.theta, rtol=1.e-10, atol=1.e-15):
+        if self.theta is None or not np.allclose(theta, self.theta.get_data(), rtol=1.e-10, atol=1.e-15):
             self.fit(theta)
 
         result = np.zeros(self.n_params)
@@ -476,7 +586,7 @@ def predict(self, testing, unc=True, deriv=True, include_nugget=True):
         :func:`mogp_emulator.GaussianProcess.GaussianProcess.predict`
         """
 
-        if self.theta is None:
+        if not self.theta.data_has_been_set():
             raise ValueError("hyperparameters have not been fit for this Gaussian Process")
 
         testing = ndarray_coerce_type_and_flags(testing)
diff --git a/mogp_emulator/HistoryMatching.py b/mogp_emulator/HistoryMatching.py
index 66fc1e07..6260fa01 100644
--- a/mogp_emulator/HistoryMatching.py
+++ b/mogp_emulator/HistoryMatching.py
@@ -1,5 +1,6 @@
 import numpy as np
 from mogp_emulator.GaussianProcess import GaussianProcessBase, PredictResult
+from mogp_emulator.MultiOutputGP import MultiOutputGPBase
 
 class HistoryMatching(object):
     r"""
@@ -80,10 +81,11 @@ def __init__(self, gp=None, obs=None, coords=None, expectations=None, threshold=
         numbers of these to allow the number of dimensions and/or number of coordinates to
         be computed also causes these to be set at this stage.
 
-        :param gp: (Optional) ``GaussianProcessBase`` object used to make predictions in
-                   history matching. Optional, can instead provide predictions directly
+        :param gp: (Optional) ``GaussianProcessBase`` or ``MultiOutputGPBase``
+                   object used to make predictions in history matching.
+                   Optional, can instead provide predictions directly
                    via the ``expectations`` argument. Default is ``None``
-        :type gp: GaussianProcessBase or None
+        :type gp: GaussianProcessBase, MultiOutputGP, or None
         :param obs: (Optional) Observations against which the predictions will be
                     compared. If provided, must either be a float (assumes no
                     uncertainty in the observations) or a list of two floats
@@ -140,6 +142,15 @@ def __init__(self, gp=None, obs=None, coords=None, expectations=None, threshold=
 
         self.update()
 
+    def get_n_obs(self):
+        """
+        Returns the number of observations
+        
+        Determines the number of observations. This is the length
+        of the first argument of the ``obs`` attribute (which will
+        be a numpy array)
+        """
+        return len(self.obs[0])
 
     def _select_expectations(self):
         r"""
@@ -183,7 +194,7 @@ def _select_expectations(self):
 
         return expectations
 
-    def get_implausibility(self, discrepancy = 0.):
+    def get_implausibility(self, discrepancy=0., rank=1):
         r"""
         Compute Implausibility measure for all query points
 
@@ -212,14 +223,29 @@ def get_implausibility(self, discrepancy = 0.):
         matches reality. In practice, the model discrepancy is essential to have
         predictions that are not overly confident, however it can be hard to estimate
         (see further discussion in Brynjarsdottir and O\'Hagan, 2014).
+        
+        If dealing with multiple outputs, the ``rank`` argument allows the user
+        to specify how to score different predictions by giving an ordering.
+        If the value is ``0``, then the maximum implausibility score will be
+        used. If the value is ``1``, then the second largest implausibility
+        is used (i.e. the scoring ignores the largest component). By default,
+        uses the second largest score, which is in general more forgiving for
+        a badly performing emulator that is overconfident in its predictions.
+        Must be a non-negative integer.
 
         As the implausibility calculation linearly sums variances, the result is
         agnostic to the precise provenance of any of the included variances
 
         :param discrepancy: Additional variance to be included in the implausibility
-                            calculation. Must be a non-negative float. Optional, default
-                            is ``0.``
-        :type discrepancy: float
+                            calculation. Must be a non-negative array or float.
+                            Optional, default is ``0.``
+        :type discrepancy: float or ndarray
+        :param rank: Scoring method for multiple outputs. Must be a non-negative
+                    integer less than the number of observations, which denotes
+                    the location in the rank ordering of implausibility values
+                    where the score is evaluated (i.e. the default value of ``1``
+                    indicates that the second largest implausibility will be used).
+        :type rank: int
         :returns: Array holding implausibility metric for all query points accounting
                   for all variances, ndarray of length ``(ncoords,)``
         :rtype: ndarray
@@ -230,21 +256,39 @@ def get_implausibility(self, discrepancy = 0.):
             raise ValueError("implausibility calculation requires that the observation value is " +
                              "set. This can be done using the set_obs method.")
 
-        assert discrepancy >= 0., "Model discrepancy variance cannot be negative"
+        assert np.all(discrepancy >= 0.), "Model discrepancy variance cannot be negative"
+        discrepancy = np.atleast_1d(discrepancy)
 
         expectations = self._select_expectations()
 
-        # Compute implausibility for each expectation value
-        Vs = np.zeros(self.ncoords)
-        Vs += expectations[1]                         # variance on expectation
-        Vs += discrepancy                             # model discrepancy
-        Vs += self.obs[1]                             # variance on observation
-        self.I = (np.abs(self.obs[0] - expectations[0]) / np.sqrt(Vs))
+        n_obs = self.get_n_obs()
+        assert n_obs == np.atleast_2d(expectations[0]).shape[0]
+        assert n_obs == np.atleast_2d(expectations[1]).shape[0]
+        
+        if n_obs == 1:
+            rank = 0
+        
+        assert rank >= 0, "rank must be a non-negative integer"
+        assert rank < n_obs, "rank must be less than the number of observations"
 
+        # Compute implausibility for each expectation value
+        Vs = np.zeros((n_obs, self.ncoords))
+        Vs += expectations[1]                           # prediction variance
+        Vs += discrepancy[:, np.newaxis]                # model discrepancy
+        Vs += self.obs[1][:, np.newaxis]                # variance on observation
+        I = (np.abs(self.obs[0][:, np.newaxis] - expectations[0]) / np.sqrt(Vs))
+        # implausibility is the -rank index in the sorted array (i.e. if rank is 0,
+        # return the largest, if rank is 1, return the second largest, etc.),
+        # which can be correctly found in O(n) using partition.
+        # partition puts the kth largest element in its correct position along
+        # the selected axis, so we just need to select the kth row of the
+        # partitioned array along axis 0
+        self.I = np.partition(I, n_obs - rank - 1, axis=0)[n_obs - rank - 1]
+        
         return self.I
 
 
-    def get_NROY(self, discrepancy = 0.):
+    def get_NROY(self, discrepancy=0., rank=1):
         r"""
         Return set of indices that are not yet ruled out
 
@@ -263,14 +307,14 @@ def get_NROY(self, discrepancy = 0.):
         :rtype: list
         """
         if self.I is None:
-            self.get_implausibility(discrepancy)
+            self.get_implausibility(discrepancy, rank)
 
         self.NROY = list(np.where(self.I <= self.threshold)[0])
 
         return self.NROY
 
 
-    def get_RO(self, discrepancy = 0.):
+    def get_RO(self, discrepancy=0., rank=1):
         r"""
         Return set of indices that have been ruled out
 
@@ -289,7 +333,7 @@ def get_RO(self, discrepancy = 0.):
         :rtype: list
         """
         if self.I is None:
-            self.get_implausibility(discrepancy)
+            self.get_implausibility(discrepancy, rank)
 
         self.RO = list(np.where(self.I > self.threshold)[0])
 
@@ -328,13 +372,15 @@ def set_obs(self, obs):
         """
         if not self.check_obs(obs):
             raise TypeError("bad input for set_obs")
-        if isinstance(obs, list):
+        if isinstance(obs, float):
+            self.obs = [np.array([obs]), np.array([0.])]
+        elif isinstance(list(obs), list):
             if len(obs) == 1:
-                self.obs = [float(obs[0]), 0.]
+                self.obs = [np.atleast_1d(obs[0]), np.array([0.])]
             else:
-                self.obs = [float(a) for a in obs]
+                self.obs = [np.atleast_1d(a) for a in obs]
         else:
-            self.obs = [float(obs), 0.]
+            self.obs = np.array(obs)
 
 
     def set_coords(self, coords):
@@ -440,12 +486,12 @@ def check_gp(self, gp):
         Returns a boolean that is True if the provided quantity is consistent with
         the requirements for a gaussian process, i.e. is of type ``GaussianProcessBase``.
 
-        :param gp: Input GP object to be checked.
-        :type gp: GaussianProcessBase
+        :param gp: Input GP or MOGP object to be checked.
+        :type gp: GaussianProcessBase or MultiOutputGPBase
         :returns: Boolean indicating if provided object is a ``GaussianProcess``
         :rtype: bool
         """
-        return isinstance(gp, GaussianProcessBase)
+        return isinstance(gp, (GaussianProcessBase, MultiOutputGPBase))
 
 
     def check_obs(self, obs):
@@ -454,12 +500,13 @@ def check_obs(self, obs):
         observations.
 
         Returns a boolean that is ``True`` if the provided quantity is consistent with
-        the requirements for the observation quantity, i.e. is a numerical value
-        that can be converted to a float, or a list of up to two such values. Also
-        checks if the provided variance is non-negative.
+        the requirements for the observation quantity. Possible options include
+        a numpy array with the first dimension having length 2, an iterable
+        of length 1 or 2 containing floats or arrays (if 2 arrays, they must
+        have the same length), or a float for a single observation with no error.
 
         :param obs: Input for observations to be checked
-        :type obs: float or list-like
+        :type obs: float, iterable, or ndarray
         :returns: Boolean indicating if provided observations are acceptable
         :rtype: bool
         """
@@ -467,19 +514,29 @@ def check_obs(self, obs):
         if obs is None:
             return False
         if isinstance(obs, list) or isinstance(obs, tuple) or isinstance(obs, np.ndarray):
-            if len(obs) > 2:
-                raise ValueError("bad input type for HistoryMatching - the specified " +
-                                 "observation parameter cannot contain more than 2 entries "+
-                                 "(value, [variance])")
-            if len(obs) <= 2:
-                try:
-                    test = [float(a) for a in obs]
-                except TypeError:
-                    raise TypeError("bad input type for HistoryMatching - the specified " +
-                                    "observation parameter must contain only numerical " +
-                                    "values")
+            if isinstance(obs, np.ndarray):
+                if obs.ndim > 2:
+                    raise ValueError("bad input for HistoryMatching, the obs parameter" +
+                                     "cannot be an array of more than 2D")
+                assert obs.shape[0] == 2, "first dimension of observations must have length 2"
+            else:
+                if len(obs) > 2:
+                    raise ValueError("bad input type for HistoryMatching - the specified " +
+                                     "observation parameter cannot contain more than 2 entries "+
+                                     "(value, [variance])")
+                    if len(obs) <= 2:
+                        try:
+                            test = [np.array(a) for a in obs]
+                        except TypeError:
+                            raise TypeError("bad input type for HistoryMatching - the specified " +
+                                            "observation parameter must contain only numerical " +
+                                            "values")
+                        if not (len(test[0]) == len(test[1]) or len(test[1]) == 1):
+                            raise ValueError("Bad input for values to history matching -- " +
+                                             "iterable arguments must be of the format " +
+                                             "(float, float), (float, array), or (array, array)")
             if len(obs) == 2:
-                assert float(obs[1]) >= 0., "variance in observations cannot be negative"
+                assert np.all(obs[1] >= 0.), "variance in observations cannot be negative"
             return True
         else:
             try:
diff --git a/mogp_emulator/MultiOutputGP.py b/mogp_emulator/MultiOutputGP.py
index cd3deb00..8857e2cb 100644
--- a/mogp_emulator/MultiOutputGP.py
+++ b/mogp_emulator/MultiOutputGP.py
@@ -10,7 +10,13 @@
 from mogp_emulator.Priors import GPPriors
 from patsy import ModelDesc
 
-class MultiOutputGP(object):
+class MultiOutputGPBase(object):
+    """Base class for Multi-Output GPs. CPU and GPU versions derive from
+    this class.
+    """
+    pass
+
+class MultiOutputGP(MultiOutputGPBase):
     """Implementation of a multiple-output Gaussian Process Emulator.
 
     Essentially a parallelized wrapper for the predict method. To fit
@@ -35,10 +41,7 @@ def __init__(self, inputs, targets, mean=None, kernel="SquaredExponential", prio
                  nugget="adaptive", inputdict={}, use_patsy=True):
         """
         Create a new multi-output GP Emulator
-        """
-
-        self.GPClass = GaussianProcess
-            
+        """ 
         if not inputdict == {}:
             warnings.warn("The inputdict interface for mean functions has been deprecated. " +
                           "You must input your mean formulae using the x[0] format directly " +
@@ -62,9 +65,9 @@ def __init__(self, inputs, targets, mean=None, kernel="SquaredExponential", prio
         if not (inputs.shape[0] == targets.shape[1]):
             raise ValueError("the first dimension of inputs must be the same length as the second dimension of targets (or first if targets is 1D))")
 
-        self.n_emulators = targets.shape[0]
-        self.n = inputs.shape[0]
-        self.D = inputs.shape[1]
+        self._n_emulators = targets.shape[0]
+        self._n = inputs.shape[0]
+        self._D = inputs.shape[1]
 
         if not isinstance(mean, list):
             mean = self.n_emulators*[mean]
@@ -97,12 +100,87 @@ def __init__(self, inputs, targets, mean=None, kernel="SquaredExponential", prio
         assert isinstance(nugget, list), "nugget must be a string, float, or a list of strings and floats"
         assert len(nugget) == self.n_emulators
 
-        self.emulators = [ self.GPClass(inputs, single_target, m, k, p, n)
+        self.emulators = [ GaussianProcess(inputs, single_target, m, k, p, n)
                            for (single_target, m, k, p, n) in zip(targets, mean, kernel, priorslist, nugget)]
 
 
+    @property
+    def inputs(self):
+        """Full array of emulator inputs
+        
+        Returns input array (2D array of shape
+        ``(n, D)``).
+        
+        :returns: Array of input values
+        :rtype: ndarray
+        """
+        return self.emulators[0].inputs
+
+    @property
+    def targets(self):
+        """Full array of emulator targets
+        
+        Returns target array (2D array of shape
+        ``(n_emulators, n)``).
+        
+        :returns: Array of target values
+        :rtype: ndarray
+        """
+        targetlist = []
+        for em in self.emulators:
+            targetlist.append(em.targets)
+        return np.array(targetlist)
+
+    @property
+    def D(self):
+        """Number of Dimensions in inputs
+        
+        Returns number of dimentions of the input data
+        
+        :returns: Number of dimentions of inputs
+        :rtype: int
+        """
+        return self._D
+
+    @property
+    def n(self):
+        """Number of training examples in inputs
+        
+        Returns length of training data.
+        
+        :returns: Length of training inputs
+        :rtype: int
+        """
+        return self._n
+
+    @property
+    def n_params(self):
+        """Returns the number of parameters for all emulators
+        as a list of integers
+        
+        :returns: Number of parameters for each emulator
+                  as a list of integers
+        :rtype: list
+        """
+        return [em.n_params for em in self.emulators]
+
+    @property
+    def n_emulators(self):
+        return self._n_emulators
+
+    def reset_fit_status(self):
+        """Reset the fit status of all emulators
+        """
+        for em in self.emulators:
+            em.theta = None
+            
+    def _process_inputs(self, inputs):
+        "Obtain inputs that are compatible with underlying GPs"
+        
+        return self.emulators[0]._process_inputs(inputs)
+
     def predict(self, testing, unc=True, deriv=False, include_nugget=True,
-                allow_not_fit=False, processes=None):
+                full_cov = False, allow_not_fit=False, processes=None):
         """Make a prediction for a set of input vectors
 
         Makes predictions for each of the emulators on a given set of
@@ -128,6 +206,13 @@ def predict(self, testing, unc=True, deriv=False, include_nugget=True,
         uncertainties should include the nugget. By default, this is
         set to ``True``.
                 
+        If desired, the full covariance can be computed by
+        setting ``full_cov=True``. In that case, the returned
+        uncertainty will have shape
+        ``(n_emulators, n_predict, n_predict)``. This argument is
+        optional and the default is to only compute the variance,
+        not the full covariance.
+                
         Derivatives have been deprecated due to changes in how the
         mean function is computed, so setting ``deriv=True`` will
         have no effect and will raise a ``DeprecationWarning``.
@@ -159,6 +244,11 @@ def predict(self, testing, unc=True, deriv=False, include_nugget=True,
                                 predictive variance. Only relevant if
                                 ``unc = True``.  Default is ``True``.
         :type include_nugget: bool
+        :param full_cov: (optional) Flag indicating if the full
+                         predictive covariance should be computed.
+                         Only relevant if ``unc = True``.
+                         Default is ``False``.
+        :type full_cov: bool
         :param allow_not_fit: (optional) Flag that allows predictions
                               to be made even if not all emulators have
                               been fit. Default is ``False`` which
@@ -204,18 +294,18 @@ def predict(self, testing, unc=True, deriv=False, include_nugget=True,
         if allow_not_fit:
             predict_method = _gp_predict_default_NaN
         else:
-            predict_method = self.GPClass.predict
-            
+            predict_method = GaussianProcess.predict
+
         serial_predict = (platform.system() == "Windows" or
                           any([isinstance(em._mean, ModelDesc) for em in self.emulators]))
 
         if serial_predict:
-            predict_vals = [predict_method(gp, testing, unc, deriv, include_nugget)
+            predict_vals = [predict_method(gp, testing, unc, deriv, include_nugget, full_cov)
                             for gp in self.emulators]
         else:
             with Pool(processes) as p:
                 predict_vals = p.starmap(predict_method,
-                                         [(gp, testing, unc, deriv, include_nugget)
+                                         [(gp, testing, unc, deriv, include_nugget, full_cov)
                                           for gp in self.emulators])
 
         # repackage predictions into numpy arrays
@@ -238,6 +328,36 @@ def __call__(self, testing, processes=None):
 
         return self.predict(testing, unc=False, deriv=False, processes=processes)[0]
 
+    def fit(self, thetas):
+        """
+        Fit all emulators
+        
+        Fit all emulators given an 2D array of hyperparameter values
+        (if all emulators have the same number of parameters) or
+        an iterable containing 1D numpy arrays (if the emulators
+        have a variable number of hyperparameter values).
+        
+        Note that this routine does not run in parallel for the CPU
+        version.
+        
+        :param thetas: hyperparameters for all emulators. 2D array
+                       or iterable containing 1D numpy arrays, must
+                       have first dimension of size ``n_emulators``
+        :type thetas: np.array or iterable
+        """
+        for thetaval, em in zip(thetas, self.emulators):
+            em.fit(thetaval)
+
+    def fit_emulator(self, index, theta):
+        """
+        Fit a specific emulator
+        :param index: index of emulator whose hyperparameters we will set 
+        :type index: int
+        :param theta: hyperparameters for all emulators. 1D array of length n_param 
+        :type theta: np.array
+        """
+        self.emulators[index].fit(theta)
+
 
     def get_indices_fit(self):
         """Returns the indices of the emulators that have been fit
@@ -353,7 +473,7 @@ def __str__(self):
                  str(self.D)+" input variables")
 
 
-def _gp_predict_default_NaN(gp, testing, unc, deriv, include_nugget):
+def _gp_predict_default_NaN(gp, testing, unc, deriv, include_nugget, full_cov):
     """Prediction method for GPs that defaults to NaN for unfit GPs
 
     Wrapper function for the ``GaussianProcess`` predict method that
@@ -388,6 +508,11 @@ def _gp_predict_default_NaN(gp, testing, unc, deriv, include_nugget):
                            predictive variance. Only relevant if
                            ``unc = True``.  Default is ``True``.
     :type include_nugget: bool
+    :param full_cov: (optional) Flag indicating if the full
+                     predictive covariance should be computed.
+                     Only relevant if ``unc = True``.
+                     Default is ``False``.
+    :type full_cov: bool
     :returns: Tuple of numpy arrays holding the predictions,
               uncertainties, and derivatives,
               respectively. Predictions and uncertainties have
diff --git a/mogp_emulator/MultiOutputGP_GPU.py b/mogp_emulator/MultiOutputGP_GPU.py
index 7bdc3580..31b21bc2 100644
--- a/mogp_emulator/MultiOutputGP_GPU.py
+++ b/mogp_emulator/MultiOutputGP_GPU.py
@@ -5,13 +5,26 @@
     GaussianProcess,
     PredictResult
 )
-from mogp_emulator.GaussianProcessGPU import GaussianProcessGPU, parse_meanfunc_formula
+
+from mogp_emulator.GaussianProcessGPU import (
+    GaussianProcessGPU, 
+    parse_meanfunc_formula,
+    create_prior_params
+)
+from mogp_emulator.MultiOutputGP import MultiOutputGPBase
 from mogp_emulator.MeanFunction import MeanBase
 from mogp_emulator.Kernel import SquaredExponential, Matern52
+from mogp_emulator.Priors import (
+    GPPriors, 
+    WeakPrior,
+    InvGammaPrior,
+    GammaPrior,
+    LogNormalPrior
+)
 from mogp_emulator import LibGPGPU
 
 
-class MultiOutputGP_GPU(object):
+class MultiOutputGP_GPU(MultiOutputGPBase):
     """Implementation of a multiple-output Gaussian Process Emulator, using the GPU.
 
     Essentially a parallelized wrapper for the predict method. To fit
@@ -65,6 +78,7 @@ def __init__(self, inputs, targets, mean=None, kernel="SquaredExponential", prio
         elif isinstance(nugget, float):
             nugtype = LibGPGPU.nugget_type(2)
             nugsize = nugget
+            nugget = "fixed"
         else:
             raise TypeError("nugget parameter must be a string or a non-negative float")
 
@@ -98,8 +112,37 @@ def __init__(self, inputs, targets, mean=None, kernel="SquaredExponential", prio
         else:
             raise ValueError("GPU implementation requires kernel to be SquaredExponential or Matern52")
 
+        # instantiate the C++ object
         self._mogp_gpu = LibGPGPU.MultiOutputGP_GPU(inputs, targets, batch_size, meanfunc, kernel_type, nugtype, nugsize)
 
+        # deal with priors - first make sure the priors are in the form of a list of correct size
+        if isinstance(priors, (GPPriors, dict)) or priors is None:
+            priorslist = self.n_emulators*[priors]
+        else:
+            priorslist = list(priors)
+            assert isinstance(priorslist, list), ("priors must be a GPPriors object, None, or arguments to construct " +
+                                                  "a GPPriors object or a list of length n_emulators containing the above")
+
+            assert len(priorslist) == self.n_emulators, "Bad length for list provided for priors to MultiOutputGP"
+
+        self._set_priors(priorslist, nugget)
+
+
+    def _set_priors(self, priorslist, nugget_type):
+        """
+        Assign a prior to each emulator.  If an entry is None, create default prior.
+        """
+        for i in range(self.n_emulators):
+            if priorslist[i]:
+                prior_params = create_prior_params(priorslist=priorslist[i])
+            else:
+                prior_params = create_prior_params(
+                    inputs=self.inputs, n_corr=self.n_corr[i], nugget_type=nugget_type
+                )
+            self._mogp_gpu.create_priors_for_emulator(
+                i, *prior_params
+            )
+
     @property
     def inputs(self):
         return self._mogp_gpu.inputs()
@@ -116,6 +159,14 @@ def D(self):
     def n(self):
         return self._mogp_gpu.n()
 
+    @property
+    def nugget_type(self):
+        return self._mogp_gpu.get_nugget_type()
+
+    @property
+    def nugget(self):
+        return self._mogp_gpu.get_nugget_size()
+
     @property
     def n_params(self):
         return self._mogp_gpu.n_data_params()
@@ -124,6 +175,10 @@ def n_params(self):
     def n_emulators(self):
         return self._mogp_gpu.n_emulators()
 
+    @property
+    def n_corr(self):
+        return self._mogp_gpu.n_corr_params()
+
     def reset_fit_status(self):
         self._mogp_gpu.reset_fit_status()
 
@@ -226,6 +281,8 @@ def predict(self, testing, unc=True, deriv=True, include_nugget=True,
         derivs = np.zeros([self.n_emulators, n_testing, self.D])
         if unc: 
             self._mogp_gpu.predict_variance_batch(testing, means, uncs)
+            if include_nugget:
+                uncs += self.nugget
         else:
             self._mogp_gpu.predict_batch(testing, means)
         if deriv:
diff --git a/mogp_emulator/__init__.py b/mogp_emulator/__init__.py
index cfb105f0..36e682c5 100644
--- a/mogp_emulator/__init__.py
+++ b/mogp_emulator/__init__.py
@@ -10,7 +10,7 @@
 from .MultiOutputGP_GPU import MultiOutputGP_GPU
 from .fitting import fit_GP_MAP
 from .MeanFunction import MeanFunction
-from .ExperimentalDesign import MonteCarloDesign, LatinHypercubeDesign
+from .ExperimentalDesign import MonteCarloDesign, LatinHypercubeDesign, MaxiMinLHC
 from .SequentialDesign import MICEDesign
 from .HistoryMatching import HistoryMatching
 from .DimensionReduction import gKDR
diff --git a/mogp_emulator/demos/multioutput_tutorial.py b/mogp_emulator/demos/multioutput_tutorial.py
new file mode 100644
index 00000000..20eca9e3
--- /dev/null
+++ b/mogp_emulator/demos/multioutput_tutorial.py
@@ -0,0 +1,87 @@
+import numpy as np
+from projectile import simulator_multioutput, print_results
+import mogp_emulator
+
+try:
+    import matplotlib.pyplot as plt
+    makeplots = True
+except ImportError:
+    makeplots = False
+
+# An end-to-end tutorial illustrating model calibration with multiple outputs using mogp_emulator
+
+# First, we need to set up our experimental design. We would like our drag coefficient to be
+# on a logarithmic scale and initial velocity to be on a linear scale. However, our simulator
+# does the drag coefficient transformation for us, so we simply can specity the exponent on
+# a linear scale.
+
+# We will use a Latin Hypercube Design. To specify, we give the distribution that we would like
+# the parameter to take. By default, we assume a uniform distribution between two endpoints, which
+# we will use for this simulation.
+
+# Once we construct the design, can draw a specified number of samples as shown.
+
+if __name__ == "__main__": # this is required for multiprocessing to work correctly!
+
+    lhd = mogp_emulator.LatinHypercubeDesign([(-5., 1.), (0., 1000.)])
+
+    n_simulations = 50
+    simulation_points = lhd.sample(n_simulations)
+
+# Run simulator. For the multioutput simulator, returns (distance, velocity) pair
+
+    simulation_output = np.array([simulator_multioutput(p) for p in simulation_points]).T
+
+# Next, fit the surrogate MOGP model using MAP with the default priors
+# Print out hyperparameter values as correlation lengths and sigma
+
+    gp = mogp_emulator.MultiOutputGP(simulation_points, simulation_output, nugget="fit")
+    gp = mogp_emulator.fit_GP_MAP(gp, n_tries=2)
+
+    print("Correlation lengths (distance)= {}".format(gp.emulators[0].theta.corr))
+    print("Correlation lengths (velocity)= {}".format(gp.emulators[1].theta.corr))
+    print("Sigma (distance)= {}".format(np.sqrt(gp.emulators[0].theta.cov)))
+    print("Sigma (velocity)= {}".format(np.sqrt(gp.emulators[1].theta.cov)))
+    print("Nugget (distance)= {}".format(np.sqrt(gp.emulators[0].theta.nugget)))
+    print("Nugget (velocity)= {}".format(np.sqrt(gp.emulators[1].theta.nugget)))
+
+# Validate emulator by comparing to true simulated value
+# To compare with the emulator, use the predict method to get mean and variance
+# values for the emulator predictions and see how many are within 2 standard
+# deviations
+
+    n_valid = 10
+    validation_points = lhd.sample(n_valid)
+    validation_output = np.array([simulator_multioutput(p) for p in validation_points]).T
+
+    predictions = gp.predict(validation_points)
+
+    print_results(validation_points, predictions.mean)
+
+# Finally, perform history matching. Sample densely from the experimental design and
+# determine which points are consistent with the data using the GP predictions
+# We compute which points are "Not Ruled Out Yet" (NROY)
+
+# Note that our observations are now vectors, with the same ordering as the
+# simulation output
+
+    n_predict = 10000
+    prediction_points = lhd.sample(n_predict)
+
+    hm = mogp_emulator.HistoryMatching(gp=gp, coords=prediction_points, obs=[np.array([2000., 100.]),
+                                                                             np.array([100., 5.])])
+
+    nroy_points = hm.get_NROY(rank=0)
+
+    print("Ruled out {} of {} points".format(n_predict - len(nroy_points), n_predict))
+
+# If plotting enabled, visualize results
+
+    if makeplots:
+        plt.figure()
+        plt.plot(prediction_points[nroy_points,0], prediction_points[nroy_points,1], "o", label="NROY points")
+        plt.plot(simulation_points[:,0], simulation_points[:,1],"o", label="Simulation Points")
+        plt.xlabel("log Drag Coefficient")
+        plt.ylabel("Launch velocity (m/s)")
+        plt.legend()
+        plt.show()
\ No newline at end of file
diff --git a/mogp_emulator/demos/projectile.py b/mogp_emulator/demos/projectile.py
index 3798c2ca..eb3bfbbc 100644
--- a/mogp_emulator/demos/projectile.py
+++ b/mogp_emulator/demos/projectile.py
@@ -46,7 +46,7 @@ def event(t, y, c):
 
 # now can define simulator
 
-def simulator(x):
+def simulator_base(x):
     "simulator to solve ODE system for projectile motion with drag. returns distance projectile travels"
 
     # unpack values
@@ -69,21 +69,35 @@ def simulator(x):
 
     results = solve_ivp(f, (0., 1.e8), y0, events=event, args = (c,))
 
+    return results
+
+def simulator(x):
+    "simulator to solve ODE system for projectile motion with drag. returns distance projectile travels"
+
+    results = simulator_base(x)
+
     return results.y_events[0][0][2]
 
+def simulator_multioutput(x):
+    "simulator to solve ODE system with multiple outputs"
+    
+    results = simulator_base(x)
+    
+    return (results.y_events[0][0][2],
+            np.sqrt(results.y_events[0][0][0]**2 + results.y_events[0][0][1]**2))
+
 # function for printing out results
 
-def print_results(inputs, predictions):
+def print_results(inputs, errors, var):
     "convenience function for printing out results and computing mean square error"
 
-    print("Target Point                   Predicted mean            Actual Value")
-    print("------------------------------------------------------------------------------")
+    print("Target Point                Standard Error             Predictive Variance")
+    print("---------------------------------------------------------------------------------")
 
     error = 0.
 
-    for pp, m in zip(inputs, predictions):
+    for pp, e, v in zip(inputs, errors, var):
         trueval = simulator(pp)
-        print("{}      {}       {}".format(pp, m, simulator(pp)))
-        error += (trueval - m)**2
+        print("{}      {}       {}".format(pp, e, v))
 
-    print("Mean squared error: {}".format(np.sqrt(error)/len(predictions)))
\ No newline at end of file
+    print("Mean squared error: {}".format(np.sum(errors**2)/len(errors)))
\ No newline at end of file
diff --git a/mogp_emulator/demos/tutorial.py b/mogp_emulator/demos/tutorial.py
index a65313ab..e752493d 100644
--- a/mogp_emulator/demos/tutorial.py
+++ b/mogp_emulator/demos/tutorial.py
@@ -1,6 +1,7 @@
 import numpy as np
 from projectile import simulator, print_results
 import mogp_emulator
+import mogp_emulator.validation
 
 try:
     import matplotlib.pyplot as plt
@@ -37,18 +38,20 @@
 print("Sigma = {}".format(np.sqrt(gp.theta.cov)))
 print("Nugget = {}".format(np.sqrt(gp.theta.nugget)))
 
-# Validate emulator by comparing to true simulated value
-# To compare with the emulator, use the predict method to get mean and variance
-# values for the emulator predictions and see how many are within 2 standard
-# deviations
+# Validate emulator by comparing to true simulated value and compute standard
+# errors
+# Errors are sorted with the largest variance first, idx values can be used to
+# similarly sort inputs or prediction results
 
 n_valid = 10
 validation_points = lhd.sample(n_valid)
 validation_output = np.array([simulator(p) for p in validation_points])
 
-predictions = gp.predict(validation_points)
+mean, var, _ = gp.predict(validation_points)
 
-print_results(validation_points, predictions.mean)
+errors, idx = mogp_emulator.validation.standard_errors(gp, validation_points, validation_output)
+
+print_results(validation_points[idx], errors, var[idx])
 
 # Finally, perform history matching. Sample densely from the experimental design and
 # determine which points are consistent with the data using the GP predictions
@@ -69,6 +72,7 @@
     plt.figure()
     plt.plot(prediction_points[nroy_points,0], prediction_points[nroy_points,1], "o", label="NROY points")
     plt.plot(simulation_points[:,0], simulation_points[:,1],"o", label="Simulation Points")
+    plt.plot(validation_points[:,0], validation_points[:,1],"o", label="Validation Points")
     plt.xlabel("log Drag Coefficient")
     plt.ylabel("Launch velocity (m/s)")
     plt.legend()
diff --git a/mogp_emulator/fitting.py b/mogp_emulator/fitting.py
index e2eeef15..bae28c10 100644
--- a/mogp_emulator/fitting.py
+++ b/mogp_emulator/fitting.py
@@ -198,6 +198,8 @@ def _fit_single_GPGPU_MAP(gp, n_tries=15, theta0=None, method='L-BFGS-B', **kwar
     if theta0 is None or len(theta0)==0:
         theta0=np.array([])
     LibGPGPU.fit_GP_MAP(gp._densegp_gpu, n_tries, theta0)
+    if not gp.theta.data_has_been_set():
+        raise RuntimeError("Fitting did not converge")
     return gp
 
 def _fit_MOGPGPU_MAP(gp, n_tries=15, theta0=None, method='L-BFGS-B', **kwargs):
@@ -323,7 +325,7 @@ def _fit_MOGP_MAP(gp, n_tries=15, theta0=None, method='L-BFGS-B',
         emulators_to_fit = gp.get_emulators_not_fit()
         thetavals = [ theta0[idx] for idx in indices_to_fit]
 
-    if platform.system() == "Windows": 
+    if platform.system() == "Windows":
 
         fit_MOGP = [fit_GP_MAP(emulator, n_tries=n_tries, theta0=t0, method=method, **kwargs)
                     for (emulator, t0) in zip(emulators_to_fit, thetavals)]
diff --git a/mogp_emulator/linalg/cholesky.py b/mogp_emulator/linalg/cholesky.py
index ae8a2182..9d5f37e9 100644
--- a/mogp_emulator/linalg/cholesky.py
+++ b/mogp_emulator/linalg/cholesky.py
@@ -1,38 +1,88 @@
 import numpy as np
 from scipy import linalg
-from scipy.linalg import lapack, cho_solve
+from scipy.linalg import lapack, cho_solve, solve_triangular
+
 
 class ChoInv(object):
     """
     Class representing inverse of the covariance matrix
+
+    :param L: Factorized :math:`{A}` square matrix using
+              pivoting. Assumes lower triangular factorization is
+              used.
+    :type L: ndarray
     """
+
     def __init__(self, L):
         L = np.array(L)
         assert L.ndim == 2, "L must be a 2D array"
         assert L.shape[0] == L.shape[1]
         self.L = L
-        
+
     def solve(self, b):
-        if self.L.shape == (0,0):
+        """Solve a Linear System factorized using Cholesky Decomposition
+
+        Solve a system :math:`{Ax = b}` where the matrix has been factorized
+        using the `cholesky` function.
+
+        :param b: Right hand side to be solved. Can be any array that
+                  satisfies the rules of the scipy `cho_solve` routine.
+        :type b: ndarray
+        :returns: Solution to the appropriate linear system as a ndarray.
+        :rtype: ndarray
+        """
+        if self.L.shape == (0, 0):
             return np.zeros(b.shape)
         elif np.any(np.array(b.shape, dtype=int) == 0):
             return np.zeros(b.shape)
-        elif self.L.shape == (1,1):
-            assert not self.L[0,0] == 0., "L must be nonzero"
-            return b/self.L[0,0]**2
+        elif self.L.shape == (1, 1):
+            assert not self.L[0, 0] == 0.0, "L must be nonzero"
+            return b / self.L[0, 0] ** 2
         else:
             return cho_solve((self.L, True), b)
-        
+
+    def solve_L(self, b):
+        """Solve a Linear System with one component of the Cholesky Decomposition
+
+        Solve a system :math:`{Lx = b}` where the :math:`L` matrix is
+        the factorized matrix found using the `cholesky` function.
+        This is effectively a matrix sqaure root.
+
+        :param b: Right hand side to be solved. Can be any array that
+                  satisfies the rules of the scipy `solve_triangular` routine.
+        :type b: ndarray
+        :returns: Solution to the appropriate linear system as a ndarray.
+        :rtype: ndarray
+        """
+        if self.L.shape == (0, 0):
+            return np.zeros(b.shape)
+        elif np.any(np.array(b.shape, dtype=int) == 0):
+            return np.zeros(b.shape)
+        elif self.L.shape == (1, 1):
+            assert not self.L[0, 0] == 0.0, "L must be nonzero"
+            return b / self.L[0, 0]
+        else:
+            return solve_triangular(self.L, b, lower=True)
+
     def logdet(self):
-        if self.L.shape == (0,0):
-            return 0.
+        """Compute the log of the matrix determinant
+
+        Computes the log determinant of the matrix. This is simply twice
+        the sum of the log of the diagonal of the factorized matrix.
+
+        :returns: Log determinant of the factorized matrix
+        :rtype: float
+        """
+        if self.L.shape == (0, 0):
+            return 0.0
         else:
-            return 2.0*np.sum(np.log(np.diag(self.L)))
-        
+            return 2.0 * np.sum(np.log(np.diag(self.L)))
+
+
 class ChoInvPivot(ChoInv):
     """
     Class representing Pivoted Cholesky factorized matrix
-    
+
     :param L: Factorized :math:`{A}` square matrix using
               pivoting. Assumes lower triangular factorization is
               used.
@@ -43,19 +93,52 @@ class ChoInvPivot(ChoInv):
               the length of the array appearing exactly once.
     :type P: list or ndarray
     """
+
     def __init__(self, L, P):
         super().__init__(L)
-        
-        assert len(P) == self.L.shape[0], "Length of pivot matrix must match linear system"
+
+        assert (
+            len(P) == self.L.shape[0]
+        ), "Length of pivot matrix must match linear system"
 
         self.P = P
-        
+
     def solve(self, b):
         """Solve a Linear System factorized using Pivoted Cholesky Decomposition
 
-        Solve a system :math:`{Ax = b}` where the matrix has been factorized using the
-        `pivot_cholesky` function. Can also solve a system which has been
-        factorized using the regular Cholesky decomposition routine if
+        Solve a system :math:`{Ax = b}` where the matrix has been factorized
+        using the `pivot_cholesky` function. The routine rearranges the order
+        of the RHS based on the pivoting order that was used, and then
+        rearranges back to the original ordering of the RHS when returning
+        the array.
+
+        :param b: Right hand side to be solved. Can be any array that
+                  satisfies the rules of the scipy `cho_solve` routine.
+        :type b: ndarray
+        :returns: Solution to the appropriate linear system as a ndarray.
+        :rtype: ndarray
+        """
+
+        if self.L.shape == (0, 0):
+            return np.zeros(b.shape)
+        elif np.any(np.array(b.shape, dtype=int) == 0):
+            return np.zeros(b.shape)
+        elif self.L.shape == (1, 1):
+            assert not self.L[0, 0] == 0.0, "L must be nonzero"
+            return b / self.L[0, 0] ** 2
+        else:
+            try:
+                return cho_solve((self.L, True), b[self.P])[_pivot_transpose(self.P)]
+            except (IndexError, ValueError):
+                raise ValueError("Bad values for pivot matrix in pivot_cho_solve")
+
+    def solve_L(self, b):
+        """Solve a Linear System with one component of the Pivoted Cholesky
+        Decomposition
+
+        Solve a system :math:`{Lx = b}` where the matrix :math:`L` is
+        the factorized matrix found using the `pivot_cholesky` function.
+        Can also solve a system which has been factorized using the regular Cholesky decomposition routine if
         `P` is the set of integers from 0 to the length of the linear system.
         The routine rearranges the order of the RHS based on the pivoting
         order that was used, and then rearranges back to the original
@@ -68,31 +151,32 @@ def solve(self, b):
         :rtype: ndarray
         """
 
-        if self.L.shape == (0,0):
+        if self.L.shape == (0, 0):
             return np.zeros(b.shape)
         elif np.any(np.array(b.shape, dtype=int) == 0):
             return np.zeros(b.shape)
-        elif self.L.shape == (1,1):
-            assert not self.L[0,0] == 0., "L must be nonzero"
-            return b/self.L[0,0]**2
+        elif self.L.shape == (1, 1):
+            assert not self.L[0, 0] == 0.0, "L must be nonzero"
+            return b / self.L[0, 0]
         else:
             try:
-                return cho_solve((self.L, True), b[self.P])[_pivot_transpose(self.P)]
+                return solve_triangular(self.L, b[self.P], lower=True)
             except (IndexError, ValueError):
                 raise ValueError("Bad values for pivot matrix in pivot_cho_solve")
 
+
 def cholesky_factor(A, nugget, nugget_type):
     """
     Interface for Cholesky factorization
-    
+
     Calls the appropriate method given how the nugget is handled and
     returns the factorized matrix as a ChoInv class along with the
     nugget value
     """
-    
+
     if isinstance(nugget, float):
-        assert nugget >= 0., "Nugget must be non-negative in cholesky_factor"
-    
+        assert nugget >= 0.0, "Nugget must be non-negative in cholesky_factor"
+
     if nugget_type == "adaptive":
         L, nugget = jit_cholesky(A)
         Ainv = ChoInv(L)
@@ -100,14 +184,14 @@ def cholesky_factor(A, nugget, nugget_type):
         L, P = pivot_cholesky(A)
         Ainv = ChoInvPivot(L, P)
     elif nugget_type in ["fit", "fixed"]:
-        A += nugget*np.eye(A.shape[0])
+        A += nugget * np.eye(A.shape[0])
         L = fixed_cholesky(A)
         Ainv = ChoInv(L)
     else:
         raise ValueError("Bad value for nugget_type in cholesky_factor")
-        
+
     return Ainv, nugget
-        
+
 
 def _check_cholesky_inputs(A):
     """
@@ -132,20 +216,22 @@ def _check_cholesky_inputs(A):
     np.testing.assert_allclose(A.T, A)
 
     diagA = np.diag(A)
-    if np.any(diagA <= 0.):
+    if np.any(diagA <= 0.0):
         raise linalg.LinAlgError("not pd: non-positive diagonal elements")
 
     return A
 
+
 def fixed_cholesky(A):
     """
     Cholesky decomposition with fixed noise level
     """
     A = _check_cholesky_inputs(A)
-    
+
     return linalg.cholesky(A, lower=True)
 
-def jit_cholesky(A, maxtries = 5):
+
+def jit_cholesky(A, maxtries=5):
     """
     Performs Jittered Cholesky Decomposition
 
@@ -175,9 +261,9 @@ def jit_cholesky(A, maxtries = 5):
     assert int(maxtries) > 0, "maxtries must be a positive integer"
 
     A = np.ascontiguousarray(A)
-    L, info = lapack.dpotrf(A, lower = 1)
+    L, info = lapack.dpotrf(A, lower=1)
     if info == 0:
-        return L, 0.
+        return L, 0.0
     else:
         diagA = np.diag(A)
         jitter = diagA.mean() * 1e-6
@@ -194,6 +280,7 @@ def jit_cholesky(A, maxtries = 5):
 
     return L, jitter
 
+
 def pivot_cholesky(A):
     """
     Pivoted cholesky decomposition routine
@@ -219,26 +306,27 @@ def pivot_cholesky(A):
               integers indicating the pivoting order needed to produce the
               factorization.
     :rtype: tuple containing an ndarray of shape `(n,n)` of floats and a
-            ndarray of shape `(n,)` of integers. 
+            ndarray of shape `(n,)` of integers.
     """
 
     A = _check_cholesky_inputs(A)
 
     A = np.ascontiguousarray(A)
-    L, P, rank, info = lapack.dpstrf(A, lower = 1)
+    L, P, rank, info = lapack.dpstrf(A, lower=1)
     L = np.tril(L)
-    
+
     if info < 0:
         raise linalg.LinAlgError("Illegal value in covariance matrix")
-    
+
     n = A.shape[0]
 
     idx = np.arange(rank, n)
-    divs = np.cumprod(np.arange(rank+1, n+1, dtype=np.float64))
-    L[idx, idx] = L[rank-1, rank-1]/divs
-    
-    return L, P-1
-    
+    divs = np.cumprod(np.arange(rank + 1, n + 1, dtype=np.float64))
+    L[idx, idx] = L[rank - 1, rank - 1] / divs
+
+    return L, P - 1
+
+
 def _pivot_transpose(P):
     """
     Invert a pivot matrix by taking its transpose
@@ -263,9 +351,8 @@ def _pivot_transpose(P):
     assert np.array(P).ndim == 1, "pivot matrix must be a list of integers"
 
     try:
-        return np.array([np.where(P == idx)[0][0] for idx in range(len(P))], dtype=np.int32)
+        return np.array(
+            [np.where(P == idx)[0][0] for idx in range(len(P))], dtype=np.int32
+        )
     except IndexError:
         raise ValueError("Bad values for pivot matrix input to pivot_transpose")
-
-
-
diff --git a/mogp_emulator/validation.py b/mogp_emulator/validation.py
new file mode 100644
index 00000000..bc73df5b
--- /dev/null
+++ b/mogp_emulator/validation.py
@@ -0,0 +1,482 @@
+import numpy as np
+from mogp_emulator.GaussianProcess import GaussianProcessBase
+from mogp_emulator.MultiOutputGP import MultiOutputGPBase
+from mogp_emulator.linalg import cholesky_factor
+from scipy.stats import f
+
+
+def mahalanobis(gp, valid_inputs, valid_targets, scaled=False):
+    """Compute the Mahalanobis distance on a validation dataset
+
+    Given a fit GP and a set of inputs and targets for validation,
+    compute the Mahalanobis distance (the correlated equivalent
+    of the sum of the squared standard errors):
+
+    .. math::
+        M = (y_{valid} - y_{pred})^T K^{-1} (y_{valid} - y_{pred})
+
+    The Mahalanobis distance is expected to follow a scaled
+    Fisher-Snedecor distribution with ``(n_valid, n - n_mean - 2)``
+    degrees of freedom. If ``scaled=True`` is selected, then
+    the returned distance will be scaled by subtracting the
+    expected mean and dividing by the standard deviation of
+    this distribution. Note that the Fisher-Snedecor distribution
+    is not symmetric, so this cannot be interpreted in the
+    same way as standard errors, but this can nevertheless be
+    a useful heuristic. By default, the Mahalanobis distance
+    is not scaled, and a convenience function
+    ``generate_mahal_dist`` is provided to simplify comparison
+    of the Mahalanobis distance to the expected distribution.
+
+    ``gp`` must be a fit ``GaussianProcess`` or ``MultiOutputGP``
+    object, ``valid_inputs`` must be valid input data to the
+    GP/MOGP, and ``valid_targets`` must be valid target data of
+    the appropraite shape for the GP/MOGP.
+
+    :param gp: A fit ``GaussianProcess`` or ``MultiOutputGP``
+               object. If the GP/MOGP has not been fit, a
+               ``ValueError`` will be raised.
+    :type gp: ``GaussianProcess`` or ``MultiOutputGP``
+    :param valid_inputs: Input points at which the GP will
+                         be validated. Must correspond to
+                         the appropriate inputs to the
+                         provided GP.
+    :type valid_inputs: ndarray
+    :param valid_targets: Target points at which the GP will
+                          be validated. Must correspond to
+                          the appropriate target shape for
+                          the provided GP.
+    :type valid_targets: ndarray
+    :param scaled: Flag indicating if the output Mahalanobis
+                   distance should be scaled by subtracting
+                   the mean and dividing by the standard
+                   deviation of the expected Fisher-Snedecor
+                   distribution. Optional, default is
+                   ``False``.
+    :type scaled: bool
+    :returns: Mahalanobis distance computed based on the GP
+              predictions on the validation data. If a
+              multiple outputs are used, then returns a
+              numpy array of shape ``(n_emulators,)``
+              holding the Mahalanobis distance for each
+              target.
+    :rtype: ndarray
+    """
+
+    pivot_errors = pivoted_errors(gp, valid_inputs, valid_targets)
+
+    if isinstance(gp, GaussianProcessBase):
+        errors = pivot_errors[0]
+    else:
+        errors = np.array([err[0] for err in pivot_errors])
+
+    M = np.sum(errors**2, axis=-1)
+
+    if scaled:
+        expected_dists = generate_mahal_dist(gp, valid_inputs)
+
+        if isinstance(gp, GaussianProcessBase):
+            M_iter = [M]
+            expected_dists_iter = [expected_dists]
+        else:
+            M_iter = M
+            expected_dists_iter = expected_dists
+
+        M_out = []
+        for M_val, dist in zip(M_iter, expected_dists_iter):
+            mean, var = dist.stats()
+            M_out.append((M_val - mean) / np.sqrt(var))
+
+        M = np.array(M_out)
+
+        if isinstance(gp, GaussianProcessBase):
+            M = M.squeeze(axis=0)
+
+    return M
+
+
+def generate_mahal_dist(gp, valid_inputs):
+    """Generate the Expected Distribution for the Mahalanobis Distance
+
+    Convenience function for generating a ``scipy.stats.f`` object
+    appropriate for the expected Mahalanobis distribution. If a
+    ``MultiOutputGP`` object is provided, then a list of distributions
+    will be returned. In all cases, the parameters will be "frozen"
+    as appropriate for the data.
+
+    :param gp: A fit ``GaussianProcess`` or ``MultiOutputGP``
+               object.
+    :type gp: ``GaussianProcess`` or ``MultiOutputGP``
+    :param valid_inputs: Input points at which the GP will
+                         be validated. Must correspond to
+                         the appropriate inputs to the
+                         provided GP.
+    :type valid_inputs: ndarray
+    :returns: ``scipy.stats`` distribution or list of distributions.
+    :rtype: scipy.stats.rv_continuous or list
+    """
+
+    if isinstance(gp, GaussianProcessBase):
+        emulators = [gp]
+    elif isinstance(gp, MultiOutputGPBase):
+        emulators = gp.emulators
+    else:
+        raise TypeError("Provided GP is not a GaussianProcess or MultiOutputGP")
+
+    n_valid = len(gp._process_inputs(valid_inputs))
+
+    outdists = []
+    for em in emulators:
+        outdists.append(f(dfn=n_valid, dfd=em.n - em.n_mean - 2, scale=n_valid))
+
+    if len(outdists) == 1:
+        return outdists[0]
+    else:
+        return outdists
+
+
+def compute_errors(gp, valid_inputs, valid_targets, method):
+    """General pattern for computing GP validation errors
+
+    Implements the general pattern of computing errors. User
+    must provide a GP to be validated, the validation inputs,
+    and the validation targets. Additionally, a class
+    must be provided in the ``method`` argument that contains
+    the information needed to compute the ordering of the
+    errors and the errors themselves. This class must derive
+    from the ``Errors`` class and provide the following: it
+    must have a boolean class attribute ``full_cov`` that
+    determines if the full covariance or the variances are
+    needed to compute the error, and a ``__call__`` method
+    that accepts three arguments (the target values, the
+    mean predicted value, and the variance/covariance of the
+    predictions). This function must return a tuple containing
+    two numpy arrays: the first contains the error values and
+    the second containing the integer indices that indicate
+    ordering the validation errors. See the provided classes
+    ``StandardErrors`` and ``PivotErrors`` for examples.
+
+    Alternatively, this function can be called using any
+    of the following strings for the method argument
+    (all strings will be transformed to lower case):
+    ``'standard'``, ``standarderrors``, ``'pivot'``,
+    ``'pivoterrors``, or the ``StandardErrors`` or
+    ``PivotErrors`` classes.
+
+    See also the convenience functions implementing
+    standard and pivoted errors.
+
+    ``gp`` must be a fit ``GaussianProcess`` or ``MultiOutputGP``
+    object, ``valid_inputs`` must be valid input data to the
+    GP/MOGP, and ``valid_targets`` must be valid target data of
+    the appropraite shape for the GP/MOGP.
+
+    Returns a tuple (GP) or a list of tuples (MOGP). Each tuple
+    applies to a single output and contains two 1D numpy arrays.
+    The first array holds the errors, and the second holds
+    integer indices indicating the order of the errors (to
+    unscramble the inputs, index the inputs using this array of
+    integers).
+
+    :param gp: A fit ``GaussianProcess`` or ``MultiOutputGP``
+               object. If the GP/MOGP has not been fit, a
+               ``ValueError`` will be raised.
+    :type gp: ``GaussianProcess`` or ``MultiOutputGP``
+    :param valid_inputs: Input points at which the GP will
+                         be validated. Must correspond to
+                         the appropriate inputs to the
+                         provided GP.
+    :type valid_inputs: ndarray
+    :param valid_targets: Target points at which the GP will
+                          be validated. Must correspond to
+                          the appropriate target shape for
+                          the provided GP.
+    :type valid_targets: ndarray
+    :param method: Class implementing the error computation
+                   method (see above) or a string indicating
+                   the method of computing the errors.
+    :type method:
+    :returns: A tuple holding two 1D numpy arrays of length
+              ``n_valid`` or a list of such tuples. The
+              first array holds the correlated errors. The
+              second array holds the integer index values that
+              indicate the ordering of the errors. If a
+              ``GaussianProcess`` is provided, a single tuple
+              will be returned, while if a ``MultiOutputGP``
+              is provided, the return value will be a list
+              of length ``n_emulators``.
+    :rtype: tuple or list of tuples
+    """
+
+    if isinstance(method, str):
+        if method.lower in ["standard", "standarderrors"]:
+            methodclass = StandardErrors()
+        elif method.lower in ["pivot", "pivoterrors"]:
+            methodclass = PivotErrors()
+        else:
+            raise ValueError("Bad value for error method in compute_errors")
+    else:
+        methodclass = method
+
+    assert issubclass(type(methodclass), Errors), "method must be a subclass of Errors"
+
+    _check_valid_data(gp, valid_inputs, valid_targets)
+
+    pred_iter = _create_iterable_preds(
+        gp, valid_inputs, valid_targets, full_cov=methodclass.full_cov
+    )
+
+    errors = []
+
+    for target, meanval, covval in pred_iter:
+        errors.append(methodclass(target, meanval, covval))
+
+    if isinstance(gp, GaussianProcessBase):
+        return errors[0]
+    else:
+        return errors
+
+
+def standard_errors(gp, valid_inputs, valid_targets):
+    """Compute standard errors on a validation dataset
+
+    Given a fit GP and a set of inputs and targets for validation,
+    compute the standard errors (number of standard devations
+    between the true and predicted values). Numbers are left
+    signed to designate the direction of the discrepancy
+    (positive values indicate the emulator predictions are
+    larger than the true values).
+
+    The standard errors are re-ordered based on the size
+    of the predictive variance. This is done to be
+    consistent with the interface for the pivoted errors.
+    This can also be useful as a heuristic to indicate where
+    the emulator predictions are most uncertain.
+
+    ``gp`` must be a fit ``GaussianProcess`` or ``MultiOutputGP``
+    object, ``valid_inputs`` must be valid input data to the
+    GP/MOGP, and ``valid_targets`` must be valid target data of
+    the appropraite shape for the GP/MOGP.
+
+    Returns a tuple (GP) or a list of tuples (MOGP). Each tuple
+    applies to a single output and contains two 1D numpy arrays.
+    The first array holds the errors, and the second holds
+    integer indices indicating the order of the errors (to
+    unscramble the inputs, index the inputs using this array of
+    integers).
+
+    :param gp: A fit ``GaussianProcess`` or ``MultiOutputGP``
+               object. If the GP/MOGP has not been fit, a
+               ``ValueError`` will be raised.
+    :type gp: ``GaussianProcess`` or ``MultiOutputGP``
+    :param valid_inputs: Input points at which the GP will
+                         be validated. Must correspond to
+                         the appropriate inputs to the
+                         provided GP.
+    :type valid_inputs: ndarray
+    :param valid_targets: Target points at which the GP will
+                          be validated. Must correspond to
+                          the appropriate target shape for
+                          the provided GP.
+    :type valid_targets: ndarray
+    :returns: A tuple holding two 1D numpy arrays of length
+              ``n_valid`` or a list of such tuples. The
+              first array holds the correlated errors. The
+              second array holds the integer index values that
+              indicate the ordering of the errors. If a
+              ``GaussianProcess`` is provided, a single tuple
+              will be returned, while if a ``MultiOutputGP``
+              is provided, the return value will be a list
+              of length ``n_emulators``.
+    :rtype: tuple or list of tuples
+    """
+    return compute_errors(gp, valid_inputs, valid_targets, method=StandardErrors())
+
+
+def pivoted_errors(gp, valid_inputs, valid_targets):
+    """Compute correlated errors on a validation dataset
+
+    Given a fit GP and a set of inputs and targets for validation,
+    compute the correlated errors (number of standard devations
+    between the true and predicted values, conditional on the
+    errors in decreasing order). Note that because the errors are
+    conditional, order matters and thus the errors are treated
+    with respect to the largest one. The routine returns both
+    the correlated errors and the index ordering of the validation
+    points (if a ``GaussianProcess`` is provided) or a list of
+    tuples containing the errors and indices indicating the
+    ordering of the errors for each target (if a ``MultiOutputGP``
+    is provided).
+
+    ``gp`` must be a fit ``GaussianProcess`` or ``MultiOutputGP``
+    object, ``valid_inputs`` must be valid input data to the
+    GP/MOGP, and ``valid_targets`` must be valid target data of
+    the appropraite shape for the GP/MOGP.
+
+    Returns a tuple (GP) or a list of tuples (MOGP). Each tuple
+    applies to a single output and contains two 1D numpy arrays.
+    The first array holds the errors, and the second holds
+    integer indices indicating the order of the errors (to
+    unscramble the inputs, index the inputs using this array of
+    integers).
+
+    :param gp: A fit ``GaussianProcess`` or ``MultiOutputGP``
+               object. If the GP/MOGP has not been fit, a
+               ``ValueError`` will be raised.
+    :type gp: ``GaussianProcess`` or ``MultiOutputGP``
+    :param valid_inputs: Input points at which the GP will
+                         be validated. Must correspond to
+                         the appropriate inputs to the
+                         provided GP.
+    :type valid_inputs: ndarray
+    :param valid_targets: Target points at which the GP will
+                          be validated. Must correspond to
+                          the appropriate target shape for
+                          the provided GP.
+    :type valid_targets: ndarray
+    :returns: Tuples holding two 1D numpy arrays of length
+              ``n_valid`` or a list of such tuples. The
+              first array holds the correlated errors. The
+              second array holds the integer index values that
+              indicate the ordering of the errors. If a
+              ``GaussianProcess`` is provided, a single tuple
+              will be returned, while if a ``MultiOutputGP``
+              is provided, the return value will be a list
+              of length ``n_emulators``.
+    :rtype: tuple or list of tuples
+    """
+
+    return compute_errors(gp, valid_inputs, valid_targets, method=PivotErrors())
+
+
+class Errors(object):
+    """
+    Base class implementing a method for computing errors
+    """
+
+    full_cov = False
+
+    def __call__(self, target, mean, cov):
+        "compute ordering array and errors"
+        raise NotImplementedError
+
+
+class StandardErrors(Errors):
+    """
+    Class implementing standard errors
+
+    This class implements the required functionality for computing
+    standard errors. This includes setting the class attribute
+    ``full_cov=False`` and implementing the ``__call__`` method
+    to compute the standard errors and their ordering given target
+    values and predicted mean/variance
+    """
+
+    full_cov = False
+
+    def __call__(self, target, mean, cov):
+        """
+        Compute standard errors and ordering array
+
+        Returns the standard errors (in decreasing order) and the
+        index ordering of the errors.
+
+        :param target: Validation target values as a 1D numpy array
+        :type target: ndarray
+        :param mean: Predicted mean values as a 1D numpy array
+        :type mean: ndarray
+        :param cov: Predicted variance values as a 1D numpy array
+        :type cov: ndarray
+        :returns: Tuple containing two 1D numpy arrays. The first
+                  is the standard errors sorted in descending
+                  order, the second is the integer indices
+                  indicating this ordering.
+        :rtype: tuple
+        """
+
+        P = np.argsort(cov)[::-1]
+        error = ((mean - target) / np.sqrt(cov))[P]
+
+        return error, P
+
+
+class PivotErrors(Errors):
+    """
+    Class implementing pivoted errors
+
+    This class implements the required functionality for computing
+    pivoted errors. This includes setting the class attribute
+    ``full_cov=True`` and implementing the ``__call__`` method
+    to compute the pivoted errors and their ordering given target
+    values and predicted mean/variance
+    """
+
+    full_cov = True
+
+    def __call__(self, target, mean, cov):
+        """
+        Compute correlated (pivoted) errors and ordering array
+
+        Returns the correlated pivoted errors (sorted in order
+        of decreasing variance conditional on all previous
+        errors) and the indices indicating this ordering.
+
+        :param target: Validation target values as a 1D numpy array
+        :type target: ndarray
+        :param mean: Predicted mean values as a 1D numpy array
+        :type mean: ndarray
+        :param cov: Predicted covariance values as a 2D numpy array
+        :type cov: ndarray
+        :returns: Tuple containing two 1D numpy arrays. The first
+                  is the correlated errors sorted in descending
+                  order conditional on the previous errors, the
+                  second is the integer indices indicating this
+                  ordering.
+        :rtype: tuple
+        """
+
+        cov_inv, _ = cholesky_factor(cov, 0.0, "pivot")
+        error = cov_inv.solve_L(mean - target)
+
+        return error, cov_inv.P
+
+
+def _create_iterable_preds(gp, valid_inputs, valid_targets, full_cov=False):
+    "Ensure that predictions and targets can be iterated over"
+
+    _check_valid_data(gp, valid_inputs, valid_targets)
+
+    mean, cov, _ = gp.predict(valid_inputs, full_cov=full_cov)
+
+    if isinstance(gp, GaussianProcessBase):
+        valid_targets_iter = [valid_targets]
+        mean_iter = [mean]
+        cov_iter = [cov]
+    else:
+        valid_targets_iter = valid_targets
+        mean_iter = mean
+        cov_iter = cov
+
+    return zip(valid_targets_iter, mean_iter, cov_iter)
+
+
+def _check_valid_data(gp, valid_inputs, valid_targets):
+    "Perform some checks on the validation data"
+
+    assert isinstance(
+        gp, (GaussianProcessBase, MultiOutputGPBase)
+    ), "Must provide a GP to validate"
+
+    valid_inputs = gp._process_inputs(valid_inputs)
+    valid_targets = np.array(valid_targets)
+
+    if isinstance(gp, GaussianProcessBase):
+        assert valid_targets.ndim == 1, "Targets for a GP must be a 1D array"
+        assert (
+            valid_targets.shape[0] == valid_inputs.shape[0]
+        ), "Bad length for validation targets"
+    else:
+        assert valid_targets.ndim == 2, "Targets for a MultiOutputGP must be a 2D array"
+        assert (
+            valid_targets.shape[1] == valid_inputs.shape[0]
+        ), "Bad shape for validation targets"
diff --git a/mogp_gpu/Makefile b/mogp_gpu/Makefile
index 9ec42e89..9f09adec 100644
--- a/mogp_gpu/Makefile
+++ b/mogp_gpu/Makefile
@@ -74,4 +74,13 @@ bin/test/test_mogp_gpu: test/test_mogp_gpu.cu obj/util.o obj/kernel.o
 	mkdir -p bin/test
 	nvcc $(CXXFLAGS) $(NVCCFLAGS) -lcublas -lcusolver  $<  obj/kernel.o obj/util.o -o $@
 
+
+bin/test/test_gpparams: test/test_gpparams.cu 
+	mkdir -p bin/test
+	nvcc $(CXXFLAGS) $(NVCCFLAGS) $<  -o $@
+
+bin/test/test_gppriors: test/test_gppriors.cu 
+	mkdir -p bin/test
+	nvcc $(CXXFLAGS) $(NVCCFLAGS) $<  -o $@
+
 .PHONY: all clean
diff --git a/mogp_gpu/src/bindings.cu b/mogp_gpu/src/bindings.cu
index 942269eb..e5244a0f 100644
--- a/mogp_gpu/src/bindings.cu
+++ b/mogp_gpu/src/bindings.cu
@@ -1,4 +1,6 @@
 #include "fitting.hpp"
+#include "gpparams.hpp"
+#include "gppriors.hpp"
 #include "multioutputgp_gpu.hpp"
 
 #include "pybind11/pybind11.h"
@@ -21,6 +23,10 @@ PYBIND11_MODULE(libgpgpu, m) {
         .def("D", &DenseGP_GPU::get_D,
              "Returns the dimension of the emulator inputs")
 
+////////////////////////////////////////
+        .def("n_corr", &DenseGP_GPU::get_n_corr,
+             "Returns the number of correlation length parameters")
+
 ////////////////////////////////////////
         .def("inputs", &DenseGP_GPU::get_inputs,
              "Returns the array of inputs points used by the emulator, with "
@@ -49,6 +55,18 @@ PYBIND11_MODULE(libgpgpu, m) {
         .def("get_theta", &DenseGP_GPU::get_theta,
              "Returns the emulator hyperparameters")
 
+////////////////////////////////////////
+        .def("get_gppriors", &DenseGP_GPU::get_gppriors,
+             py::return_value_policy::reference,
+             "Returns the GPPriors object")
+////////////////////////////////////////
+        .def("set_gppriors", &DenseGP_GPU::set_gppriors,
+             "Set the GPPriors object")
+
+///////////////////////////////////////
+        .def("create_gppriors", &DenseGP_GPU::create_gppriors,  
+             "Instantiate and configure the GPPriors object")
+
 ////////////////////////////////////////
         .def("predict", &DenseGP_GPU::predict,
              R"(Single point predictive mean.
@@ -126,10 +144,20 @@ Currently unused by :class`GaussianProcessGPU`, but useful for testing.
 
 ////////////////////////////////////////
 
-        .def("fit", &DenseGP_GPU::fit,
+        .def("fit", py::overload_cast<GPParams&>(&DenseGP_GPU::fit),
+             R"(Sets theta, and updates the internal state of the emulator accordingly.
+
+:param theta: Values for the hyperparameters (length ``n_data``)
+
+:returns: None
+)",
+             py::arg("theta"))
+////////////////////////////////////////
+
+.def("fit", py::overload_cast<vec>(&DenseGP_GPU::fit),
              R"(Sets theta, and updates the internal state of the emulator accordingly.
 
-:param theta: Values for the hyperparameters (length ``n_params``)
+:param theta: Values for the hyperparameters (length ``n_data``)
 
 :returns: None
 )",
@@ -166,7 +194,11 @@ Currently unused by :class`GaussianProcessGPU`, but useful for testing.
              py::arg("invQt_h"))
 
 ////////////////////////////////////////
-        .def("get_logpost", &DenseGP_GPU::get_logpost,
+        .def("get_logpost", py::overload_cast<vec>(&DenseGP_GPU::get_logpost),
+	     "Return the log posterior"
+            )
+////////////////////////////////////////
+        .def("get_logpost", py::overload_cast<GPParams&>(&DenseGP_GPU::get_logpost),
 	     "Return the log posterior"
             )
 
@@ -248,14 +280,22 @@ likelihood of the hyperparameters, from the current state of the emulator.)
                      "Return the number of inputs to the GP")
        .def("D", &MultiOutputGP_GPU::get_D,
                      "Return the number of dimensions of the GP")
-       .def("n_data_params", &MultiOutputGP_GPU::get_n_data_params,
-                     "Return the number of fitting parameters of the GP")
+       .def("n_data_params", &MultiOutputGP_GPU::n_data_params,
+                     "Return the number of data parameters of the GP")
+       .def("n_corr_params", &MultiOutputGP_GPU::n_corr_params,
+                     "Return the number of correlation length parameters of the GP")
+       .def("get_nugget_type", &MultiOutputGP_GPU::get_nugget_type,
+                     "Return the nugget type")
+       .def("get_nugget_size", &MultiOutputGP_GPU::get_nugget_size,
+                     "Return the nugget size")
        .def("get_fitted_indices", &MultiOutputGP_GPU::get_fitted_indices,
                      "Return a vector of indices of emulators where fitting succeeded")
        .def("get_unfitted_indices", &MultiOutputGP_GPU::get_unfitted_indices,
                      "Return a vector of indices of emulators where fitting failed or not done yet")
        .def("n_emulators", &MultiOutputGP_GPU::n_emulators,
                      "Return the number of GP emulators")
+       .def("create_priors_for_emulator", &MultiOutputGP_GPU::create_priors_for_emulator,
+                     "Instantiate GPPriors object for specified emulator index")
        .def("reset_fit_status", &MultiOutputGP_GPU::reset_fit_status,
                      "reset the fit status of all emulators")
        .def("targets", &MultiOutputGP_GPU::get_targets,
@@ -263,11 +303,18 @@ likelihood of the hyperparameters, from the current state of the emulator.)
        .def("targets_at_index", &MultiOutputGP_GPU::get_targets_at_index,
                          "Return the targets of the GP at specified index",
                          py::arg("index"))
-       .def("fit_emulator", &MultiOutputGP_GPU::fit_emulator,
+       .def("fit_emulator", py::overload_cast<unsigned int, GPParams&>(&MultiOutputGP_GPU::fit_emulator),
+                         "Set hyperparameters of the GP at specified index",
+                         py::arg("index"),
+                         py::arg("theta"))  
+       .def("fit_emulator", py::overload_cast<unsigned int, vec&>(&MultiOutputGP_GPU::fit_emulator),
                          "Set hyperparameters of the GP at specified index",
                          py::arg("index"),
                          py::arg("theta"))  
-       .def("fit", &MultiOutputGP_GPU::fit,
+       .def("fit", py::overload_cast<mat_ref>(&MultiOutputGP_GPU::fit),
+                         "Set hyperparameters of all emulators",
+                         py::arg("thetas"))  
+       .def("fit", py::overload_cast<std::vector<GPParams>>(&MultiOutputGP_GPU::fit),
                          "Set hyperparameters of all emulators",
                          py::arg("thetas"))  
        .def("predict", &MultiOutputGP_GPU::predict,
@@ -297,7 +344,9 @@ likelihood of the hyperparameters, from the current state of the emulator.)
         .def("kernel_deriv", &SquaredExponentialKernel::kernel_deriv,
              "Calculate the derivative of the covariance matrix wrt hyperparameters")
         .def("kernel_inputderiv", &SquaredExponentialKernel::kernel_inputderiv,
-	     "Derivative of covariance matrix wrt inputs");
+	     "Derivative of covariance matrix wrt inputs")
+        .def("get_n_params", &SquaredExponentialKernel::get_n_params,
+	     "Number of correlation length parameters");
 
 ////////////////////////////////////////
     py::class_<Matern52Kernel>(m, "Matern52Kernel")
@@ -307,7 +356,10 @@ likelihood of the hyperparameters, from the current state of the emulator.)
         .def("kernel_deriv", &Matern52Kernel::kernel_deriv,
              "Calculate the derivative of the covariance matrix wrt hyperparameters")
         .def("kernel_inputderiv", &Matern52Kernel::kernel_inputderiv,
-	     "Derivative of covariance matrix wrt inputs");
+	     "Derivative of covariance matrix wrt inputs")
+        .def("get_n_params", &SquaredExponentialKernel::get_n_params,
+	     "Number of correlation length parameters");
+;
 
 ////////////////////////////////////////
     py::class_<BaseMeanFunc>(m, "BaseMeanFunc");
@@ -360,23 +412,199 @@ likelihood of the hyperparameters, from the current state of the emulator.)
         .def("get_n_params", &ZeroMeanFunc::get_n_params,
 	     "Number of parameters of mean function");
 
-    py::enum_<kernel_type>(m, "kernel_type")
+     py::class_<GPParams>(m, "GPParameters")
+      .def(py::init<int, int, nugget_type>())
+        .def("get_data", &GPParams::get_data,
+             "get the raw data")
+        .def("set_data", &GPParams::set_data,
+             "set the raw data")
+        .def("get_n_data", &GPParams::get_n_data,
+             "get size of data array")
+        .def("get_mean", &GPParams::get_mean,
+	     "get the mean function parameters")
+        .def("set_mean", &GPParams::set_mean,
+	     "set the mean function parameters")
+        .def("get_n_mean", &GPParams::get_n_mean,
+	     "get the number of mean function parameters")
+        .def("get_nugget_type", &GPParams::get_nugget_type,
+	     "get the nugget_type")
+        .def("set_nugget_type", &GPParams::set_nugget_type,
+	     "set nugget_type")
+        .def("get_nugget_size", &GPParams::get_nugget_size,
+	     "get the nugget_size")
+        .def("set_nugget_size", &GPParams::set_nugget_size,
+	     "set nugget_size")
+        .def("get_corr", &GPParams::get_corr,
+	     "get the transformed correlation parameters")
+        .def("set_corr", &GPParams::set_corr,
+	     "set correlation parameters using transformed values")
+        .def("get_n_corr", &GPParams::get_n_corr,
+	     "get the number of correlation parameters")
+        .def("get_cov", &GPParams::get_cov,
+	     "get the transformed covariance parameter")
+        .def("set_cov", &GPParams::set_cov,
+	     "set covariance parameter using transformed values")
+        .def("get_corr_raw", &GPParams::get_corr_raw,
+	     "get the raw correlation parameters")
+        .def("unset_data", &GPParams::unset_data,
+	     "unset the flag to say params have been set")
+        .def("data_has_been_set", &GPParams::data_has_been_set,
+	     "get the flag to say whether or not params have been set")
+        .def("test_same_shape", py::overload_cast<GPParams&>(&GPParams::test_same_shape, py::const_),
+	     "test whether two GPParams objects are the same shape")
+        .def("test_same_shape", py::overload_cast<vec&>(&GPParams::test_same_shape, py::const_),
+	     "test whether two vectors of params are the same shape");
+     ////////////////////////////////////////
+     py::class_<WeakPrior>(m, "WeakPrior")
+      .def(py::init<>())
+          .def("sample", py::overload_cast<>(&WeakPrior::sample),
+          "sample from the prior")
+          .def("sample", py::overload_cast<const BaseTransform&>(&WeakPrior::sample),
+          "sample from the prior");
+     ////////////////////////////////////////
+     py::class_<InvGammaPrior, WeakPrior>(m, "InvGammaPrior")
+      .def(py::init<double, double>())
+        .def("logp", &InvGammaPrior::logp,
+             "Calculate the log probability")
+        .def("dlogpdx", &InvGammaPrior::dlogpdx,
+             "Calculate the derivative of the log probability")
+        .def("d2logpdx2", &InvGammaPrior::d2logpdx2,
+	     "Second Derivative of log prob wrt input")
+        .def("sample_x", &InvGammaPrior::sample_x,
+	     "Sample from the distribution");
+     //////////////////////////////////////
+     py::class_<GammaPrior, WeakPrior>(m, "GammaPrior")
+          .def(py::init<double, double>())
+            .def("logp", &GammaPrior::logp,
+                 "Calculate the log probability")
+            .def("dlogpdx", &GammaPrior::dlogpdx,
+                 "Calculate the derivative of the log probability")
+            .def("d2logpdx2", &GammaPrior::d2logpdx2,
+              "Second Derivative of log prob wrt input")
+            .def("sample_x", &GammaPrior::sample_x,
+              "Sample from the distribution");
+     //////////////////////////////////////
+     py::class_<BaseTransform>(m, "BaseTransform");
+     ////////////////////////////////////////
+     py::class_<CovTransform, BaseTransform>(m, "CovTransform")
+          .def(py::init<>()) 
+          .def("raw_to_scaled", py::overload_cast<REAL>(&CovTransform::raw_to_scaled, py::const_),
+          "transform from raw values to scaled ones")
+          .def("raw_to_scaled", py::overload_cast<vec>(&CovTransform::raw_to_scaled, py::const_),
+          "transform from raw values to scaled ones")
+          .def("scaled_to_raw", py::overload_cast<REAL>(&CovTransform::scaled_to_raw, py::const_),
+          "transform from scaled values to raw ones")
+          .def("scaled_to_raw", py::overload_cast<vec>(&CovTransform::scaled_to_raw, py::const_),
+          "transform from scaled values to raw ones")
+          .def("dscaled_draw", py::overload_cast<REAL>(&CovTransform::dscaled_draw, py::const_),
+          "derivative of scaled values wrt raw ones")
+          .def("dscaled_draw", py::overload_cast<vec>(&CovTransform::dscaled_draw, py::const_),
+          "derivative of scaled values wrt raw ones")
+          .def("2scaled_draw2", py::overload_cast<REAL>(&CovTransform::d2scaled_draw2, py::const_),
+          "second derivative of scaled values wrt raw ones")
+          .def("d2scaled_draw2", py::overload_cast<vec>(&CovTransform::d2scaled_draw2, py::const_),
+          "second derivative of scaled values wrt raw ones");
+        
+     ////////////////////////////////////////
+     py::class_<CorrTransform, BaseTransform>(m, "CorrTransform")
+          .def(py::init<>()) 
+          .def("raw_to_scaled", py::overload_cast<REAL>(&CorrTransform::raw_to_scaled, py::const_),
+          "transform from raw values to scaled ones")
+          .def("raw_to_scaled", py::overload_cast<vec>(&CorrTransform::raw_to_scaled, py::const_),
+          "transform from raw values to scaled ones")
+          .def("scaled_to_raw", py::overload_cast<REAL>(&CorrTransform::scaled_to_raw, py::const_),
+          "transform from scaled values to raw ones")
+          .def("scaled_to_raw", py::overload_cast<vec>(&CorrTransform::scaled_to_raw, py::const_),
+          "transform from scaled values to raw ones")
+          .def("dscaled_draw", py::overload_cast<REAL>(&CorrTransform::dscaled_draw, py::const_),
+          "derivative of scaled values wrt raw ones")
+          .def("dscaled_draw", py::overload_cast<vec>(&CorrTransform::dscaled_draw, py::const_),
+          "derivative of scaled values wrt raw ones")
+          .def("2scaled_draw2", py::overload_cast<REAL>(&CorrTransform::d2scaled_draw2, py::const_),
+          "second derivative of scaled values wrt raw ones")
+          .def("d2scaled_draw2", py::overload_cast<vec>(&CorrTransform::d2scaled_draw2, py::const_),
+          "second derivative of scaled values wrt raw ones");
+     ////////////////////////////////////////
+     py::class_<GPPriors>(m, "GPPriors")
+          .def(py::init<int, nugget_type>())
+          .def("set_corr", py::overload_cast<std::vector<WeakPrior*>>(&GPPriors::set_corr),
+               "set correlation length priors")
+          .def("set_corr", py::overload_cast<>(&GPPriors::set_corr),
+               "set correlation length priors")
+          .def("get_corr", &GPPriors::get_corr,
+               "get correlation length priors")
+          .def("set_cov", py::overload_cast<>(&GPPriors::set_cov),
+               "set covariance priors")
+          .def("set_cov", py::overload_cast<WeakPrior*>(&GPPriors::set_cov),
+               "set covariance priors")
+          .def("get_cov", &GPPriors::get_cov,
+               "get covariance priors")
+          .def("get_logp", &GPPriors::logp,
+               "get log posterior ")
+          .def("get_dlogpdtheta", &GPPriors::dlogpdtheta,
+               "get derivative of log posterior wrt theta")
+          .def("get_d2logpdtheta2", &GPPriors::d2logpdtheta2,
+               "get 2nd derivative of log posterior wrt theta")
+          .def("create_corr_priors", &GPPriors::create_corr_priors,
+               "create default correlation length priors")
+          .def("create_cov_prior", &GPPriors::create_cov_prior,
+               "create default covariance parameter prior")
+          .def("make_prior", &GPPriors::make_prior,
+               "instantiate prior of chosen type")
+          .def("sample", &GPPriors::sample,
+               "sample from the priors");
+
+     ////////////////////////////////////////
+     py::class_<MeanPriors>(m, "MeanPriors")
+          .def(py::init<vec, mat>())
+          .def("get_mean", &MeanPriors::get_mean,
+               "return the meanfunction vector")
+          .def("get_cov", &MeanPriors::get_cov,
+               "return the covariance matrix")
+          .def("get_n_params", &MeanPriors::get_n_params,
+               "how many parameters?")
+          .def("has_weak_priors", &MeanPriors::has_weak_priors,
+               "do we only have non-informative priors?")
+          .def("dm_dot_b", &MeanPriors::dm_dot_b,
+               "dot product of design matrix with mean")
+          .def("inv_cov", &MeanPriors::get_inv_cov,
+               "inverse covariance matrix")
+          .def("inv_cov_b", &MeanPriors::get_inv_cov_b,
+               "dot product of inverse covariance matrix with mean")
+          .def("logdet_cov", &MeanPriors::logdet_cov,
+               "log of the determinant of the covariance")
+          .def("set_prior_dists", py::overload_cast<>(&MeanPriors::set_prior_dists),
+               "specify default distribution for meanfunc priors")
+          .def("set_prior_dists", py::overload_cast<std::vector<prior_type>, 
+                                  std::vector<std::vector<REAL>>>(&MeanPriors::set_prior_dists),
+               "specify the distribution for meanfunc priors")
+          .def("sample", &MeanPriors::sample,
+               "sample from the meanfunc prior distributions");
+
+     ////////////////  enums   /////////////////////
+     py::enum_<kernel_type>(m, "kernel_type")
         .value("SquaredExponential", SQUARED_EXPONENTIAL)
         .value("Matern52", MATERN52);
 
-    py::enum_<nugget_type>(m, "nugget_type")
+     py::enum_<nugget_type>(m, "nugget_type")
         .value("adaptive", NUG_ADAPTIVE)
         .value("fixed", NUG_FIXED)
         .value("fit", NUG_FIT);
 
-    m.def("have_compatible_device", &have_compatible_device);
+     py::enum_<prior_type>(m, "prior_type")
+        .value("LogNormal", LOGNORMAL)
+        .value("Gamma", GAMMA)
+        .value("InvGamma", INVGAMMA)
+        .value("Weak", WEAK);
+
+     m.def("have_compatible_device", &have_compatible_device);
 
-    m.def("fit_GP_MAP", py::overload_cast<DenseGP_GPU&, const int, const std::vector<double>>(&fit_single_GP_MAP), 
+     m.def("fit_GP_MAP", py::overload_cast<DenseGP_GPU&, const int, const std::vector<double>>(&fit_single_GP_MAP), 
           py::return_value_policy::reference);
-    m.def("fit_GP_MAP", py::overload_cast<MultiOutputGP_GPU&, const int, const std::vector<double>>(&fit_GP_MAP), 
+     m.def("fit_GP_MAP", py::overload_cast<MultiOutputGP_GPU&, const int, const std::vector<double>>(&fit_GP_MAP), 
           py::return_value_policy::reference);
 
-    m.doc() = R"(
+     m.doc() = R"(
 The libgpgpu library
 --------------------
 
diff --git a/mogp_gpu/src/densegp_gpu.hpp b/mogp_gpu/src/densegp_gpu.hpp
index b0f58a25..43852899 100644
--- a/mogp_gpu/src/densegp_gpu.hpp
+++ b/mogp_gpu/src/densegp_gpu.hpp
@@ -25,9 +25,8 @@ These in turn use CUDA kernels defined in the file cov_gpu.cu
 #include "util.hpp"
 #include "kernel.hpp"
 #include "meanfunc.hpp"
-
-
-
+#include "gpparams.hpp"
+#include "gppriors.hpp"
 
 //////////////////////////////////////////////////////////////////////////////////////
 ////////////////////////// The Gaussian Process  class ///////////////////////////////
@@ -39,7 +38,7 @@ class DenseGP_GPU {
     unsigned int testing_size;
 
     // Number of training points, dimension of input
-    unsigned int n, D;
+    unsigned int n, D, n_corr;
 
     // inputs
     mat inputs;
@@ -62,9 +61,6 @@ class DenseGP_GPU {
     // hyperparameters
     vec current_theta;
 
-    // flag for whether we have fit hyperparameters
-    bool theta_fitted;
-
     // handle for CUBLAS calls
     cublasHandle_t cublasHandle;
 
@@ -89,12 +85,15 @@ class DenseGP_GPU {
     // adaptive nugget jitter, or specified fixed, or fitted, nugget size
     double nug_size;
 
+    // object holding all our hyperparameters
+    GPParams gptheta;
+
+    // object holding all our priors
+    GPPriors* priors;
+
     // kernel hyperparameters on the device
     thrust::device_vector<REAL> theta_d;
 
-    // mean function hyperparameters
-    vec meanfunc_params;
-
     // inputs, on the device, row major order
     thrust::device_vector<REAL> inputs_d;
 
@@ -143,51 +142,95 @@ class DenseGP_GPU {
         return targets;
     }
 
+    int get_n_kernel_params(void) const
+    {
+        return D + 1 + int(get_nugget_type()==NUG_FIT);
+    }
+
+    int get_n_corr(void) const
+    {
+        return n_corr;
+    }
+
     int get_n_params(void) const
     {
         return D + 1 + int(nug_type == NUG_FIT);
     }
 
-    vec get_theta(void)
+    GPParams get_theta(void) const
     {
-        return current_theta;
+        return gptheta;
     }
 
     double get_nugget_size(void) const
     {
-        return nug_size;
+        return gptheta.get_nugget_size();
     }
 
     void set_nugget_size(double nugget_size) 
     {
-        nug_size = nugget_size;
+        gptheta.set_nugget_size(nugget_size);
     }
 
     void set_nugget_type(nugget_type nugtype)
     {
-        nug_type = nugtype;
+        gptheta.set_nugget_type(nugtype);
     }
 
     nugget_type get_nugget_type(void) const
     {
-        return nug_type;
+        return gptheta.get_nugget_type();
     }
 
-    BaseMeanFunc* get_meanfunc(void)
+    BaseMeanFunc* get_meanfunc(void) const
     {
         return meanfunc;
     }
 
-    BaseKernel* get_kernel(void)
+    BaseKernel* get_kernel(void) const
     {
         return kernel;
     }
 
-    kernel_type get_kernel_type(void)
+    kernel_type get_kernel_type(void) const
     {
         return kern_type;
     }
 
+    GPPriors* get_gppriors(void) const
+    {
+        // This is the version of the function to be bound to the Python class
+        // WARNING - this will cause a (small) memory leak!
+        // but the alternative is that the destructor of the 
+        // GPPriors class will be called when Python does its garbage collection
+        return priors;
+    }
+
+    void set_gppriors(GPPriors* gppriors_)
+    {
+        priors = gppriors_;
+    }
+
+    void create_gppriors(int n_corr,
+                         std::vector< std::pair< prior_type, std::vector<REAL> > > corr_params,
+                         std::pair<prior_type, std::vector<REAL> > cov_params,
+                         std::pair<prior_type, std::vector<REAL> > nug_params)
+//                         prior_type ptype_cov_, REAL cov_p1_, REAL cov_p2_,
+  //                       prior_type ptype_nug_, REAL nug_p1_, REAL nug_p2_)
+    {
+        priors = new GPPriors(n_corr, nug_type=get_nugget_type());
+        priors->create_corr_priors(corr_params);
+        priors->create_cov_prior(cov_params);
+        priors->set_nugget(nug_params);
+        // default MeanPrior if required
+        if (gptheta.get_n_mean() > 0) {
+            int n_mean = gptheta.get_n_mean();
+            MeanPriors* mean_priors = new MeanPriors(vec::Zero(n_mean),mat::Identity(n_mean,n_mean));
+            priors->set_mean(mean_priors);
+        }
+
+    }
+
     // make a single prediction (mainly for testing - most use-cases will use predict_batch or predict_deriv_batch)
     double predict(mat_ref testing)
     {
@@ -202,6 +245,7 @@ class DenseGP_GPU {
         CUBLASDOT(cublasHandle, n, dev_ptr(work_d), 1, dev_ptr(invQt_d), 1,
                   &result);
 	// evaluate the mean function and add to the result
+        vec meanfunc_params = gptheta.get_mean();
         vec meanfunc_vals = meanfunc->mean_f(testing, meanfunc_params);
         return double(result + meanfunc_vals(0));
     }
@@ -247,6 +291,7 @@ class DenseGP_GPU {
 
         var = kappa - var.array();
 	    // evaluate the mean function and add to the result
+        vec meanfunc_params = gptheta.get_mean();
 	    vec meanfunc_vals = meanfunc->mean_f(testing, meanfunc_params);
         return REAL(result + meanfunc_vals(0));
     }
@@ -287,6 +332,7 @@ class DenseGP_GPU {
         thrust::copy(result_d.begin(), result_d.end(), result.data());
 
 	// evaluate the mean function and add to the result
+    vec meanfunc_params = gptheta.get_mean();
 	vec meanfunc_vals = meanfunc->mean_f(testing, meanfunc_params);
 	result += meanfunc_vals;
     }
@@ -309,7 +355,6 @@ class DenseGP_GPU {
                                      "than the maximum batch size");
         }
 
-
         thrust::device_vector<REAL> testing_d(testing.data(), testing.data() + Nbatch * D);
         thrust::device_vector<REAL> mean_d(Nbatch);
 
@@ -355,6 +400,7 @@ class DenseGP_GPU {
         // copy back means
         thrust::copy(mean_d.begin(), mean_d.end(), mean.data());
 	// evaluate the mean function and add to the mean
+    vec meanfunc_params = gptheta.get_mean();
 	vec meanfunc_vals = meanfunc->mean_f(testing, meanfunc_params);
 	mean += meanfunc_vals;
         // copy back variances
@@ -396,6 +442,7 @@ class DenseGP_GPU {
         thrust::copy(result_d.begin(), result_d.end(), result.data());
 
 	// evaluate deriv of meanfunc wrt test points, and add to result
+    vec meanfunc_params = gptheta.get_mean();
 	mat meanfunc_inputderiv = meanfunc->mean_inputderiv(testing, meanfunc_params);
 	result += meanfunc_inputderiv.transpose();
     }
@@ -433,36 +480,50 @@ class DenseGP_GPU {
         thrust::copy(chol_lower_d.begin(), chol_lower_d.end(), result.data());
     }
 
-    // Update the hyperparameters, and K, invQ and invQt which depend on them
-    void fit(vec_ref theta)
+    // Update the hyperparameters, given a new GPParams object
+    void fit(GPParams& theta_)
     {
+        if (! gptheta.test_same_shape(theta_)) 
+            throw std::runtime_error("Shape of new GPParams object does not match existing one");
+        gptheta = theta_;
+        vec new_params = gptheta.get_data();
+        fit(new_params);
+    }
 
-	    int param_switch_index = meanfunc->get_n_params();
-	    meanfunc_params = theta.block(0,0,param_switch_index,1);
-
-	    // evaluate the mean function at the input values and subtract from targets
+    void fit(vec params) {
+        // test that the new parameters are the expected size
+        if (! gptheta.test_same_shape(params)) 
+            throw std::runtime_error("Shape of new GPParams object does not match existing one");
+        int switch_index = gptheta.get_n_mean();
+        gptheta.set_data(params.block(switch_index,0,gptheta.get_n_data(), 1));
+        gptheta.set_mean(params.block(0,0, switch_index, 1));
+        if (nug_type == NUG_FIT) {
+            // set the nugget size to be the last element of params
+            gptheta.set_nugget_size(params[params.size()-1]);
+        }
+        vec meanfunc_params = gptheta.get_mean();
+        vec kernel_params = gptheta.get_data();
+        // evaluate the mean function at the input values and subtract from targets
 	    vec new_targets = targets - meanfunc->mean_f(inputs, meanfunc_params);
-	    thrust::copy(new_targets.data(), new_targets.data() + n, targets_d.begin());
-
-	    // copy all the parameters _after_ those for the mean function to the device vector theta_d
-        thrust::copy(theta.data() + param_switch_index,
-		         theta.data() + D + 2 + param_switch_index,
-		         theta_d.begin());
-
-	    // calculate covariance matrix, put result into K_d
+        thrust::copy(new_targets.data(), new_targets.data() + n, targets_d.begin());
+	    // copy all the non-meanfunc parameters to the device vector theta_d
+        thrust::copy(kernel_params.data(),
+		             kernel_params.data() + gptheta.get_n_data(),
+		             theta_d.begin());
+        // calculate covariance matrix, put result into K_d
         kernel->cov_batch_gpu(dev_ptr(K_d), n, n, D, dev_ptr(inputs_d),
 		    	      dev_ptr(inputs_d), dev_ptr(theta_d));
 	    /// copy the covariance matrix K_d into work_mat_d
 	    thrust::copy(K_d.begin(), K_d.end(), work_mat_d.begin());
-
         thrust::device_vector<int> info_d(1);
         int info_h;
         cusolverStatus_t status;
 	    int factorisation_status;
 	    // for adaptive nugget start with a nugget of zero and increase by small amount
 	    // until we find a value where factorization succeeds.
-	    if (nug_type == NUG_ADAPTIVE) {
+	    if (gptheta.get_nugget_type() == NUG_ADAPTIVE) {
 	        double mean_diag;
+            double tmp_nug_size = 0.;
 	        const int max_tries = 5;
 	        int itry = 0;
 	        while (itry < max_tries) {
@@ -479,19 +540,18 @@ class DenseGP_GPU {
 			            // elements should all be positive)
 		    	        cublasDasum(cublasHandle, n, dev_ptr(K_d), n+1, &mean_diag);
 			            mean_diag /= n;
-			            nug_size = 1e-6 * mean_diag;
+			            tmp_nug_size = 1e-6 * mean_diag;
 		            }
-		            add_diagonal(n, nug_size, dev_ptr(work_mat_d));
+		            add_diagonal(n, tmp_nug_size, dev_ptr(work_mat_d));
 		        }
 
 		        factorisation_status = calc_cholesky_factors();
 		        if (factorisation_status == 0) {
                     break;
 		        }
-		        nug_size *= 10;
+		        tmp_nug_size *= 10;
 		        itry++;
 	        }
-
 	        // if none of the factorization attempts succeeded:
 	        if (itry == max_tries) {
 		        std::string msg;
@@ -499,40 +559,28 @@ class DenseGP_GPU {
 		        smsg << "All attempts at factorization failed. Last return code " << factorisation_status;
 		        throw std::runtime_error(smsg.str());
 	        }
-	    // for fixed nugget, add "nugget_size" to the diagonal of the matrix.
-	    } else if (nug_type == NUG_FIXED) {
-	        add_diagonal(n, nug_size, dev_ptr(work_mat_d));
-	        factorisation_status = calc_cholesky_factors();
-	        if (factorisation_status != 0) {
-                throw std::runtime_error("Unable to factorize matrix using fixed nugget");
-	        }
-
-	    } else if (nug_type == NUG_FIT) {
-	        // set to exp(last-element-of-theta)
-	        nug_size = exp( theta(theta.size()-1) );
-
-	        add_diagonal(n, nug_size, dev_ptr(work_mat_d));
+            gptheta.set_nugget_size(tmp_nug_size);
+	    // for fixed or fitted nugget, add "nugget_size" to the diagonal of the matrix.
+	    } else if ((gptheta.get_nugget_type() == NUG_FIXED) || 
+                    (gptheta.get_nugget_type() == NUG_FIT)) {
+	        add_diagonal(n, gptheta.get_nugget_size(), dev_ptr(work_mat_d));
 	        factorisation_status = calc_cholesky_factors();
 	        if (factorisation_status != 0) {
-                throw std::runtime_error("Unable to factorize matrix using fitted nugget");
+                throw std::runtime_error("Unable to factorize matrix using selected nugget type");
 	        }
 	    } else throw std::runtime_error("Unrecognized nugget_type");
 
-
         // get the inverse covariance matrix invQ by solving the system of linear eqns
 	    //    work_mat_d . invQ_d = I
 	    // where work_mat_d is holding the current covariance matrix.
-
         identity_device(n, dev_ptr(invQ_d));
 
         status = cusolverDnDpotrs(cusolverHandle, CUBLAS_FILL_MODE_LOWER, n, n,
                                   dev_ptr(work_mat_d), n, dev_ptr(invQ_d), n,
                                   dev_ptr(info_d));
-
         thrust::copy(info_d.begin(), info_d.end(), &info_h);
         check_cusolver_status(status, info_h);
 
-
         // invQt - product of inverse covariance matrix with the target values
         thrust::copy(targets_d.begin(), targets_d.end(), invQt_d.begin());
         status = cusolverDnDpotrs(cusolverHandle, CUBLAS_FILL_MODE_LOWER, n, 1,
@@ -541,42 +589,36 @@ class DenseGP_GPU {
 
         thrust::copy(info_d.begin(), info_d.end(), &info_h);
         check_cusolver_status(status, info_h);
-
         // logdetC - sum the log of the diagonal elements of the Cholesky factor of covariance matrix (in work_mat_d)
         thrust::device_vector<double> logdetC_d(1);
 
         sum_log_diag(n, dev_ptr(work_mat_d), dev_ptr(logdetC_d), dev_ptr(sum_buffer_d), sum_buffer_size_bytes);
 
         thrust::copy(logdetC_d.begin(), logdetC_d.end(), &logdetC);
-
 	    //copy work_mat_d into the lower triangular Cholesky factor
 	    thrust::copy(work_mat_d.begin(), work_mat_d.begin()+n*n, chol_lower_d.begin());
 
 	    //set the flag to say we have fitted theta
-	    theta_fitted = true;
-        // copy mean function params then kernel params into current_theta
-	    current_theta.block(0,0,param_switch_index,1) = meanfunc_params;
-        thrust::copy(theta_d.begin(), theta_d.end(), current_theta.data()+param_switch_index);
-        
+	    gptheta.set_fitted_ok();
         // update the current_logpost
         double logpost;
         CUBLASDOT(cublasHandle, n, dev_ptr(targets_d), 1, dev_ptr(invQt_d), 1,
                   &logpost);
-
         logpost += logdetC + n * log(2.0 * M_PI);
 
         logpost = 0.5 * logpost;
-        current_logpost = logpost;
+        // subtract the priors log posterior
+        current_logpost = logpost - priors->logp(gptheta);
     }
 
     bool get_theta_fit_status(void)
     {
-        return theta_fitted;
+        return gptheta.data_has_been_set(); 
     }
 
     void reset_theta_fit_status(void)
     {
-        theta_fitted = false;
+        gptheta.unset_data();
     }
 
     void get_K(mat_ref K_h)
@@ -594,11 +636,21 @@ class DenseGP_GPU {
         thrust::copy(invQt_d.begin(), invQt_d.end(), invQt_h.data());
     }
 
-    double get_logpost(vec_ref new_theta)
+    double get_logpost(vec new_theta)
     {
+	    if (gptheta.test_close(new_theta)) {
+	        return current_logpost;
+	    }
+	    
+        //theta has changed - refit
+	    fit(new_theta);
+
+	    return current_logpost;
+    }
 
-        bool theta_close = (new_theta - current_theta).norm() < 1e-8;
-	    if (theta_close) {
+    double get_logpost(GPParams& new_theta)
+    {
+	    if (gptheta.test_close(new_theta)) {
 	        return current_logpost;
 	    }
 	    
@@ -679,6 +731,7 @@ class DenseGP_GPU {
 
 	    // first elements of result (up to meanfunc_nparam) will be from meanfunc,
 	    // then the next (D+1) from the Kernel.
+        vec meanfunc_params = gptheta.get_mean();
 	    int meanfunc_nparam = meanfunc_params.rows();
 	    if (meanfunc_nparam > 0) { // only copy data to device and do calculation if we need to.
 	        mat meanfunc_deriv = meanfunc->mean_deriv(inputs, meanfunc_params);
@@ -706,17 +759,18 @@ class DenseGP_GPU {
 	        REAL invQtSq = std::numeric_limits<double>::quiet_NaN();
 	        CUBLASDOT(cublasHandle, n, dev_ptr(invQt_d), 1, dev_ptr(invQt_d), 1,
 		        &invQtSq);
-
 	        // set the last element of result, putting it all together
 	        result(result.size()-1) = 0.5 * nug_size * (tr_invQ - invQtSq);
 	    }
+        // subtract the values from the priors
+        result = result - priors->dlogpdtheta(gptheta);  //priors_dlogpdtheta;
     }
 
     // destructor - just to see when it is being called
     ~DenseGP_GPU() {
-       // std::cout<<"In destructor of DenseGP_GPU"<<std::endl;
         delete meanfunc;
         delete kernel;
+        delete priors;
     }
 
     // constructor
@@ -732,20 +786,18 @@ class DenseGP_GPU {
         : testing_size(testing_size_)
         , n(inputs_.rows())
         , D(inputs_.cols())
+        , n_corr(0)
         , K_d(n * n, 0.0)
         , invQ_d(n * n, 0.0)
         , chol_lower_d(n * n, 0.0)
-        , theta_d(D + 2)
-	    , meanfunc_params(1) // resize later if necessary
         , inputs(inputs_)
         , targets(targets_)
 	    , kern_type(kern_)
         , kernel(0)
-	 //   , meanfunc(mean_)
-	    , nug_type(nugtype_)
+        , priors(0)
+	    , meanfunc(mean_)
+        , nug_type(nugtype_)
         , nug_size(nugsize_)
-	    , current_theta(1) // resize later
-        , theta_fitted(false)
         , inputs_d(inputs_.data(), inputs_.data() + D * n)
         , targets_d(targets_.data(), targets_.data() + n)
         , logdetC(0.0)
@@ -759,7 +811,6 @@ class DenseGP_GPU {
         , kappa_d(testing_size, 0.0)
         , invQk_d(testing_size * n, 0.0)
     {
-        //std::cout<<"in constructor of DenseGP_GPU"<<std::endl;
         cublasStatus_t cublas_status = cublasCreate(&cublasHandle);
         if (cublas_status != CUBLAS_STATUS_SUCCESS)
         {
@@ -787,22 +838,11 @@ class DenseGP_GPU {
 
         // if mean function is not provided, assume zero everywhere.
         if (!mean_) {
-           // std::cout<<"Creating a new zero meanfunction"<<std::endl;
             meanfunc = new ZeroMeanFunc();
         } else {
             // clone the meanfunction that we were given
-           // std::cout<<"Cloning the mean function"<<std::endl;
             meanfunc = mean_->clone();
         }
-         // make a polynomial mean function 
-        //std::vector< std::pair<int, int> > dims_powers;
-        //dims_powers.push_back(std::make_pair<int, int>(0,1));
-        //dims_powers.push_back(std::make_pair<int, int>(1,1));
-        // meanfunc = new PolyMeanFunc(dims_powers);   
-       //REAL mfval = 0.5;
-       //meanfunc = new FixedMeanFunc(mfval);
-       //    meanfunc = new ConstMeanFunc();
-                   
         
 	    if (kern_type == SQUARED_EXPONENTIAL) {
 	        kernel = new SquaredExponentialKernel();
@@ -810,10 +850,16 @@ class DenseGP_GPU {
 	        kernel = new Matern52Kernel();
 	    } else throw std::runtime_error("Unrecognized kernel type\n");
 
-	    // resize meanfunc_params vector here
-	    meanfunc_params.resize(meanfunc->get_n_params(),1);
+        int n_mean = meanfunc->get_n_params();
+        n_corr = kernel->get_n_params(inputs);
+        gptheta = GPParams(n_mean, n_corr, nugtype_, nugsize_);
+        // resize the device vector that will store kernel hyperparameters
+        theta_d.resize(gptheta.get_n_data());
 	    // resize the device vector that will store derivative of mean function
 	    meanfunc_deriv_d.resize(meanfunc->get_n_params() * inputs.rows());
+
+        gptheta.set_nugget_type(nugtype_);
+        gptheta.set_nugget_size(nugsize_);
         // resize current_theta vector
         current_theta.resize(meanfunc->get_n_params() + get_n_params(),1);
     }
diff --git a/mogp_gpu/src/fitting.hpp b/mogp_gpu/src/fitting.hpp
index e8574d4e..01f23b07 100644
--- a/mogp_gpu/src/fitting.hpp
+++ b/mogp_gpu/src/fitting.hpp
@@ -52,14 +52,12 @@ class GPWrapper {
     return logpostderiv;
   }
 
-  vec theta() {
+  GPParams theta() {
     return gp.get_theta();
   }
-
 };
 
 
-
 void fit_single_GP_MAP(DenseGP_GPU& gp, const int n_tries=15, const std::vector<double> theta0=std::vector<double>()) {
   // Fit the hyperparameters of a Gaussian Process by minimizing the
   // negative log-posterior.
@@ -68,23 +66,17 @@ void fit_single_GP_MAP(DenseGP_GPU& gp, const int n_tries=15, const std::vector<
   std::vector< std::vector<double> > all_params;  
   // if we're given an initial set of theta values, put this at the start of all_params
   if (theta0.size() > 0) {
-    if (theta0.size() != gp.get_n_params() ) 
+    if (theta0.size() != gp.get_n_params() )  // note this include meanfunc params
       throw std::runtime_error("length of theta0 must equal n_params of GP.");
     all_params.push_back(theta0);
   }
 
-  // generate random starting values for theta
-  std::random_device rd;
-  std::mt19937 e2(rd());
-
-  std::uniform_real_distribution<> dist(-2.5, 2.5);
-
+  // use priors to generate starting values for theta
   
   for (int itry=all_params.size(); itry<n_tries; ++itry) {
-    std::vector<REAL> v(gp.get_n_params());
-    std::generate(v.begin(), v.end(), [&dist, &e2](){ return dist(e2);});
-   
-    all_params.push_back(v);
+    //std::vector<REAL> meanfunc_params(gp.get_theta().get_n_mean(),0.);
+    std::vector<REAL> prior_samples = gp.get_gppriors()->sample();
+    all_params.push_back(prior_samples);
   }
   
   std::vector<double> minvals;
@@ -106,7 +98,11 @@ void fit_single_GP_MAP(DenseGP_GPU& gp, const int n_tries=15, const std::vector<
             },
             theta, -1);
       minvals.push_back(minf);
-      thetavals.push_back(gpw.theta());    
+      vec meanfunc_params = gpw.theta().get_mean();
+      vec kernel_params = gpw.theta().get_data();
+      vec fitted_params(meanfunc_params.size() + kernel_params.size() );
+      fitted_params << meanfunc_params, kernel_params;
+      thetavals.push_back(fitted_params);    
     } catch(std::exception &e) {
       std::cout << "dlib optimization failed: " << e.what() << std::endl;
     }
@@ -117,11 +113,9 @@ void fit_single_GP_MAP(DenseGP_GPU& gp, const int n_tries=15, const std::vector<
     gp.reset_theta_fit_status();
   } else {
       int minvalIndex = std::min_element(minvals.begin(),minvals.end()) - minvals.begin();
-      vec best_theta = thetavals.at(minvalIndex);
-      gp.fit(best_theta);
-
+      vec best_params = thetavals.at(minvalIndex);
+      gp.fit(best_params);
   }
-  
   return;
 }
 
diff --git a/mogp_gpu/src/gpparams.hpp b/mogp_gpu/src/gpparams.hpp
new file mode 100644
index 00000000..c476b538
--- /dev/null
+++ b/mogp_gpu/src/gpparams.hpp
@@ -0,0 +1,240 @@
+#ifndef GPPARAMS_HPP
+#define GPPARAMS_HPP
+
+#include <iostream>
+#include <math.h>
+#include "types.hpp"
+
+/* GPParams class, and helper functions
+Following the convention in GPParams.py we store the raw 
+data values in the "data" vector, and provide easy access to
+transformed subsets of this.
+    The transformations between the raw parameters :math:`\theta` and the
+    transformed ones :math:`(\beta, l, \sigma^2, \eta^2)` are as follows:
+    
+    1. **Mean:** No transformation; :math:`{\beta = \theta}`
+    2. **Correlation:** The raw values are transformed via
+       :math:`{l = \exp(-0.5\theta)}` such that the transformed values are the
+       correlation length associated with the given input.
+    3. **Covariance:** The raw value is transformed via
+       :math:`{\sigma^2 = \exp(\theta)}` so that the transformed value is
+       the covariance.
+    4. **Nugget:** The raw value is transformed via
+       :math:`{\eta^2 = \exp(\theta)}` so that the transformed value is
+       the variance associated with the nugget noise.
+*/
+
+class BaseTransform {
+    public:
+        virtual REAL raw_to_scaled(REAL r) const = 0;
+        virtual REAL scaled_to_raw(REAL s) const = 0;
+        virtual vec raw_to_scaled(vec r) const = 0;
+        virtual vec scaled_to_raw(vec s) const= 0;
+        virtual REAL dscaled_draw(REAL s) const= 0;
+        virtual REAL d2scaled_draw2(REAL s) const= 0;
+        virtual vec dscaled_draw(vec s) const = 0;
+        virtual vec d2scaled_draw2(vec s) const= 0;
+
+};
+
+
+class CorrTransform : public BaseTransform {
+
+public:
+    virtual REAL raw_to_scaled(REAL r) const{
+        return exp(-0.5 * r);
+    } 
+    virtual  REAL scaled_to_raw(REAL s) const{
+        return -2.0*log(s);
+    }   
+    virtual vec raw_to_scaled(vec r) const{
+        return Eigen::exp(-0.5*r.array());//.exp();    
+    }
+    virtual vec scaled_to_raw(vec s) const{
+        return -2.0*s.array().log();
+    } 
+    virtual REAL dscaled_draw(REAL s) const{
+        return -0.5*s;
+    }
+    virtual REAL d2scaled_draw2(REAL s) const{
+        return -0.25*s;
+    }
+    virtual vec dscaled_draw(vec s) const{
+        return -0.5*s;
+    }
+    virtual vec d2scaled_draw2(vec s) const{
+        return -0.25*s;
+    }
+};
+
+
+class CovTransform : public BaseTransform {
+
+public:
+
+    virtual REAL raw_to_scaled(REAL r) const {
+        return exp(r);
+    } 
+    virtual REAL scaled_to_raw(REAL s) const {
+        return log(s);
+    }
+    virtual vec raw_to_scaled(vec r) const { 
+        return r.array().exp();    
+    }
+    virtual vec scaled_to_raw(vec s) const {
+        return s.array().log();
+    }
+    virtual REAL dscaled_draw(REAL s) const {
+        return s;
+    }
+    virtual REAL d2scaled_draw2(REAL s) const {
+        return s;
+    }
+    virtual vec dscaled_draw(vec s) const {
+        return s;
+    }
+    virtual vec d2scaled_draw2(vec s) const {
+        return s;
+    }
+};
+
+class GPParams {
+    
+public:
+    // constructor, taking a string for the nugget
+    GPParams(int _n_mean=0, int _n_corr=1, nugget_type _nugget=NUG_FIT, REAL _nugget_size=0.0) 
+    : n_mean(_n_mean)
+    , n_corr(_n_corr)
+    , n_data(_n_corr + 1 + int(_nugget==NUG_FIT))
+    , nugget_size(_nugget_size)
+    , nug_type(_nugget)
+    , has_data(false)
+    {
+        mean = vec::Zero(n_mean);
+        data = vec::Zero(n_data);
+    }
+    // alternative constructor, for the fixed nugget case
+    GPParams(int _n_mean, int _n_corr, REAL _nugget) 
+    : GPParams(_n_mean, _n_corr, NUG_FIXED, _nugget ) {}
+
+    // getter and setter methods
+
+    inline int get_n_data() const { return n_corr + 1 + int(nug_type==NUG_FIT);}
+
+    void set_mean(vec _new_mean) {
+        mean = _new_mean;
+        n_mean = _new_mean.size();
+    }
+
+    inline vec get_mean() const { return mean; }
+
+    inline int get_n_mean() const { return n_mean; }  
+
+    inline vec get_data() const { return data; }
+
+    void set_data(vec _new_data) {
+        /// also resets Mean, covariance, nugget
+        if (_new_data.size() != n_data) 
+            throw std::runtime_error("New data not correct shape");
+        data = _new_data;
+        has_data = true;
+        mean = vec::Zero(n_mean);
+        if (nug_type == NUG_ADAPTIVE) nugget_size = 0.;
+    }
+
+    vec get_corr_raw() const {
+        return data.block(0, 0, n_corr, 1);
+    }
+
+    vec get_corr() const {
+        return CorrTransform().raw_to_scaled(get_corr_raw());
+    }
+
+    void set_corr(vec _new_corr) {
+        data.block(0,0,n_corr,1) = CorrTransform().scaled_to_raw(_new_corr);
+    }
+
+    inline int get_n_corr() { return n_corr; }
+
+    int get_cov_index() const {
+        // determine where in the data array the covariance parameter is
+        if (nug_type == NUG_FIT) return n_data - 2;
+        else return n_data - 1;
+    }
+
+    // transformed covariance parameter
+    REAL get_cov() const {
+       
+        if (! has_data) return 0.;
+        else return CovTransform().raw_to_scaled(data[get_cov_index()]);
+        
+    }
+
+    void set_cov(REAL _cov) {
+        if (! has_data) throw  std::runtime_error("Need to set data before setting covariance parameter");
+        data[get_cov_index()] = CovTransform().scaled_to_raw(_cov);
+    }
+
+    inline void set_nugget_type(nugget_type _nug_type) { nug_type = _nug_type; }
+
+    inline nugget_type get_nugget_type() const { return nug_type; }
+
+    void set_nugget_size(REAL _nug_size) {
+        nugget_size = _nug_size;
+        if ((nug_type == NUG_FIT) && (data.size()>0)) {
+            data[n_data -1] = nugget_size;
+        }
+    }
+
+    REAL get_nugget_size() const {    
+        if (nug_type != NUG_FIT) {
+            return nugget_size; 
+        } else {
+            return CovTransform().raw_to_scaled(data[n_data-1]);
+        }
+    }   
+
+    bool test_same_shape(GPParams& other) const {
+        if (n_mean != other.get_n_mean()) return false;
+        if (n_corr != other.get_n_corr()) return false;
+        if (nug_type != other.get_nugget_type()) return false;
+        return true;
+    }
+
+    bool test_same_shape(vec& other) const {
+        if ((get_n_mean() + get_n_data()) != other.size()) return false;
+        return true;
+    }
+
+    bool test_close(GPParams& other, REAL threshold=1e-8) const {
+        return (get_data() - other.get_data()).norm() < threshold;
+    }
+
+    bool test_close(vec& other, REAL threshold=1e-8) const {
+        return (get_data() - other).norm() < threshold;
+    }
+
+    void unset_data() { 
+        mean = vec::Zero(n_mean);
+        data = vec::Zero(n_data);
+        has_data = false;
+    }
+
+    void set_fitted_ok() {
+        has_data = true;
+    }
+
+    inline bool data_has_been_set() const { return has_data; }
+
+private:
+    REAL nugget_size;
+    nugget_type nug_type;
+    int n_mean;
+    int n_corr;
+    int n_data;
+    bool has_data;
+    vec data;
+    vec mean;
+};
+
+#endif
\ No newline at end of file
diff --git a/mogp_gpu/src/gppriors.hpp b/mogp_gpu/src/gppriors.hpp
new file mode 100644
index 00000000..06585f8d
--- /dev/null
+++ b/mogp_gpu/src/gppriors.hpp
@@ -0,0 +1,514 @@
+#ifndef GPPRIORS_HPP
+#define GPPRIORS_HPP
+
+#include <iostream>
+#include <vector>
+#include <math.h>
+#include <random>
+#include <boost/math/distributions/inverse_gamma.hpp>
+#include <boost/math/distributions/gamma.hpp>
+#include <boost/math/distributions/lognormal.hpp>
+#include "types.hpp"
+#include "util.hpp"
+#include "gpparams.hpp"
+
+
+class WeakPrior {
+    public:
+
+        virtual ~WeakPrior() {
+        }
+        inline virtual REAL logp(REAL x) {return 0.;}
+
+        inline virtual REAL dlogpdx(REAL x) {return 0.;}
+
+        inline virtual REAL d2logpdx2(REAL x) {return 0.;}
+
+        virtual REAL dlogpdtheta(REAL x, const BaseTransform& transform) {
+            return dlogpdx(x)* transform.dscaled_draw(x);
+        }
+
+        virtual REAL d2logpdtheta2(REAL x, const BaseTransform& transform) {
+            REAL term1 = d2logpdx2(x) * pow(transform.dscaled_draw(x),2); 
+            REAL term2 = dlogpdx(x) * (transform.d2scaled_draw2(x));
+            return term1 + term2;
+        }
+
+        virtual REAL sample(const BaseTransform& transform) {
+            return sample();
+        }
+
+        virtual REAL sample() { 
+            // see https://en.cppreference.com/w/cpp/numeric/random/random_device
+            // for details on this random number generator, and 
+            // https://stackoverflow.com/questions/36663027/how-to-seed-random-number-generator
+            // for suggestions on how to seed e.g. with current time.
+            std::random_device rd;
+            std::mt19937 e2(rd());
+            std::uniform_real_distribution<> dist(0.,5.);
+            return dist(e2) - 0.5;
+        }
+
+};
+/// introduce intermediate pure virtual base class between WeakPrior
+/// and the concrete Prior distributions, to implement the sample(transform) method
+class PriorDist : public WeakPrior {
+    public:
+
+        virtual REAL sample_x() = 0;
+
+        virtual REAL sample(const BaseTransform& transform) {
+            return transform.scaled_to_raw(sample_x());
+        }
+};
+
+class NormalPrior : public PriorDist {
+    public:
+        NormalPrior(REAL _mean, REAL _std)
+        : mean(_mean)
+        , std(_std) {}
+
+        REAL logp(REAL x) {
+            // computes log probability at given value
+            return -0.5 * pow((x - mean)/std, 2) - log(std) -0.5* log(2*M_PI);
+        }
+
+        REAL dlogpdx(REAL x) {
+            // derivative of log probability wrt scaled parameter
+            return -1.0*(x - mean)/pow(std,2) ;
+        }
+        REAL d2logpdx2(REAL x) {
+            // second derivative of log prob
+            return -1.0*pow(std, -2) ;
+
+        }
+        REAL sample_x() {
+            std::random_device rd;
+            std::mt19937 e2(rd());
+            std::normal_distribution<> dist(mean, std);
+            return dist(e2);
+        }
+
+    private:
+        REAL mean;
+        REAL std;
+};
+
+class LogNormalPrior : public PriorDist {
+    public:
+        LogNormalPrior(REAL _shape, REAL _scale)
+        : shape(_shape)
+        , scale(_scale) {
+            assert(shape > 0.);
+            assert(scale > 0.);
+        }
+
+        REAL logp(REAL x) {
+            // computes log probability at given value
+            assert(x> 0.);
+            return -0.5 * pow(log(x/scale)/shape,2) - 0.5*log(2*M_PI) - log(x) -log(shape);
+        }
+
+        REAL dlogpdx(REAL x) {
+            // derivative of log probability wrt scaled parameter
+            assert(x > 0.);
+            return -1.0*log(x/scale)/pow(shape,2)/x - 1./x;
+        }
+
+        REAL d2logpdx2(REAL x) {
+            // second derivative of log prob
+            assert(x>0.);
+            return (-1./pow(shape,2) + log(x/scale)/pow(shape,2) + 1.)/pow(x,2);
+        }
+
+        REAL sample_x() {
+            std::random_device rd;
+            std::mt19937 e2(rd());
+            std::lognormal_distribution<> dist(scale, shape);
+            return dist(e2);
+        }
+
+    private:
+        REAL shape;
+        REAL scale;
+};
+
+class GammaPrior : public PriorDist {
+    public:
+        GammaPrior(REAL _shape, REAL _scale)
+        : shape(_shape)
+        , scale(_scale) {
+            assert(shape > 0.);
+            assert(scale > 0.);
+        }
+
+        REAL logp(REAL x) {
+            // computes log probability at given value
+            assert(x> 0.);
+            return (-1.0*shape*log(scale) - lgamma(shape) +
+                (shape - 1.)*log(x) - x/scale);
+            
+        }
+
+        REAL dlogpdx(REAL x) {
+            // derivative of log probability wrt scaled parameter
+            assert(x > 0.);
+            return (shape - 1.)/x - 1./scale ;
+        }
+
+        REAL d2logpdx2(REAL x) {
+            // second derivative of log prob
+            assert(x>0.);
+            return -1.0*(shape - 1.)/pow(x,2);
+        }
+
+        REAL sample_x() {
+            std::random_device rd;
+            std::mt19937 e2(rd());
+            std::gamma_distribution<> dist(scale, shape);
+            return dist(e2);
+        }
+
+    private:
+        REAL shape;
+        REAL scale;
+};
+
+class InvGammaPrior : public PriorDist {
+    public:
+        InvGammaPrior(REAL _shape, REAL _scale)
+        : shape(_shape)
+        , scale(_scale) {
+            assert(shape > 0.);
+            assert(scale > 0.);
+        }
+
+        virtual ~InvGammaPrior() {
+        }
+
+        REAL logp(REAL x) {
+            return (shape*log(scale) - lgamma(shape) -
+                (shape + 1.)*log(x) - scale/x) ;       
+        }
+
+        REAL dlogpdx(REAL x) {
+            return -1.0*(shape + 1.)/x + scale/pow(x,2);
+        }
+
+        REAL d2logpdx2(REAL x) {
+            return (shape + 1)/pow(x,2) - 2.*scale/pow(x,3);
+        }
+
+        REAL sample_x() {
+            std::random_device rd;
+            std::mt19937 e2(rd());
+            std::gamma_distribution<> dist(shape, scale);
+            // inverse gamma is 1/Gamma, but multiply by scale^2 for 
+            // consistency with Python version
+            return scale*scale/dist(e2);      
+        }
+
+    private:
+        REAL shape;
+        REAL scale;
+};
+
+class MeanPriors {
+// prior distributions on mean function parameters.
+// Default instantiation will result in weak priors for these parameters.
+// Code is included that makes it possible to set other prior distributions,
+// but currently those would not be included in the logpost calculation, and 
+// are here as a placeholder for future work.
+public: 
+    MeanPriors(vec mean_, mat cov_)
+    : mean(mean_)
+    , cov(cov_) {
+        set_prior_dists();
+        weak_priors_only = true;
+    }
+
+    MeanPriors() {
+        vec default_mean;
+        mat default_cov;
+        MeanPriors(default_mean, default_cov);
+    }
+
+    inline vec get_mean() const {return mean;}
+
+    inline mat get_cov() const { return cov; }
+
+    inline int get_n_params() const { return mean.size();}
+
+    inline bool has_weak_priors() const { return weak_priors_only; }
+
+    void set_prior_dists() {
+        // This will create a set of weak priors for the meanfunc parameters
+        prior_dists.clear();
+        for (int i=0; i< get_n_params(); ++i ) {
+            prior_dists.push_back(new WeakPrior());
+        }
+    }
+
+    void set_prior_dists(std::vector<prior_type> ptype, 
+                         std::vector<std::vector<REAL> > priorparams) {
+        // This function can be called to set non-weak priors on meanfunc params.
+        // HOWEVER, these would not currently be included in the logposterior calculation.
+        // The code is therefore here as a placeholder for future work.
+        prior_dists.clear();
+        if ((ptype.size() != get_n_params()) || (priorparams.size() != get_n_params()))
+            throw std::runtime_error("Number of prior types and parameter sets must equal number of meanfunc parameters.");
+        for (int i=0; i< get_n_params(); ++i ) {
+            if (ptype[i] == INVGAMMA && priorparams[i].size()==2) {
+                prior_dists.push_back(new InvGammaPrior(priorparams[i][0],priorparams[i][1])); // shape and scale
+                weak_priors_only = false;
+            } else if (ptype[i] == GAMMA && priorparams[i].size()==2) {
+                prior_dists.push_back(new GammaPrior(priorparams[i][0],priorparams[i][1])); // shape and scale
+                weak_priors_only = false;
+            } else if (ptype[i] == LOGNORMAL && priorparams[i].size()==2) {
+                prior_dists.push_back(new LogNormalPrior(priorparams[i][0],priorparams[i][1])); // shape and scale
+                weak_priors_only = false;
+            } else {
+                std::cout<<"Unknown prior dist requested for MeanFunc param. Creating weak prior"<<std::endl;
+                prior_dists.push_back(new WeakPrior());
+            }
+        }
+    }
+    
+    std::vector<REAL> sample(const BaseTransform& transform) {
+        std::vector<REAL> vals;
+        for (int i=0; i< get_n_params(); ++i) {
+            vals.push_back(prior_dists[i]->sample(transform));
+        }
+        return vals;
+    }
+
+    // Function to multiply means of meanfunc parameter prior distributions with 
+    // design matrix.  This is not currently called in the C++ / GPU implementation
+    // but could be used in the future when the meanfunc is fit analytically.
+    vec dm_dot_b(mat_ref dm) const {
+        if (has_weak_priors()) return vec::Zero(dm.cols());
+        if (dm.cols() != mean.size()) 
+            throw std::runtime_error("Number of columns in design matrix doesn't match number of meanfunc parameters");
+        return dm * mean;
+    }
+
+    mat get_inv_cov() const {
+        if (has_weak_priors()) return mat::Zero(cov.rows(), cov.cols());
+        return cov.inverse();
+    }
+
+    vec get_inv_cov_b() const {
+       if (cov.size() == 0) return vec::Zero(1);
+       else  return get_inv_cov() * mean;
+    }
+
+    REAL logdet_cov() const {
+        if (has_weak_priors()) return 0.;
+        return log(cov.determinant());
+    }
+    
+private:
+    vec mean;
+    mat cov;
+    bool weak_priors_only;
+    std::vector<WeakPrior*> prior_dists;
+
+};
+
+class GPPriors {
+// The GPPriors object is the container that holds prior distributions for the 
+// correlation length and covariance parameters, the nugget (if fitted), and the
+// mean function parameters.  The current default implementation will only include 
+// weak priors for the mean function parameters.
+public:
+
+    GPPriors(int n_corr_, nugget_type nug_type_=NUG_FIT, MeanPriors* mean_=NULL, WeakPrior* cov_=NULL, WeakPrior* nug_prior_=NULL)
+    : n_corr(n_corr_)
+    , nug_type(nug_type_)
+    , mean_prior(mean_)
+    , cov_prior(cov_)
+    , nug_prior(nug_prior_) {
+    }
+
+
+    ~GPPriors() {
+        if (mean_prior != NULL) delete mean_prior;
+        if (cov_prior != NULL) delete cov_prior; 
+        if (nug_prior != NULL) delete nug_prior;          
+        corr_priors.clear(); 
+    }
+
+    inline MeanPriors* get_mean() { return mean_prior; }
+
+    inline void set_mean() { mean_prior = new MeanPriors();}
+
+    inline void set_mean(MeanPriors* _mean) {mean_prior = _mean;}
+
+    inline int n_mean() { return mean_prior->get_n_params();}
+
+    inline std::vector<WeakPrior*> get_corr() { return corr_priors; }
+
+    void set_corr(std::vector<WeakPrior*> _newcorr) {
+        corr_priors = _newcorr;
+        n_corr = corr_priors.size();
+    }
+
+    void set_corr() {
+        for (int i=0; i< n_corr; ++i) {
+            corr_priors.push_back(new WeakPrior());
+        }   
+    }
+
+    inline WeakPrior* get_cov() { return cov_prior;}
+
+    void set_cov(WeakPrior* _newcov) {
+        cov_prior = _newcov;
+    }
+
+    void set_cov() {
+        WeakPrior* wp = new WeakPrior();
+        set_cov(wp);
+    }
+
+    inline nugget_type get_nugget_type() { return nug_type;}
+
+    void set_nugget(WeakPrior* _newnugget) {
+        if (nug_type == NUG_FIT) nug_prior = _newnugget;
+    }
+
+    void set_nugget() {
+        if (nug_type == NUG_FIT) {
+            WeakPrior* wp = new WeakPrior();
+            set_nugget(wp);
+        }
+    }
+
+    void set_nugget(std::pair< prior_type, std::vector<REAL> > params_) {
+        if (nug_type == NUG_FIT) {
+            WeakPrior* wp = make_prior(params_.first, params_.second);
+            set_nugget(wp);
+        }
+    }
+
+    void check_theta(GPParams theta) {
+        assert(n_corr == theta.get_n_corr());
+        assert(nug_type == theta.get_nugget_type());
+        assert(theta.get_data().size() > 0);
+    }
+
+    REAL logp(GPParams theta) {
+        // Note that currently this calculation does not include any 
+        // non-weak MeanFunc priors
+        check_theta(theta);
+
+        REAL logposterior = 0.;
+        for (int i=0; i< corr_priors.size(); ++i) {
+            logposterior += corr_priors[i]->logp(theta.get_corr()[i]);
+        }
+        
+        logposterior += cov_prior->logp(theta.get_cov());
+        
+        if (nug_type == NUG_FIT) {
+            logposterior += nug_prior->logp(theta.get_nugget_size());
+        }
+
+
+        return logposterior;
+    }
+
+    vec dlogpdtheta(GPParams theta) {
+        // Note that currently this calculation does not include any 
+        // non-weak MeanFunc priors
+        check_theta(theta);
+
+        vec partials(theta.get_n_data());
+        //correlation length parameters
+        for (int i=0; i< corr_priors.size(); ++i) {
+            partials[i] = corr_priors[i]->dlogpdtheta(theta.get_corr()[i],CorrTransform());   
+        }
+        // covariance parameter
+        partials[theta.get_n_corr()] = cov_prior->dlogpdtheta(theta.get_cov(), CovTransform());
+        // nugget parameter (if fitted nugget)
+        if (nug_type == NUG_FIT) {
+        //    partials.push_back(nug_prior->dlogpdtheta(theta.get_nugget_size(), CovTransform()));
+            partials[theta.get_n_data()-1] = nug_prior->dlogpdtheta(theta.get_nugget_size(), CovTransform());  
+        }
+        return partials;
+    }
+
+    vec d2logpdtheta2(GPParams theta) {
+        // Note that currently this calculation does not include any 
+        // non-weak MeanFunc priors
+        check_theta(theta);
+
+        vec hessian(theta.get_n_data());
+        //correlation length parameters
+        for (int i=0; i< corr_priors.size(); ++i) {
+            hessian[i] = corr_priors[i]->d2logpdtheta2(theta.get_corr()[i],CorrTransform());
+        }
+        // covariance parameter
+        hessian[theta.get_n_corr()] = cov_prior->d2logpdtheta2(theta.get_cov(), CovTransform());
+        // nugget parameter (if fitted nugget)
+        if (nug_type == NUG_FIT) {
+            hessian[theta.get_n_data()-1] = nug_prior->d2logpdtheta2(theta.get_nugget_size(), CovTransform()); 
+        }
+        return hessian;
+    }
+
+    std::vector<REAL> sample() {
+        std::vector<REAL> samples; 
+        if ( mean_prior != NULL) {
+            samples = mean_prior->sample(CorrTransform());
+        } 
+        for (int i=0; i< corr_priors.size(); ++i) {
+            samples.push_back(corr_priors[i]->sample(CorrTransform()));
+        }  
+        samples.push_back(cov_prior->sample(CovTransform()));
+        if (nug_type == NUG_FIT) {
+            samples.push_back(nug_prior->sample(CovTransform())); 
+        }
+        return samples;
+    }
+
+    WeakPrior* make_prior(prior_type ptype, std::vector<REAL> priorparams) { 
+
+        if (ptype == INVGAMMA && priorparams.size()==2) {
+            return new InvGammaPrior(priorparams[0], priorparams[1]); // shape and scale
+        } else if (ptype == GAMMA && priorparams.size()==2) {
+            return new GammaPrior(priorparams[0], priorparams[1]); // shape and scale
+        } else if (ptype == LOGNORMAL && priorparams.size()==2) {
+            return new LogNormalPrior(priorparams[0], priorparams[1]); // shape and scale
+        } else if (ptype == WEAK) {
+            return new WeakPrior();
+        } else {
+            std::cout<<"Non-weak Prior must be Inverse Gamma, Gamma, or LogNormal"<<std::endl;  
+            return new WeakPrior();
+        }
+    }
+
+    void create_corr_priors(std::vector< std::pair< prior_type, std::vector<REAL> > > all_params) {
+        for (int i=0; i< all_params.size(); ++i) {
+            prior_type ptype = all_params[i].first;
+            std::vector<REAL> priorparams = all_params[i].second;
+            WeakPrior* new_prior = make_prior(ptype, priorparams);
+            corr_priors.push_back(new_prior);
+        }
+    }
+
+    void create_cov_prior(std::pair< prior_type, std::vector<REAL> > params_) {
+       cov_prior = make_prior(params_.first, params_.second);
+    }
+
+private:
+    MeanPriors* mean_prior;
+    std::vector<WeakPrior*> corr_priors;
+    int n_corr;
+    WeakPrior* cov_prior;
+    WeakPrior* nug_prior;
+    nugget_type nug_type;
+};
+
+
+
+
+#endif
\ No newline at end of file
diff --git a/mogp_gpu/src/kernel.hpp b/mogp_gpu/src/kernel.hpp
index db714a00..74b911dc 100644
--- a/mogp_gpu/src/kernel.hpp
+++ b/mogp_gpu/src/kernel.hpp
@@ -106,6 +106,9 @@ class BaseKernel {
     return result;
   }
 
+  int get_n_params(mat inputs) {
+    return inputs.cols();
+  }
 };
 
 
diff --git a/mogp_gpu/src/multioutputgp_gpu.hpp b/mogp_gpu/src/multioutputgp_gpu.hpp
index 4b4de187..d5321239 100644
--- a/mogp_gpu/src/multioutputgp_gpu.hpp
+++ b/mogp_gpu/src/multioutputgp_gpu.hpp
@@ -92,19 +92,45 @@ class MultiOutputGP_GPU {
         return emulators.size();
     }
 
-    int get_n_data_params(void) const
+
+    std::vector<unsigned int> n_data_params(void) const
     {   
-        if (emulators.size() > 0)
-            return emulators[0]->get_n_params();
-        return 0;
+        std::vector<unsigned int> n_params;
+        for (unsigned int i=0; i< emulators.size(); ++i) 
+            n_params.push_back(emulators[i]->get_theta().get_n_data());
+        return n_params;
     }
 
+    std::vector<unsigned int> n_corr_params(void) const
+    {   
+        std::vector<unsigned int> n_params;
+        for (unsigned int i=0; i< emulators.size(); ++i) 
+            n_params.push_back(emulators[i]->get_theta().get_n_corr());
+        return n_params;
+    }
+
+    inline nugget_type get_nugget_type(void) const { return nug_type;}
+
+    inline REAL get_nugget_size(void) const { return nug_size;}
+
     void reset_fit_status(void) {
         for (unsigned int idx=0; idx < emulators.size(); ++idx) {
             emulators[idx]->reset_theta_fit_status();
         }   
     }
 
+    void create_priors_for_emulator(unsigned int emulator_index,
+                                    int n_corr,
+                                    std::vector< std::pair< prior_type, std::vector<REAL> > > corr_params,
+                                    std::pair<prior_type, std::vector<REAL> > cov_params,
+                                    std::pair<prior_type, std::vector<REAL> > nug_params)
+                                   
+    {   
+        if (emulators.size() <= emulator_index)
+            throw std::runtime_error("Invalid emulator index for setting priors");
+        emulators[emulator_index]->create_gppriors(n_corr, corr_params, cov_params, nug_params);
+    }
+
     std::vector<unsigned int> get_fitted_indices(void) const
     {
         std::vector<unsigned int> fitted_indices;
@@ -187,18 +213,28 @@ class MultiOutputGP_GPU {
         }
     }
 
-    void fit(mat_ref thetas) {
+    void fit(std::vector<GPParams> thetas) {
         #pragma omp parallel for
         for (unsigned int i=0; i< emulators.size(); ++i) {      
-            emulators[i]->fit(thetas.row(i));
+            emulators[i]->fit(thetas[i]);
         }   
     }
 
+    void fit(mat_ref thetas) {
+        #pragma omp parallel for
+        for (unsigned int i=0; i< emulators.size(); ++i) {      
+            emulators[i]->fit(thetas.row(i));
+        }   
+    }    
 
-    void fit_emulator(unsigned int index, vec_ref theta) {
+    void fit_emulator(unsigned int index, GPParams& theta) {
         emulators.at(index)->fit(theta);
     }
 
+    void fit_emulator(unsigned int index, vec& theta) {
+        emulators.at(index)->fit(theta);
+    }    
+
     void create_emulators() {
         unsigned int testing_size_per_emulator = testing_size / targets.size();
         for (auto targ : targets) {
diff --git a/mogp_gpu/src/types.hpp b/mogp_gpu/src/types.hpp
index 27db46d1..0ecf0aac 100644
--- a/mogp_gpu/src/types.hpp
+++ b/mogp_gpu/src/types.hpp
@@ -31,4 +31,7 @@ enum nugget_type {NUG_ADAPTIVE, NUG_FIT, NUG_FIXED};
 // enum to allow python code to select Kernel function
 enum kernel_type {SQUARED_EXPONENTIAL, MATERN52};
 
+// enum for Prior types
+enum prior_type {INVGAMMA, GAMMA, LOGNORMAL, WEAK};
+
 #endif
diff --git a/mogp_gpu/src/util.hpp b/mogp_gpu/src/util.hpp
index e6753598..b5936a20 100644
--- a/mogp_gpu/src/util.hpp
+++ b/mogp_gpu/src/util.hpp
@@ -53,5 +53,9 @@ void add_diagonal(int N, double b, double *A);
 void sum_log_diag(int N, double *A, double *result, double *work, size_t work_size);
 void trace(int N, double *A, double *result, double *work, size_t work_size);
 
+template<typename Base, typename T>
+inline bool instanceof(const T *ptr) {
+    return dynamic_cast<const Base*>(ptr) != nullptr;
+}
 
 #endif
diff --git a/setup.py b/setup.py
index 78badfce..b199e4ec 100644
--- a/setup.py
+++ b/setup.py
@@ -8,8 +8,8 @@
 
 # version information
 MAJOR = 0
-MINOR = 6
-MICRO = 1
+MINOR = 7
+MICRO = 0
 PRERELEASE = 0
 ISRELEASED = True
 version = "{}.{}.{}".format(MAJOR, MINOR, MICRO)


Test patch:
diff --git a/mogp_emulator/tests/test_ExperimentalDesign.py b/mogp_emulator/tests/test_ExperimentalDesign.py
index f12d29bd..5c7e585f 100644
--- a/mogp_emulator/tests/test_ExperimentalDesign.py
+++ b/mogp_emulator/tests/test_ExperimentalDesign.py
@@ -4,7 +4,12 @@
 import pytest
 from inspect import signature
 import types
-from ..ExperimentalDesign import ExperimentalDesign, MonteCarloDesign, LatinHypercubeDesign
+from ..ExperimentalDesign import (
+    ExperimentalDesign,
+    MonteCarloDesign,
+    LatinHypercubeDesign,
+    MaxiMinLHC,
+)
 
 def test_ExperimentalDesign_init():
     "test the init method of ExperimentalDesign"
@@ -237,4 +242,62 @@ def test_LatinHypercubeDesign_sample():
     sample_expected = np.array([[3.0319074209362724, 1.4378127211566698],
                                 [3.8539110146362914, 0.186886121478102 ],
                                 [0.3690344832370781, 1.1036076298948434]])
-    assert_allclose(sample, sample_expected)
\ No newline at end of file
+    assert_allclose(sample, sample_expected)
+
+def test_MaxiMinLHC(monkeypatch):
+    "test the MaxiMin version of an LHC"
+
+    class MockLHCSample:
+        def __init__(self):
+            self.i = 0
+
+        def __call__(self, n_samples):
+            if self.i == 0:
+                self.i = 1
+                return np.array([[0.3], [0.7]])
+            elif self.i == 1:
+                self.i = 2
+                return np.array([[0.25], [0.75]])
+            else:
+                return np.array([[0.0], [1.0]])
+
+    def reset_lhc_draw_samples():
+        monkeypatch.setattr(
+            "mogp_emulator.ExperimentalDesign.LatinHypercubeDesign._draw_samples",
+            MockLHCSample(),
+        )
+
+    reset_lhc_draw_samples()
+
+    ed = MaxiMinLHC(1)
+    sample = ed._draw_samples(2, n_tries=1)
+    assert_allclose(sample, np.array([[0.3], [0.7]]))
+
+    reset_lhc_draw_samples()
+
+    ed = MaxiMinLHC(1)
+    sample = ed._draw_samples(2, n_tries=2)
+    assert_allclose(sample, np.array([[0.25], [0.75]]))
+
+    reset_lhc_draw_samples()
+
+    ed = MaxiMinLHC(1)
+    sample = ed._draw_samples(2, n_tries=3)
+    assert_allclose(sample, np.array([[0.0], [1.0]]))
+
+    reset_lhc_draw_samples()
+
+    ed = MaxiMinLHC(1)
+    sample = ed._draw_samples(2)
+    assert_allclose(sample, np.array([[0.0], [1.0]]))
+
+    with pytest.raises(AssertionError):
+        ed = MaxiMinLHC(1)
+        sample = ed._draw_samples(-1)
+
+    with pytest.raises(AssertionError):
+        ed = MaxiMinLHC(1)
+        sample = ed._draw_samples(1, n_tries=-1)
+
+    ed = MaxiMinLHC(1)
+    sample = ed._draw_samples(1, metric="cityblock")
diff --git a/mogp_emulator/tests/test_GaussianProcess.py b/mogp_emulator/tests/test_GaussianProcess.py
index 02a1e97d..a62fbb7d 100644
--- a/mogp_emulator/tests/test_GaussianProcess.py
+++ b/mogp_emulator/tests/test_GaussianProcess.py
@@ -194,7 +194,7 @@ def test_GaussianProcessGPU_nugget(x, y):
     assert gp.nugget_type == "adaptive"
 
     gp.nugget = "fit"
-    assert gp.nugget == 0.
+    assert gp.nugget == 1.
     assert gp.nugget_type == "fit"
 
     gp.nugget = 1.
@@ -321,16 +321,23 @@ def test_GaussianProcess_theta_GPParams(x, y):
         gp.theta = GPParams(n_corr=1, nugget="adaptive")
 
 @pytest.mark.skipif(not gpu_usable(), reason=GPU_NOT_FOUND_MSG)
-@pytest.mark.parametrize("nugget,sn", [(0., 1.), ("adaptive", 0.)]) # ("fit", np.log(1.e-6))])
-def test_GaussianProcessGPU_theta(x, y, nugget, sn):
+@pytest.mark.parametrize("mean,nugget,sn", [(None, 0., 1.), (None, "adaptive", 0.),
+                                           ("x[0]", "fit", np.log(1.e-6))])
+def test_GaussianProcessGPU_theta(x, y, mean, nugget, sn):
     "test the theta property of GaussianProcess (effectively the same as fit)"
 
-    # zero mean, zero nugget
+    if isinstance(nugget, float):
+        nugget_type = "fixed"
+    else:
+        nugget_type = nugget
 
-    gp = GaussianProcessGPU(x, y, nugget=nugget)
+    gp = GaussianProcessGPU(x, y, mean=mean, nugget=nugget, priors=GPPriors(n_corr=3, nugget_type=nugget_type))
 
+    with pytest.raises(RuntimeError):
+        gp.theta = np.ones(gp.n_params + 1)
+    
     theta = np.ones(gp.n_params)
-    if gp.nugget_type == "fit":
+    if nugget == "fit":
         theta[-1] = sn
 
     gp.theta = theta
@@ -341,17 +348,41 @@ def test_GaussianProcessGPU_theta(x, y, nugget, sn):
     else:
         assert_allclose(gp.nugget, np.exp(sn))
         noise = np.exp(sn)*np.eye(x.shape[0])
+
+ #   Q = np.exp(theta[switch + gp.D])*gp.kernel.kernel_f(x, x, theta[switch:(switch + gp.D)]) + noise
+#    ym = y - np.dot(gp._dm, theta[:switch]) TBD for GPU
+
     Q = gp.kernel.kernel_f(x, x, theta[:gp.D]) + noise
 
     L_expect = np.linalg.cholesky(Q)
-    invQt_expect = np.linalg.solve(Q, y)
-    logpost_expect = 0.5*(np.log(np.linalg.det(Q)) +
-                          np.dot(y, invQt_expect) +
-                          gp.n*np.log(2.*np.pi))
+#    invQt_expect = np.linalg.solve(Q, ym)  TBD for GPU
+#    logpost_expect = 0.5*(np.log(np.linalg.det(Q)) +
+#                          np.dot(ym, invQt_expect) +
+#                          gp.n*np.log(2.*np.pi))
 
-    assert_allclose(L_expect, gp.L)
-    assert_allclose(invQt_expect, gp.invQt)
-    assert_allclose(logpost_expect, gp.current_logpost)
+ #   assert_allclose(L_expect, gp.L)
+ #   assert_allclose(invQt_expect, gp.Kinv_t)
+ #   assert_allclose(logpost_expect, gp.current_logpost)
+
+@pytest.mark.skipif(not gpu_usable(), reason=GPU_NOT_FOUND_MSG)
+def test_GaussianProcess_theta_GPParams(x, y):
+    "test that we can set parameters using a GPParams object"
+
+    gp = GaussianProcess(x, y)
+
+    gpp = GPParams(n_corr=3, nugget="adaptive")
+    gpp.set_data(np.ones(4))
+
+    gp.theta = gpp
+
+    assert_allclose(gp.theta.get_data(), np.ones(4))
+    assert_allclose(gp.theta.nugget, 0.)
+
+    with pytest.raises(AssertionError):
+        gp.theta = GPParams(n_corr=3, nugget="fit")
+        
+    with pytest.raises(AssertionError):
+        gp.theta = GPParams(n_corr=1, nugget="adaptive")
 
 def test_GaussianProcess_theta_pivot():
     "test that pivoting works as expected"
@@ -500,11 +531,8 @@ def test_GaussianProcess_fit_logposterior(x, y, mean, nugget, sn):
                           
     mean_expect = np.linalg.solve(A, np.dot(gp._dm.T, Kinv_t_expect))
 
-    Kinv_t_mean_expect = np.linalg.solve(K, y - np.dot(gp._dm, mean_expect))
-
     assert_allclose(L_expect, gp.Kinv.L)
     assert_allclose(Kinv_t_expect, gp.Kinv_t)
-    assert_allclose(Kinv_t_mean_expect, gp.Kinv_t_mean, atol=1.e-10)
     assert_allclose(LA_expect, gp.Ainv.L)
     assert_allclose(mean_expect, gp.theta.mean)
     assert_allclose(logpost_expect, gp.current_logpost)
@@ -550,27 +578,37 @@ def test_GaussianProcessGPU_logposterior(x, y):
     "test logposterior method of GaussianProcessGPU"
 
     # logposterior already tested, but check that parameters are re-fit if changed
-
-    gp = GaussianProcessGPU(x, y, nugget = 0.)
+    gp = GaussianProcessGPU(x, y, nugget = 0., priors=GPPriors(n_corr=3, nugget_type="fixed"))
 
     theta = np.ones(gp.n_params)
     gp.fit(theta)
 
     theta = np.zeros(gp.n_params)
 
-    Q = gp.kernel.kernel_f(x, x, theta[:gp.D])
 
-    L_expect = np.linalg.cholesky(Q)
-    invQt_expect = np.linalg.solve(Q, y)
-    logpost_expect = 0.5*(np.log(np.linalg.det(Q)) +
-                          np.dot(y, invQt_expect) +
+    K = np.exp(theta[-2])*gp.kernel.kernel_f(x, x, theta[:-1])
+
+#    Q = gp.kernel.kernel_f(x, x, theta[:gp.D])
+
+    L_expect = np.linalg.cholesky(K)
+    Kinv_t_expect = np.linalg.solve(K, y)
+    logpost_expect = 0.5*(np.log(np.linalg.det(K)) +
+                          np.dot(y, Kinv_t_expect) +
                           gp.n*np.log(2.*np.pi))
 
     assert_allclose(logpost_expect, gp.logposterior(theta))
     assert_allclose(gp.L, L_expect)
-    assert_allclose(invQt_expect, gp.invQt)
+    assert_allclose(Kinv_t_expect, gp.Kinv_t)
     assert_allclose(logpost_expect, gp.current_logpost)
 
+    # check we can set theta back to none correctly
+
+    gp.theta = None
+    assert_allclose(gp.theta.get_data(), np.zeros(gp.n_params)) #GPU implementation resets to zero
+    assert gp.Kinv_t is None
+    assert gp.current_logpost is None
+  
+
 @pytest.fixture
 def dx():
     return 1.e-6
@@ -621,8 +659,11 @@ def test_GaussianProcessGPU_logpost_deriv(x, y, dx, nugget, sn):
     gp = GaussianProcessGPU(x, y, nugget=nugget)
 
     n = gp.n_params
-    theta = np.ones(n)
-    theta[-1] = sn
+    theta = np.zeros(n)
+    theta[:2] = -np.ones(2)
+    theta[2] = -2.
+    if gp.nugget_type == "fit":
+        theta[-1] = sn
 
     deriv = np.zeros(n)
 
@@ -631,7 +672,7 @@ def test_GaussianProcessGPU_logpost_deriv(x, y, dx, nugget, sn):
         dx_array[i] = dx
         deriv[i] = (gp.logposterior(theta) - gp.logposterior(theta - dx_array))/dx
 
-    assert_allclose(deriv, gp.logpost_deriv(theta), atol=1.e-7, rtol=1.e-5)
+    assert_allclose(deriv, gp.logpost_deriv(theta), atol=1.e-4, rtol=1.e-4)
 
 @pytest.mark.parametrize("mean,nugget,sn", [(None, 0., 1.),
                                             (None, "pivot", 1.), ("x[0]", "fit", np.log(1.e-6))])
@@ -793,7 +834,7 @@ def test_GaussianProcess_priors(x, y, dx, priors, nugget, sn):
     #
     # assert_allclose(hess, gp.logpost_hessian(theta), rtol=1.e-5, atol=1.e-5)
 
-def test_GaussianProcess_predict(x, y, dx):
+def test_GaussianProcess_predict(x, y):
     "test the predict method of GaussianProcess"
 
     # zero mean
@@ -857,6 +898,7 @@ def test_GaussianProcess_predict(x, y, dx):
     R = dm_test.T - np.dot(gp._dm.T, np.linalg.solve(K, Ktest.T))
 
     mu_expect = m + np.dot(Ktest, gp.Kinv_t_mean)
+    var_expect = np.exp(theta[-1]) - np.diag(np.dot(Ktest, np.linalg.solve(K, Ktest.T)))
     var_expect += np.diag(np.dot(R.T, np.linalg.solve(np.dot(gp._dm.T, np.linalg.solve(K, gp._dm)),
                                                       R)))
 
@@ -927,29 +969,24 @@ def test_GaussianProcess_predict(x, y, dx):
 
 @pytest.mark.skipif(not gpu_usable(), reason=GPU_NOT_FOUND_MSG)
 def test_GaussianProcessGPU_predict(x, y, dx):
-    "test the predict method of GaussianProcessGPU"
+    "test the predict method of GaussianProcessGPU vs an equivalent GaussianProcess"
 
     # zero mean
-
     gp = GaussianProcessGPU(x, y, nugget=0.)
     theta = np.ones(gp.n_params)
-
     gp.fit(theta)
 
     x_test = np.array([[2., 3., 4.]])
 
     mu, var, deriv = gp.predict(x_test)
 
-    K = gp.kernel.kernel_f(x, x, theta[:gp.D])
-    Ktest = gp.kernel.kernel_f(x_test, x, theta[:gp.D])
-
-    mu_expect = np.dot(Ktest, gp.invQt)
-    var_expect = np.exp(theta[-2]) - np.diag(np.dot(Ktest, np.linalg.solve(K, Ktest.T)))
+    # GP (CPU implementation) to compare it to
+    gpcpu = GaussianProcess(x, y, nugget=0.)
+    gpcpu.fit(theta)
+    mu_expect, var_expect, _ = gpcpu.predict(x_test)
 
     D = gp.D
-
     deriv_expect = np.zeros((1, D))
-
     for i in range(D):
         dx_array = np.zeros(D)
         dx_array[i] = dx
@@ -957,8 +994,7 @@ def test_GaussianProcessGPU_predict(x, y, dx):
 
     assert_allclose(mu, mu_expect)
     assert_allclose(var, var_expect)
-    assert_allclose(deriv, deriv_expect, atol=1.e-7, rtol=1.e-7)
-
+   
     # check that reshaping works as expected
 
     x_test = np.array([2., 3., 4.])
@@ -967,7 +1003,6 @@ def test_GaussianProcessGPU_predict(x, y, dx):
 
     assert_allclose(mu, mu_expect)
     assert_allclose(var, var_expect)
-    assert_allclose(deriv, deriv_expect, atol=1.e-7, rtol=1.e-7)
 
     # check that with 1D input data can handle 1D prediction data correctly
 
@@ -1046,22 +1081,19 @@ def test_GaussianProcessGPU_predict_nugget(x, y):
 
     gp = GaussianProcessGPU(x, y, nugget=nugget)
     theta = np.ones(gp.n_params)
-
     gp.fit(theta)
-
     preds = gp.predict(x)
 
-    K = gp.kernel.kernel_f(x, x, theta[:gp.D])
+    gpcpu = GaussianProcess(x, y, nugget=nugget)
+    gpcpu.fit(theta)
+    preds_expect = gpcpu.predict(x)
 
-    var_expect = np.exp(theta[-2]) + nugget - np.diag(np.dot(K, np.linalg.solve(K + np.eye(gp.n)*nugget, K)))
-
-    assert_allclose(preds.unc, var_expect, atol=1.e-7)
+    assert_allclose(preds.unc, preds_expect.unc, atol=1.e-7)
 
     preds = gp.predict(x, include_nugget=False)
-
-    var_expect = np.exp(theta[-2]) - np.diag(np.dot(K, np.linalg.solve(K + np.eye(gp.n)*nugget, K)))
-
-    assert_allclose(preds.unc, var_expect, atol=1.e-7)
+    preds_expect = gpcpu.predict(x, include_nugget=False)
+    
+    assert_allclose(preds.unc, preds_expect.unc, atol=1.e-7)
 
 def test_GaussianProcess_predict_pivot():
     "test that pivoting gives same predictions as standard version"
@@ -1106,6 +1138,59 @@ def test_GaussianProcess_predict_variance():
 
     assert_allclose(np.zeros(101), var, atol = 1.e-3)
 
+def test_GaussianProcess_predict_full_cov(x, y):
+    "Test the predictions with full covariance"
+
+    # zero mean
+
+    gp = GaussianProcess(x, y, nugget=0.0)
+    theta = np.ones(gp.n_params)
+
+    gp.fit(theta)
+
+    x_test = np.array([[2.0, 3.0, 4.0]])
+
+    mu, var, deriv = gp.predict(x_test, full_cov=True)
+
+    K = np.exp(theta[-1]) * gp.kernel.kernel_f(x, x, theta[:-1])
+    Kpredict = np.exp(theta[-1]) * gp.kernel.kernel_f(x_test, x_test, theta[:-1])
+    Ktest = np.exp(theta[-1]) * gp.kernel.kernel_f(x_test, x, theta[:-1])
+
+    mu_expect = np.dot(Ktest, gp.Kinv_t_mean)
+    var_expect = Kpredict - np.dot(Ktest, np.linalg.solve(K, Ktest.T))
+
+    assert_allclose(mu, mu_expect)
+    assert_allclose(var, var_expect)
+
+    # nonzero mean function
+
+    gp = GaussianProcess(x, y, mean="x[0]", nugget=0.0)
+
+    theta = np.ones(gp.n_params)
+
+    gp.fit(theta)
+
+    x_test = np.array([[2.0, 3.0, 4.0]])
+
+    mu, var, deriv = gp.predict(x_test, full_cov=True)
+
+    dm_test = gp.get_design_matrix(x_test)
+
+    m = np.dot(dm_test, gp.theta.mean)
+    K = np.exp(theta[-1]) * gp.kernel.kernel_f(x, x, theta[:-1])
+    Kpredict = np.exp(theta[-1]) * gp.kernel.kernel_f(x_test, x_test, theta[:-1])
+    Ktest = np.exp(theta[-1]) * gp.kernel.kernel_f(x_test, x, theta[:-1])
+    R = dm_test.T - np.dot(gp._dm.T, np.linalg.solve(K, Ktest.T))
+
+    mu_expect = m + np.dot(Ktest, gp.Kinv_t_mean)
+    var_expect = Kpredict - np.diag(np.dot(Ktest, np.linalg.solve(K, Ktest.T)))
+    var_expect += np.dot(
+        R.T, np.linalg.solve(np.dot(gp._dm.T, np.linalg.solve(K, gp._dm)), R)
+    )
+
+    assert_allclose(mu, mu_expect)
+    assert_allclose(var, var_expect)
+
 def test_GaussianProcess_predict_failures(x, y):
     "test situations where predict method of GaussianProcess should fail"
 
diff --git a/mogp_emulator/tests/test_HistoryMatching.py b/mogp_emulator/tests/test_HistoryMatching.py
index 4c8f212b..3a8c3b7d 100644
--- a/mogp_emulator/tests/test_HistoryMatching.py
+++ b/mogp_emulator/tests/test_HistoryMatching.py
@@ -4,6 +4,7 @@
 from ..HistoryMatching import HistoryMatching
 from ..fitting import fit_GP_MAP
 from ..GaussianProcess import GaussianProcess, PredictResult
+from ..MultiOutputGP import MultiOutputGP
 from ..LibGPGPU import gpu_usable
 from ..GaussianProcessGPU import GaussianProcessGPU
 
@@ -283,7 +284,7 @@ def test_HistoryMatching_init():
     hm = HistoryMatching(gp=gp, obs=1., coords=coords, expectations=None, threshold=5.)
 
     assert hm.gp == gp
-    assert_allclose(hm.obs, [1., 0.])
+    assert_allclose(hm.obs, [[1.], [0.]])
     assert_allclose(hm.coords, coords)
     assert hm.expectations is None
     assert hm.ndim == 1
@@ -300,7 +301,7 @@ def test_HistoryMatching_init():
                          threshold=5.)
 
     assert hm.gp == None
-    assert_allclose(hm.obs, [1., 0.1])
+    assert_allclose(hm.obs, [[1.], [0.1]])
     assert hm.coords == None
     for a in range(3):
         assert_allclose(hm.expectations[a], expectations[a])
@@ -403,6 +404,59 @@ def test_HistoryMatching_get_implausibility():
     with pytest.raises(AssertionError):
         hm.get_implausibility(-1.)
 
+def test_HistoryMatching_get_implausibility_mogp():
+    "Test implausibility with multiple outputs"
+    
+    # correct functionality
+
+    expectations = PredictResult(mean=np.array([[2., 10.], [4., 6.]]), unc=np.array([[0.5, 0.5], [0.5, 0.5]]),
+                                 deriv=None)
+    hm = HistoryMatching(obs=[[1., 5.], 0.5], expectations=expectations)
+    I = hm.get_implausibility()
+
+    assert_allclose(I, [1., 1.])
+    assert_allclose(hm.I, [1., 1.])
+
+    I = hm.get_implausibility(1.)
+
+    assert_allclose(I, [1./np.sqrt(2.), 1./np.sqrt(2.)])
+    assert_allclose(hm.I, [1./np.sqrt(2.), 1./np.sqrt(2.)])
+    
+    I = hm.get_implausibility(np.array([1., 1.]))
+
+    assert_allclose(I, [1./np.sqrt(2.), 1./np.sqrt(2.)])
+    assert_allclose(hm.I, [1./np.sqrt(2.), 1./np.sqrt(2.)])
+
+    np.random.seed(57483)
+    gp = fit_GP_MAP(np.reshape(np.linspace(0., 1.), (-1, 1)), [np.linspace(0., 1.), np.linspace(0., 1.)**2])
+    coords = np.array([[0.1], [1.]])
+    obs = [np.array([1., 5.]), np.array([0.01, 0.01])]
+    mean, unc, _ = gp.predict(coords)
+    I_exp = np.abs(mean - obs[0][:, np.newaxis])/np.sqrt(unc + obs[1][:,np.newaxis])
+    I_exp2 = [np.max(x) for x in I_exp.T]
+
+    hm = HistoryMatching(gp=gp, obs=obs, coords=coords)
+    I = hm.get_implausibility(rank=0)
+
+    assert_allclose(I, I_exp2)
+    assert_allclose(hm.I, I_exp2)
+    
+    np.random.seed(57483)
+    gp = fit_GP_MAP(np.reshape(np.linspace(0., 1.), (-1, 1)), [np.linspace(0., 1.),
+                                                               np.linspace(0., 1.)**2,
+                                                               np.linspace(0., 1.)**3])
+    coords = np.array([[0.1], [1.]])
+    obs = [np.array([1., 5., 4.]), np.array([0.01, 0.01, 0.01])]
+    mean, unc, _ = gp.predict(coords)
+    I_exp = np.abs(mean - obs[0][:, np.newaxis])/np.sqrt(unc + obs[1][:,np.newaxis])
+    I_exp2 = [np.min(x) for x in I_exp.T]
+
+    hm = HistoryMatching(gp=gp, obs=obs, coords=coords)
+    I = hm.get_implausibility(rank=2)
+
+    assert_allclose(I, I_exp2)
+    assert_allclose(hm.I, I_exp2)
+
 def test_HistoryMatching_get_NROY():
     "test the get_NROY method of HistoryMatching"
 
@@ -455,6 +509,14 @@ def test_HistoryMatching_set_gp():
     hm = HistoryMatching()
     hm.set_gp(gp)
 
+    assert hm.gp == gp
+    
+    gp = MultiOutputGP(np.reshape(np.linspace(0., 1.), (-1, 1)), [np.linspace(0., 1.),
+                                                                  np.linspace(0., 1.)**2])
+    
+    hm = HistoryMatching()
+    hm.set_gp(gp)
+
     assert hm.gp == gp
 
     # bad type for GP
@@ -491,12 +553,28 @@ def test_HistoryMatching_set_obs():
     hm = HistoryMatching()
     hm.set_obs(obs)
 
-    assert_allclose(hm.obs, obs)
+    assert_allclose(hm.obs, np.reshape(obs, (-1, 1)))
 
     obs = 1.
     hm.set_obs(obs)
 
-    assert_allclose(hm.obs, [obs, 0.])
+    assert_allclose(hm.obs, [[obs], [0.]])
+    
+    obs = [np.array([1., 2.]), 2.]
+
+    hm = HistoryMatching()
+    hm.set_obs(obs)
+
+    assert_allclose(hm.obs[0], obs[0])
+    assert_allclose(hm.obs[1], obs[1])
+    
+    obs = [np.array([1., 2.]), np.array([2., 4.])]
+
+    hm = HistoryMatching()
+    hm.set_obs(obs)
+
+    assert_allclose(hm.obs[0], obs[0])
+    assert_allclose(hm.obs[1], obs[1])
 
     # wrong number of values for obs
 
@@ -512,6 +590,11 @@ def test_HistoryMatching_set_obs():
 
     with pytest.raises(TypeError):
         hm.set_obs("abc")
+        
+    # array with wrong first dimension size
+    
+    with pytest.raises(AssertionError):
+        hm.set_obs(np.ones((3, 2)))
 
 def test_HistoryMatching_set_coords():
     "test the set_coords method of HistoryMatching"
diff --git a/mogp_emulator/tests/test_MultiOutputGP.py b/mogp_emulator/tests/test_MultiOutputGP.py
index 4d8a653a..eaec0a16 100644
--- a/mogp_emulator/tests/test_MultiOutputGP.py
+++ b/mogp_emulator/tests/test_MultiOutputGP.py
@@ -24,13 +24,12 @@ def test_MultiOutputGP_init(x, y):
     "Test function for correct functioning of the init method of GaussianProcess"
 
     gp = MultiOutputGP(x, y)
-    assert_allclose(x, gp.emulators[0].inputs)
-    assert_allclose(x, gp.emulators[1].inputs)
-    assert_allclose(y[0], gp.emulators[0].targets)
-    assert_allclose(y[1], gp.emulators[1].targets)
+    assert_allclose(x, gp.inputs)
+    assert_allclose(y, gp.targets)
     assert gp.D == 3
     assert gp.n == 2
     assert gp.n_emulators == 2
+    assert gp.n_params == [4, 4]
     
     gp = MultiOutputGP(x, y, mean="1")
     gp = MultiOutputGP(x, y, mean=["1", "x[0]"])
@@ -63,10 +62,9 @@ def test_MultiOutputGP_predict(x, y, dx):
     "test the predict method of GaussianProcess"
 
     gp = MultiOutputGP(x, y, nugget=0.)
-    theta = np.ones(gp.emulators[0].n_params)
+    thetas = [np.ones((n,)) for n in gp.n_params]
 
-    gp.emulators[0].fit(theta)
-    gp.emulators[1].fit(theta)
+    gp.fit(thetas)
 
     x_test = np.array([[2., 3., 4.]])
 
@@ -74,27 +72,30 @@ def test_MultiOutputGP_predict(x, y, dx):
 
     for i in range(2):
 
-        K = np.exp(theta[-1])*gp.emulators[i].kernel.kernel_f(x, x, theta[:-1])
-        Ktest = np.exp(theta[-1])*gp.emulators[i].kernel.kernel_f(x_test, x, theta[:-1])
+        K = np.exp(thetas[i][-1])*gp.emulators[i].kernel.kernel_f(x, x, thetas[i][:-1])
+        Ktest = np.exp(thetas[i][-1])*gp.emulators[i].kernel.kernel_f(x_test, x, thetas[i][:-1])
 
         mu_expect = np.dot(Ktest, gp.emulators[i].Kinv_t)
-        var_expect = np.exp(theta[-1]) - np.diag(np.dot(Ktest, np.linalg.solve(K, Ktest.T)))
+        var_expect = np.exp(thetas[i][-1]) - np.diag(np.dot(Ktest, np.linalg.solve(K, Ktest.T)))
 
         assert_allclose(mu[i], mu_expect)
         assert_allclose(var[i], var_expect)
+        
+    gp.reset_fit_status()
 
     nugget = 1.
     gp = MultiOutputGP(x, y, nugget=nugget)
-    theta = np.ones(gp.emulators[0].n_params)
+    theta0 = np.ones(gp.emulators[0].n_params)
+    theta1 = np.ones(gp.emulators[1].n_params)
 
-    gp.emulators[0].fit(theta)
-    gp.emulators[1].fit(theta)
+    gp.fit_emulator(0, theta0)
+    gp.fit_emulator(1, theta1)
 
     x_test = np.array([[2., 3., 4.]])
 
     mu, var, deriv = gp.predict(x_test)
 
-    for i in range(2):
+    for i, theta in zip(range(2), [theta0, theta1]):
 
         K = np.exp(theta[-1])*gp.emulators[i].kernel.kernel_f(x, x, theta[:-1]) + np.eye(gp.emulators[i].n)*nugget
         Ktest = np.exp(theta[-1])*gp.emulators[i].kernel.kernel_f(x_test, x, theta[:-1])
@@ -137,7 +138,7 @@ def test_MultiOutputGP_GPU_predict(x, y, dx):
     "test the predict method of GaussianProcess"
 
     gp = MultiOutputGP_GPU(x, y, nugget=0.)
-    theta = np.ones(gp.n_params)
+    theta = np.ones(gp.n_params[0])
 
     gp.fit_emulator(0,theta)
     gp.fit_emulator(1,theta)
@@ -159,7 +160,7 @@ def test_MultiOutputGP_GPU_predict(x, y, dx):
 
     nugget = 1.
     gp = MultiOutputGP_GPU(x, y, nugget=nugget)
-    theta = np.ones(gp.n_params)
+    theta = np.ones(gp.n_params[0])
 
     gp.fit_emulator(0,theta)
     gp.fit_emulator(1,theta)
@@ -227,7 +228,7 @@ def test_MultiOutputGP_GPU_check(x, y):
     """
 
     gp = MultiOutputGP_GPU(x, y, nugget=0.)
-    theta = np.ones(gp.n_params)
+    theta = np.ones(gp.n_params[0])
 
     assert gp.get_indices_fit() == []
     assert gp.get_indices_not_fit() == [0, 1]
@@ -239,4 +240,3 @@ def test_MultiOutputGP_GPU_check(x, y):
     gp.fit_emulator(1,theta) 
     assert gp.get_indices_fit() == [0, 1]
     assert gp.get_indices_not_fit() == []
-    
\ No newline at end of file
diff --git a/mogp_emulator/tests/test_fitting.py b/mogp_emulator/tests/test_fitting.py
index 6f43ffa1..11447d81 100644
--- a/mogp_emulator/tests/test_fitting.py
+++ b/mogp_emulator/tests/test_fitting.py
@@ -58,9 +58,11 @@ def test_fit_GP_MAP_GPU():
 
     gp = fit_GP_MAP(gp, theta0=theta_exp)
 
+    # remove value checks as machine precision means they might fail
+    # instead test basic functionality
     assert isinstance(gp, GaussianProcessGPU)
-    assert_allclose(gp.theta.data, theta_exp)
-    assert_allclose(gp.current_logpost, logpost_exp)
+    assert gp.theta.data_has_been_set()
+    assert gp.theta.get_data().shape ==  theta_exp.shape
 
 def test_fit_GP_MAP_failures():
     "test failures of fit_GP_MAP"
@@ -109,7 +111,7 @@ def test_fit_GP_MAP_GPU_failures():
     # minimization fails
 
     with pytest.raises(RuntimeError):
-        fit_GP_MAP(gp, n_tries=1, theta0=-10000.*np.ones(3))
+        fit_GP_MAP(gp, n_tries=1, theta0=-1000000.*np.ones(3))
 
     gp = GaussianProcessGPU(x, y, nugget=0.)
 
diff --git a/mogp_emulator/tests/test_linalg.py b/mogp_emulator/tests/test_linalg.py
index cd0e1b85..162a44ab 100644
--- a/mogp_emulator/tests/test_linalg.py
+++ b/mogp_emulator/tests/test_linalg.py
@@ -28,6 +28,10 @@ def test_ChoInv(A, b):
     
     assert_allclose(Ainv.solve(b), x)
     
+    x = np.linalg.solve(L, b)
+    
+    assert_allclose(Ainv.solve_L(b), x)
+    
     assert_allclose(np.log(np.linalg.det(A)), Ainv.logdet())
     
     assert Ainv.solve(np.zeros((3,0))).shape == (3,0)
@@ -54,6 +58,17 @@ def test_ChoInvPivot(A, b):
     x_pivot = Ainv.solve(b)
 
     assert_allclose(x, x_pivot)
+    
+    x = np.linalg.solve(Ainv.L, b[Ainv.P])
+    
+    x_pivot = Ainv.solve_L(b)
+    
+    assert_allclose(x, x_pivot)
+    
+    prod = np.dot(b, np.linalg.solve(A, b))
+    prod_pivot = np.dot(x_pivot, x_pivot)
+    
+    assert_allclose(prod, prod_pivot)
 
     with pytest.raises(AssertionError):
         ChoInvPivot(L_pivot, np.array([0, 2, 1, 1], dtype=np.int32)).solve(b)
diff --git a/mogp_emulator/tests/test_validation.py b/mogp_emulator/tests/test_validation.py
new file mode 100644
index 00000000..2402eee6
--- /dev/null
+++ b/mogp_emulator/tests/test_validation.py
@@ -0,0 +1,169 @@
+import numpy as np
+from numpy.testing import assert_allclose, assert_equal
+import pytest
+from ..GaussianProcess import GaussianProcess, PredictResult
+from ..MultiOutputGP import MultiOutputGP
+from ..validation import mahalanobis, standard_errors, pivoted_errors
+
+targets1 = np.array([0.5, 2.1, 2.8])
+targets2 = np.array([2.5, 2.9, 3.5])
+mean1 = np.array([1.0, 2.0, 3.0])
+mean2 = np.array([2.0, 3.0, 4.0])
+err1 = mean1 - targets1
+err2 = mean2 - targets2
+unc1 = np.array([[0.1, 0.05, 0.02], [0.05, 0.2, 0.01], [0.02, 0.01, 0.15]])
+
+
+def mock_predict(testing, unc=True, deriv=False, include_nugget=True, full_cov=False):
+    mean = mean1
+    if full_cov:
+        unc = unc1
+    else:
+        unc = np.diag(unc1)
+    return PredictResult(mean=mean, unc=unc, deriv=None)
+
+
+def mock_predict_mogp(
+    testing, unc=True, deriv=False, include_nugget=True, full_cov=False, processes=None
+):
+    mean = np.vstack([mean1, mean2])
+    if full_cov:
+        unc = np.stack([unc1] * 2)
+    else:
+        unc = np.vstack([np.diag(unc1)] * 2)
+    return PredictResult(
+        mean=mean,
+        unc=unc,
+        deriv=None,
+    )
+
+
+def mock_stats(*args, **kwds):
+    return 2.0, 3.0
+
+
+@pytest.fixture
+def valid_inputs():
+    return np.reshape(mean1, (-1, 1))
+
+
+@pytest.fixture
+def valid_targets():
+    return targets1
+
+
+@pytest.fixture
+def valid_targets_mogp():
+    return np.vstack([targets1, targets2])
+
+
+def test_standard_errors_GP(valid_inputs, valid_targets, monkeypatch):
+    "test standard errors for a GP"
+
+    monkeypatch.setattr("mogp_emulator.GaussianProcess.predict", mock_predict)
+
+    gp = GaussianProcess(valid_inputs, valid_targets, nugget=0.0)
+
+    errors = standard_errors(gp, valid_inputs, valid_targets)
+
+    idx1 = np.array([1, 2, 0])
+
+    assert_allclose(errors[0], err1[idx1] / np.sqrt(np.diag(unc1)[idx1]))
+    assert_equal(errors[1], idx1)
+
+
+def test_standard_errors_MOGP(valid_inputs, valid_targets_mogp, monkeypatch):
+    "test standard errors for a MOGP"
+
+    monkeypatch.setattr("mogp_emulator.MultiOutputGP.predict", mock_predict_mogp)
+
+    gp = MultiOutputGP(valid_inputs, valid_targets_mogp, nugget=0.0)
+
+    errors = standard_errors(gp, valid_inputs, valid_targets_mogp)
+
+    idx1 = np.array([1, 2, 0])
+    errors_expect = [
+        err1[idx1] / np.sqrt(np.diag(unc1))[idx1],
+        err2[idx1] / np.sqrt(np.diag(unc1))[idx1],
+    ]
+
+    for e, e_expect in zip(errors, errors_expect):
+        assert_allclose(e[0], e_expect)
+        assert_equal(e[1], idx1)
+
+
+def test_pivoted_errors_GP(valid_inputs, valid_targets, monkeypatch):
+    "test correlated errors for a MOGP"
+
+    monkeypatch.setattr("mogp_emulator.GaussianProcess.predict", mock_predict)
+
+    gp = GaussianProcess(valid_inputs, valid_targets, nugget=0.0)
+
+    errors = pivoted_errors(gp, valid_inputs, valid_targets)
+
+    idx1 = np.array([1, 2, 0])
+    A = np.linalg.cholesky(unc1[idx1][:, idx1])
+    b = np.linalg.solve(A, err1[idx1])
+
+    assert_allclose(errors[0], b)
+    assert_equal(errors[1], idx1)
+
+
+def test_pivoted_errors_MOGP(valid_inputs, valid_targets_mogp, monkeypatch):
+    "test correlated errors for a MOGP"
+
+    monkeypatch.setattr("mogp_emulator.MultiOutputGP.predict", mock_predict_mogp)
+
+    gp = MultiOutputGP(valid_inputs, valid_targets_mogp, nugget=0.0)
+
+    errors = pivoted_errors(gp, valid_inputs, valid_targets_mogp)
+
+    idx1 = np.array([1, 2, 0])
+    A = np.linalg.cholesky(unc1[idx1][:, idx1])
+    b = [np.linalg.solve(A, e[idx1]) for e in [err1, err2]]
+
+    for e, bval in zip(errors, b):
+        assert_allclose(e[0], bval)
+        assert_equal(e[1], idx1)
+
+
+def test_mahalanobis_GP(valid_inputs, valid_targets, monkeypatch):
+    "test correlated errors for a MOGP"
+
+    monkeypatch.setattr("mogp_emulator.GaussianProcess.predict", mock_predict)
+    monkeypatch.setattr(
+        "scipy.stats._distn_infrastructure.rv_generic.stats", mock_stats
+    )
+
+    gp = GaussianProcess(valid_inputs, valid_targets, nugget=0.0)
+
+    M = mahalanobis(gp, valid_inputs, valid_targets)
+
+    M_expect = np.dot(err1, np.linalg.solve(unc1, err1))
+
+    assert_allclose(M, M_expect)
+
+    M = mahalanobis(gp, valid_inputs, valid_targets, scaled=True)
+
+    assert_allclose(M, (M_expect - 2.0) / np.sqrt(3.0))
+
+
+def test_mahalanobis_MOGP(valid_inputs, valid_targets_mogp, monkeypatch):
+    "test correlated errors for a MOGP"
+
+    monkeypatch.setattr("mogp_emulator.MultiOutputGP.predict", mock_predict_mogp)
+    monkeypatch.setattr(
+        "scipy.stats._distn_infrastructure.rv_generic.stats", mock_stats
+    )
+
+    gp = MultiOutputGP(valid_inputs, valid_targets_mogp, nugget=0.0)
+
+    M = mahalanobis(gp, valid_inputs, valid_targets_mogp)
+
+    M_expect = np.array([np.dot(e, np.linalg.solve(unc1, e)) for e in [err1, err2]])
+
+    assert_allclose(M, M_expect)
+
+    M = mahalanobis(gp, valid_inputs, valid_targets_mogp, scaled=True)
+
+    assert_allclose(M, (M_expect - 2.0) / np.sqrt(3.0))
diff --git a/mogp_gpu/test/test_densegp_gpu.cu b/mogp_gpu/test/test_densegp_gpu.cu
index 0e9a148c..cb266c10 100644
--- a/mogp_gpu/test/test_densegp_gpu.cu
+++ b/mogp_gpu/test/test_densegp_gpu.cu
@@ -10,7 +10,7 @@
 
 #include <math.h>
 
-#include "../src/mog_gpu.hpp"
+#include "../src/densegp_gpu.hpp"
 
 typedef double REAL;
 
@@ -30,7 +30,8 @@ void testGP() {
   kernel_type kernel = SQUARED_EXPONENTIAL;
   // instantiate the GP
   DenseGP_GPU gp(inputs, targets, max_batch_size, meanfunc);
-  vec theta = vec::Constant(gp.get_n_params(),1, -1.0);
+  vec theta = vec::Constant(gp.get_n_data(),1, -1.0);
+  
   gp.fit(theta);
   mat x_predict(2,3);
   x_predict << 2., 3., 4.,
@@ -42,7 +43,28 @@ void testGP() {
 
 }
 
+void test_set_params() {
+  mat inputs(2,3);
+
+  inputs << 1., 2., 3.,
+            4., 5., 6.;
+  vec targets(2);
+  targets << 4., 6.;
+         
+  unsigned int max_batch_size = 2000;
+  ZeroMeanFunc* meanfunc = new ZeroMeanFunc();
+
+  kernel_type kernel = SQUARED_EXPONENTIAL;
+  // instantiate the GP
+  DenseGP_GPU gp(inputs, targets, max_batch_size, meanfunc);
+  GPParams theta = gp.get_theta();
+  assert(theta.get_nugget_type() == NUG_ADAPTIVE);
+  assert(theta.get_nugget_size() == 0.0);
+  gp.fit(theta);
+}
+
 int main() {
   testGP();
+  test_set_params();
   return 0;
 }
\ No newline at end of file
diff --git a/mogp_gpu/test/test_gpparams.cu b/mogp_gpu/test/test_gpparams.cu
new file mode 100644
index 00000000..3e9744e8
--- /dev/null
+++ b/mogp_gpu/test/test_gpparams.cu
@@ -0,0 +1,98 @@
+#include <iostream>
+
+#include <vector>
+#include <algorithm>
+#include <string>
+#include <sstream>
+#include <assert.h>
+#include <math.h>
+#include <stdexcept>
+#include "../src/types.hpp"
+#include "../src/gpparams.hpp"
+
+typedef double REAL;
+
+void test_corr_transforms() 
+{
+    std::cout<<"testing corr transforms"<<std::endl;
+    vec x_raw(3);
+    x_raw << 1.0, 2.0, 3.0 ;
+    CorrTransform corr_trans;
+    vec x_scaled = corr_trans.raw_to_scaled(x_raw);
+    assert (x_scaled[0] == exp(-0.5));
+    assert (x_scaled[1] == exp(-1.0));
+    assert (x_scaled[2] == exp(-1.5));
+    vec x_raw_new = corr_trans.scaled_to_raw(x_scaled);
+    assert (x_raw_new[0] == x_raw[0]);
+    assert (x_raw_new[1] == x_raw[1]);
+    assert (x_raw_new[2] == x_raw[2]);
+}
+
+void test_cov_transforms() 
+{
+    std::cout<<"testing cov transforms"<<std::endl;
+    vec x_raw(3);
+    x_raw << 1.0, 2.0, 3.0 ;
+    CovTransform cov_trans;
+    vec x_scaled = cov_trans.raw_to_scaled(x_raw);
+    assert (x_scaled[0] == exp(1.0));
+    assert (x_scaled[1] == exp(2.0));
+    assert (x_scaled[2] == exp(3.0));
+    vec x_raw_new = cov_trans.scaled_to_raw(x_scaled);
+    assert (x_raw_new[0] == x_raw[0]);
+    assert (x_raw_new[1] == x_raw[1]);
+    assert (x_raw_new[2] == x_raw[2]);
+}
+
+void test_construct_with_nug_type()
+{
+    std::cout<<"construct GPParams with nug_type for nugget"<<std::endl;
+    GPParams gpp(4,3, NUG_ADAPTIVE);
+    //std::cout<<" size is "<<gpp.get_corr_raw().size()<<std::endl;
+    assert (gpp.get_n_data() == 4);
+    assert (gpp.get_corr_raw().size()==3);
+    GPParams gpp_fit(4,3, NUG_FIT);
+    assert (gpp_fit.get_n_data()==5);
+
+}
+
+void test_construct_with_float()
+{
+    std::cout<<"construct GPParams with float for nugget"<<std::endl;
+    GPParams gpp(4,3, 0.5);
+    assert (gpp.get_n_mean() == 4);
+    assert (gpp.get_n_data() == 4);
+    assert (gpp.get_corr_raw().size()==3);
+    assert (gpp.get_nugget_type()==NUG_FIXED);
+}
+
+void test_eigen_block() {
+    vec x(9);
+    x << 1., 2., 3., 4., 5., 6., 7., 8., 9.;
+    vec x2 = x.block(1,0,4,1);
+    vec x3 = x.block(4,0,4,1);
+    std::cout<<"x2 "<<std::endl<<x2<<std::endl;
+    std::cout<<"x3 "<<std::endl<<x3<<std::endl;
+}
+
+void test_concatenate_vecs() {
+    GPParams gpp(4,3, 0.5);   
+    vec mean(4);
+    mean << 2., 3., 4., 5.;
+    gpp.set_mean(mean);
+    vec all_params(7); 
+    all_params << gpp.get_mean(), gpp.get_data();
+    assert (all_params.size() == 7);
+    std::cout<<"all params "<<all_params<<std::endl;
+}
+
+int main(void)
+{
+    test_construct_with_nug_type();
+    test_construct_with_float();
+    test_corr_transforms();
+    test_cov_transforms();
+    test_eigen_block();
+    test_concatenate_vecs();
+    return 0;
+}
diff --git a/mogp_gpu/test/test_gppriors.cu b/mogp_gpu/test/test_gppriors.cu
new file mode 100644
index 00000000..ab7705e5
--- /dev/null
+++ b/mogp_gpu/test/test_gppriors.cu
@@ -0,0 +1,120 @@
+#include <iostream>
+
+#include <vector>
+#include <set>
+#include <algorithm>
+#include <string>
+#include <sstream>
+#include <assert.h>
+#include <math.h>
+#include <stdexcept>
+#include <random>
+#include <boost/math/distributions/inverse_gamma.hpp>
+#include <boost/math/tools/roots.hpp>
+
+#include "../src/types.hpp"
+#include "../src/util.hpp"
+#include "../src/gpparams.hpp"
+#include "../src/gppriors.hpp"
+
+typedef double REAL;
+
+void test_normal_prior() {
+    NormalPrior np(2.,2.);
+    assert(abs(np.logp(3.) - (-1.73709)) < 0.001);
+    assert(np.dlogpdx(3.) == -0.25);
+    assert(np.d2logpdx2(3.) == -0.25);
+    std::cout<<"Normal Prior OK"<<std::endl;
+}
+
+void test_lognormal_prior() {
+    LogNormalPrior lnp(2.,2.);
+    std::cout<<" logp(3) is "<<lnp.logp(3.)<<std::endl;
+    std::cout<<" dlogpdx(3) is "<<lnp.dlogpdx(3.)<<std::endl;
+    std::cout<<" d2logpdx2(3) is "<<lnp.d2logpdx2(3.)<<std::endl;
+    assert(abs(lnp.logp(3.) - (-2.7312)) < 0.001);
+    assert(abs(lnp.dlogpdx(3.) - (-0.36712)) < 0.001);
+    assert(abs(lnp.d2logpdx2(3.) - (0.09459)) < 0.001);
+    std::cout<<"LogNormal Prior OK"<<std::endl;
+}
+
+void test_gamma_prior() {
+    GammaPrior gamp(2.,2.);
+    std::cout<<" logp(3) is "<<gamp.logp(3.)<<std::endl;
+    std::cout<<" dlogpdx(3) is "<<gamp.dlogpdx(3.)<<std::endl;
+    std::cout<<" d2logpdx2(3) is "<<gamp.d2logpdx2(3.)<<std::endl;
+    assert(abs(gamp.logp(3.) - (-1.78768)) < 0.001);
+    assert(abs(gamp.dlogpdx(3.) - (-0.16667)) < 0.001);
+    assert(abs(gamp.d2logpdx2(3.) - (-0.11111)) < 0.001);
+    std::cout<<"Gamma Prior OK"<<std::endl;
+}
+
+void test_invgamma_prior() {
+    InvGammaPrior igamp(2.,2.);
+    std::cout<<" logp(3) is "<<igamp.logp(3.)<<std::endl;
+    std::cout<<" dlogpdx(3) is "<<igamp.dlogpdx(3.)<<std::endl;
+    std::cout<<" d2logpdx2(3) is "<<igamp.d2logpdx2(3.)<<std::endl;
+    assert(abs(igamp.logp(3.) - (-2.5762)) < 0.001);
+    assert(abs(igamp.dlogpdx(3.) - (-0.77777)) < 0.001);
+    assert(abs(igamp.d2logpdx2(3.) - (0.185185)) < 0.001);
+    std::cout<<"InvGamma Prior OK"<<std::endl;
+   // WeakPrior wp = gamp.default_prior(1.,5.);
+}
+
+/*
+void test_spacings() {
+    vec x(5);
+    x << 1., 4., 3., 6., 3.;
+    assert(max_spacing(x) == 5.);
+    assert(median_spacing(x)==2.);
+    vec x2(9);
+    x2 << 1.2, 7.3, 4.6, 5.5, 5.7, 3.1, 6.5 ,3.0, 4.2;
+    std::cout<<" spacings "<<max_spacing(x2)<<" "<<median_spacing(x2)<<std::endl;
+    assert(max_spacing(x2) == 6.1);
+    assert(abs(median_spacing(x2)-0.8)<0.0001);
+}
+*/
+
+void test_isinstance() {
+    LogNormalPrior lnp(2.,2.);
+    bool is_lnp = instanceof<GammaPrior>(&lnp);
+    std::cout<<" is_lnp? "<<is_lnp<<std::endl;
+    
+}
+
+void test_root_finding() {
+    std::cout<<" testing root finding"<<std::endl;
+    boost::math::inverse_gamma_distribution<> dist(2., 2.);
+    REAL cdfval = boost::math::cdf(dist, 17.);
+    std::cout<<"cdfval "<<cdfval<<std::endl;
+    REAL min_val = 1.;
+    REAL max_val= 10.;
+    REAL mode = sqrt(min_val*max_val);
+    
+    auto f = [mode, max_val](REAL x) {
+        REAL a = exp(x);
+        return boost::math::cdf(boost::math::inverse_gamma_distribution<>(a,(1. + a)*mode), max_val) - 0.995;
+    };
+    std::cout<<" value of f(5.) is "<<f(5.)<<std::endl;
+    std::cout<<" value of f(8.) is "<<f(8.)<<std::endl;
+
+    int get_digits = std::numeric_limits<REAL>::digits - 3;
+    boost::math::tools::eps_tolerance<REAL> tol(get_digits); 
+    boost::uintmax_t it = 20;
+//    std::pair<REAL, REAL> r = boost::math::tools::bracket_and_solve_root( f, 5., 2., true, tol, it);
+    std::pair<REAL, REAL> r = boost::math::tools::bisect(f, -10., 10., tol, it);
+    REAL answer = r.first + (r.second - r.first)/2;
+    std::cout<<" value of f(5.) is "<<f(5.)<<" answer is "<<answer<<std::endl;
+}
+
+
+int main(void)
+{
+    test_normal_prior();
+    test_lognormal_prior();
+    test_gamma_prior();
+    test_invgamma_prior();
+    test_root_finding();
+ 
+    return 0;
+}