diff --git a/.gitignore b/.gitignore
index d494630f..2b6a92cd 100644
--- a/.gitignore
+++ b/.gitignore
@@ -45,4 +45,7 @@ iterate.dat
 
 # git merge files #
 ###################
-*.orig
\ No newline at end of file
+*.orig
+
+# pycharm IDE stuff
+.idea/
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 167b4ddc..88284db3 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,105 @@
 # Changelog
 
+## v1.6.1 (2017-02-28)
+
+### Fix
+
+* Beiwang will add GMM in full. [mzwiessele]
+
+### Other
+
+* Bump version: 1.6.0 → 1.6.1. [mzwiessele]
+
+
+## v1.6.0 (2017-02-28)
+
+### Fix
+
+* Kernel tests and variational tests. [mzwiessele]
+
+* Plotting tests for new matplotlib. [mzwiessele]
+
+* Model tests numpy integer error. [mzwiessele]
+
+* Replot with new matplotlib. [mzwiessele]
+
+* Offline plotting workaround with squeezing arrays. [mzwiessele]
+
+* Fixed numpy 1.12 indexing and shape preservation. [mzwiessele]
+
+### Other
+
+* Bump version: 1.5.9 → 1.6.0. [mzwiessele]
+
+* Merge branch 'devel' into alexfeld-offline_plotly. [mzwiessele]
+
+* Merge branch 'devel' into alexfeld-offline_plotly. [mzwiessele]
+
+* Merge branch 'offline_plotly' of git://github.com/alexfeld/GPy into alexfeld-offline_plotly. [mzwiessele]
+
+* Provide two classes for plotly plots to remove global variable. [Alex Feldstein]
+
+* Add offline plotting for plotly. [Alex Feldstein]
+
+
+## v1.5.9 (2017-02-23)
+
+### Other
+
+* Bump version: 1.5.8 → 1.5.9. [mzwiessele]
+
+* Merge remote-tracking branch 'origin/deploy' into devel. [mzwiessele]
+
+* Merge pull request #455 from SheffieldML/devel. [Max Zwiessele]
+
+  1.5.6
+
+
+## v1.5.8 (2017-02-23)
+
+### Fix
+
+* Predictive_gradients for new posterior class. [mzwiessele]
+
+* Removed additional dict line. [mzwiessele]
+
+* Plotting also allows 3D (capitals) [mzwiessele]
+
+* Fallback for when no environment variables are set (#467) [Safrone]
+
+  * fix: dev: add or in home directory getting
+
+  adds another or when getting the home directory with os.getenv() so that if neither $HOME nor $USERPROFILE environment variable is set, os.path.join() will not fail by getting a None and the config will revert to the default configuration file.
+
+  * fix: remove extra statement
+
+### Other
+
+* Bump version: 1.5.7 → 1.5.8. [mzwiessele]
+
+* Update ss_gplvm.py. [Zhenwen Dai]
+
+  resolve the future warning: FutureWarning:comparison to `None` will result in an elementwise object comparison in the future.
+
+* Merge pull request #472 from SheffieldML/predictive_gradients. [Max Zwiessele]
+
+  fix: predictive_gradients for new posterior class
+
+* Merge pull request #470 from SheffieldML/plotting_fix. [Max Zwiessele]
+
+  Plotting fix
+
+* Bump version: 1.5.6 → 1.5.7. [mzwiessele]
+
+* Changed the order of the operations, ensuring that the covariance matrix is symmetric despite numerical precision issues. Suggested by Alan. [Teo de Campos]
+
+* Delete gmm_bayesian_gplvm.py. [beiwang]
+
+* Gmm_creation. [beiwang]
+
+* Gmm_creation. [beiwang]
+
+
 ## v1.5.6 (2016-11-07)
 
 ### New
diff --git a/GPy/__version__.py b/GPy/__version__.py
index d0113b35..f49459c7 100644
--- a/GPy/__version__.py
+++ b/GPy/__version__.py
@@ -1 +1 @@
-__version__ = "1.5.6"
+__version__ = "1.6.1"
diff --git a/GPy/core/gp.py b/GPy/core/gp.py
index f0910d9d..b90c95c1 100644
--- a/GPy/core/gp.py
+++ b/GPy/core/gp.py
@@ -339,9 +339,17 @@ class GP(Model):
         # 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}
-        alpha = -2.*np.dot(kern.K(Xnew, self._predictive_variable), self.posterior.woodbury_inv)
-        dv_dX += kern.gradients_X(alpha, Xnew, self._predictive_variable)
-        return mean_jac, dv_dX
+        if self.posterior.woodbury_inv.ndim == 3:
+            tmp = np.empty(dv_dX.shape + (self.posterior.woodbury_inv.shape[2],))
+            tmp[:] = 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)
+        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
 
     def predict_jacobian(self, Xnew, kern=None, full_cov=False):
         """
diff --git a/GPy/core/gp_grid.py b/GPy/core/gp_grid.py
index bdb1b614..3cd7d8cb 100644
--- a/GPy/core/gp_grid.py
+++ b/GPy/core/gp_grid.py
@@ -68,12 +68,12 @@ class GpGrid(GP):
         for b in (B.T):
             x = b
             N = 1
-            G = np.zeros(D)
+            G = np.zeros(D, dtype=np.int_)
             for d in range(D):
                 G[d] = len(A[d])
             N = np.prod(G)
             for d in range(D-1, -1, -1):
-                X = np.reshape(x, (G[d], np.round(N/G[d])), order='F')
+                X = np.reshape(x, (G[d], int(np.round(N/G[d]))), order='F')
                 Z = np.dot(A[d], X)
                 Z = Z.T
                 x = np.reshape(Z, (-1, 1), order='F')
diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py
index 81e1b773..fe84df72 100644
--- a/GPy/examples/dimensionality_reduction.py
+++ b/GPy/examples/dimensionality_reduction.py
@@ -260,7 +260,7 @@ def _simulate_sincos(D1, D2, D3, N, num_inducing, plot_sim=False):
 
     x = _np.linspace(0, 4 * _np.pi, N)[:, None]
     s1 = _np.vectorize(lambda x: _np.sin(x))
-    s2 = _np.vectorize(lambda x: _np.cos(x) ** 2)
+    s2 = _np.vectorize(lambda x: _np.cos(x))
     s3 = _np.vectorize(lambda x:-_np.exp(-_np.cos(2 * x)))
     sS = _np.vectorize(lambda x: _np.cos(x))
 
@@ -302,8 +302,8 @@ def _simulate_sincos(D1, D2, D3, N, num_inducing, plot_sim=False):
 
 def _generate_high_dimensional_output(D1, D2, D3, s1, s2, s3, sS):
     S1 = _np.hstack([s1, sS])
-    S2 = _np.hstack([s2, s3, sS])
-    S3 = _np.hstack([s3, sS])
+    S2 = _np.hstack([sS])
+    S3 = _np.hstack([s1, s3, sS])
     Y1 = S1.dot(_np.random.randn(S1.shape[1], D1))
     Y2 = S2.dot(_np.random.randn(S2.shape[1], D2))
     Y3 = S3.dot(_np.random.randn(S3.shape[1], D3))
diff --git a/GPy/inference/latent_function_inference/gaussian_grid_inference.py b/GPy/inference/latent_function_inference/gaussian_grid_inference.py
index abdd78a3..aec52a46 100644
--- a/GPy/inference/latent_function_inference/gaussian_grid_inference.py
+++ b/GPy/inference/latent_function_inference/gaussian_grid_inference.py
@@ -36,12 +36,12 @@ class GaussianGridInference(LatentFunctionInference):
         x = b
         N = 1
         D = len(A)
-        G = np.zeros((D,1))
+        G = np.zeros((D), dtype=np.int_)
         for d in range(0, D):
             G[d] = len(A[d])
         N = np.prod(G)
         for d in range(D-1, -1, -1):
-            X = np.reshape(x, (G[d], np.round(N/G[d])), order='F')
+            X = np.reshape(x, (G[d], int(np.round(N/G[d]))), order='F')
             Z = np.dot(A[d], X)
             Z = Z.T
             x = np.reshape(Z, (-1, 1), order='F')
diff --git a/GPy/kern/src/kern.py b/GPy/kern/src/kern.py
index 931591fd..198ebd13 100644
--- a/GPy/kern/src/kern.py
+++ b/GPy/kern/src/kern.py
@@ -46,11 +46,11 @@ class Kern(Parameterized):
         self.input_dim = int(input_dim)
 
         if active_dims is None:
-            active_dims = np.arange(input_dim)
+            active_dims = np.arange(input_dim, dtype=np.int_)
 
         self.active_dims = np.atleast_1d(np.asarray(active_dims, np.int_))
 
-        self._all_dims_active = np.atleast_1d(self.active_dims).astype(int)
+        self._all_dims_active = np.atleast_1d(self.active_dims).astype(np.int_)
 
         assert self.active_dims.size == self.input_dim, "input_dim={} does not match len(active_dim)={}".format(self.input_dim, self._all_dims_active.size)
 
diff --git a/GPy/kern/src/standard_periodic.py b/GPy/kern/src/standard_periodic.py
index d45e8d16..974bce6a 100644
--- a/GPy/kern/src/standard_periodic.py
+++ b/GPy/kern/src/standard_periodic.py
@@ -55,7 +55,6 @@ class StdPeriodic(Kern):
 
     def __init__(self, input_dim, variance=1., period=None, lengthscale=None, ARD1=False, ARD2=False, active_dims=None, name='std_periodic',useGPU=False):
         super(StdPeriodic, self).__init__(input_dim, active_dims, name, useGPU=useGPU)
-        self.input_dim = input_dim
         self.ARD1 = ARD1 # correspond to periods
         self.ARD2 = ARD2 # correspond to lengthscales
 
@@ -66,7 +65,7 @@ class StdPeriodic(Kern):
                 period = np.asarray(period)
                 assert period.size == 1, "Only one period needed for non-ARD kernel"
             else:
-                period = np.ones(1.0)
+                period = np.ones(1)
         else:
             if period is not None:
                 period = np.asarray(period)
diff --git a/GPy/kern/src/stationary.py b/GPy/kern/src/stationary.py
index c10384bb..79d4b954 100644
--- a/GPy/kern/src/stationary.py
+++ b/GPy/kern/src/stationary.py
@@ -135,7 +135,7 @@ class Stationary(Kern):
             #X2, = self._slice_X(X2)
             X1sq = np.sum(np.square(X),1)
             X2sq = np.sum(np.square(X2),1)
-            r2 = -2.*np.dot(X, X2.T) + X1sq[:,None] + X2sq[None,:]
+            r2 = -2.*np.dot(X, X2.T) + (X1sq[:,None] + X2sq[None,:])
             r2 = np.clip(r2, 0, np.inf)
             return np.sqrt(r2)
 
diff --git a/GPy/likelihoods/binomial.py b/GPy/likelihoods/binomial.py
index 306f0a1a..e63c009a 100644
--- a/GPy/likelihoods/binomial.py
+++ b/GPy/likelihoods/binomial.py
@@ -63,7 +63,8 @@ class Binomial(Likelihood):
         :rtype: float
         """
         N = Y_metadata['trials']
-        assert N.shape == y.shape
+        np.testing.assert_array_equal(N.shape, y.shape)
+
         nchoosey = special.gammaln(N+1) - special.gammaln(y+1) - special.gammaln(N-y+1)
         return nchoosey + y*np.log(inv_link_f) + (N-y)*np.log(1.-inv_link_f)
 
@@ -83,7 +84,8 @@ class Binomial(Likelihood):
         :rtype: Nx1 array
         """
         N = Y_metadata['trials']
-        assert N.shape == y.shape
+        np.testing.assert_array_equal(N.shape, y.shape)
+
         return y/inv_link_f - (N-y)/(1.-inv_link_f)
 
     def d2logpdf_dlink2(self, inv_link_f, y, Y_metadata=None):
@@ -108,7 +110,7 @@ class Binomial(Likelihood):
             (the distribution for y_i depends only on inverse link of f_i not on inverse link of f_(j!=i)
         """
         N = Y_metadata['trials']
-        assert N.shape == y.shape
+        np.testing.assert_array_equal(N.shape, y.shape)
         return -y/np.square(inv_link_f) - (N-y)/np.square(1.-inv_link_f)
 
     def d3logpdf_dlink3(self, inv_link_f, y, Y_metadata=None):
@@ -131,7 +133,8 @@ class Binomial(Likelihood):
             (the distribution for y_i depends only on inverse link of f_i not on inverse link of f_(j!=i)
         """
         N = Y_metadata['trials']
-        assert N.shape == y.shape
+        np.testing.assert_array_equal(N.shape, y.shape)
+
         inv_link_f2 = np.square(inv_link_f)
         return 2*y/inv_link_f**3 - 2*(N-y)/(1.-inv_link_f)**3
 
diff --git a/GPy/models/ss_gplvm.py b/GPy/models/ss_gplvm.py
index 7d10bab6..95f58cad 100644
--- a/GPy/models/ss_gplvm.py
+++ b/GPy/models/ss_gplvm.py
@@ -193,7 +193,7 @@ class SSGPLVM(SparseGP_MPI):
         self.init = init
         self.sharedX = sharedX
         
-        if X == None:
+        if X is None:
             from ..util.initialization import initialize_latent
             X, fracs = initialize_latent(init, input_dim, Y)
         else:
diff --git a/GPy/models/state_space.py b/GPy/models/state_space.py
index 0e3498d8..3ef1023f 100644
--- a/GPy/models/state_space.py
+++ b/GPy/models/state_space.py
@@ -293,13 +293,13 @@ class StateSpace(Model):
 
             # Update step (only if there is data)
             if not np.isnan(Y[:,k]):
-                 if Y.shape[0]==1:
-                     K = PF[:,:,k].dot(H.T)/(H.dot(PF[:,:,k]).dot(H.T) + R)
-                 else:
-                     LL = linalg.cho_factor(H.dot(PF[:,:,k]).dot(H.T) + R)
-                     K = linalg.cho_solve(LL, H.dot(PF[:,:,k].T)).T
-                 MF[:,k] += K.dot(Y[:,k]-H.dot(MF[:,k]))
-                 PF[:,:,k] -= K.dot(H).dot(PF[:,:,k])
+                if Y.shape[0]==1:
+                    K = PF[:,:,k].dot(H.T)/(H.dot(PF[:,:,k]).dot(H.T) + R)
+                else:
+                    LL = linalg.cho_factor(H.dot(PF[:,:,k]).dot(H.T) + R)
+                    K = linalg.cho_solve(LL, H.dot(PF[:,:,k].T)).T
+                MF[:,k] += K.dot(Y[:,k]-H.dot(MF[:,k]))
+                PF[:,:,k] -= K.dot(H).dot(PF[:,:,k])
 
         # Return values
         return (MF, PF)
diff --git a/GPy/models/state_space_main.py b/GPy/models/state_space_main.py
index d0406e96..65763a05 100644
--- a/GPy/models/state_space_main.py
+++ b/GPy/models/state_space_main.py
@@ -215,7 +215,7 @@ class R_handling_Python(Measurement_Callables_Class):
             inv_R_square_root(k)
         """
         self.R = R
-        self.index = index
+        self.index = np.asarray(index, np.int_)
         self.R_time_var_index = int(R_time_var_index)
         self.dR = dR
 
@@ -350,7 +350,7 @@ class Q_handling_Python(Dynamic_Callables_Class):
         """
 
         self.Q = Q
-        self.index = index
+        self.index = np.asarray(index, np.int_)
         self.Q_time_var_index = Q_time_var_index
         self.dQ = dQ
 
@@ -427,7 +427,7 @@ class Std_Dynamic_Callables_Python(Q_handling_Class):
               self).__init__(Q, index, Q_time_var_index, unique_Q_number, dQ)
 
         self.A = A
-        self.A_time_var_index = A_time_var_index
+        self.A_time_var_index = np.asarray(A_time_var_index, np.int_)
         self.dA = dA
 
     def f_a(self, k, m, A):
diff --git a/GPy/plotting/__init__.py b/GPy/plotting/__init__.py
index 359a841a..ad62a00f 100644
--- a/GPy/plotting/__init__.py
+++ b/GPy/plotting/__init__.py
@@ -2,10 +2,10 @@
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 current_lib = [None]
 
-supported_libraries = ['matplotlib', 'plotly', 'none']
+supported_libraries = ['matplotlib', 'plotly', 'plotly_online', 'plotly_offline', 'none']
 error_suggestion = "Please make sure you specify your plotting library in your configuration file (<User>/.config/GPy/user.cfg).\n\n[plotting]\nlibrary = <library>\n\nCurrently supported libraries: {}".format(", ".join(supported_libraries))
 
-def change_plotting_library(lib):
+def change_plotting_library(lib, **kwargs):
     try:
         #===========================================================================
         # Load in your plotting library here and
@@ -19,10 +19,14 @@ def change_plotting_library(lib):
             from .matplot_dep.plot_definitions import MatplotlibPlots
             from .matplot_dep import visualize, mapping_plots, priors_plots, ssgplvm, svig_plots, variational_plots, img_plots
             current_lib[0] = MatplotlibPlots()
-        if lib == 'plotly':
+        if lib in ['plotly', 'plotly_online']:
             import plotly
-            from .plotly_dep.plot_definitions import PlotlyPlots
-            current_lib[0] = PlotlyPlots()
+            from .plotly_dep.plot_definitions import PlotlyPlotsOnline
+            current_lib[0] = PlotlyPlotsOnline(**kwargs)
+        if lib == 'plotly_offline':
+            import plotly
+            from .plotly_dep.plot_definitions import PlotlyPlotsOffline
+            current_lib[0] = PlotlyPlotsOffline(**kwargs)
         if lib == 'none':
             current_lib[0] = None
         inject_plotting()
diff --git a/GPy/plotting/gpy_plot/gp_plots.py b/GPy/plotting/gpy_plot/gp_plots.py
index 702aeb7b..230d47f0 100644
--- a/GPy/plotting/gpy_plot/gp_plots.py
+++ b/GPy/plotting/gpy_plot/gp_plots.py
@@ -88,12 +88,12 @@ def _plot_mean(self, canvas, helper_data, helper_prediction,
             update_not_existing_kwargs(kwargs, pl().defaults.meanplot_1d)  # @UndefinedVariable
             plots = dict(gpmean=[pl().plot(canvas, Xgrid[:, free_dims], mu, label=label, **kwargs)])
         else:
-            if projection == '2d':
+            if projection.lower() in '2d':
                 update_not_existing_kwargs(kwargs, pl().defaults.meanplot_2d)  # @UndefinedVariable
                 plots = dict(gpmean=[pl().contour(canvas, x[:,0], y[0,:],
                                                mu.reshape(resolution, resolution).T,
                                                levels=levels, label=label, **kwargs)])
-            elif projection == '3d':
+            elif projection.lower() in '3d':
                 update_not_existing_kwargs(kwargs, pl().defaults.meanplot_3d)  # @UndefinedVariable
                 plots = dict(gpmean=[pl().surface(canvas, x, y,
                                                mu.reshape(resolution, resolution),
diff --git a/GPy/plotting/matplot_dep/plot_definitions.py b/GPy/plotting/matplot_dep/plot_definitions.py
index 0e3bc32d..9aa72583 100644
--- a/GPy/plotting/matplot_dep/plot_definitions.py
+++ b/GPy/plotting/matplot_dep/plot_definitions.py
@@ -92,7 +92,7 @@ class MatplotlibPlots(AbstractPlottingLibrary):
         if title is not None: ax.figure.suptitle(title)
         return ax
 
-    def show_canvas(self, ax):
+    def show_canvas(self, ax, **kwargs):
         ax.figure.canvas.draw()
         return ax.figure
 
diff --git a/GPy/plotting/plotly_dep/plot_definitions.py b/GPy/plotting/plotly_dep/plot_definitions.py
index 9e021fd8..9b659551 100644
--- a/GPy/plotting/plotly_dep/plot_definitions.py
+++ b/GPy/plotting/plotly_dep/plot_definitions.py
@@ -31,8 +31,8 @@ import numpy as np
 from ..abstract_plotting_library import AbstractPlottingLibrary
 from .. import Tango
 from . import defaults
+import plotly
 from plotly import tools
-from plotly import plotly as py
 from plotly.graph_objs import Scatter, Scatter3d, Line,\
     Marker, ErrorX, ErrorY, Bar, Heatmap, Trace,\
     Annotations, Annotation, Contour, Font, Surface
@@ -52,9 +52,9 @@ SYMBOL_MAP = {
     'd': 'diamond',
 }
 
-class PlotlyPlots(AbstractPlottingLibrary):
+class PlotlyPlotsBase(AbstractPlottingLibrary):
     def __init__(self):
-        super(PlotlyPlots, self).__init__()
+        super(PlotlyPlotsBase, self).__init__()
         self._defaults = defaults.__dict__
         self.current_states = dict()
 
@@ -96,7 +96,10 @@ class PlotlyPlots(AbstractPlottingLibrary):
                 xref, yref = figure._grid_ref[row-1][col-1]
                 for a in traces:
                     append_annotation(a, xref, yref)
-            elif isinstance(traces, (Trace)):
+            # elif isinstance(traces, (Trace)):  # doesn't work
+            # elif type(traces) in [v for k,v in go.__dict__.iteritems()]:
+            elif isinstance(traces, (Scatter, Scatter3d, ErrorX,
+                        ErrorY, Bar, Heatmap, Trace, Contour, Surface)):
                 try:
                     append_trace(traces, row, col)
                 except PlotlyDictKeyError:
@@ -114,15 +117,7 @@ class PlotlyPlots(AbstractPlottingLibrary):
         return canvas
 
     def show_canvas(self, canvas, filename=None, **kwargs):
-        figure, _, _ = canvas
-        if len(figure.data) == 0:
-            # add mock data
-            figure.append_trace(Scatter(x=[], y=[], name='', showlegend=False), 1, 1)
-        from ..gpy_plot.plot_util import in_ipynb
-        if in_ipynb():
-            return py.iplot(figure, filename=filename)#self.current_states[hex(id(figure))]['filename'])
-        else:
-            return py.plot(figure, filename=filename)#self.current_states[hex(id(figure))]['filename'])
+        return NotImplementedError
 
     def scatter(self, ax, X, Y, Z=None, color=Tango.colorsHex['mediumBlue'], cmap=None, label=None, marker='o', marker_kwargs=None, **kwargs):
         try:
@@ -133,7 +128,9 @@ class PlotlyPlots(AbstractPlottingLibrary):
         marker_kwargs = marker_kwargs or {}
         if 'symbol' not in marker_kwargs:
             marker_kwargs['symbol'] = marker
+        X, Y = np.squeeze(X), np.squeeze(Y)
         if Z is not None:
+            Z = np.squeeze(Z)
             return Scatter3d(x=X, y=Y, z=Z, mode='markers',
                              showlegend=label is not None,
                              marker=Marker(color=color, colorscale=cmap, **marker_kwargs),
@@ -145,7 +142,9 @@ class PlotlyPlots(AbstractPlottingLibrary):
     def plot(self, ax, X, Y, Z=None, color=None, label=None, line_kwargs=None, **kwargs):
         if 'mode' not in kwargs:
             kwargs['mode'] = 'lines'
+        X, Y = np.squeeze(X), np.squeeze(Y)
         if Z is not None:
+            Z = np.squeeze(Z)
             return Scatter3d(x=X, y=Y, z=Z, showlegend=label is not None, line=Line(color=color, **line_kwargs or {}), name=label, **kwargs)
         return Scatter(x=X, y=Y, showlegend=label is not None, line=Line(color=color, **line_kwargs or {}), name=label, **kwargs)
 
@@ -191,7 +190,9 @@ class PlotlyPlots(AbstractPlottingLibrary):
             error_kwargs.update(dict(array=error[1], arrayminus=error[0], symmetric=False))
         else:
             error_kwargs.update(dict(array=error, symmetric=True))
+        X, Y = np.squeeze(X), np.squeeze(Y)
         if Z is not None:
+            Z = np.squeeze(Z)
             return Scatter3d(x=X, y=Y, z=Z, mode='markers',
                              error_x=ErrorX(color=color, **error_kwargs or {}),
                              marker=Marker(size='0'), name=label,
@@ -208,7 +209,9 @@ class PlotlyPlots(AbstractPlottingLibrary):
             error_kwargs.update(dict(array=error[1], arrayminus=error[0], symmetric=False))
         else:
             error_kwargs.update(dict(array=error, symmetric=True))
+        X, Y = np.squeeze(X), np.squeeze(Y)
         if Z is not None:
+            Z = np.squeeze(Z)
             return Scatter3d(x=X, y=Y, z=Z, mode='markers',
                              error_y=ErrorY(color=color, **error_kwargs or {}),
                              marker=Marker(size='0'), name=label,
@@ -232,7 +235,7 @@ class PlotlyPlots(AbstractPlottingLibrary):
 
     def imshow_interact(self, ax, plot_function, extent=None, label=None, resolution=None, vmin=None, vmax=None, **imshow_kwargs):
         # TODO stream interaction?
-        super(PlotlyPlots, self).imshow_interact(ax, plot_function)
+        super(PlotlyPlotsBase, self).imshow_interact(ax, plot_function)
 
     def annotation_heatmap(self, ax, X, annotation, extent=None, label='Gradient', imshow_kwargs=None, **annotation_kwargs):
         imshow_kwargs.setdefault('label', label)
@@ -259,7 +262,7 @@ class PlotlyPlots(AbstractPlottingLibrary):
         return imshow, annotations
 
     def annotation_heatmap_interact(self, ax, plot_function, extent, label=None, resolution=15, imshow_kwargs=None, **annotation_kwargs):
-        super(PlotlyPlots, self).annotation_heatmap_interact(ax, plot_function, extent)
+        super(PlotlyPlotsBase, self).annotation_heatmap_interact(ax, plot_function, extent)
 
     def contour(self, ax, X, Y, C, levels=20, label=None, **kwargs):
         return Contour(x=X, y=Y, z=C,
@@ -312,3 +315,35 @@ class PlotlyPlots(AbstractPlottingLibrary):
                                        name=None, line=Line(width=1, smoothing=0, color=fcolor), mode='none', fill='tonextx',
                                        legendgroup='density', hoverinfo='none', **kwargs))
         return polycol
+
+
+class PlotlyPlotsOnline(PlotlyPlotsBase):
+    def __init__(self):
+        super(PlotlyPlotsOnline, self).__init__()
+
+    def show_canvas(self, canvas, filename=None, **kwargs):
+        figure, _, _ = canvas
+        if len(figure.data) == 0:
+            # add mock data
+            figure.append_trace(Scatter(x=[], y=[], name='', showlegend=False), 1, 1)
+        from ..gpy_plot.plot_util import in_ipynb
+        if in_ipynb():
+            return plotly.plotly.iplot(figure, filename=filename, **kwargs)
+        else:
+            return plotly.plotly.plot(figure, filename=filename, **kwargs)#self.current_states[hex(id(figure))]['filename'])
+
+class PlotlyPlotsOffline(PlotlyPlotsBase):
+    def __init__(self):
+        super(PlotlyPlotsOffline, self).__init__()
+
+    def show_canvas(self, canvas, filename=None, **kwargs):
+        figure, _, _ = canvas
+        if len(figure.data) == 0:
+            # add mock data
+            figure.append_trace(Scatter(x=[], y=[], name='', showlegend=False), 1, 1)
+        from ..gpy_plot.plot_util import in_ipynb
+        if in_ipynb():
+            plotly.offline.init_notebook_mode(connected=True)
+            return plotly.offline.iplot(figure, filename=filename, **kwargs)#self.current_states[hex(id(figure))]['filename'])
+        else:
+            return plotly.offline.plot(figure, filename=filename, **kwargs)
diff --git a/GPy/testing/baseline/bayesian_gplvm_gradient.npz b/GPy/testing/baseline/bayesian_gplvm_gradient.npz
index b48db887..37c46b7a 100644
Binary files a/GPy/testing/baseline/bayesian_gplvm_gradient.npz and b/GPy/testing/baseline/bayesian_gplvm_gradient.npz differ
diff --git a/GPy/testing/baseline/bayesian_gplvm_inducing.npz b/GPy/testing/baseline/bayesian_gplvm_inducing.npz
index 82e44f54..7ac5c1c9 100644
Binary files a/GPy/testing/baseline/bayesian_gplvm_inducing.npz and b/GPy/testing/baseline/bayesian_gplvm_inducing.npz differ
diff --git a/GPy/testing/baseline/bayesian_gplvm_inducing_3d.npz b/GPy/testing/baseline/bayesian_gplvm_inducing_3d.npz
index e55b4b01..d5a3b391 100644
Binary files a/GPy/testing/baseline/bayesian_gplvm_inducing_3d.npz and b/GPy/testing/baseline/bayesian_gplvm_inducing_3d.npz differ
diff --git a/GPy/testing/baseline/bayesian_gplvm_latent.npz b/GPy/testing/baseline/bayesian_gplvm_latent.npz
index 3e27bc64..710ab07e 100644
Binary files a/GPy/testing/baseline/bayesian_gplvm_latent.npz and b/GPy/testing/baseline/bayesian_gplvm_latent.npz differ
diff --git a/GPy/testing/baseline/bayesian_gplvm_latent_3d.npz b/GPy/testing/baseline/bayesian_gplvm_latent_3d.npz
index b4da4571..0d50befe 100644
Binary files a/GPy/testing/baseline/bayesian_gplvm_latent_3d.npz and b/GPy/testing/baseline/bayesian_gplvm_latent_3d.npz differ
diff --git a/GPy/testing/baseline/bayesian_gplvm_magnification.npz b/GPy/testing/baseline/bayesian_gplvm_magnification.npz
index dd83f22d..411ab821 100644
Binary files a/GPy/testing/baseline/bayesian_gplvm_magnification.npz and b/GPy/testing/baseline/bayesian_gplvm_magnification.npz differ
diff --git a/GPy/testing/baseline/coverage_3d_plot.npz b/GPy/testing/baseline/coverage_3d_plot.npz
index 6372cb73..4df41bde 100644
Binary files a/GPy/testing/baseline/coverage_3d_plot.npz and b/GPy/testing/baseline/coverage_3d_plot.npz differ
diff --git a/GPy/testing/baseline/coverage_annotation_interact.npz b/GPy/testing/baseline/coverage_annotation_interact.npz
index 16fd4e38..ae102c5e 100644
Binary files a/GPy/testing/baseline/coverage_annotation_interact.npz and b/GPy/testing/baseline/coverage_annotation_interact.npz differ
diff --git a/GPy/testing/baseline/coverage_gradient.npz b/GPy/testing/baseline/coverage_gradient.npz
index a55dd97b..54d289e8 100644
Binary files a/GPy/testing/baseline/coverage_gradient.npz and b/GPy/testing/baseline/coverage_gradient.npz differ
diff --git a/GPy/testing/baseline/coverage_imshow_interact.npz b/GPy/testing/baseline/coverage_imshow_interact.npz
index 55dd013c..b270db72 100644
Binary files a/GPy/testing/baseline/coverage_imshow_interact.npz and b/GPy/testing/baseline/coverage_imshow_interact.npz differ
diff --git a/GPy/testing/baseline/gp_2d_data.npz b/GPy/testing/baseline/gp_2d_data.npz
index 0c703f69..6cfde842 100644
Binary files a/GPy/testing/baseline/gp_2d_data.npz and b/GPy/testing/baseline/gp_2d_data.npz differ
diff --git a/GPy/testing/baseline/gp_2d_in_error.npz b/GPy/testing/baseline/gp_2d_in_error.npz
index beb00720..3bbfc02f 100644
Binary files a/GPy/testing/baseline/gp_2d_in_error.npz and b/GPy/testing/baseline/gp_2d_in_error.npz differ
diff --git a/GPy/testing/baseline/gp_2d_inducing.npz b/GPy/testing/baseline/gp_2d_inducing.npz
index 2a249e27..931e2cd4 100644
Binary files a/GPy/testing/baseline/gp_2d_inducing.npz and b/GPy/testing/baseline/gp_2d_inducing.npz differ
diff --git a/GPy/testing/baseline/gp_2d_mean.npz b/GPy/testing/baseline/gp_2d_mean.npz
index dd9ef64d..f4619bfd 100644
Binary files a/GPy/testing/baseline/gp_2d_mean.npz and b/GPy/testing/baseline/gp_2d_mean.npz differ
diff --git a/GPy/testing/baseline/gp_3d_data.npz b/GPy/testing/baseline/gp_3d_data.npz
index e190118a..eae953dd 100644
Binary files a/GPy/testing/baseline/gp_3d_data.npz and b/GPy/testing/baseline/gp_3d_data.npz differ
diff --git a/GPy/testing/baseline/gp_3d_inducing.npz b/GPy/testing/baseline/gp_3d_inducing.npz
index b6316363..c0bd33b7 100644
Binary files a/GPy/testing/baseline/gp_3d_inducing.npz and b/GPy/testing/baseline/gp_3d_inducing.npz differ
diff --git a/GPy/testing/baseline/gp_3d_mean.npz b/GPy/testing/baseline/gp_3d_mean.npz
index ae0a36ca..09a35cd0 100644
Binary files a/GPy/testing/baseline/gp_3d_mean.npz and b/GPy/testing/baseline/gp_3d_mean.npz differ
diff --git a/GPy/testing/baseline/gp_class_likelihood.npz b/GPy/testing/baseline/gp_class_likelihood.npz
index 288b15ef..40c95423 100644
Binary files a/GPy/testing/baseline/gp_class_likelihood.npz and b/GPy/testing/baseline/gp_class_likelihood.npz differ
diff --git a/GPy/testing/baseline/gp_class_raw.npz b/GPy/testing/baseline/gp_class_raw.npz
index 08bb56df..207196bd 100644
Binary files a/GPy/testing/baseline/gp_class_raw.npz and b/GPy/testing/baseline/gp_class_raw.npz differ
diff --git a/GPy/testing/baseline/gp_class_raw_link.npz b/GPy/testing/baseline/gp_class_raw_link.npz
index 01490f4c..5ff15ea2 100644
Binary files a/GPy/testing/baseline/gp_class_raw_link.npz and b/GPy/testing/baseline/gp_class_raw_link.npz differ
diff --git a/GPy/testing/baseline/gp_conf.npz b/GPy/testing/baseline/gp_conf.npz
index 3eaf087e..801d4c1c 100644
Binary files a/GPy/testing/baseline/gp_conf.npz and b/GPy/testing/baseline/gp_conf.npz differ
diff --git a/GPy/testing/baseline/gp_data.npz b/GPy/testing/baseline/gp_data.npz
index 2a4606df..d25fa6cf 100644
Binary files a/GPy/testing/baseline/gp_data.npz and b/GPy/testing/baseline/gp_data.npz differ
diff --git a/GPy/testing/baseline/gp_density.npz b/GPy/testing/baseline/gp_density.npz
index 6a3235fa..e3fed380 100644
Binary files a/GPy/testing/baseline/gp_density.npz and b/GPy/testing/baseline/gp_density.npz differ
diff --git a/GPy/testing/baseline/gp_in_error.npz b/GPy/testing/baseline/gp_in_error.npz
index 072fbc12..971372ad 100644
Binary files a/GPy/testing/baseline/gp_in_error.npz and b/GPy/testing/baseline/gp_in_error.npz differ
diff --git a/GPy/testing/baseline/gp_mean.npz b/GPy/testing/baseline/gp_mean.npz
index 432b2c25..34553aa2 100644
Binary files a/GPy/testing/baseline/gp_mean.npz and b/GPy/testing/baseline/gp_mean.npz differ
diff --git a/GPy/testing/baseline/gp_out_error.npz b/GPy/testing/baseline/gp_out_error.npz
index 4e86a06d..1f7d9a9d 100644
Binary files a/GPy/testing/baseline/gp_out_error.npz and b/GPy/testing/baseline/gp_out_error.npz differ
diff --git a/GPy/testing/baseline/gp_samples.npz b/GPy/testing/baseline/gp_samples.npz
index 774a576b..36668c74 100644
Binary files a/GPy/testing/baseline/gp_samples.npz and b/GPy/testing/baseline/gp_samples.npz differ
diff --git a/GPy/testing/baseline/gplvm_gradient.npz b/GPy/testing/baseline/gplvm_gradient.npz
index fb6a1658..affb3339 100644
Binary files a/GPy/testing/baseline/gplvm_gradient.npz and b/GPy/testing/baseline/gplvm_gradient.npz differ
diff --git a/GPy/testing/baseline/gplvm_latent.npz b/GPy/testing/baseline/gplvm_latent.npz
index fac4cca9..2ce41632 100644
Binary files a/GPy/testing/baseline/gplvm_latent.npz and b/GPy/testing/baseline/gplvm_latent.npz differ
diff --git a/GPy/testing/baseline/gplvm_latent_3d.npz b/GPy/testing/baseline/gplvm_latent_3d.npz
index 2fe42704..e772b631 100644
Binary files a/GPy/testing/baseline/gplvm_latent_3d.npz and b/GPy/testing/baseline/gplvm_latent_3d.npz differ
diff --git a/GPy/testing/baseline/gplvm_magnification.npz b/GPy/testing/baseline/gplvm_magnification.npz
index fb84f382..9909c5d5 100644
Binary files a/GPy/testing/baseline/gplvm_magnification.npz and b/GPy/testing/baseline/gplvm_magnification.npz differ
diff --git a/GPy/testing/baseline/kern_ARD.npz b/GPy/testing/baseline/kern_ARD.npz
index 7735251e..c2bed84e 100644
Binary files a/GPy/testing/baseline/kern_ARD.npz and b/GPy/testing/baseline/kern_ARD.npz differ
diff --git a/GPy/testing/baseline/kern_cov_1d.npz b/GPy/testing/baseline/kern_cov_1d.npz
index 3bb0ac83..be277cb3 100644
Binary files a/GPy/testing/baseline/kern_cov_1d.npz and b/GPy/testing/baseline/kern_cov_1d.npz differ
diff --git a/GPy/testing/baseline/kern_cov_2d.npz b/GPy/testing/baseline/kern_cov_2d.npz
index 78a5930b..12697c91 100644
Binary files a/GPy/testing/baseline/kern_cov_2d.npz and b/GPy/testing/baseline/kern_cov_2d.npz differ
diff --git a/GPy/testing/baseline/kern_cov_3d.npz b/GPy/testing/baseline/kern_cov_3d.npz
index a76f302f..e07a8154 100644
Binary files a/GPy/testing/baseline/kern_cov_3d.npz and b/GPy/testing/baseline/kern_cov_3d.npz differ
diff --git a/GPy/testing/baseline/kern_cov_no_lim.npz b/GPy/testing/baseline/kern_cov_no_lim.npz
index 019db745..5c9171ca 100644
Binary files a/GPy/testing/baseline/kern_cov_no_lim.npz and b/GPy/testing/baseline/kern_cov_no_lim.npz differ
diff --git a/GPy/testing/baseline/sparse_gp_class_likelihood.npz b/GPy/testing/baseline/sparse_gp_class_likelihood.npz
index 64eb8299..aab87332 100644
Binary files a/GPy/testing/baseline/sparse_gp_class_likelihood.npz and b/GPy/testing/baseline/sparse_gp_class_likelihood.npz differ
diff --git a/GPy/testing/baseline/sparse_gp_class_raw.npz b/GPy/testing/baseline/sparse_gp_class_raw.npz
index c3876b70..1a8d3dd2 100644
Binary files a/GPy/testing/baseline/sparse_gp_class_raw.npz and b/GPy/testing/baseline/sparse_gp_class_raw.npz differ
diff --git a/GPy/testing/baseline/sparse_gp_class_raw_link.npz b/GPy/testing/baseline/sparse_gp_class_raw_link.npz
index 957003e9..d0962d67 100644
Binary files a/GPy/testing/baseline/sparse_gp_class_raw_link.npz and b/GPy/testing/baseline/sparse_gp_class_raw_link.npz differ
diff --git a/GPy/testing/baseline/sparse_gp_data_error.npz b/GPy/testing/baseline/sparse_gp_data_error.npz
index 04a05f55..e9574ce1 100644
Binary files a/GPy/testing/baseline/sparse_gp_data_error.npz and b/GPy/testing/baseline/sparse_gp_data_error.npz differ
diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py
index c2725c8e..053fce35 100644
--- a/GPy/testing/kernel_tests.py
+++ b/GPy/testing/kernel_tests.py
@@ -421,6 +421,21 @@ class KernelGradientTestsContinuous(unittest.TestCase):
         k.randomize()
         self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose))
 
+    def test_OU(self):
+        k = GPy.kern.OU(self.D-1, ARD=True)
+        k.randomize()
+        self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose))
+
+    def test_RatQuad(self):
+        k = GPy.kern.RatQuad(self.D-1, ARD=True)
+        k.randomize()
+        self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose))
+
+    def test_ExpQuad(self):
+        k = GPy.kern.ExpQuad(self.D-1, ARD=True)
+        k.randomize()
+        self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose))
+
     def test_integral(self):
         k = GPy.kern.Integral(1)
         k.randomize()
diff --git a/GPy/testing/likelihood_tests.py b/GPy/testing/likelihood_tests.py
index 346b4b83..629705ab 100644
--- a/GPy/testing/likelihood_tests.py
+++ b/GPy/testing/likelihood_tests.py
@@ -119,7 +119,7 @@ class TestNoiseModels(object):
         self.integer_Y = np.where(tmp > 0, tmp, 0)
         self.ns = np.random.poisson(50, size=self.N)[:, None]
         p = np.abs(np.cos(2*np.pi*self.X + np.random.normal(scale=.2, size=(self.N, self.D)))).mean(1)
-        self.binomial_Y = np.array([np.random.binomial(self.ns[i], p[i]) for i in range(p.shape[0])])[:, None]
+        self.binomial_Y = np.array([np.random.binomial(int(self.ns[i]), p[i]) for i in range(p.shape[0])])[:, None]
         
         self.var = 0.2
         self.deg_free = 4.0
@@ -570,7 +570,6 @@ class TestNoiseModels(object):
         white_var = 1e-4
         kernel = GPy.kern.RBF(X.shape[1]) + GPy.kern.White(X.shape[1])
         laplace_likelihood = GPy.inference.latent_function_inference.Laplace()
-
         m = GPy.core.GP(X.copy(), Y.copy(), kernel, likelihood=model, Y_metadata=Y_metadata, inference_method=laplace_likelihood)
         m.kern.white.constrain_fixed(white_var)
 
diff --git a/GPy/testing/model_tests.py b/GPy/testing/model_tests.py
index 2b969f3e..c6dc50f1 100644
--- a/GPy/testing/model_tests.py
+++ b/GPy/testing/model_tests.py
@@ -1,6 +1,6 @@
 # Copyright (c) 2012, GPy authors (see AUTHORS.txt).
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
-
+from __future__ import division
 
 import unittest
 import numpy as np
@@ -54,7 +54,7 @@ class MiscTests(unittest.TestCase):
         m.randomize()
         m.optimize()
         # Compute the mean of model prediction on 1e5 Monte Carlo samples
-        Xp = np.random.uniform(size=(1e5,2))
+        Xp = np.random.uniform(size=(int(1e5),2))
         Xp[:,0] = Xp[:,0]*15-5
         Xp[:,1] = Xp[:,1]*15
         _, var = m.predict(Xp)
@@ -759,16 +759,18 @@ class GradientTests(np.testing.TestCase):
 
     def test_GP_EP_probit(self):
         N = 20
-        X = np.hstack([np.random.normal(5, 2, N / 2), np.random.normal(10, 2, N / 2)])[:, None]
-        Y = np.hstack([np.ones(N / 2), np.zeros(N / 2)])[:, None]
+        Nhalf = int(N/2)
+        X = np.hstack([np.random.normal(5, 2, Nhalf), np.random.normal(10, 2, Nhalf)])[:, None]
+        Y = np.hstack([np.ones(Nhalf), np.zeros(Nhalf)])[:, None]
         kernel = GPy.kern.RBF(1)
         m = GPy.models.GPClassification(X, Y, kernel=kernel)
         self.assertTrue(m.checkgrad())
 
     def test_sparse_EP_DTC_probit(self):
         N = 20
-        X = np.hstack([np.random.normal(5, 2, N / 2), np.random.normal(10, 2, N / 2)])[:, None]
-        Y = np.hstack([np.ones(N / 2), np.zeros(N / 2)])[:, None]
+        Nhalf = int(N/2)
+        X = np.hstack([np.random.normal(5, 2, Nhalf), np.random.normal(10, 2, Nhalf)])[:, None]
+        Y = np.hstack([np.ones(Nhalf), np.zeros(Nhalf)])[:, None]
         Z = np.linspace(0, 15, 4)[:, None]
         kernel = GPy.kern.RBF(1)
         m = GPy.models.SparseGPClassification(X, Y, kernel=kernel, Z=Z)
@@ -776,8 +778,9 @@ class GradientTests(np.testing.TestCase):
 
     def test_sparse_EP_DTC_probit_uncertain_inputs(self):
         N = 20
-        X = np.hstack([np.random.normal(5, 2, N / 2), np.random.normal(10, 2, N / 2)])[:, None]
-        Y = np.hstack([np.ones(N / 2), np.zeros(N / 2)])[:, None]
+        Nhalf = int(N/2)
+        X = np.hstack([np.random.normal(5, 2, Nhalf), np.random.normal(10, 2, Nhalf)])[:, None]
+        Y = np.hstack([np.ones(Nhalf), np.zeros(Nhalf)])[:, None]
         Z = np.linspace(0, 15, 4)[:, None]
         X_var = np.random.uniform(0.1, 0.2, X.shape)
         kernel = GPy.kern.RBF(1)
diff --git a/GPy/testing/plotting_tests.py b/GPy/testing/plotting_tests.py
index f115c801..dce8b91a 100644
--- a/GPy/testing/plotting_tests.py
+++ b/GPy/testing/plotting_tests.py
@@ -206,7 +206,7 @@ def test_figure():
     with warnings.catch_warnings():
         warnings.simplefilter("ignore")
 
-        ax, _ = pl().new_canvas(num=1)
+        ax, _ = pl().new_canvas(num="imshow_interact")
         def test_func(x):
             return x[:, 0].reshape(3,3)
         pl().imshow_interact(ax, test_func, extent=(-1,1,-1,1), resolution=3)
@@ -228,7 +228,7 @@ def test_figure():
 
         pl().fill_gradient(ax, x, y, facecolors=['r', 'g'], array=array, cmap=cmap)
 
-        ax, _ = pl().new_canvas(num=4, figsize=(4,3), projection='3d', xlabel='x', ylabel='y', zlabel='z', title='awsome title', xlim=(-1,1), ylim=(-1,1), zlim=(-3,3))
+        ax, _ = pl().new_canvas(num="3d_plot", figsize=(4,3), projection='3d', xlabel='x', ylabel='y', zlabel='z', title='awsome title', xlim=(-1,1), ylim=(-1,1), zlim=(-3,3))
         z = 2-np.abs(np.linspace(-2,2,(100)))+1
         x, y = z*np.sin(np.linspace(-2*np.pi,2*np.pi,(100))), z*np.cos(np.linspace(-np.pi,np.pi,(100)))
 
diff --git a/GPy/testing/variational_tests.py b/GPy/testing/variational_tests.py
new file mode 100644
index 00000000..89053b81
--- /dev/null
+++ b/GPy/testing/variational_tests.py
@@ -0,0 +1,63 @@
+'''
+Copyright (c) 2015, Max Zwiessele
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+* Redistributions of source code must retain the above copyright notice, this
+  list of conditions and the following disclaimer.
+
+* Redistributions in binary form must reproduce the above copyright notice,
+  this list of conditions and the following disclaimer in the documentation
+  and/or other materials provided with the distribution.
+
+* Neither the name of paramax nor the names of its
+  contributors may be used to endorse or promote products derived from
+  this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+'''
+import unittest
+import GPy, numpy as np
+
+class KLGrad(GPy.core.Model):
+            def __init__(self, Xvar, kl):   
+                super(KLGrad, self).__init__(name="klgrad")     
+                self.kl = kl
+                self.link_parameter(Xvar)
+                self.Xvar = Xvar
+                self._obj = 0
+            def parameters_changed(self):
+                self.Xvar.gradient[:] = 0
+                self.kl.update_gradients_KL(self.Xvar)
+                self._obj = self.kl.KL_divergence(self.Xvar)
+            def objective_function(self):
+                return self._obj
+        
+class Test(unittest.TestCase):
+
+    def setUp(self):
+        np.random.seed(12345)
+        self.Xvar = GPy.core.parameterization.variational.NormalPosterior(
+            np.random.uniform(0,1,(10,3)), 
+            np.random.uniform(1e-5,.01, (10,3))
+            )
+
+
+    def testNormal(self):
+        klgrad = KLGrad(self.Xvar, GPy.core.parameterization.variational.NormalPrior())
+        np.testing.assert_(klgrad.checkgrad())
+
+if __name__ == "__main__":
+    #import sys;sys.argv = ['', 'Test.testNormal']
+    unittest.main()
\ No newline at end of file
diff --git a/GPy/util/config.py b/GPy/util/config.py
index 53675efe..c6e09e1a 100644
--- a/GPy/util/config.py
+++ b/GPy/util/config.py
@@ -13,7 +13,6 @@ except ImportError:
     config = configparser.ConfigParser()
     from configparser import NoOptionError
 
-
 # This is the default configuration file that always needs to be present.
 default_file = os.path.abspath(os.path.join(os.path.dirname( __file__ ), '..', 'defaults.cfg'))
 
@@ -22,7 +21,7 @@ default_file = os.path.abspath(os.path.join(os.path.dirname( __file__ ), '..', '
 local_file = os.path.abspath(os.path.join(os.path.dirname( __file__ ), '..', 'installation.cfg'))
 
 # This specifies configurations specific to the user (it is found in the user home directory)
-home = os.getenv('HOME') or os.getenv('USERPROFILE')
+home = os.getenv('HOME') or os.getenv('USERPROFILE') or ''
 user_file = os.path.join(home,'.config','GPy', 'user.cfg')
 
 # Read in the given files.
diff --git a/appveyor.yml b/appveyor.yml
index 3738ea1d..4ffda8f9 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.5.6
+  gpy_version: 1.6.1
   matrix:
     - PYTHON_VERSION: 2.7
       MINICONDA: C:\Miniconda-x64
diff --git a/setup.cfg b/setup.cfg
index 77a46f72..514f36eb 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 1.5.6
+current_version = 1.6.1
 tag = True
 commit = True