diff --git a/GPy/core/model.py b/GPy/core/model.py
index c067d51d..21bcf0c7 100644
--- a/GPy/core/model.py
+++ b/GPy/core/model.py
@@ -485,20 +485,17 @@ class Model(Parameterized):
         if not hasattr(self, 'kern'):
             raise ValueError, "this model has no kernel"
 
-        k = [p for p in self.kern._parameters_ if hasattr(p, "ARD") and p.ARD]
-        if (not len(k) == 1):
-            raise ValueError, "cannot determine sensitivity for this kernel"
-        k = k[0]
-        from ..kern.parts.rbf import RBF
-        from ..kern.parts.rbf_inv import RBFInv
-        from ..kern.parts.linear import Linear
+        k = self.kern#[p for p in self.kern._parameters_ if hasattr(p, "ARD") and p.ARD]
+        from ..kern import RBF, Linear#, RBFInv
+
         if isinstance(k, RBF):
             return 1. / k.lengthscale
-        elif isinstance(k, RBFInv):
-            return k.inv_lengthscale
+        #elif isinstance(k, RBFInv):
+        #    return k.inv_lengthscale
         elif isinstance(k, Linear):
             return k.variances
-
+        else:
+            raise ValueError, "cannot determine sensitivity for this kernel"
 
     def pseudo_EM(self, stop_crit=.1, **kwargs):
         """
diff --git a/GPy/core/parameterization/index_operations.py b/GPy/core/parameterization/index_operations.py
index bfd0bf21..b5399741 100644
--- a/GPy/core/parameterization/index_operations.py
+++ b/GPy/core/parameterization/index_operations.py
@@ -83,11 +83,21 @@ class ParameterIndexOperations(object):
     def iterproperties(self):
         return self._properties.iterkeys()
     
-    def shift(self, start, size):
+    def shift_right(self, start, size):
         for ind in self.iterindices():
             toshift = ind>=start
-            if toshift.size > 0:
-                ind[toshift] += size
+            ind[toshift] += size
+
+    def shift_left(self, start, size):
+        for v, ind in self.items():
+            todelete = (ind>=start) * (ind<start+size)
+            if todelete.size != 0:
+                ind = ind[~todelete]
+            toshift = ind>=start
+            if toshift.size != 0:
+                ind[toshift] -= size
+            if ind.size != 0: self._properties[v] = ind
+            else: del self._properties[v]
     
     def clear(self):
         self._properties.clear()
@@ -183,7 +193,7 @@ class ParameterIndexOperationsView(object):
             yield i 
 
 
-    def shift(self, start, size):
+    def shift_right(self, start, size):
         raise NotImplementedError, 'Shifting only supported in original ParamIndexOperations'
     
 
diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py
index 45b57eab..c2c8a05a 100644
--- a/GPy/core/parameterization/parameter_core.py
+++ b/GPy/core/parameterization/parameter_core.py
@@ -390,6 +390,7 @@ class Parameterizable(Constrainable):
         import copy
         from .index_operations import ParameterIndexOperations, ParameterIndexOperationsView
         from .array_core import ParamList
+
         dc = dict()
         for k, v in self.__dict__.iteritems():
             if k not in ['_direct_parent_', '_parameters_', '_parent_index_'] + self.parameter_names():
@@ -399,18 +400,21 @@ class Parameterizable(Constrainable):
                     dc[k] = copy.deepcopy(v)
             if k == '_parameters_':
                 params = [p.copy() for p in v]
-        # dc = copy.deepcopy(self.__dict__)
+        
         dc['_direct_parent_'] = None
         dc['_parent_index_'] = None
         dc['_parameters_'] = ParamList()
+        dc['constraints'].clear()
+        dc['priors'].clear()
+        dc['size'] = 0
+
         s = self.__new__(self.__class__)
         s.__dict__ = dc
-        # import ipdb;ipdb.set_trace()
+        
         for p in params:
             s.add_parameter(p)
-        # dc._notify_parent_change()
+        
         return s
-        # return copy.deepcopy(self)
 
     def _notify_parameters_changed(self):
         self.parameters_changed()
diff --git a/GPy/core/parameterization/parameterized.py b/GPy/core/parameterization/parameterized.py
index 177cc217..d463ed43 100644
--- a/GPy/core/parameterization/parameterized.py
+++ b/GPy/core/parameterization/parameterized.py
@@ -87,8 +87,8 @@ class Parameterized(Parameterizable, Pickleable, Observable, Gradcheckable):
                 self._parameters_.append(param)
             else:
                 start = sum(p.size for p in self._parameters_[:index])
-                self.constraints.shift(start, param.size)
-                self.priors.shift(start, param.size)
+                self.constraints.shift_right(start, param.size)
+                self.priors.shift_right(start, param.size)
                 self.constraints.update(param.constraints, start)
                 self.priors.update(param.priors, start)
                 self._parameters_.insert(index, param)
@@ -113,15 +113,19 @@ class Parameterized(Parameterizable, Pickleable, Observable, Gradcheckable):
         """
         if not param in self._parameters_:
             raise RuntimeError, "Parameter {} does not belong to this object, remove parameters directly from their respective parents".format(param._short())
-        del self._parameters_[param._parent_index_]
+        
+        start = sum([p.size for p in self._parameters_[:param._parent_index_]])
+        self._remove_parameter_name(param)
         self.size -= param.size
+        del self._parameters_[param._parent_index_]
         
         param._disconnect_parent()
-        self._remove_parameter_name(param)
-        
-        #self._notify_parent_change()
+        self.constraints.shift_left(start, param.size)
         self._connect_fixes()
-
+        self._connect_parameters()
+        self._notify_parent_change()
+        
+        
     def _connect_parameters(self):
         # connect parameterlist to this parameterized object
         # This just sets up the right connection for the params objects
diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py
index c8e79e6c..3ba54d34 100644
--- a/GPy/examples/dimensionality_reduction.py
+++ b/GPy/examples/dimensionality_reduction.py
@@ -74,7 +74,7 @@ def gplvm_oil_100(optimize=True, verbose=1, plot=True):
     data = GPy.util.datasets.oil_100()
     Y = data['X']
     # create simple GP model
-    kernel = GPy.kern.RBF(6, ARD=True) + GPy.kern.bias(6)
+    kernel = GPy.kern.RBF(6, ARD=True) + GPy.kern.Bias(6)
     m = GPy.models.GPLVM(Y, 6, kernel=kernel)
     m.data_labels = data['Y'].argmax(axis=1)
     if optimize: m.optimize('scg', messages=verbose)
@@ -190,17 +190,22 @@ def _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim=False):
     _np.random.seed(1234)
     
     x = _np.linspace(0, 4 * _np.pi, N)[:, None]
-    s1 = _np.vectorize(lambda x: _np.sin(x))
+    s1 = _np.vectorize(lambda x: -_np.sin(x))
     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.sin(2 * x))
+    sS = _np.vectorize(lambda x: x*_np.sin(x))
 
     s1 = s1(x)
     s2 = s2(x)
     s3 = s3(x)
     sS = sS(x)
 
-    S1 = _np.hstack([s1, sS])
+    s1 -= s1.mean(); s1 /= s1.std(0)
+    s2 -= s2.mean(); s2 /= s2.std(0)
+    s3 -= s3.mean(); s3 /= s3.std(0)
+    sS -= sS.mean(); sS /= sS.std(0)
+
+    S1 = _np.hstack([s1, s2, sS])
     S2 = _np.hstack([s2, s3, sS])
     S3 = _np.hstack([s3, sS])
 
@@ -271,7 +276,7 @@ def bgplvm_simulation(optimize=True, verbose=1,
     D1, D2, D3, N, num_inducing, Q = 15, 5, 8, 30, 3, 10
     _, _, Ylist = _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim)
     Y = Ylist[0]
-    k = kern.linear(Q, ARD=True)# + kern.white(Q, _np.exp(-2)) # + kern.bias(Q)
+    k = kern.Linear(Q, ARD=True)# + kern.white(Q, _np.exp(-2)) # + kern.bias(Q)
     m = BayesianGPLVM(Y, Q, init="PCA", num_inducing=num_inducing, kernel=k)
     
     if optimize:
@@ -291,10 +296,10 @@ def bgplvm_simulation_missing_data(optimize=True, verbose=1,
     from GPy.models import BayesianGPLVM
     from GPy.inference.latent_function_inference.var_dtc import VarDTCMissingData
 
-    D1, D2, D3, N, num_inducing, Q = 15, 5, 8, 30, 3, 10
+    D1, D2, D3, N, num_inducing, Q = 15, 5, 8, 30, 5, 9
     _, _, Ylist = _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim)
     Y = Ylist[0]
-    k = kern.linear(Q, ARD=True)# + kern.white(Q, _np.exp(-2)) # + kern.bias(Q)
+    k = kern.Linear(Q, ARD=True)# + kern.white(Q, _np.exp(-2)) # + kern.bias(Q)
     
     inan = _np.random.binomial(1, .6, size=Y.shape).astype(bool)
     m = BayesianGPLVM(Y.copy(), Q, init="random", num_inducing=num_inducing, kernel=k)
diff --git a/GPy/inference/latent_function_inference/var_dtc.py b/GPy/inference/latent_function_inference/var_dtc.py
index c2f179ac..349cd72d 100644
--- a/GPy/inference/latent_function_inference/var_dtc.py
+++ b/GPy/inference/latent_function_inference/var_dtc.py
@@ -326,14 +326,14 @@ class VarDTCMissingData(object):
         # gradients:
         if uncertain_inputs:
             grad_dict = {'dL_dKmm': dL_dKmm, 
-                         'dL_dpsi0':dL_dpsi0, 
-                         'dL_dpsi1':dL_dpsi1, 
-                         'dL_dpsi2':dL_dpsi2, 
+                         'dL_dpsi0':dL_dpsi0_all, 
+                         'dL_dpsi1':dL_dpsi1_all, 
+                         'dL_dpsi2':dL_dpsi2_all, 
                          'partial_for_likelihood':partial_for_likelihood}
         else:
             grad_dict = {'dL_dKmm': dL_dKmm, 
-                         'dL_dKdiag':dL_dpsi0, 
-                         'dL_dKnm':dL_dpsi1, 
+                         'dL_dKdiag':dL_dpsi0_all, 
+                         'dL_dKnm':dL_dpsi1_all, 
                          'partial_for_likelihood':partial_for_likelihood}
 
         #get sufficient things for posterior prediction
diff --git a/GPy/kern/__init__.py b/GPy/kern/__init__.py
index 16c13066..e5dc6d35 100644
--- a/GPy/kern/__init__.py
+++ b/GPy/kern/__init__.py
@@ -3,6 +3,7 @@ from _src.white import White
 from _src.kern import Kern
 from _src.linear import Linear
 from _src.brownian import Brownian
+from _src.stationary import Exponential, Matern32, Matern52, ExpQuad
 #from _src.bias import Bias
 #import coregionalize
 #import exponential
diff --git a/GPy/kern/_src/Matern32.py b/GPy/kern/_src/Matern32.py
deleted file mode 100644
index 08fa452c..00000000
--- a/GPy/kern/_src/Matern32.py
+++ /dev/null
@@ -1,139 +0,0 @@
-# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
-# Licensed under the BSD 3-clause license (see LICENSE.txt)
-
-
-from kernpart import Kernpart
-import numpy as np
-from scipy import integrate
-
-class Matern32(Kernpart):
-    """
-    Matern 3/2 kernel:
-
-    .. math::
-
-       k(r) = \\sigma^2 (1 + \\sqrt{3} r) \exp(- \sqrt{3} r) \\ \\ \\ \\  \\text{ where  } r = \sqrt{\sum_{i=1}^input_dim \\frac{(x_i-y_i)^2}{\ell_i^2} }
-
-    :param input_dim: the number of input dimensions
-    :type input_dim: int
-    :param variance: the variance :math:`\sigma^2`
-    :type variance: float
-    :param lengthscale: the vector of lengthscale :math:`\ell_i`
-    :type lengthscale: array or list of the appropriate size (or float if there is only one lengthscale parameter)
-    :param ARD: Auto Relevance Determination. If equal to "False", the kernel is isotropic (ie. one single lengthscale parameter \ell), otherwise there is one lengthscale parameter per dimension.
-    :type ARD: Boolean
-    :rtype: kernel object
-
-    """
-
-    def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False):
-        self.input_dim = input_dim
-        self.ARD = ARD
-        if ARD == False:
-            self.num_params = 2
-            self.name = 'Mat32'
-            if lengthscale is not None:
-                lengthscale = np.asarray(lengthscale)
-                assert lengthscale.size == 1, "Only one lengthscale needed for non-ARD kernel"
-            else:
-                lengthscale = np.ones(1)
-        else:
-            self.num_params = self.input_dim + 1
-            self.name = 'Mat32'
-            if lengthscale is not None:
-                lengthscale = np.asarray(lengthscale)
-                assert lengthscale.size == self.input_dim, "bad number of lengthscales"
-            else:
-                lengthscale = np.ones(self.input_dim)
-        self._set_params(np.hstack((variance, lengthscale.flatten())))
-
-    def _get_params(self):
-        """return the value of the parameters."""
-        return np.hstack((self.variance, self.lengthscale))
-
-    def _set_params(self, x):
-        """set the value of the parameters."""
-        assert x.size == self.num_params
-        self.variance = x[0]
-        self.lengthscale = x[1:]
-
-    def _get_param_names(self):
-        """return parameter names."""
-        if self.num_params == 2:
-            return ['variance', 'lengthscale']
-        else:
-            return ['variance'] + ['lengthscale_%i' % i for i in range(self.lengthscale.size)]
-
-    def K(self, X, X2, target):
-        """Compute the covariance matrix between X and X2."""
-        if X2 is None: X2 = X
-        dist = np.sqrt(np.sum(np.square((X[:, None, :] - X2[None, :, :]) / self.lengthscale), -1))
-        np.add(self.variance * (1 + np.sqrt(3.) * dist) * np.exp(-np.sqrt(3.) * dist), target, target)
-
-    def Kdiag(self, X, target):
-        """Compute the diagonal of the covariance matrix associated to X."""
-        np.add(target, self.variance, target)
-
-    def _param_grad_helper(self, dL_dK, X, X2, target):
-        """derivative of the covariance matrix with respect to the parameters."""
-        if X2 is None: X2 = X
-        dist = np.sqrt(np.sum(np.square((X[:, None, :] - X2[None, :, :]) / self.lengthscale), -1))
-        dvar = (1 + np.sqrt(3.) * dist) * np.exp(-np.sqrt(3.) * dist)
-        invdist = 1. / np.where(dist != 0., dist, np.inf)
-        dist2M = np.square(X[:, None, :] - X2[None, :, :]) / self.lengthscale ** 3
-        # dl = (self.variance* 3 * dist * np.exp(-np.sqrt(3.)*dist))[:,:,np.newaxis] * dist2M*invdist[:,:,np.newaxis]
-        target[0] += np.sum(dvar * dL_dK)
-        if self.ARD == True:
-            dl = (self.variance * 3 * dist * np.exp(-np.sqrt(3.) * dist))[:, :, np.newaxis] * dist2M * invdist[:, :, np.newaxis]
-            # dl = self.variance*dvar[:,:,None]*dist2M*invdist[:,:,None]
-            target[1:] += (dl * dL_dK[:, :, None]).sum(0).sum(0)
-        else:
-            dl = (self.variance * 3 * dist * np.exp(-np.sqrt(3.) * dist)) * dist2M.sum(-1) * invdist
-            # dl = self.variance*dvar*dist2M.sum(-1)*invdist
-            target[1] += np.sum(dl * dL_dK)
-
-    def dKdiag_dtheta(self, dL_dKdiag, X, target):
-        """derivative of the diagonal of the covariance matrix with respect to the parameters."""
-        target[0] += np.sum(dL_dKdiag)
-
-    def gradients_X(self, dL_dK, X, X2, target):
-        """derivative of the covariance matrix with respect to X."""
-        if X2 is None:
-            dist = np.sqrt(np.sum(np.square((X[:, None, :] - X[None, :, :]) / self.lengthscale), -1))[:, :, None]
-            ddist_dX = 2*(X[:, None, :] - X[None, :, :]) / self.lengthscale ** 2 / np.where(dist != 0., dist, np.inf)
-
-        else:
-            dist = np.sqrt(np.sum(np.square((X[:, None, :] - X2[None, :, :]) / self.lengthscale), -1))[:, :, None]
-            ddist_dX = (X[:, None, :] - X2[None, :, :]) / self.lengthscale ** 2 / np.where(dist != 0., dist, np.inf)
-        gradients_X = -np.transpose(3 * self.variance * dist * np.exp(-np.sqrt(3) * dist) * ddist_dX, (1, 0, 2))
-        target += np.sum(gradients_X * dL_dK.T[:, :, None], 0)
-
-    def dKdiag_dX(self, dL_dKdiag, X, target):
-        pass
-
-    def Gram_matrix(self, F, F1, F2, lower, upper):
-        """
-        Return the Gram matrix of the vector of functions F with respect to the RKHS norm. The use of this function is limited to input_dim=1.
-
-        :param F: vector of functions
-        :type F: np.array
-        :param F1: vector of derivatives of F
-        :type F1: np.array
-        :param F2: vector of second derivatives of F
-        :type F2: np.array
-        :param lower,upper: boundaries of the input domain
-        :type lower,upper: floats
-        """
-        assert self.input_dim == 1
-        def L(x, i):
-            return(3. / self.lengthscale ** 2 * F[i](x) + 2 * np.sqrt(3) / self.lengthscale * F1[i](x) + F2[i](x))
-        n = F.shape[0]
-        G = np.zeros((n, n))
-        for i in range(n):
-            for j in range(i, n):
-                G[i, j] = G[j, i] = integrate.quad(lambda x : L(x, i) * L(x, j), lower, upper)[0]
-        Flower = np.array([f(lower) for f in F])[:, None]
-        F1lower = np.array([f(lower) for f in F1])[:, None]
-        # print "OLD \n", np.dot(F1lower,F1lower.T), "\n \n"
-        # return(G)
-        return(self.lengthscale ** 3 / (12.*np.sqrt(3) * self.variance) * G + 1. / self.variance * np.dot(Flower, Flower.T) + self.lengthscale ** 2 / (3.*self.variance) * np.dot(F1lower, F1lower.T))
diff --git a/GPy/kern/_src/Matern52.py b/GPy/kern/_src/Matern52.py
deleted file mode 100644
index 7d36254c..00000000
--- a/GPy/kern/_src/Matern52.py
+++ /dev/null
@@ -1,145 +0,0 @@
-# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
-# Licensed under the BSD 3-clause license (see LICENSE.txt)
-
-
-from kernpart import Kernpart
-import numpy as np
-import hashlib
-from scipy import integrate
-
-class Matern52(Kernpart):
-    """
-    Matern 5/2 kernel:
-
-    .. math::
-
-       k(r) = \sigma^2 (1 + \sqrt{5} r + \\frac53 r^2) \exp(- \sqrt{5} r) \ \ \ \ \  \\text{ where  } r = \sqrt{\sum_{i=1}^input_dim \\frac{(x_i-y_i)^2}{\ell_i^2} }
-
-    :param input_dim: the number of input dimensions
-    :type input_dim: int
-    :param variance: the variance :math:`\sigma^2`
-    :type variance: float
-    :param lengthscale: the vector of lengthscale :math:`\ell_i`
-    :type lengthscale: array or list of the appropriate size (or float if there is only one lengthscale parameter)
-    :param ARD: Auto Relevance Determination. If equal to "False", the kernel is isotropic (ie. one single lengthscale parameter \ell), otherwise there is one lengthscale parameter per dimension.
-    :type ARD: Boolean
-    :rtype: kernel object
-
-    """
-    def __init__(self,input_dim,variance=1.,lengthscale=None,ARD=False):
-        self.input_dim = input_dim
-        self.ARD = ARD
-        if ARD == False:
-            self.num_params = 2
-            self.name = 'Mat52'
-            if lengthscale is not None:
-                lengthscale = np.asarray(lengthscale)
-                assert lengthscale.size == 1, "Only one lengthscale needed for non-ARD kernel"
-            else:
-                lengthscale = np.ones(1)
-        else:
-            self.num_params = self.input_dim + 1
-            self.name = 'Mat52'
-            if lengthscale is not None:
-                lengthscale = np.asarray(lengthscale)
-                assert lengthscale.size == self.input_dim, "bad number of lengthscales"
-            else:
-                lengthscale = np.ones(self.input_dim)
-        self._set_params(np.hstack((variance,lengthscale.flatten())))
-
-    def _get_params(self):
-        """return the value of the parameters."""
-        return np.hstack((self.variance,self.lengthscale))
-
-    def _set_params(self,x):
-        """set the value of the parameters."""
-        assert x.size == self.num_params
-        self.variance = x[0]
-        self.lengthscale = x[1:]
-
-    def _get_param_names(self):
-        """return parameter names."""
-        if self.num_params == 2:
-            return ['variance','lengthscale']
-        else:
-            return ['variance']+['lengthscale_%i'%i for i in range(self.lengthscale.size)]
-
-    def K(self,X,X2,target):
-        """Compute the covariance matrix between X and X2."""
-        if X2 is None: X2 = X
-        dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1))
-        np.add(self.variance*(1+np.sqrt(5.)*dist+5./3*dist**2)*np.exp(-np.sqrt(5.)*dist), target,target)
-
-    def Kdiag(self,X,target):
-        """Compute the diagonal of the covariance matrix associated to X."""
-        np.add(target,self.variance,target)
-
-    def _param_grad_helper(self,dL_dK,X,X2,target):
-        """derivative of the covariance matrix with respect to the parameters."""
-        if X2 is None: X2 = X
-        dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1))
-        invdist = 1./np.where(dist!=0.,dist,np.inf)
-        dist2M = np.square(X[:,None,:]-X2[None,:,:])/self.lengthscale**3
-        dvar = (1+np.sqrt(5.)*dist+5./3*dist**2)*np.exp(-np.sqrt(5.)*dist)
-        dl = (self.variance * 5./3 * dist * (1 + np.sqrt(5.)*dist ) * np.exp(-np.sqrt(5.)*dist))[:,:,np.newaxis] * dist2M*invdist[:,:,np.newaxis]
-        target[0] += np.sum(dvar*dL_dK)
-        if self.ARD:
-            dl = (self.variance * 5./3 * dist * (1 + np.sqrt(5.)*dist ) * np.exp(-np.sqrt(5.)*dist))[:,:,np.newaxis] * dist2M*invdist[:,:,np.newaxis]
-            #dl = (self.variance* 3 * dist * np.exp(-np.sqrt(3.)*dist))[:,:,np.newaxis] * dist2M*invdist[:,:,np.newaxis]
-            target[1:] += (dl*dL_dK[:,:,None]).sum(0).sum(0)
-        else:
-            dl = (self.variance * 5./3 * dist * (1 + np.sqrt(5.)*dist ) * np.exp(-np.sqrt(5.)*dist)) * dist2M.sum(-1)*invdist
-            #dl = (self.variance* 3 * dist * np.exp(-np.sqrt(3.)*dist)) * dist2M.sum(-1)*invdist
-            target[1] += np.sum(dl*dL_dK)
-
-    def dKdiag_dtheta(self,dL_dKdiag,X,target):
-        """derivative of the diagonal of the covariance matrix with respect to the parameters."""
-        target[0] += np.sum(dL_dKdiag)
-
-    def gradients_X(self,dL_dK,X,X2,target):
-        """derivative of the covariance matrix with respect to X."""
-        if X2 is None:
-            dist = np.sqrt(np.sum(np.square((X[:,None,:]-X[None,:,:])/self.lengthscale),-1))[:,:,None]
-            ddist_dX = 2*(X[:,None,:]-X[None,:,:])/self.lengthscale**2/np.where(dist!=0.,dist,np.inf)
-        else:
-            dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1))[:,:,None]
-            ddist_dX = (X[:,None,:]-X2[None,:,:])/self.lengthscale**2/np.where(dist!=0.,dist,np.inf)
-        gradients_X = -  np.transpose(self.variance*5./3*dist*(1+np.sqrt(5)*dist)*np.exp(-np.sqrt(5)*dist)*ddist_dX,(1,0,2))
-        target += np.sum(gradients_X*dL_dK.T[:,:,None],0)
-
-    def dKdiag_dX(self,dL_dKdiag,X,target):
-        pass
-
-    def Gram_matrix(self,F,F1,F2,F3,lower,upper):
-        """
-        Return the Gram matrix of the vector of functions F with respect to the RKHS norm. The use of this function is limited to input_dim=1.
-
-        :param F: vector of functions
-        :type F: np.array
-        :param F1: vector of derivatives of F
-        :type F1: np.array
-        :param F2: vector of second derivatives of F
-        :type F2: np.array
-        :param F3: vector of third derivatives of F
-        :type F3: np.array
-        :param lower,upper: boundaries of the input domain
-        :type lower,upper: floats
-        """
-        assert self.input_dim == 1
-        def L(x,i):
-            return(5*np.sqrt(5)/self.lengthscale**3*F[i](x) + 15./self.lengthscale**2*F1[i](x)+ 3*np.sqrt(5)/self.lengthscale*F2[i](x) + F3[i](x))
-        n = F.shape[0]
-        G = np.zeros((n,n))
-        for i in range(n):
-            for j in range(i,n):
-                G[i,j] = G[j,i] = integrate.quad(lambda x : L(x,i)*L(x,j),lower,upper)[0]
-        G_coef = 3.*self.lengthscale**5/(400*np.sqrt(5))
-        Flower = np.array([f(lower) for f in F])[:,None]
-        F1lower = np.array([f(lower) for f in F1])[:,None]
-        F2lower = np.array([f(lower) for f in F2])[:,None]
-        orig = 9./8*np.dot(Flower,Flower.T) + 9.*self.lengthscale**4/200*np.dot(F2lower,F2lower.T)
-        orig2 = 3./5*self.lengthscale**2 * ( np.dot(F1lower,F1lower.T) + 1./8*np.dot(Flower,F2lower.T) + 1./8*np.dot(F2lower,Flower.T))
-        return(1./self.variance* (G_coef*G + orig + orig2))
-
-
-
diff --git a/GPy/kern/_src/coregionalize.py b/GPy/kern/_src/coregionalize.py
index 69fc27ef..74cd2a1d 100644
--- a/GPy/kern/_src/coregionalize.py
+++ b/GPy/kern/_src/coregionalize.py
@@ -5,6 +5,7 @@ from kern import Kern
 import numpy as np
 from scipy import weave
 from ...core.parameterization import Param
+from ...core.parameterization.transformations import Logexp
 
 class Coregionalize(Kern):
     """
@@ -20,7 +21,7 @@ class Coregionalize(Kern):
        k_2(x, y)=\mathbf{B} k(x, y)
 
     it is obtained as the tensor product between a covariance function
-    k(x,y) and B.
+    k(x, y) and B.
 
     :param output_dim: number of outputs to coregionalize
     :type output_dim: int
@@ -29,7 +30,7 @@ class Coregionalize(Kern):
     :param W: a low rank matrix that determines the correlations between the different outputs, together with kappa it forms the coregionalization matrix B
     :type W: numpy array of dimensionality (num_outpus, W_columns)
     :param kappa: a vector which allows the outputs to behave independently
-    :type kappa: numpy array of dimensionality  (output_dim,)
+    :type kappa: numpy array of dimensionality  (output_dim, )
 
     .. note: see coregionalization examples in GPy.examples.regression for some usage.
     """
@@ -37,18 +38,18 @@ class Coregionalize(Kern):
         super(Coregionalize, self).__init__(input_dim=1, name=name)
         self.output_dim = output_dim
         self.rank = rank
-        if self.rank>output_dim-1:
+        if self.rank>output_dim:
             print("Warning: Unusual choice of rank, it should normally be less than the output_dim.")
         if W is None:
-            W = 0.5*np.random.randn(self.output_dim,self.rank)/np.sqrt(self.rank)
+            W = 0.5*np.random.randn(self.output_dim, self.rank)/np.sqrt(self.rank)
         else:
-            assert W.shape==(self.output_dim,self.rank)
-        self.W = Param('W',W)
+            assert W.shape==(self.output_dim, self.rank)
+        self.W = Param('W', W)
         if kappa is None:
             kappa = 0.5*np.ones(self.output_dim)
         else:
-            assert kappa.shape==(self.output_dim,)
-        self.kappa = Param('kappa', kappa)
+            assert kappa.shape==(self.output_dim, )
+        self.kappa = Param('kappa', kappa, Logexp())
         self.add_parameters(self.W, self.kappa)
         self.parameters_changed()
 
@@ -56,54 +57,58 @@ class Coregionalize(Kern):
     def parameters_changed(self):
         self.B = np.dot(self.W, self.W.T) + np.diag(self.kappa)
 
-    def K(self,index,index2,target):
-        index = np.asarray(index,dtype=np.int)
+    def K(self, X, X2=None):
+        index = np.asarray(X, dtype=np.int)
 
         #here's the old code (numpy)
         #if index2 is None:
             #index2 = index
         #else:
-            #index2 = np.asarray(index2,dtype=np.int)
+            #index2 = np.asarray(index2, dtype=np.int)
         #false_target = target.copy()
-        #ii,jj = np.meshgrid(index,index2)
-        #ii,jj = ii.T, jj.T
-        #false_target += self.B[ii,jj]
+        #ii, jj = np.meshgrid(index, index2)
+        #ii, jj = ii.T, jj.T
+        #false_target += self.B[ii, jj]
 
-        if index2 is None:
+
+        if X2 is None:
+            target = np.empty((X.shape[0], X.shape[0]), dtype=np.float64)
             code="""
             for(int i=0;i<N; i++){
-              target[i+i*N] += B[index[i]+output_dim*index[i]];
+              target[i+i*N] = B[index[i]+output_dim*index[i]];
               for(int j=0; j<i; j++){
-                  target[j+i*N] += B[index[i]+output_dim*index[j]];
-                  target[i+j*N] += target[j+i*N];
+                  target[j+i*N] = B[index[i]+output_dim*index[j]];
+                  target[i+j*N] = target[j+i*N];
                 }
               }
             """
-            N,B,output_dim = index.size, self.B, self.output_dim
-            weave.inline(code,['target','index','N','B','output_dim'])
+            N, B, output_dim = index.size, self.B, self.output_dim
+            weave.inline(code, ['target', 'index', 'N', 'B', 'output_dim'])
         else:
-            index2 = np.asarray(index2,dtype=np.int)
+            index2 = np.asarray(X2, dtype=np.int)
+            target = np.empty((X.shape[0], X2.shape[0]), dtype=np.float64)
             code="""
             for(int i=0;i<num_inducing; i++){
               for(int j=0; j<N; j++){
-                  target[i+j*num_inducing] += B[output_dim*index[j]+index2[i]];
+                  target[i+j*num_inducing] = B[output_dim*index[j]+index2[i]];
                 }
               }
             """
-            N,num_inducing,B,output_dim = index.size,index2.size, self.B, self.output_dim
-            weave.inline(code,['target','index','index2','N','num_inducing','B','output_dim'])
+            N, num_inducing, B, output_dim = index.size, index2.size, self.B, self.output_dim
+            weave.inline(code, ['target', 'index', 'index2', 'N', 'num_inducing', 'B', 'output_dim'])
+        return target
 
 
-    def Kdiag(self,index,target):
-        target += np.diag(self.B)[np.asarray(index,dtype=np.int).flatten()]
+    def Kdiag(self, X):
+        return np.diag(self.B)[np.asarray(X, dtype=np.int).flatten()]
 
-    def update_gradients_full(self,dL_dK, index, index2=None):
-        index = np.asarray(index,dtype=np.int)
+    def update_gradients_full(self, dL_dK, X, X2=None):
+        index = np.asarray(X, dtype=np.int)
         dL_dK_small = np.zeros_like(self.B)
-        if index2 is None:
+        if X2 is None:
             index2 = index
         else:
-            index2 = np.asarray(index2,dtype=np.int)
+            index2 = np.asarray(X2, dtype=np.int)
 
         code="""
         for(int i=0; i<num_inducing; i++){
@@ -113,28 +118,24 @@ class Coregionalize(Kern):
         }
         """
         N, num_inducing, output_dim = index.size, index2.size, self.output_dim
-        weave.inline(code, ['N','num_inducing','output_dim','dL_dK','dL_dK_small','index','index2'])
+        weave.inline(code, ['N', 'num_inducing', 'output_dim', 'dL_dK', 'dL_dK_small', 'index', 'index2'])
 
         dkappa = np.diag(dL_dK_small)
         dL_dK_small += dL_dK_small.T
-        dW = (self.W[:,None,:]*dL_dK_small[:,:,None]).sum(0)
+        dW = (self.W[:, None, :]*dL_dK_small[:, :, None]).sum(0)
 
         self.W.gradient = dW
         self.kappa.gradient = dkappa
 
-    def update_gradients_sparse(self, dL_dKmm, dL_dKnm, dL_dKdiag, X, Z):
-        raise NotImplementedError, "some code below"
-    #def dKdiag_dtheta(self,dL_dKdiag,index,target):
-        #index = np.asarray(index,dtype=np.int).flatten()
-        #dL_dKdiag_small = np.zeros(self.output_dim)
-        #for i in range(self.output_dim):
-            #dL_dKdiag_small[i] += np.sum(dL_dKdiag[index==i])
-        #dW = 2.*self.W*dL_dKdiag_small[:,None]
-        #dkappa = dL_dKdiag_small
-        #target += np.hstack([dW.flatten(),dkappa])
+    def update_gradients_diag(self, dL_dKdiag, X):
+        index = np.asarray(X, dtype=np.int).flatten()
+        dL_dKdiag_small = np.array([dL_dKdiag[index==i] for i in xrange(output_dim)])
+        self.W.gradient = 2.*self.W*dL_dKdiag_small[:, None]
+        self.kappa.gradient = dL_dKdiag_small
+
+    def gradients_X(self, dL_dK, X, X2=None):
+        return np.zeros(X.shape)
+
+    def gradients_X_diag(self, dL_dKdiag, X):
+        return np.zeros(X.shape)
 
-    def gradients_X(self,dL_dK,X,X2):
-        if X2 is None:
-            return np.zeros((X.shape[0], X.shape[0]))
-        else:
-            return np.zeros((X.shape[0], X2.shape[0]))
diff --git a/GPy/kern/_src/exponential.py b/GPy/kern/_src/exponential.py
deleted file mode 100644
index 372d4d9b..00000000
--- a/GPy/kern/_src/exponential.py
+++ /dev/null
@@ -1,129 +0,0 @@
-# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
-# Licensed under the BSD 3-clause license (see LICENSE.txt)
-
-
-from kernpart import Kernpart
-import numpy as np
-from scipy import integrate
-
-class Exponential(Kernpart):
-    """
-    Exponential kernel (aka Ornstein-Uhlenbeck or Matern 1/2)
-
-    .. math::
-
-       k(r) = \sigma^2 \exp(- r) \ \ \ \ \  \\text{ where  } r = \sqrt{\sum_{i=1}^input_dim \\frac{(x_i-y_i)^2}{\ell_i^2} }
-
-    :param input_dim: the number of input dimensions
-    :type input_dim: int
-    :param variance: the variance :math:`\sigma^2`
-    :type variance: float
-    :param lengthscale: the vector of lengthscale :math:`\ell_i`
-    :type lengthscale: array or list of the appropriate size (or float if there is only one lengthscale parameter)
-    :param ARD: Auto Relevance Determination. If equal to "False", the kernel is isotropic (ie. one single lengthscale parameter \ell), otherwise there is one lengthscale parameter per dimension.
-    :type ARD: Boolean
-    :param name: the name of the kernel
-    :rtype: kernel object
-
-    """
-    def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, name='exp'):
-        self.input_dim = input_dim
-        self.ARD = ARD
-        self.variance = variance
-        self.name = name
-        if ARD == False:
-            self.num_params = 2
-            if lengthscale is not None:
-                lengthscale = np.asarray(lengthscale)
-                assert lengthscale.size == 1, "Only one lengthscale needed for non-ARD kernel"
-            else:
-                lengthscale = np.ones(1)
-        else:
-            self.num_params = self.input_dim + 1
-            if lengthscale is not None:
-                lengthscale = np.asarray(lengthscale)
-                assert lengthscale.size == self.input_dim, "bad number of lengthscales"
-            else:
-                lengthscale = np.ones(self.input_dim)
-        #self._set_params(np.hstack((variance, lengthscale.flatten())))
-        self.set_as_parameter('variance', 'lengthscale')
-
-#     def _get_params(self):
-#         """return the value of the parameters."""
-#         return np.hstack((self.variance, self.lengthscale))
-# 
-#     def _set_params(self, x):
-#         """set the value of the parameters."""
-#         assert x.size == self.num_params
-#         self.variance = x[0]
-#         self.lengthscale = x[1:]
-# 
-#     def _get_param_names(self):
-#         """return parameter names."""
-#         if self.num_params == 2:
-#             return ['variance', 'lengthscale']
-#         else:
-#             return ['variance'] + ['lengthscale_%i' % i for i in range(self.lengthscale.size)]
-
-    def K(self, X, X2, target):
-        """Compute the covariance matrix between X and X2."""
-        if X2 is None: X2 = X
-        dist = np.sqrt(np.sum(np.square((X[:, None, :] - X2[None, :, :]) / self.lengthscale), -1))
-        np.add(self.variance * np.exp(-dist), target, target)
-
-    def Kdiag(self, X, target):
-        """Compute the diagonal of the covariance matrix associated to X."""
-        np.add(target, self.variance, target)
-
-    def _param_grad_helper(self, dL_dK, X, X2, target):
-        """derivative of the covariance matrix with respect to the parameters."""
-        if X2 is None: X2 = X
-        dist = np.sqrt(np.sum(np.square((X[:, None, :] - X2[None, :, :]) / self.lengthscale), -1))
-        invdist = 1. / np.where(dist != 0., dist, np.inf)
-        dist2M = np.square(X[:, None, :] - X2[None, :, :]) / self.lengthscale ** 3
-        dvar = np.exp(-dist)
-        target[0] += np.sum(dvar * dL_dK)
-        if self.ARD == True:
-            dl = self.variance * dvar[:, :, None] * dist2M * invdist[:, :, None]
-            target[1:] += (dl * dL_dK[:, :, None]).sum(0).sum(0)
-        else:
-            dl = self.variance * dvar * dist2M.sum(-1) * invdist
-            target[1] += np.sum(dl * dL_dK)
-
-    def dKdiag_dtheta(self, dL_dKdiag, X, target):
-        """derivative of the diagonal of the covariance matrix with respect to the parameters."""
-        # NB: derivative of diagonal elements wrt lengthscale is 0
-        target[0] += np.sum(dL_dKdiag)
-
-    def gradients_X(self, dL_dK, X, X2, target):
-        """derivative of the covariance matrix with respect to X."""
-        if X2 is None: X2 = X
-        dist = np.sqrt(np.sum(np.square((X[:, None, :] - X2[None, :, :]) / self.lengthscale), -1))[:, :, None]
-        ddist_dX = (X[:, None, :] - X2[None, :, :]) / self.lengthscale ** 2 / np.where(dist != 0., dist, np.inf)
-        gradients_X = -np.transpose(self.variance * np.exp(-dist) * ddist_dX, (1, 0, 2))
-        target += np.sum(gradients_X * dL_dK.T[:, :, None], 0)
-
-    def dKdiag_dX(self, dL_dKdiag, X, target):
-        pass
-
-    def Gram_matrix(self, F, F1, lower, upper):
-        """
-        Return the Gram matrix of the vector of functions F with respect to the RKHS norm. The use of this function is limited to input_dim=1.
-
-        :param F: vector of functions
-        :type F: np.array
-        :param F1: vector of derivatives of F
-        :type F1: np.array
-        :param lower,upper: boundaries of the input domain
-        :type lower,upper: floats
-        """
-        assert self.input_dim == 1
-        def L(x, i):
-            return(1. / self.lengthscale * F[i](x) + F1[i](x))
-        n = F.shape[0]
-        G = np.zeros((n, n))
-        for i in range(n):
-            for j in range(i, n):
-                G[i, j] = G[j, i] = integrate.quad(lambda x : L(x, i) * L(x, j), lower, upper)[0]
-        Flower = np.array([f(lower) for f in F])[:, None]
-        return(self.lengthscale / 2. / self.variance * G + 1. / self.variance * np.dot(Flower, Flower.T))
diff --git a/GPy/kern/_src/kern.py b/GPy/kern/_src/kern.py
index 5fe29d51..3ef231b3 100644
--- a/GPy/kern/_src/kern.py
+++ b/GPy/kern/_src/kern.py
@@ -67,7 +67,7 @@ class Kern(Parameterized):
         See GPy.plotting.matplot_dep.plot_ARD
         """
         assert "matplotlib" in sys.modules, "matplotlib package has not been imported."
-        from ..plotting.matplot_dep import kernel_plots
+        from ...plotting.matplot_dep import kernel_plots
         return kernel_plots.plot_ARD(self,*args)
 
 
diff --git a/GPy/kern/_src/stationary.py b/GPy/kern/_src/stationary.py
new file mode 100644
index 00000000..7cc2e695
--- /dev/null
+++ b/GPy/kern/_src/stationary.py
@@ -0,0 +1,211 @@
+# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
+# Licensed under the BSD 3-clause license (see LICENSE.txt)
+
+
+from kern import Kern
+from ...core.parameterization import Param
+from ...core.parameterization.transformations import Logexp
+from ... import util
+import numpy as np
+from scipy import integrate
+
+class Stationary(Kern):
+    def __init__(self, input_dim, variance, lengthscale, ARD, name):
+        super(Stationary, self).__init__(input_dim, name)
+        self.ARD = ARD
+        if not ARD:
+            if lengthscale is None:
+                lengthscale = np.ones(1)
+            else:
+                lengthscale = np.asarray(lengthscale)
+                assert lengthscale.size == 1 "Only  lengthscale needed for non-ARD kernel"
+        else:
+            if lengthscale is not None:
+                lengthscale = np.asarray(lengthscale)
+                assert lengthscale.size in [1, input_dim], "Bad lengthscales"
+                if lengthscale.size != input_dim:
+                    lengthscale = np.ones(input_dim)*lengthscale
+            else:
+                lengthscale = np.ones(self.input_dim)
+        self.lengthscale = Param('lengthscale', lengthscale, Logexp())
+        self.variance = Param('variance', variance, Logexp())
+        assert self.variance.size==1
+        self.add_parameters(self.variance, self.lengthscale)
+
+    def _dist(self, X, X2):
+        if X2 is None:
+            X2 = X
+        return X[:, None, :] - X2[None, :, :]
+
+    def _scaled_dist(self, X, X2=None):
+        return np.sqrt(np.sum(np.square(self._dist(X, X2) / self.lengthscale), -1))
+
+    def Kdiag(self, X):
+        ret = np.empty(X.shape[0])
+        ret[:] = self.variance
+        return ret
+
+    def update_gradients_diag(self, dL_dKdiag, X):
+        self.variance.gradient = np.sum(dL_dKdiag)
+        self.lengthscale.gradient = 0.
+
+    def update_gradients_full(self, dL_dK, X, X2=None):
+        K = self.K(X, X2)
+        self.variance.gradient = np.sum(K * dL_dK)/self.variance
+
+        rinv = self._inv_dist(X, X2)
+        dL_dr = self.dK_dr(X, X2) * dL_dK
+        x_xl3 = np.square(self._dist(X, X2)) / self.lengthscale**3
+
+        if self.ARD:
+            self.lengthscale.gradient = -((dL_dr*rinv)[:,:,None]*x_xl3).sum(0).sum(0)
+        else:
+            self.lengthscale.gradient = -((dL_dr*rinv)[:,:,None]*x_xl3).sum()
+
+    def _inv_dist(self, X, X2=None):
+        dist = self._scaled_dist(X, X2)
+        if X2 is None:
+            nondiag = util.diag.offdiag_view(dist)
+            nondiag[:] = 1./nondiag
+            return dist
+        else:
+            return 1./np.where(dist != 0., dist, np.inf)
+
+    def gradients_X(self, dL_dK, X, X2=None):
+        dL_dr = self.dK_dr(X, X2) * dL_dK
+        invdist = self._inv_dist(X, X2)
+        ret = np.sum((invdist*dL_dr)[:,:,None]*self._dist(X, X2),1)/self.lengthscale**2
+        if X2 is None:
+            ret *= 2.
+        return ret
+
+    def gradients_X_diag(self, dL_dKdiag, X):
+        return np.zeros(X.shape)
+
+
+
+
+class Exponential(Stationary):
+    def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, name='Exponential'):
+        super(Exponential, self).__init__(input_dim, variance, lengthscale, ARD, name)
+
+    def K(self, X, X2=None):
+        dist = self._scaled_dist(X, X2)
+        return self.variance * np.exp(-0.5 * dist)
+
+    def dK_dr(self, X, X2):
+        return -0.5*self.K(X, X2)
+
+class Matern32(Stationary):
+    """
+    Matern 3/2 kernel:
+
+    .. math::
+
+       k(r) = \\sigma^2 (1 + \\sqrt{3} r) \exp(- \sqrt{3} r) \\ \\ \\ \\  \\text{ where  } r = \sqrt{\sum_{i=1}^input_dim \\frac{(x_i-y_i)^2}{\ell_i^2} }
+
+    """
+
+    def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, name='Mat32'):
+        super(Matern32, self).__init__(input_dim, variance, lengthscale, ARD, name)
+
+    def K(self, X, X2=None):
+        dist = self._scaled_dist(X, X2)
+        return self.variance * (1. + np.sqrt(3.) * dist) * np.exp(-np.sqrt(3.) * dist)
+
+    def dK_dr(self, X, X2):
+        dist = self._scaled_dist(X, X2)
+        return -3.*self.variance*dist*np.exp(-np.sqrt(3.)*dist)
+
+    def Gram_matrix(self, F, F1, F2, lower, upper):
+        """
+        Return the Gram matrix of the vector of functions F with respect to the
+        RKHS norm. The use of this function is limited to input_dim=1.
+
+        :param F: vector of functions
+        :type F: np.array
+        :param F1: vector of derivatives of F
+        :type F1: np.array
+        :param F2: vector of second derivatives of F
+        :type F2: np.array
+        :param lower,upper: boundaries of the input domain
+        :type lower,upper: floats
+        """
+        assert self.input_dim == 1
+        def L(x, i):
+            return(3. / self.lengthscale ** 2 * F[i](x) + 2 * np.sqrt(3) / self.lengthscale * F1[i](x) + F2[i](x))
+        n = F.shape[0]
+        G = np.zeros((n, n))
+        for i in range(n):
+            for j in range(i, n):
+                G[i, j] = G[j, i] = integrate.quad(lambda x : L(x, i) * L(x, j), lower, upper)[0]
+        Flower = np.array([f(lower) for f in F])[:, None]
+        F1lower = np.array([f(lower) for f in F1])[:, None]
+        return(self.lengthscale ** 3 / (12.*np.sqrt(3) * self.variance) * G + 1. / self.variance * np.dot(Flower, Flower.T) + self.lengthscale ** 2 / (3.*self.variance) * np.dot(F1lower, F1lower.T))
+
+
+class Matern52(Stationary):
+    """
+    Matern 5/2 kernel:
+
+    .. math::
+
+       k(r) = \sigma^2 (1 + \sqrt{5} r + \\frac53 r^2) \exp(- \sqrt{5} r) \ \ \ \ \  \\text{ where  } r = \sqrt{\sum_{i=1}^input_dim \\frac{(x_i-y_i)^2}{\ell_i^2} }
+       """
+
+    def K(self, X, X2=None):
+        r = self._scaled_dist(X, X2)
+        return self.variance*(1+np.sqrt(5.)*r+5./3*r**2)*np.exp(-np.sqrt(5.)*r)
+
+    def dK_dr(self, X, X2):
+        r = self._scaled_dist(X, X2)
+        return self.variance*(10./3*r -5.*r -5.*np.sqrt(5.)/3*r**2)*np.exp(-np.sqrt(5.)*r)
+
+    def Gram_matrix(self,F,F1,F2,F3,lower,upper):
+        """
+        Return the Gram matrix of the vector of functions F with respect to the RKHS norm. The use of this function is limited to input_dim=1.
+
+        :param F: vector of functions
+        :type F: np.array
+        :param F1: vector of derivatives of F
+        :type F1: np.array
+        :param F2: vector of second derivatives of F
+        :type F2: np.array
+        :param F3: vector of third derivatives of F
+        :type F3: np.array
+        :param lower,upper: boundaries of the input domain
+        :type lower,upper: floats
+        """
+        assert self.input_dim == 1
+        def L(x,i):
+            return(5*np.sqrt(5)/self.lengthscale**3*F[i](x) + 15./self.lengthscale**2*F1[i](x)+ 3*np.sqrt(5)/self.lengthscale*F2[i](x) + F3[i](x))
+        n = F.shape[0]
+        G = np.zeros((n,n))
+        for i in range(n):
+            for j in range(i,n):
+                G[i,j] = G[j,i] = integrate.quad(lambda x : L(x,i)*L(x,j),lower,upper)[0]
+        G_coef = 3.*self.lengthscale**5/(400*np.sqrt(5))
+        Flower = np.array([f(lower) for f in F])[:,None]
+        F1lower = np.array([f(lower) for f in F1])[:,None]
+        F2lower = np.array([f(lower) for f in F2])[:,None]
+        orig = 9./8*np.dot(Flower,Flower.T) + 9.*self.lengthscale**4/200*np.dot(F2lower,F2lower.T)
+        orig2 = 3./5*self.lengthscale**2 * ( np.dot(F1lower,F1lower.T) + 1./8*np.dot(Flower,F2lower.T) + 1./8*np.dot(F2lower,Flower.T))
+        return(1./self.variance* (G_coef*G + orig + orig2))
+
+
+
+
+class ExpQuad(Stationary):
+    def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, name='ExpQuad'):
+        super(ExpQuad, self).__init__(input_dim, variance, lengthscale, ARD, name)
+
+    def K(self, X, X2=None):
+        r = self._scaled_dist(X, X2)
+        return self.variance * np.exp(-0.5 * r**2)
+
+    def dK_dr(self, X, X2):
+        dist = self._scaled_dist(X, X2)
+        return -dist*self.K(X, X2)
+
+
+
diff --git a/GPy/plotting/matplot_dep/dim_reduction_plots.py b/GPy/plotting/matplot_dep/dim_reduction_plots.py
index 74292c05..3f4ea9b0 100644
--- a/GPy/plotting/matplot_dep/dim_reduction_plots.py
+++ b/GPy/plotting/matplot_dep/dim_reduction_plots.py
@@ -1,8 +1,8 @@
 import pylab as pb
 import numpy as np
-from ... import util
 from latent_space_visualizations.controllers.imshow_controller import ImshowController,ImAnnotateController
-from GPy.util.misc import param_to_array
+from ...util.misc import param_to_array
+from .base_plots import x_frame2D
 import itertools
 import Tango
 from matplotlib.cm import get_cmap
@@ -37,7 +37,7 @@ def plot_latent(model, labels=None, which_indices=None,
     if ax is None:
         fig = pb.figure(num=fignum)
         ax = fig.add_subplot(111)
-    util.plot.Tango.reset()
+    Tango.reset()
 
     if labels is None:
         labels = np.ones(model.num_data)
@@ -46,7 +46,7 @@ def plot_latent(model, labels=None, which_indices=None,
     X = param_to_array(model.X)
 
     # first, plot the output variance as a function of the latent space
-    Xtest, xx, yy, xmin, xmax = util.plot.x_frame2D(X[:, [input_1, input_2]], resolution=resolution)
+    Xtest, xx, yy, xmin, xmax = x_frame2D(X[:, [input_1, input_2]], resolution=resolution)
     Xtest_full = np.zeros((Xtest.shape[0], model.X.shape[1]))
 
     def plot_function(x):
@@ -87,7 +87,7 @@ def plot_latent(model, labels=None, which_indices=None,
         else:
             x = X[index, input_1]
             y = X[index, input_2]
-        ax.scatter(x, y, marker=m, s=s, color=util.plot.Tango.nextMedium(), label=this_label)
+        ax.scatter(x, y, marker=m, s=s, color=Tango.nextMedium(), label=this_label)
 
     ax.set_xlabel('latent dimension %i' % input_1)
     ax.set_ylabel('latent dimension %i' % input_2)
@@ -120,7 +120,7 @@ def plot_magnification(model, labels=None, which_indices=None,
     if ax is None:
         fig = pb.figure(num=fignum)
         ax = fig.add_subplot(111)
-    util.plot.Tango.reset()
+    Tango.reset()
 
     if labels is None:
         labels = np.ones(model.num_data)
@@ -128,7 +128,7 @@ def plot_magnification(model, labels=None, which_indices=None,
     input_1, input_2 = most_significant_input_dimensions(model, which_indices)
 
     # first, plot the output variance as a function of the latent space
-    Xtest, xx, yy, xmin, xmax = util.plot.x_frame2D(model.X[:, [input_1, input_2]], resolution=resolution)
+    Xtest, xx, yy, xmin, xmax = x_frame2D(model.X[:, [input_1, input_2]], resolution=resolution)
     Xtest_full = np.zeros((Xtest.shape[0], model.X.shape[1]))
 
     def plot_function(x):
@@ -165,7 +165,7 @@ def plot_magnification(model, labels=None, which_indices=None,
         else:
             x = model.X[index, input_1]
             y = model.X[index, input_2]
-        ax.scatter(x, y, marker=m, s=s, color=util.plot.Tango.nextMedium(), label=this_label)
+        ax.scatter(x, y, marker=m, s=s, color=Tango.nextMedium(), label=this_label)
 
     ax.set_xlabel('latent dimension %i' % input_1)
     ax.set_ylabel('latent dimension %i' % input_2)
@@ -205,7 +205,7 @@ def plot_steepest_gradient_map(model, fignum=None, ax=None, which_indices=None,
         return dmu_dX[indices, argmax], np.array(labels)[argmax]
 
     if ax is None:
-        fig = pyplot.figure(num=fignum)
+        fig = pb.figure(num=fignum)
         ax = fig.add_subplot(111)
 
     if data_labels is None:
@@ -241,7 +241,7 @@ def plot_steepest_gradient_map(model, fignum=None, ax=None, which_indices=None,
     ax.legend()
     ax.figure.tight_layout()
     if updates:
-        pyplot.show()
+        pb.show()
         clear = raw_input('Enter to continue')
         if clear.lower() in 'yes' or clear == '':
             controller.deactivate()
diff --git a/GPy/plotting/matplot_dep/kernel_plots.py b/GPy/plotting/matplot_dep/kernel_plots.py
index 30157294..3436c4ff 100644
--- a/GPy/plotting/matplot_dep/kernel_plots.py
+++ b/GPy/plotting/matplot_dep/kernel_plots.py
@@ -1,7 +1,6 @@
 # Copyright (c) 2012, GPy authors (see AUTHORS.txt).
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
-import sys
 import numpy as np
 import pylab as pb
 import Tango
@@ -29,22 +28,23 @@ def plot_ARD(kernel, fignum=None, ax=None, title='', legend=False):
     xticklabels = []
     bars = []
     x0 = 0
-    for p in kernel._parameters_:
-        c = Tango.nextMedium()
-        if hasattr(p, 'ARD') and p.ARD:
-            if title is None:
-                ax.set_title('ARD parameters, %s kernel' % p.name)
-            else:
-                ax.set_title(title)
-            if isinstance(p, Linear):
-                ard_params = p.variances
-            else:
-                ard_params = 1. / p.lengthscale
-
-            x = np.arange(x0, x0 + len(ard_params))
-            bars.append(ax.bar(x, ard_params, align='center', color=c, edgecolor='k', linewidth=1.2, label=p.name.replace("_"," ")))
-            xticklabels.extend([r"$\mathrm{{{name}}}\ {x}$".format(name=p.name, x=i) for i in np.arange(len(ard_params))])
-            x0 += len(ard_params)
+    #for p in kernel._parameters_:
+    p = kernel
+    c = Tango.nextMedium()
+    if hasattr(p, 'ARD') and p.ARD:
+        if title is None:
+            ax.set_title('ARD parameters, %s kernel' % p.name)
+        else:
+            ax.set_title(title)
+        if isinstance(p, Linear):
+            ard_params = p.variances
+        else:
+            ard_params = 1. / p.lengthscale
+        x = np.arange(x0, x0 + len(ard_params))
+        from ...util.misc import param_to_array
+        bars.append(ax.bar(x, param_to_array(ard_params), align='center', color=c, edgecolor='k', linewidth=1.2, label=p.name.replace("_"," ")))
+        xticklabels.extend([r"$\mathrm{{{name}}}\ {x}$".format(name=p.name, x=i) for i in np.arange(len(ard_params))])
+        x0 += len(ard_params)
     x = np.arange(x0)
     transOffset = offset_copy(ax.transData, fig=fig,
                               x=0., y= -2., units='points')
diff --git a/GPy/plotting/matplot_dep/models_plots.py b/GPy/plotting/matplot_dep/models_plots.py
index 47c8642e..59c32775 100644
--- a/GPy/plotting/matplot_dep/models_plots.py
+++ b/GPy/plotting/matplot_dep/models_plots.py
@@ -56,7 +56,10 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
     if ax is None:
         fig = pb.figure(num=fignum)
         ax = fig.add_subplot(111)
-
+    
+    X, Y = param_to_array(model.X, model.Y)
+    if model.has_uncertain_inputs(): X_variance = model.X_variance
+    
     #work out what the inputs are for plotting (1D or 2D)
     fixed_dims = np.array([i for i,v in fixed_inputs])
     free_dims = np.setdiff1d(np.arange(model.input_dim),fixed_dims)
@@ -66,7 +69,7 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
 
         #define the frame on which to plot
         resolution = resolution or 200
-        Xnew, xmin, xmax = x_frame1D(model.X[:,free_dims], plot_limits=plot_limits)
+        Xnew, xmin, xmax = x_frame1D(X[:,free_dims], plot_limits=plot_limits)
         Xgrid = np.empty((Xnew.shape[0],model.input_dim))
         Xgrid[:,free_dims] = Xnew
         for i,v in fixed_inputs:
@@ -77,13 +80,13 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
             m, v = model._raw_predict(Xgrid)
             lower = m - 2*np.sqrt(v)
             upper = m + 2*np.sqrt(v)
-            Y = model.Y
+            Y = Y
         else:
             m, v, lower, upper = model.predict(Xgrid)
-            Y = model.Y
+            Y = Y
         for d in which_data_ycols:
             gpplot(Xnew, m[:, d], lower[:, d], upper[:, d], axes=ax, edgecol=linecol, fillcol=fillcol)
-            ax.plot(model.X[which_data_rows,free_dims], Y[which_data_rows, d], 'kx', mew=1.5)
+            ax.plot(X[which_data_rows,free_dims], Y[which_data_rows, d], 'kx', mew=1.5)
 
         #optionally plot some samples
         if samples: #NOTE not tested with fixed_inputs
@@ -95,8 +98,8 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
         
         #add error bars for uncertain (if input uncertainty is being modelled)
         if hasattr(model,"has_uncertain_inputs") and model.has_uncertain_inputs():
-            ax.errorbar(model.X[which_data_rows, free_dims], model.Y[which_data_rows, which_data_ycols],
-                        xerr=2 * np.sqrt(model.X_variance[which_data_rows, free_dims]),
+            ax.errorbar(X[which_data_rows, free_dims].flatten(), Y[which_data_rows, which_data_ycols].flatten(),
+                        xerr=2 * np.sqrt(X_variance[which_data_rows, free_dims].flatten()),
                         ecolor='k', fmt=None, elinewidth=.5, alpha=.5)
 
 
@@ -120,7 +123,7 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
 
         #define the frame for plotting on
         resolution = resolution or 50
-        Xnew, _, _, xmin, xmax = x_frame2D(model.X[:,free_dims], plot_limits, resolution)
+        Xnew, _, _, xmin, xmax = x_frame2D(X[:,free_dims], plot_limits, resolution)
         Xgrid = np.empty((Xnew.shape[0],model.input_dim))
         Xgrid[:,free_dims] = Xnew
         for i,v in fixed_inputs:
@@ -130,14 +133,14 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
         #predict on the frame and plot
         if plot_raw:
             m, _ = model._raw_predict(Xgrid)
-            Y = model.Y
+            Y = Y
         else:
             m, _, _, _ = model.predict(Xgrid)
             Y = model.data
         for d in which_data_ycols:
             m_d = m[:,d].reshape(resolution, resolution).T
             ax.contour(x, y, m_d, levels, vmin=m.min(), vmax=m.max(), cmap=pb.cm.jet)
-            ax.scatter(model.X[which_data_rows, free_dims[0]], model.X[which_data_rows, free_dims[1]], 40, Y[which_data_rows, d], cmap=pb.cm.jet, vmin=m.min(), vmax=m.max(), linewidth=0.)
+            ax.scatter(X[which_data_rows, free_dims[0]], X[which_data_rows, free_dims[1]], 40, Y[which_data_rows, d], cmap=pb.cm.jet, vmin=m.min(), vmax=m.max(), linewidth=0.)
 
         #set the limits of the plot to some sensible values
         ax.set_xlim(xmin[0], xmax[0])
diff --git a/GPy/testing/index_operations_tests.py b/GPy/testing/index_operations_tests.py
index d5ef7007..171db5cc 100644
--- a/GPy/testing/index_operations_tests.py
+++ b/GPy/testing/index_operations_tests.py
@@ -24,6 +24,13 @@ class Test(unittest.TestCase):
         self.param_index.remove(one, [1])
         self.assertListEqual(self.param_index[one].tolist(), [3])        
 
+    def test_shift_left(self):
+        self.param_index.shift_left(1, 2)
+        self.assertListEqual(self.param_index[three].tolist(), [2,5])
+        self.assertListEqual(self.param_index[two].tolist(), [0,3])
+        self.assertListEqual(self.param_index[one].tolist(), [1])        
+
+
     def test_index_view(self):
         #=======================================================================
         #          0    1    2    3    4    5    6    7    8    9
diff --git a/GPy/testing/parameterized_tests.py b/GPy/testing/parameterized_tests.py
index ff57606a..6f13d294 100644
--- a/GPy/testing/parameterized_tests.py
+++ b/GPy/testing/parameterized_tests.py
@@ -10,8 +10,8 @@ import numpy as np
 class Test(unittest.TestCase):
 
     def setUp(self):
-        self.rbf = GPy.kern.rbf(1)
-        self.white = GPy.kern.white(1)
+        self.rbf = GPy.kern.RBF(1)
+        self.white = GPy.kern.White(1)
         from GPy.core.parameterization import Param
         from GPy.core.parameterization.transformations import Logistic
         self.param = Param('param', np.random.rand(25,2), Logistic(0, 1))
@@ -39,14 +39,13 @@ class Test(unittest.TestCase):
         
         
     def test_remove_parameter(self):
-        from GPy.core.parameterization.transformations import FIXED, UNFIXED, __fixed__
+        from GPy.core.parameterization.transformations import FIXED, UNFIXED, __fixed__, Logexp
         self.white.fix()
         self.test1.remove_parameter(self.white)
         self.assertIs(self.test1._fixes_,None)
         
         self.assertListEqual(self.white._fixes_.tolist(), [FIXED])
-        self.assertIs(self.white.constraints,self.white.white.constraints._param_index_ops)
-        self.assertEquals(self.white.white.constraints._offset, 0)
+        self.assertEquals(self.white.constraints._offset, 0)
         self.assertIs(self.test1.constraints, self.rbf.constraints._param_index_ops)
         self.assertIs(self.test1.constraints, self.param.constraints._param_index_ops)        
         
@@ -57,18 +56,19 @@ class Test(unittest.TestCase):
         self.assertListEqual(self.test1.constraints[__fixed__].tolist(), [0])
         self.assertIs(self.white._fixes_,None)
         self.assertListEqual(self.test1._fixes_.tolist(),[FIXED] + [UNFIXED] * 52)
+        
         self.test1.remove_parameter(self.white)
         self.assertIs(self.test1._fixes_,None)
         self.assertListEqual(self.white._fixes_.tolist(), [FIXED])
-        self.assertIs(self.white.constraints,self.white.white.constraints._param_index_ops)
         self.assertIs(self.test1.constraints, self.rbf.constraints._param_index_ops)
-        self.assertIs(self.test1.constraints, self.param.constraints._param_index_ops)        
+        self.assertIs(self.test1.constraints, self.param.constraints._param_index_ops)
+        self.assertListEqual(self.test1.constraints[Logexp()].tolist(), [0,1])
         
     def test_add_parameter_already_in_hirarchy(self):
         self.test1.add_parameter(self.white._parameters_[0])
         
     def test_default_constraints(self):
-        self.assertIs(self.rbf.rbf.variance.constraints._param_index_ops, self.rbf.constraints._param_index_ops)
+        self.assertIs(self.rbf.variance.constraints._param_index_ops, self.rbf.constraints._param_index_ops)
         self.assertIs(self.test1.constraints, self.rbf.constraints._param_index_ops)
         self.assertListEqual(self.rbf.constraints.indices()[0].tolist(), range(2))
         from GPy.core.parameterization.transformations import Logexp
diff --git a/GPy/util/__init__.py b/GPy/util/__init__.py
index c10fea4c..f93bb0ec 100644
--- a/GPy/util/__init__.py
+++ b/GPy/util/__init__.py
@@ -12,6 +12,7 @@ import decorators
 import classification
 import subarray_and_sorting
 import caching
+import diag
 
 try:
     import sympy
diff --git a/GPy/util/diag.py b/GPy/util/diag.py
index 3d6b4dc9..3044ed54 100644
--- a/GPy/util/diag.py
+++ b/GPy/util/diag.py
@@ -11,14 +11,14 @@ import numpy as np
 def view(A, offset=0):
     """
     Get a view on the diagonal elements of a 2D array.
-    
-    This is actually a view (!) on the diagonal of the array, so you can 
+
+    This is actually a view (!) on the diagonal of the array, so you can
     in-place adjust the view.
-    
+
     :param :class:`ndarray` A: 2 dimensional numpy array
     :param int offset: view offset to give back (negative entries allowed)
     :rtype: :class:`ndarray` view of diag(A)
-    
+
     >>> import numpy as np
     >>> X = np.arange(9).reshape(3,3)
     >>> view(X)
@@ -36,7 +36,7 @@ def view(A, offset=0):
     """
     from numpy.lib.stride_tricks import as_strided
     assert A.ndim == 2, "only implemented for 2 dimensions"
-    assert A.shape[0] == A.shape[1], "attempting to get the view of non-square matrix?!" 
+    assert A.shape[0] == A.shape[1], "attempting to get the view of non-square matrix?!"
     if offset > 0:
         return as_strided(A[0, offset:], shape=(A.shape[0] - offset, ), strides=((A.shape[0]+1)*A.itemsize, ))
     elif offset < 0:
@@ -44,6 +44,12 @@ def view(A, offset=0):
     else:
         return as_strided(A, shape=(A.shape[0], ), strides=((A.shape[0]+1)*A.itemsize, ))
 
+def offdiag_view(A, offset=0):
+    from numpy.lib.stride_tricks import as_strided
+    assert A.ndim == 2, "only implemented for 2 dimensions"
+    Af = as_strided(A, shape=(A.size,), strides=(A.itemsize,))
+    return as_strided(Af[(1+offset):], shape=(A.shape[0]-1, A.shape[1]), strides=(A.strides[0] + A.itemsize, A.strides[1]))
+
 def _diag_ufunc(A,b,offset,func):
     dA = view(A, offset); func(dA,b,dA)
     return A
@@ -51,11 +57,11 @@ def _diag_ufunc(A,b,offset,func):
 def times(A, b, offset=0):
     """
     Times the view of A with b in place (!).
-    Returns modified A 
+    Returns modified A
     Broadcasting is allowed, thus b can be scalar.
-    
+
     if offset is not zero, make sure b is of right shape!
-    
+
     :param ndarray A: 2 dimensional array
     :param ndarray-like b: either one dimensional or scalar
     :param int offset: same as in view.
@@ -67,11 +73,11 @@ multiply = times
 def divide(A, b, offset=0):
     """
     Divide the view of A by b in place (!).
-    Returns modified A 
+    Returns modified A
     Broadcasting is allowed, thus b can be scalar.
-    
+
     if offset is not zero, make sure b is of right shape!
-    
+
     :param ndarray A: 2 dimensional array
     :param ndarray-like b: either one dimensional or scalar
     :param int offset: same as in view.
@@ -84,9 +90,9 @@ def add(A, b, offset=0):
     Add b to the view of A in place (!).
     Returns modified A.
     Broadcasting is allowed, thus b can be scalar.
-    
+
     if offset is not zero, make sure b is of right shape!
-    
+
     :param ndarray A: 2 dimensional array
     :param ndarray-like b: either one dimensional or scalar
     :param int offset: same as in view.
@@ -99,16 +105,16 @@ def subtract(A, b, offset=0):
     Subtract b from the view of A in place (!).
     Returns modified A.
     Broadcasting is allowed, thus b can be scalar.
-    
+
     if offset is not zero, make sure b is of right shape!
-    
+
     :param ndarray A: 2 dimensional array
     :param ndarray-like b: either one dimensional or scalar
     :param int offset: same as in view.
     :rtype: view of A, which is adjusted inplace
     """
     return _diag_ufunc(A, b, offset, np.subtract)
-        
+
 if __name__ == '__main__':
     import doctest
-    doctest.testmod()
\ No newline at end of file
+    doctest.testmod()