diff --git a/.travis.yml b/.travis.yml
index 63fa1c5e..f526bd13 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -49,11 +49,6 @@ after_success:
   - coveralls
 
 before_deploy:
-  - cd doc
-  - pip install sphinx_rtd_theme
-  - sphinx-apidoc -o source/ ../GPy
-  - make html
-  - cd ../
   - if [[ "$TRAVIS_OS_NAME" == "linux" ]];
     then
       export DIST='sdist';
@@ -69,5 +64,9 @@ deploy:
     secure: "vMEOlP7DQhFJ7hQAKtKC5hrJXFl5BkUt4nXdosWWiw//Kg8E+PPLg88XPI2gqIosir9wwgtbSBBbbwCxkM6uxRNMpoNR8Ixyv9fmSXp4rLl7bbBY768W7IRXKIBjpuEy2brQjoT+CwDDSzUkckHvuUjJDNRvUv8ab4P/qYO1LG4="
   on:
     branch: deploy
+  edge:
+    branch: v1.8.45
   distributions: $DIST
+  skip_existing: true
   skip_cleanup: true
+  skip_upload_docs: false
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 7e0d5c81..46cb6f69 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,22 +1,311 @@
 # Changelog
 
-## v1.8.5 (2017-12-01)
-
-### New Features
-
-* Implement [Latent Variable Multiple Output Gaussian Processes (LVMOGP)](https://arxiv.org/abs/1705.09862) [Zhenwen Dai]
-
-* Add mean function functionality to dtc inference method [Mark Pullin]
-
-* Allow non-zero mean GP prior for EP [Pablo Moreno]
+## v1.9.2 (2018-02-22)
 
 ### Fix
 
-* Fix DSYR function interface (to support SciPy 1.0) [Pablo Moreno]
+* Rtd. [mzwiessele]
 
-* Fix scipy=1.0.0 incompatibility of lyapunov [Alan Saul]
+* Rtd. [mzwiessele]
+
+* Rtd. [mzwiessele]
+
+* Rtd. [mzwiessele]
+
+* Rtd. [mzwiessele]
+
+* Rtd. [mzwiessele]
+
+* Rtd. [mzwiessele]
+
+* Rtd. [mzwiessele]
+
+* Rtd. [mzwiessele]
+
+* Rtd. [mzwiessele]
+
+* Rtd. [mzwiessele]
+
+* Rtd. [mzwiessele]
+
+* Rtd. [mzwiessele]
+
+* Rtd. [mzwiessele]
+
+* Rtd. [mzwiessele]
+
+* Rtd. [mzwiessele]
+
+* Rtd. [mzwiessele]
+
+### Other
+
+* Bump version: 1.9.1 → 1.9.2. [mzwiessele]
+
+
+## v1.9.1 (2018-02-22)
+
+### Fix
+
+* Paramz newest version. [mzwiessele]
+
+### Other
+
+* Bump version: 1.9.0 → 1.9.1. [mzwiessele]
+
+
+## v1.9.0 (2018-02-22)
+
+### Other
+
+* Bump version: 1.8.7 → 1.9.0. [mzwiessele]
+
+
+## v1.8.7 (2018-02-22)
+
+### Fix
+
+* Merge deploy back into devel. [mzwiessele]
+
+### Other
+
+* Bump version: 1.8.6 → 1.8.7. [mzwiessele]
+
+* Deploy version 1.8.5. [Zhenwen Dai]
+
+  * added extended version of MLP function with multiple hidden layers and different activation functions
+
+  * Update mapping_tests.py
+
+  Make output of gradient check verbose to diagnose error
+
+  * Update mapping_tests.py
+
+  Remove verbosity again after gradient checks passed without problem with verbosity
+
+  * the implementation of SVI-MOGP
+
+  * Try to fix the issue with model_tests
+
+  * updated mapping test to pass gradient checks
+
+  * Fix random seed for reproducible results in tests
+
+  * Add mean function functionality to dtc inference method
+
+  * Fix DSYR function (See https://github.com/scipy/scipy/issues/8155)
+
+  * Updated sde_kern to work with scipy=1.0.0
+
+  * Trying to fix tests for Matplotlib plotting issue
+
+  * Testing Again #575
+
+  * Figured it must be a matplotlib import error #575
+
+  New import matplotlib must be missing a package
+
+  * Removed ImageComparisonFailure #575
+
+  ImageComparisonFailure no longer exists which causes issues with travis testing using the most recent matplotlib
+
+  * Fix EP for non-zero mean GP priors
+
+  * improve the documentation for LVMOGP
+
+  * remove non-ascii characters
+
+  * Small correction to doc
+
+  * add type into docstring
+
+  * update changelog for 1.8.5
+
+  * bump the version: 1.8.4 -> 1.8.5
+
+
+## v1.8.6 (2018-02-22)
+
+### Fix
+
+* Gamma prior no assignment after init. [mzwiessele]
+
+* #568, product kernel resolution. [mzwiessele]
+
+* #590. [Max Zwiessele]
+
+  Y_normalized was not used for running optimization
+
+* Appveyor comment missing. [mzwiessele]
+
+### Other
+
+* Bump version: 1.8.5 → 1.8.6. [mzwiessele]
+
+* Merge pull request #597 from marpulli/devel. [Max Zwiessele]
+
+  Allow calculation of full predictive covariance matrices with multipl…
+
+* Allow calculation of full predictive covariance matrices with multiple outputs and normalization. [Mark Pullin]
+
+* Merge pull request #600 from marpulli/plotting_fix. [Max Zwiessele]
+
+  Fix visible dimensions for plotting inducing points
+
+* Fix visible dimensions for plotting inducing points. [Mark Pullin]
+
+* Merge pull request #599 from marpulli/grads_efficiency. [Zhenwen Dai]
+
+  Make predictive_gradients more efficient
+
+* Make predictive_gradients more efficient. [Mark Pullin]
+
+* Merge pull request #587 from esiivola/feature-multioutput. [Zhenwen Dai]
+
+  Merge the implementation of Multioutput kernel
+
+* Changed two function names so that they follow the python naming convention. [Siivola Eero]
+
+* Merge remote-tracking branch 'origin' into feature-multioutput. [Eero Siivola]
+
+* Merge pull request #592 from SheffieldML/sparsegp-normalization. [Zhenwen Dai]
+
+  fix: #590
+
+* Merge pull request #589 from apaleyes/devel. [Zhenwen Dai]
+
+  Implemented utility function to compute covariance between points in GP Model
+
+* Moved posterior_covariance to Posterior class. [Andrei Paleyes]
+
+* Implemented utility function to compute covariance between points in GP Model. [Andrei Paleyes]
+
+* Changed the structure of multioutput kernel so that it doesn't change the API of Kernels + documented the class. [Eero Siivola]
+
+* Merge remote-tracking branch 'origin/devel' into feature-multioutput. [Eero Siivola]
+
+* Merge pull request #585 from YoshikawaMasashi/devel. [Zhenwen Dai]
+
+  modify the MLP kernel equation
+
+* Modify the MLP kernel equation. [masashi yoshikawa]
+
+* Added multioutput kern and tests. [Eero Siivola]
+
+* Multioutput kernel + initial test. [Siivola Eero]
+
+* Multioutput kernel + initial test. [Siivola Eero]
+
+* Change dtype for Python 3 in robot_wirelss. [Neil Lawrence]
+
+* Bump the version: 1.8.4 -> 1.8.5. [Zhenwen Dai]
+
+* Update changelog for 1.8.5. [Zhenwen Dai]
+
+* Merge pull request #579 from SheffieldML/multi_out_doc. [Zhenwen Dai]
+
+  Improve the documentation for LVMOGP
+
+* Add type into docstring. [Zhenwen Dai]
+
+* Merge branch 'devel' of github.com:SheffieldML/GPy into multi_out_doc. [Zhenwen Dai]
+
+* Remove non-ascii characters. [Zhenwen Dai]
+
+* Improve the documentation for LVMOGP. [Zhenwen Dai]
+
+* Merge pull request #580 from marpulli/devel. [Zhenwen Dai]
+
+  Small correction to doc
+
+* Small correction to doc. [Mark Pullin]
+
+* Merge pull request #578 from pgmoren/devel. [Zhenwen Dai]
+
+  Fix EP for non-zero mean GP priors (binary classification)
+
+* Fix EP for non-zero mean GP priors. [Moreno]
+
+* Merge pull request #572 from marpulli/devel. [Alan Saul]
+
+  Add mean function functionality to dtc inference method
+
+* Add mean function functionality to dtc inference method. [Mark Pullin]
+
+* Merge pull request #573 from pgmoren/devel. [Zhenwen Dai]
+
+  Fix DSYR function (See https://github.com/scipy/scipy/issues/8155)
+
+* Fix DSYR function (See https://github.com/scipy/scipy/issues/8155) [Moreno]
+
+* Merge pull request #574 from alansaul/lyapunov_fix. [Alan Saul]
+
+  Fixing scipy=1.0.0 incompatibility of lyapunov discovered in PR #573. Coverage issue should be resolved by PR #575.
+
+* Updated sde_kern to work with scipy=1.0.0. [Alan Saul]
+
+* Merge pull request #575 from SheffieldML/matplotlib_testing. [Alan Saul]
+
+  Fixing tests for Matplotlib plotting issue
+
+* Removed ImageComparisonFailure #575. [Alan Saul]
+
+  ImageComparisonFailure no longer exists which causes issues with travis testing using the most recent matplotlib
+
+* Figured it must be a matplotlib import error #575. [Alan Saul]
+
+  New import matplotlib must be missing a package
+
+* Testing Again #575. [Alan Saul]
+
+* Trying to fix tests for Matplotlib plotting issue. [Alan Saul]
+
+* Merge pull request #526 from msbauer/mlp_extended. [Zhenwen Dai]
+
+  added extended version of MLP function
+
+* Fix random seed for reproducible results in tests. [msbauer]
+
+* Updated mapping test to pass gradient checks. [msbauer]
+
+* Update mapping_tests.py. [msbauer]
+
+  Remove verbosity again after gradient checks passed without problem with verbosity
+
+* Update mapping_tests.py. [msbauer]
+
+  Make output of gradient check verbose to diagnose error
+
+* Added extended version of MLP function with multiple hidden layers and different activation functions. [Bauer]
+
+* Merge pull request #562 from SheffieldML/external-mo. [Zhenwen Dai]
+
+  Release the implementation of LVMOGP
+
+* Try to fix the issue with model_tests. [Zhenwen Dai]
+
+* Merge with new changes from devel. [Zhenwen Dai]
+
+* Merge pull request #561 from SheffieldML/deploy. [Max Zwiessele]
+
+  Deploy
+
+* Merge pull request #560 from SheffieldML/devel. [Max Zwiessele]
+
+  appveyor twine upload error fix
+
+* Merge branch 'deploy' into devel. [Max Zwiessele]
+
+* Merge pull request #558 from SheffieldML/devel. [Max Zwiessele]
+
+  Uniform prior fix for other domains
+
+* Merge pull request #559 from SheffieldML/PS-upload-error. [Max Zwiessele]
+
+  Update appveyor.yml
+
+* The implementation of SVI-MOGP. [Zhenwen Dai]
 
-* Fix tests for Matplotlib plotting issue [Alan Saul]
 
 ## v1.8.4 (2017-10-06)
 
diff --git a/GPy/__version__.py b/GPy/__version__.py
index 89c6ad8e..2cbc28c3 100644
--- a/GPy/__version__.py
+++ b/GPy/__version__.py
@@ -1 +1 @@
-__version__ = "1.8.5"
+__version__ = "1.9.2"
diff --git a/GPy/core/gp.py b/GPy/core/gp.py
index d6f42c1c..453de407 100644
--- a/GPy/core/gp.py
+++ b/GPy/core/gp.py
@@ -282,10 +282,12 @@ class GP(Model):
             mu += self.mean_function.f(Xnew)
         return mu, var
 
-    def predict(self, Xnew, full_cov=False, Y_metadata=None, kern=None, likelihood=None, include_likelihood=True):
+    def predict(self, Xnew, full_cov=False, Y_metadata=None, kern=None, 
+                likelihood=None, include_likelihood=True):
         """
-        Predict the function(s) at the new point(s) Xnew. This includes the likelihood
-        variance added to the predicted underlying function (usually referred to as f).
+        Predict the function(s) at the new point(s) Xnew. This includes the
+        likelihood variance added to the predicted underlying function
+        (usually referred to as f).
 
         In order to predict without adding in the likelihood give
         `include_likelihood=False`, or refer to self.predict_noiseless().
@@ -295,33 +297,49 @@ class GP(Model):
         :param full_cov: whether to return the full covariance matrix, or just
                          the diagonal
         :type full_cov: bool
-        :param Y_metadata: metadata about the predicting point to pass to the likelihood
+        :param Y_metadata: metadata about the predicting point to pass to the
+                           likelihood
         :param kern: The kernel to use for prediction (defaults to the model
                      kern). this is useful for examining e.g. subprocesses.
-        :param bool include_likelihood: Whether or not to add likelihood noise to the predicted underlying latent function f.
+        :param include_likelihood: Whether or not to add likelihood noise to
+                                   the predicted underlying latent function f.
+        :type include_likelihood: bool
 
         :returns: (mean, var):
             mean: posterior mean, a Numpy array, Nnew x self.input_dim
-            var: posterior variance, a Numpy array, Nnew x 1 if full_cov=False, Nnew x Nnew otherwise
+            var: posterior variance, a Numpy array, Nnew x 1 if full_cov=False,
+                 Nnew x Nnew otherwise
 
-           If full_cov and self.input_dim > 1, the return shape of var is Nnew x Nnew x self.input_dim. If self.input_dim == 1, the return shape is Nnew x Nnew.
-           This is to allow for different normalizations of the output dimensions.
+            If full_cov and self.input_dim > 1, the return shape of var is
+            Nnew x Nnew x self.input_dim. If self.input_dim == 1, the return
+            shape is Nnew x Nnew. This is to allow for different normalizations
+            of the output dimensions.
 
-        Note: If you want the predictive quantiles (e.g. 95% confidence interval) use :py:func:"~GPy.core.gp.GP.predict_quantiles".
+        Note: If you want the predictive quantiles (e.g. 95% confidence
+        interval) use :py:func:"~GPy.core.gp.GP.predict_quantiles".
         """
-        #predict the latent function values
-        mu, var = self._raw_predict(Xnew, full_cov=full_cov, kern=kern)
+
+        # Predict the latent function values
+        mean, var = self._raw_predict(Xnew, full_cov=full_cov, kern=kern)
 
         if include_likelihood:
             # now push through likelihood
             if likelihood is None:
                 likelihood = self.likelihood
-            mu, var = likelihood.predictive_values(mu, var, full_cov, Y_metadata=Y_metadata)
+            mean, var = likelihood.predictive_values(mean, var, full_cov,
+                                                     Y_metadata=Y_metadata)
 
         if self.normalizer is not None:
-            mu, var = self.normalizer.inverse_mean(mu), self.normalizer.inverse_variance(var)
+            mean = self.normalizer.inverse_mean(mean)
 
-        return mu, var
+            # We need to create 3d array for the full covariance matrix with
+            # multiple outputs.
+            if full_cov & (mean.shape[1] > 1):
+                var = self.normalizer.inverse_covariance(var)
+            else:
+                var = self.normalizer.inverse_variance(var)
+
+        return mean, var
 
     def predict_noiseless(self,  Xnew, full_cov=False, Y_metadata=None, kern=None):
         """
@@ -376,13 +394,16 @@ class GP(Model):
 
     def predictive_gradients(self, Xnew, kern=None):
         """
-        Compute the derivatives of the predicted latent function with respect to X*
+        Compute the derivatives of the predicted latent function with respect
+        to X*
 
         Given a set of points at which to predict X* (size [N*,Q]), compute the
         derivatives of the mean and variance. Resulting arrays are sized:
-         dmu_dX* -- [N*, Q ,D], where D is the number of output in this GP (usually one).
+            dmu_dX* -- [N*, Q ,D], where D is the number of output in this GP
+            (usually one).
 
-        Note that this is not the same as computing the mean and variance of the derivative of the function!
+        Note that this is not the same as computing the mean and variance of
+        the derivative of the function!
 
          dv_dX*  -- [N*, Q],    (since all outputs have the same variance)
         :param X: The points at which to get the predictive gradients
@@ -393,25 +414,32 @@ class GP(Model):
         """
         if kern is None:
             kern = self.kern
-        mean_jac = np.empty((Xnew.shape[0],Xnew.shape[1],self.output_dim))
+        mean_jac = np.empty((Xnew.shape[0], Xnew.shape[1], self.output_dim))
 
         for i in range(self.output_dim):
-            mean_jac[:,:,i] = kern.gradients_X(self.posterior.woodbury_vector[:,i:i+1].T, Xnew, self._predictive_variable)
+            mean_jac[:, :, i] = kern.gradients_X(
+                self.posterior.woodbury_vector[:, i:i+1].T, Xnew,
+                self._predictive_variable)
 
-        # gradients wrt the diagonal part k_{xx}
-        dv_dX = kern.gradients_X(np.eye(Xnew.shape[0]), Xnew)
-        #grads wrt 'Schur' part K_{xf}K_{ff}^{-1}K_{fx}
+        # Gradients wrt the diagonal part k_{xx}
+        dv_dX = kern.gradients_X_diag(np.ones(Xnew.shape[0]), Xnew)
+
+        # Grads wrt 'Schur' part K_{xf}K_{ff}^{-1}K_{fx}
         if self.posterior.woodbury_inv.ndim == 3:
-            tmp = np.empty(dv_dX.shape + (self.posterior.woodbury_inv.shape[2],))
-            tmp[:] = dv_dX[:,:,None]
+            var_jac = np.empty(dv_dX.shape +
+                               (self.posterior.woodbury_inv.shape[2],))
+            var_jac[:] = dv_dX[:, :, None]
             for i in range(self.posterior.woodbury_inv.shape[2]):
-                alpha = -2.*np.dot(kern.K(Xnew, self._predictive_variable), self.posterior.woodbury_inv[:, :, i])
-                tmp[:, :, i] += kern.gradients_X(alpha, Xnew, self._predictive_variable)
+                alpha = -2.*np.dot(kern.K(Xnew, self._predictive_variable),
+                                   self.posterior.woodbury_inv[:, :, i])
+                var_jac[:, :, i] += kern.gradients_X(alpha, Xnew,
+                                                     self._predictive_variable)
         else:
-            tmp = dv_dX
-            alpha = -2.*np.dot(kern.K(Xnew, self._predictive_variable), self.posterior.woodbury_inv)
-            tmp += kern.gradients_X(alpha, Xnew, self._predictive_variable)
-        return mean_jac, tmp
+            var_jac = dv_dX
+            alpha = -2.*np.dot(kern.K(Xnew, self._predictive_variable),
+                               self.posterior.woodbury_inv)
+            var_jac += kern.gradients_X(alpha, Xnew, self._predictive_variable)
+        return mean_jac, var_jac
 
     def predict_jacobian(self, Xnew, kern=None, full_cov=False):
         """
@@ -564,7 +592,7 @@ class GP(Model):
         if self.output_dim == 1:
             return sim_one_dim(m, v)
         else:
-            fsim = np.empty((self.output_dim, self.num_data, size))
+            fsim = np.empty((self.output_dim, X.shape[0], size))
             for d in range(self.output_dim):
                 if full_cov and v.ndim == 3:
                     fsim[d] = sim_one_dim(m[:, d], v[:, :, d])
diff --git a/GPy/core/parameterization/priors.py b/GPy/core/parameterization/priors.py
index 3d69f39e..cbff4ca0 100644
--- a/GPy/core/parameterization/priors.py
+++ b/GPy/core/parameterization/priors.py
@@ -288,9 +288,17 @@ class Gamma(Prior):
         cls._instances.append(weakref.ref(o))
         return cls._instances[-1]()
 
+    @property
+    def a(self):
+        return self._a
+
+    @property
+    def b(self):
+        return self._b
+
     def __init__(self, a, b):
-        self.a = float(a)
-        self.b = float(b)
+        self._a = float(a)
+        self._b = float(b)
         self.constant = -gammaln(self.a) + a * np.log(b)
 
     def __str__(self):
@@ -333,8 +341,8 @@ class Gamma(Prior):
         return self.a, self.b
 
     def __setstate__(self, state):
-        self.a = state[0]
-        self.b = state[1]
+        self._a = state[0]
+        self._b = state[1]
         self.constant = -gammaln(self.a) + self.a * np.log(self.b)
 
 class InverseGamma(Gamma):
@@ -360,8 +368,8 @@ class InverseGamma(Gamma):
         return cls._instances[-1]()
 
     def __init__(self, a, b):
-        self.a = float(a)
-        self.b = float(b)
+        self._a = float(a)
+        self._b = float(b)
         self.constant = -gammaln(self.a) + a * np.log(b)
 
     def __str__(self):
diff --git a/GPy/examples/state_space.py b/GPy/examples/state_space.py
index 5a213f45..898d1676 100644
--- a/GPy/examples/state_space.py
+++ b/GPy/examples/state_space.py
@@ -4,23 +4,26 @@ import matplotlib.pyplot as plt
 
 import GPy.models.state_space_model as SS_model
 
-X = np.linspace(0, 10, 2000)[:, None]
-Y = np.sin(X) + np.random.randn(*X.shape)*0.1
+def state_space_example():
+    X = np.linspace(0, 10, 2000)[:, None]
+    Y = np.sin(X) + np.random.randn(*X.shape)*0.1
 
-kernel1 = GPy.kern.Matern32(X.shape[1])
-m1  = GPy.models.GPRegression(X,Y, kernel1)
+    kernel1 = GPy.kern.Matern32(X.shape[1])
+    m1  = GPy.models.GPRegression(X,Y, kernel1)
 
-print(m1)
-m1.optimize(optimizer='bfgs',messages=True)
+    print(m1)
+    m1.optimize(optimizer='bfgs',messages=True)
 
-print(m1)
+    print(m1)
 
-kernel2 = GPy.kern.sde_Matern32(X.shape[1])
-#m2  = SS_model.StateSpace(X,Y, kernel2)
-m2 = GPy.models.StateSpace(X,Y, kernel2)
-print(m2)
+    kernel2 = GPy.kern.sde_Matern32(X.shape[1])
+    #m2  = SS_model.StateSpace(X,Y, kernel2)
+    m2 = GPy.models.StateSpace(X,Y, kernel2)
+    print(m2)
 
-m2.optimize(optimizer='bfgs',messages=True)
+    m2.optimize(optimizer='bfgs',messages=True)
 
-print(m2)
+    print(m2)
+
+    return m1, m2
 
diff --git a/GPy/inference/latent_function_inference/expectation_propagation.py b/GPy/inference/latent_function_inference/expectation_propagation.py
index e92b58cb..61d3feff 100644
--- a/GPy/inference/latent_function_inference/expectation_propagation.py
+++ b/GPy/inference/latent_function_inference/expectation_propagation.py
@@ -132,6 +132,13 @@ class posteriorParamsDTC(posteriorParamsBase):
         self.mu += (delta_v-delta_tau*self.mu[i])*si
         #mu = np.dot(Sigma, v_tilde)
 
+    def to_dict(self):
+        return { "mu": self.mu.tolist(), "Sigma_diag": self.Sigma_diag.tolist()}
+
+    @staticmethod
+    def from_dict(input_dict):
+        return posteriorParamsDTC(np.array(input_dict["mu"]), np.array(input_dict["Sigma_diag"]))
+
     @staticmethod
     def _recompute(LLT0, Kmn, ga_approx):
         LLT = LLT0 + np.dot(Kmn*ga_approx.tau[None,:],Kmn.T)
@@ -533,3 +540,35 @@ class EPDTC(EPBase, VarDTC):
             #Posterior distribution parameters update
             if self.parallel_updates == False:
                 post_params._update_rank1(LLT, Kmn, delta_v, delta_tau, i)
+
+
+    def to_dict(self):
+        input_dict = super(EPDTC, self)._to_dict()
+        input_dict["class"] = "GPy.inference.latent_function_inference.expectation_propagation.EPDTC"
+        if self.ga_approx_old is not  None:
+            input_dict["ga_approx_old"] = self.ga_approx_old.to_dict()
+        if self._ep_approximation is not  None:
+            input_dict["_ep_approximation"] = {}
+            input_dict["_ep_approximation"]["post_params"] = self._ep_approximation[0].to_dict()
+            input_dict["_ep_approximation"]["ga_approx"] = self._ep_approximation[1].to_dict()
+            input_dict["_ep_approximation"]["cav_params"] = self._ep_approximation[2].to_dict()
+            input_dict["_ep_approximation"]["log_Z_tilde"] = self._ep_approximation[3].tolist()
+
+        return input_dict
+
+    @staticmethod
+    def _from_dict(inference_class, input_dict):
+        ga_approx_old = input_dict.pop('ga_approx_old', None)
+        if ga_approx_old is not None:
+            ga_approx_old = gaussianApproximation.from_dict(ga_approx_old)
+        _ep_approximation_dict = input_dict.pop('_ep_approximation', None)
+        _ep_approximation = []
+        if _ep_approximation is not None:
+            _ep_approximation.append(posteriorParamsDTC.from_dict(_ep_approximation_dict["post_params"]))
+            _ep_approximation.append(gaussianApproximation.from_dict(_ep_approximation_dict["ga_approx"]))
+            _ep_approximation.append(cavityParams.from_dict(_ep_approximation_dict["cav_params"]))
+            _ep_approximation.append(np.array(_ep_approximation_dict["log_Z_tilde"]))
+        ee = EPDTC(**input_dict)
+        ee.ga_approx_old = ga_approx_old
+        ee._ep_approximation = _ep_approximation
+        return ee
diff --git a/GPy/kern/src/prod.py b/GPy/kern/src/prod.py
index 43314e7a..31e62392 100644
--- a/GPy/kern/src/prod.py
+++ b/GPy/kern/src/prod.py
@@ -31,13 +31,16 @@ class Prod(CombinationKernel):
 
     """
     def __init__(self, kernels, name='mul'):
-        for i, kern in enumerate(kernels[:]):
+        _newkerns = []
+        for kern in kernels:
             if isinstance(kern, Prod):
-                del kernels[i]
-                for part in kern.parts[::-1]:
-                    kern.unlink_parameter(part)
-                    kernels.insert(i, part)
-        super(Prod, self).__init__(kernels, name)
+                for part in kern.parts:
+                    #kern.unlink_parameter(part)
+                    _newkerns.append(part.copy())
+            else:
+                _newkerns.append(kern.copy())
+
+        super(Prod, self).__init__(_newkerns, name)
 
     def to_dict(self):
         input_dict = super(Prod, self)._to_dict()
diff --git a/GPy/plotting/gpy_plot/gp_plots.py b/GPy/plotting/gpy_plot/gp_plots.py
index 230d47f0..a12fc858 100644
--- a/GPy/plotting/gpy_plot/gp_plots.py
+++ b/GPy/plotting/gpy_plot/gp_plots.py
@@ -337,7 +337,7 @@ def plot(self, plot_limits=None, fixed_inputs=None,
         plot_data = False
     plots = {}
     if hasattr(self, 'Z') and plot_inducing:
-        plots.update(_plot_inducing(self, canvas, visible_dims, projection, 'Inducing'))
+        plots.update(_plot_inducing(self, canvas, free_dims, projection, 'Inducing'))
     if plot_data:
         plots.update(_plot_data(self, canvas, which_data_rows, which_data_ycols, free_dims, projection, "Data"))
         plots.update(_plot_data_error(self, canvas, which_data_rows, which_data_ycols, free_dims, projection, "Data Error"))
diff --git a/GPy/testing/model_tests.py b/GPy/testing/model_tests.py
index c8d097a3..3558b9bb 100644
--- a/GPy/testing/model_tests.py
+++ b/GPy/testing/model_tests.py
@@ -118,6 +118,51 @@ class MiscTests(unittest.TestCase):
         from scipy.stats import norm
         np.testing.assert_allclose((mu+(norm.ppf(qs/100.)*np.sqrt(var))).flatten(), np.array(q95).flatten())
 
+    def test_multioutput_regression_with_normalizer(self):
+        """
+        Test that normalizing works in multi-output case
+        """
+
+        # Create test inputs
+        X = self.X
+        Y1 = np.sin(X) + np.random.randn(*X.shape) * 0.2
+        Y2 = -np.sin(X) + np.random.randn(*X.shape) * 0.05
+        Y = np.hstack((Y1, Y2))
+
+        mu, std = Y.mean(0), Y.std(0)
+        m = GPy.models.GPRegression(X, Y, normalizer=True)
+        m.optimize(messages=True)
+        assert(m.checkgrad())
+        k = GPy.kern.RBF(1)
+        m2 = GPy.models.GPRegression(X, (Y-mu)/std, normalizer=False)
+        m2[:] = m[:]
+
+        mu1, var1 = m.predict(m.X, full_cov=True)
+        mu2, var2 = m2.predict(m2.X, full_cov=True)
+        np.testing.assert_allclose(mu1, (mu2*std)+mu)
+        np.testing.assert_allclose(var1, var2[:, :, None]*std[None, None, :]**2)
+
+        mu1, var1 = m.predict(m.X, full_cov=False)
+        mu2, var2 = m2.predict(m2.X, full_cov=False)
+
+        np.testing.assert_allclose(mu1, (mu2*std)+mu)
+        np.testing.assert_allclose(var1, var2*std[None, :]**2)
+
+        q50n = m.predict_quantiles(m.X, (50,))
+        q50 = m2.predict_quantiles(m2.X, (50,))
+
+        np.testing.assert_allclose(q50n[0], (q50[0]*std)+mu)
+
+        # Test variance component:
+        qs = np.array([2.5, 97.5])
+        # The quantiles get computed before unormalization
+        # And transformed using the mean transformation:
+        c = np.random.choice(X.shape[0])
+        q95 = m2.predict_quantiles(X[[c]], qs)
+        mu, var = m2.predict(X[[c]])
+        from scipy.stats import norm
+        np.testing.assert_allclose((mu.T+(norm.ppf(qs/100.)*np.sqrt(var))).T.flatten(), np.array(q95).flatten())
+
     def check_jacobian(self):
         try:
             import autograd.numpy as np, autograd as ag, GPy, matplotlib.pyplot as plt
diff --git a/GPy/testing/serialization_tests.py b/GPy/testing/serialization_tests.py
index 80dfd219..7eb3fe5c 100644
--- a/GPy/testing/serialization_tests.py
+++ b/GPy/testing/serialization_tests.py
@@ -116,11 +116,45 @@ class Test(unittest.TestCase):
         np.testing.assert_array_equal(e1._ep_approximation[2].v[:], e1_r._ep_approximation[2].v[:])
         np.testing.assert_array_equal(e1._ep_approximation[3][:], e1_r._ep_approximation[3][:])
 
+
+        e1 = GPy.inference.latent_function_inference.expectation_propagation.EPDTC(ep_mode="nested")
+        e1.ga_approx_old = GPy.inference.latent_function_inference.expectation_propagation.gaussianApproximation(np.random.rand(10),np.random.rand(10))
+        e1._ep_approximation = []
+        e1._ep_approximation.append(GPy.inference.latent_function_inference.expectation_propagation.posteriorParamsDTC(np.random.rand(10),np.random.rand(10)))
+        e1._ep_approximation.append(GPy.inference.latent_function_inference.expectation_propagation.gaussianApproximation(np.random.rand(10),np.random.rand(10)))
+        e1._ep_approximation.append(GPy.inference.latent_function_inference.expectation_propagation.cavityParams(10))
+        e1._ep_approximation[-1].v = np.random.rand(10)
+        e1._ep_approximation[-1].tau = np.random.rand(10)
+        e1._ep_approximation.append(np.random.rand(10))
+        e1_r = GPy.inference.latent_function_inference.LatentFunctionInference.from_dict(e1.to_dict())
+
+
+        assert type(e1) == type(e1_r)
+        assert e1.epsilon==e1_r.epsilon
+        assert e1.eta==e1_r.eta
+        assert e1.delta==e1_r.delta
+        assert e1.always_reset==e1_r.always_reset
+        assert e1.max_iters==e1_r.max_iters
+        assert e1.ep_mode==e1_r.ep_mode
+        assert e1.parallel_updates==e1_r.parallel_updates
+
+        np.testing.assert_array_equal(e1.ga_approx_old.tau[:], e1_r.ga_approx_old.tau[:])
+        np.testing.assert_array_equal(e1.ga_approx_old.v[:], e1_r.ga_approx_old.v[:])
+        np.testing.assert_array_equal(e1._ep_approximation[0].mu[:], e1_r._ep_approximation[0].mu[:])
+        np.testing.assert_array_equal(e1._ep_approximation[0].Sigma_diag[:], e1_r._ep_approximation[0].Sigma_diag[:])
+        np.testing.assert_array_equal(e1._ep_approximation[1].tau[:], e1_r._ep_approximation[1].tau[:])
+        np.testing.assert_array_equal(e1._ep_approximation[1].v[:], e1_r._ep_approximation[1].v[:])
+        np.testing.assert_array_equal(e1._ep_approximation[2].tau[:], e1_r._ep_approximation[2].tau[:])
+        np.testing.assert_array_equal(e1._ep_approximation[2].v[:], e1_r._ep_approximation[2].v[:])
+        np.testing.assert_array_equal(e1._ep_approximation[3][:], e1_r._ep_approximation[3][:])
+
+
         e2 = GPy.inference.latent_function_inference.exact_gaussian_inference.ExactGaussianInference()
         e2_r = GPy.inference.latent_function_inference.LatentFunctionInference.from_dict(e2.to_dict())
 
         assert type(e2) == type(e2_r)
 
+
     def test_serialize_deserialize_model(self):
         np.random.seed(fixed_seed)
         N = 20
diff --git a/GPy/testing/util_tests.py b/GPy/testing/util_tests.py
index 5cd275c2..bdab63e8 100644
--- a/GPy/testing/util_tests.py
+++ b/GPy/testing/util_tests.py
@@ -28,7 +28,8 @@
 # OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #===============================================================================
 
-import unittest, numpy as np
+import unittest
+import numpy as np
 import GPy
 
 class TestDebug(unittest.TestCase):
@@ -225,3 +226,17 @@ class TestUnivariateGaussian(unittest.TestCase):
         for i in range(len(pySols)):
           diff += abs(derivLogCdfNormal(self.zz[i]) - pySols[i])
         self.assertTrue(diff  < 1e-8)
+
+class TestStandardize(unittest.TestCase):
+    def setUp(self):
+        self.normalizer = GPy.util.normalizer.Standardize()
+        y = np.stack([np.random.randn(10), 2*np.random.randn(10)], axis=1)
+        self.normalizer.scale_by(y)
+    
+    def test_inverse_covariance(self):
+        """
+        Test inverse covariance outputs correct size
+        """
+        covariance = np.random.rand(100, 100)
+        output = self.normalizer.inverse_covariance(covariance)
+        self.assertTrue(output.shape == (100, 100, 2))
\ No newline at end of file
diff --git a/GPy/util/normalizer.py b/GPy/util/normalizer.py
index b62ac35b..7a3ee020 100644
--- a/GPy/util/normalizer.py
+++ b/GPy/util/normalizer.py
@@ -3,30 +3,46 @@ Created on Aug 27, 2014
 
 @author: Max Zwiessele
 '''
-import logging
 import numpy as np
 
+
 class _Norm(object):
     def __init__(self):
         pass
+
     def scale_by(self, Y):
         """
         Use data matrix Y as normalization space to work in.
         """
         raise NotImplementedError
+
     def normalize(self, Y):
         """
         Project Y into normalized space
         """
         if not self.scaled():
             raise AttributeError("Norm object not initialized yet, try calling scale_by(data) first.")
+
     def inverse_mean(self, X):
         """
         Project the normalized object X into space of Y
         """
         raise NotImplementedError
+
     def inverse_variance(self, var):
         return var
+
+    def inverse_covariance(self, covariance):
+        """
+        Convert scaled covariance to unscaled.
+        Args:
+            covariance - numpy array of shape (n, n)
+        Returns:
+            covariance - numpy array of shape (n, n, m) where m is number of
+                         outputs
+        """
+        raise NotImplementedError
+
     def scaled(self):
         """
         Whether this Norm object has been initialized.
@@ -57,17 +73,25 @@ class _Norm(object):
 class Standardize(_Norm):
     def __init__(self):
         self.mean = None
+
     def scale_by(self, Y):
         Y = np.ma.masked_invalid(Y, copy=False)
         self.mean = Y.mean(0).view(np.ndarray)
         self.std = Y.std(0).view(np.ndarray)
+
     def normalize(self, Y):
         super(Standardize, self).normalize(Y)
         return (Y-self.mean)/self.std
+
     def inverse_mean(self, X):
         return (X*self.std)+self.mean
+
     def inverse_variance(self, var):
         return (var*(self.std**2))
+
+    def inverse_covariance(self, covariance):
+        return (covariance[..., np.newaxis]*(self.std**2))
+
     def scaled(self):
         return self.mean is not None
 
@@ -87,29 +111,3 @@ class Standardize(_Norm):
         if "std" in input_dict:
             s.std = np.array(input_dict["std"])
         return s
-
-# Inverse variance to be implemented, disabling for now
-# If someone in the future want to implement this,
-# we need to implement the inverse variance for
-# normalization. This means, we need to know the factor
-# for the variance to be multiplied to the variance in
-# normalized space. This is easy to compute for standardization
-# (see above) but gets tricky here.
-# class Normalize(_Norm):
-#     def __init__(self):
-#         self.ymin = None
-#         self.ymax = None
-#     def scale_by(self, Y):
-#         Y = np.ma.masked_invalid(Y, copy=False)
-#         self.ymin = Y.min(0).view(np.ndarray)
-#         self.ymax = Y.max(0).view(np.ndarray)
-#     def normalize(self, Y):
-#         super(Normalize, self).normalize(Y)
-#         return (Y - self.ymin) / (self.ymax - self.ymin) - .5
-#     def inverse_mean(self, X):
-#         return (X + .5) * (self.ymax - self.ymin) + self.ymin
-#     def inverse_variance(self, var):
-#
-#         return (var*(self.std**2))
-#     def scaled(self):
-#         return (self.ymin is not None) and (self.ymax is not None)
diff --git a/README.md b/README.md
index ffbf6a34..f03213b0 100644
--- a/README.md
+++ b/README.md
@@ -5,7 +5,7 @@ The Gaussian processes framework in Python.
 * GPy [homepage](http://sheffieldml.github.io/GPy/)
 * Tutorial [notebooks](http://nbviewer.ipython.org/github/SheffieldML/notebook/blob/master/GPy/index.ipynb)
 * User [mailing-list](https://lists.shef.ac.uk/sympa/subscribe/gpy-users)
-* Developer [documentation](http://pythonhosted.org/GPy/)
+* Developer [documentation](http://gpy.readthedocs.io/)
 * Travis-CI [unit-tests](https://travis-ci.org/SheffieldML/GPy)
 * [![licence](https://img.shields.io/badge/licence-BSD-blue.svg)](http://opensource.org/licenses/BSD-3-Clause)
 
diff --git a/appveyor.yml b/appveyor.yml
index dd315e9d..b736d6b4 100644
--- a/appveyor.yml
+++ b/appveyor.yml
@@ -3,7 +3,7 @@ environment:
     secure: 8/ZjXFwtd1S7ixd7PJOpptupKKEDhm2da/q3unabJ00=
   COVERALLS_REPO_TOKEN:
     secure: d3Luic/ESkGaWnZrvWZTKrzO+xaVwJWaRCEP0F+K/9DQGPSRZsJ/Du5g3s4XF+tS
-  gpy_version: 1.8.5
+  gpy_version: 1.9.2
   matrix:
     - PYTHON_VERSION: 2.7
       MINICONDA: C:\Miniconda-x64
diff --git a/doc/source/conf.py b/doc/source/conf.py
index 1f9c98b6..b968fe48 100644
--- a/doc/source/conf.py
+++ b/doc/source/conf.py
@@ -22,32 +22,52 @@ import shlex
 #for p in os.walk('../../GPy'):
 #    sys.path.append(p[0])
 sys.path.insert(0, os.path.abspath('../../'))
-#sys.path.insert(0, os.path.abspath('../../GPy/'))
+sys.path.insert(0, os.path.abspath('../../GPy/'))
 
 on_rtd = os.environ.get('READTHEDOCS', None) == 'True'
 
+import sys
+from unittest.mock import MagicMock
+
+class Mock(MagicMock):
+    @classmethod
+    def __getattr__(cls, name):
+        return MagicMock()
+
+MOCK_MODULES = [
+    "GPy.util.linalg.linalg_cython",
+    "GPy.util.linalg_cython",
+    "sympy",
+    'GPy.kern.stationary_cython',
+    "sympy.utilities",
+    "sympy.utilities.lambdify",
+]
+
+sys.modules.update((mod_name, Mock()) for mod_name in MOCK_MODULES)
+
 #on_rtd = True
 if on_rtd:
     # sys.path.append(os.path.abspath('../GPy'))
 
     import subprocess
 
+    # build extensions:
+    # proc = subprocess.Popen("cd ../../; python setup.py build_ext install", stdout=subprocess.PIPE, shell=True)
+    # (out, err) = proc.communicate()
+    # print("build_ext develop:")
+    # print(out)
+
+    # print current folder:
     proc = subprocess.Popen("pwd", stdout=subprocess.PIPE, shell=True)
     (out, err) = proc.communicate()
-    print "program output:", out
-    proc = subprocess.Popen("ls ../../", stdout=subprocess.PIPE, shell=True)
-    (out, err) = proc.communicate()
-    print "program output:", out
+    print("$ pwd: ")
+    print(out)
+
     #Lets regenerate our rst files from the source, -P adds private modules (i.e kern._src)
     proc = subprocess.Popen("sphinx-apidoc -P -f -o . ../../GPy", stdout=subprocess.PIPE, shell=True)
     (out, err) = proc.communicate()
-    print "program output:", out
-    #proc = subprocess.Popen("whereis numpy", stdout=subprocess.PIPE, shell=True)
-    #(out, err) = proc.communicate()
-    #print "program output:", out
-    #proc = subprocess.Popen("whereis matplotlib", stdout=subprocess.PIPE, shell=True)
-    #(out, err) = proc.communicate()
-    #print "program output:", out
+    print("$ Apidoc:")
+    print(out)
 
 
 # -- General configuration ------------------------------------------------
@@ -77,15 +97,6 @@ extensions = [
 #    def __getattr__(cls, name):
 #            return Mock()
 #
-MOCK_MODULES = ['scipy.linalg.blas', 'blas', 'scipy.optimize', 'scipy.optimize.linesearch', 'scipy.linalg',
-                'scipy', 'scipy.special', 'scipy.integrate', 'scipy.io', 'scipy.stats',
-                'sympy', 'sympy.utilities.iterables', 'sympy.utilities.lambdify',
-                'sympy.utilities', 'sympy.utilities.codegen', 'sympy.core.cache',
-                'sympy.core', 'sympy.parsing', 'sympy.parsing.sympy_parser',
-                'nose', 'nose.tools'
-                ]
-
-autodoc_mock_imports = MOCK_MODULES
 #
 #sys.modules.update((mod_name, Mock()) for mod_name in MOCK_MODULES)
 #
@@ -97,6 +108,7 @@ autodoc_default_flags = ['members',
                          #'special-members',
                          #'inherited-members',
                          'show-inheritance']
+
 autodoc_member_order = 'groupwise'
 add_function_parentheses = False
 add_module_names = False
@@ -144,7 +156,21 @@ print version
 #
 # This is also used if you do content translation via gettext catalogs.
 # Usually you set "language" from the command line for these cases.
-language = None
+language = 'python'
+
+# autodoc:
+autoclass_content = 'both'
+autodoc_default_flags = ['members',
+                         #'undoc-members',
+                         #'private-members',
+                         #'special-members',
+                         #'inherited-members',
+                         'show-inheritance']
+autodoc_member_order = 'groupwise'
+add_function_parentheses = False
+add_module_names = False
+modindex_common_prefix = ['paramz']
+show_authors = True
 
 # There are two options for replacing |today|: either, you set today to some
 # non-false value, then it is used:
@@ -172,7 +198,7 @@ exclude_patterns = []
 #show_authors = False
 
 # The name of the Pygments (syntax highlighting) style to use.
-#pygments_style = 'sphinx'
+pygments_style = 'sphinx'
 
 # A list of ignored prefixes for module index sorting.
 #modindex_common_prefix = []
@@ -217,7 +243,7 @@ html_theme = 'sphinx_rtd_theme'
 # Add any paths that contain custom static files (such as style sheets) here,
 # relative to this directory. They are copied after the builtin static files,
 # so a file named "default.css" will overwrite the builtin "default.css".
-#html_static_path = ['_static']
+html_static_path = ['_static']
 
 # Add any extra paths that contain custom files (such as robots.txt or
 # .htaccess) here, relative to this directory. These files are copied
@@ -242,16 +268,16 @@ html_theme = 'sphinx_rtd_theme'
 #html_additional_pages = {}
 
 # If false, no module index is generated.
-#html_domain_indices = False
+html_domain_indices = False
 
 # If false, no index is generated.
-#html_use_index = False
+html_use_index = False
 
 # If true, the index is split into individual pages for each letter.
 html_split_index = True
 
 # If true, links to the reST sources are added to the pages.
-#html_show_sourcelink = True
+html_show_sourcelink = True
 
 # If true, "Created using Sphinx" is shown in the HTML footer. Default is True.
 #html_show_sphinx = True
@@ -286,9 +312,9 @@ htmlhelp_basename = 'GPydoc'
 
 # -- Options for LaTeX output ---------------------------------------------
 
-#latex_elements = {
+latex_elements = {
 # The paper size ('letterpaper' or 'a4paper').
-#'papersize': 'letterpaper',
+'papersize': 'a4paper',
 
 # The font size ('10pt', '11pt' or '12pt').
 #'pointsize': '10pt',
@@ -297,8 +323,8 @@ htmlhelp_basename = 'GPydoc'
 #'preamble': '',
 
 # Latex figure (float) alignment
-#'figure_align': 'htbp',
-#}
+'figure_align': 'htbp',
+}
 
 # Grouping the document tree into LaTeX files. List of tuples
 # (source start file, target name, title,
diff --git a/doc/source/requirements.txt b/doc/source/requirements.txt
index dd3ba36f..5ae1e857 100644
--- a/doc/source/requirements.txt
+++ b/doc/source/requirements.txt
@@ -1 +1,10 @@
-paramz
\ No newline at end of file
+numpy
+scipy
+six
+decorator
+matplotlib
+paramz
+cython
+mock
+sympy
+nose
\ No newline at end of file
diff --git a/setup.cfg b/setup.cfg
index 21160939..0f1d4075 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 1.8.5
+current_version = 1.9.2
 tag = True
 commit = True
 
@@ -12,3 +12,4 @@ upload-dir = doc/build/html
 
 [medatdata]
 description-file = README.rst
+
diff --git a/setup.py b/setup.py
index 7d3a5355..5e4357b5 100644
--- a/setup.py
+++ b/setup.py
@@ -94,15 +94,18 @@ ext_mods = [Extension(name='GPy.kern.src.stationary_cython',
             Extension(name='GPy.util.linalg_cython',
                       sources=['GPy/util/linalg_cython.c'],
                       include_dirs=[np.get_include(),'.'],
-                      extra_compile_args=compile_flags),
+                      extra_compile_args=compile_flags,
+                      extra_link_args = link_args),
             Extension(name='GPy.kern.src.coregionalize_cython',
                       sources=['GPy/kern/src/coregionalize_cython.c'],
                       include_dirs=[np.get_include(),'.'],
-                      extra_compile_args=compile_flags),
+                      extra_compile_args=compile_flags,
+                      extra_link_args = link_args),
             Extension(name='GPy.models.state_space_cython',
                       sources=['GPy/models/state_space_cython.c'],
                       include_dirs=[np.get_include(),'.'],
-                      extra_compile_args=compile_flags)]
+                      extra_compile_args=compile_flags,
+                      extra_link_args = link_args)]
 
 setup(name = 'GPy',
       version = __version__,
@@ -150,7 +153,7 @@ setup(name = 'GPy',
       py_modules = ['GPy.__init__'],
       test_suite = 'GPy.testing',
       setup_requires = ['numpy>=1.7'],
-      install_requires = ['numpy>=1.7', 'scipy>=0.16', 'six', 'paramz>=0.8.5'],
+      install_requires = ['numpy>=1.7', 'scipy>=0.16', 'six', 'paramz>=0.9.0'],
       extras_require = {'docs':['sphinx'],
                         'optional':['mpi4py',
                                     'ipython>=4.0.0',