diff --git a/GPy/core/parameterised.py b/GPy/core/parameterised.py
index ab656f52..007f1b25 100644
--- a/GPy/core/parameterised.py
+++ b/GPy/core/parameterised.py
@@ -333,8 +333,7 @@ class parameterised(object):
         header_string = ["{h:^{col}}".format(h = header[i], col = cols[i]) for i in range(len(cols))]
         header_string = map(lambda x: '|'.join(x), [header_string])
         separator = '-'*len(header_string[0])
-        param_string = ["{n:^{c0}}|{v:^{c1}}|{c:^{c2}}|{t:^{c3}}".format(n = names[i], v = values[i], c = constraints[i], t = ties[i],
-                                                                        c0 = cols[0], c1 = cols[1], c2 = cols[2], c3 = cols[3]) for i in range(len(values))]
+        param_string = ["{n:^{c0}}|{v:^{c1}}|{c:^{c2}}|{t:^{c3}}".format(n = names[i], v = values[i], c = constraints[i], t = ties[i], c0 = cols[0], c1 = cols[1], c2 = cols[2], c3 = cols[3]) for i in range(len(values))]
 
 
         return ('\n'.join([header_string[0], separator]+param_string)) + '\n'
diff --git a/GPy/kern/__init__.py b/GPy/kern/__init__.py
index 17bc607e..008a2e1a 100644
--- a/GPy/kern/__init__.py
+++ b/GPy/kern/__init__.py
@@ -2,5 +2,5 @@
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
 
-from constructors import rbf, Matern32, Matern52, exponential, linear, white, bias, finite_dimensional, spline, Brownian, rbf_sympy, sympykern, periodic_exponential, periodic_Matern32, periodic_Matern52, product_orthogonal
+from constructors import rbf, Matern32, Matern52, exponential, linear, white, bias, finite_dimensional, spline, Brownian, rbf_sympy, sympykern, periodic_exponential, periodic_Matern32, periodic_Matern52, product, product_orthogonal
 from kern import kern
diff --git a/GPy/kern/constructors.py b/GPy/kern/constructors.py
index 1e76f559..af82bb24 100644
--- a/GPy/kern/constructors.py
+++ b/GPy/kern/constructors.py
@@ -18,6 +18,7 @@ from Brownian import Brownian as Brownianpart
 from periodic_exponential import periodic_exponential as periodic_exponentialpart
 from periodic_Matern32 import periodic_Matern32 as periodic_Matern32part
 from periodic_Matern52 import periodic_Matern52 as periodic_Matern52part
+from product import product as productpart
 from product_orthogonal import product_orthogonal as product_orthogonalpart
 #TODO these s=constructors are not as clean as we'd like. Tidy the code up
 #using meta-classes to make the objects construct properly wthout them.
@@ -242,10 +243,21 @@ def periodic_Matern52(D,variance=1., lengthscale=None, period=2*np.pi,n_freq=10,
     part = periodic_Matern52part(D,variance, lengthscale, period, n_freq, lower, upper)
     return kern(D, [part])
 
+def product(k1,k2):
+    """
+     Construct a product kernel over D from two kernels over D
+
+    :param k1, k2: the kernels to multiply
+    :type k1, k2: kernpart
+    :rtype: kernel object
+    """
+    part = productpart(k1,k2)
+    return kern(k1.D, [part])
+
 def product_orthogonal(k1,k2):
     """
-     Construct a product kernel
-     
+     Construct a product kernel over D1 x D2 from a kernel over D1 and another over D2.
+
     :param k1, k2: the kernels to multiply
     :type k1, k2: kernpart
     :rtype: kernel object
diff --git a/GPy/kern/kern.py b/GPy/kern/kern.py
index 3528dba0..de4ebc8a 100644
--- a/GPy/kern/kern.py
+++ b/GPy/kern/kern.py
@@ -8,6 +8,7 @@ from ..core.parameterised import parameterised
 from kernpart import kernpart
 import itertools
 from product_orthogonal import product_orthogonal
+from product import product
 
 class kern(parameterised):
     def __init__(self,D,parts=[], input_slices=None):
@@ -149,24 +150,38 @@ class kern(parameterised):
         """
         Shortcut for `prod_orthogonal`. Note that `+` assumes that we sum 2 kernels defines on the same space whereas `*` assumes that the kernels are defined on different subspaces.
         """
-        return self.prod_orthogonal(other)
+        return self.prod(other)
 
-    def prod_orthogonal(self,other):
+    def prod(self,other):
         """
-        multiply two kernels. Both kernels are defined on separate spaces. Note that the constrains on the parameters of the kernels to multiply will be lost.
+        multiply two kernels defined on the same spaces.
         :param other: the other kernel to be added
         :type other: GPy.kern
         """
         K1 = self.copy()
         K2 = other.copy()
-        K1.unconstrain('')
-        K2.unconstrain('')
 
-        prev_ties = K1.tied_indices + [arr + K1.Nparam for arr in K2.tied_indices]
-        K1.untie_everything()
-        K2.untie_everything()
+        newkernparts = [product(k1,k2) for k1, k2 in itertools.product(K1.parts,K2.parts)]
 
-        D = K1.D + K2.D
+        slices = []
+        for sl1, sl2 in itertools.product(K1.input_slices,K2.input_slices):
+            s1, s2 = [False]*K1.D, [False]*K2.D
+            s1[sl1], s2[sl2] = [True], [True]
+            slices += [s1+s2]
+
+        newkern = kern(K1.D, newkernparts, slices)
+        newkern._follow_constrains(K1,K2)
+
+        return newkern
+
+    def prod_orthogonal(self,other):
+        """
+        multiply two kernels. Both kernels are defined on separate spaces.
+        :param other: the other kernel to be added
+        :type other: GPy.kern
+        """
+        K1 = self.copy()
+        K2 = other.copy()
 
         newkernparts = [product_orthogonal(k1,k2) for k1, k2 in itertools.product(K1.parts,K2.parts)]
 
@@ -176,9 +191,14 @@ class kern(parameterised):
             s1[sl1], s2[sl2] = [True], [True]
             slices += [s1+s2]
 
-        newkern = kern(D, newkernparts, slices)
+        newkern = kern(K1.D + K2.D, newkernparts, slices)
+        newkern._follow_constrains(K1,K2)
 
-        # create the ties
+        return newkern
+
+    def _follow_constrains(self,K1,K2):
+
+        # Build the array that allows to go from the initial indices of the param to the new ones
         K1_param = []
         n = 0
         for k1 in K1.parts:
@@ -192,19 +212,40 @@ class kern(parameterised):
         index_param = []
         for p1 in K1_param:
             for p2 in K2_param:
-                index_param += [0] + p1[1:] + p2[1:]
+                index_param += p1 + p2
         index_param = np.array(index_param)
 
+        # Get the ties and constrains of the kernels before the multiplication
+        prev_ties = K1.tied_indices + [arr + K1.Nparam for arr in K2.tied_indices]
+
+        prev_constr_pos = np.append(K1.constrained_positive_indices, K1.Nparam + K2.constrained_positive_indices)
+        prev_constr_neg = np.append(K1.constrained_negative_indices, K1.Nparam + K2.constrained_negative_indices)
+
+        prev_constr_fix = K1.constrained_fixed_indices + [arr + K1.Nparam for arr in K2.constrained_fixed_indices]
+        prev_constr_fix_values = K1.constrained_fixed_values + K2.constrained_fixed_values
+
+        prev_constr_bou = K1.constrained_bounded_indices + [arr + K1.Nparam for arr in K2.constrained_bounded_indices]
+        prev_constr_bou_low = K1.constrained_bounded_lowers + K2.constrained_bounded_lowers
+        prev_constr_bou_upp = K1.constrained_bounded_uppers + K2.constrained_bounded_uppers
+
         # follow the previous ties
         for arr in prev_ties:
             for j in arr:
                 index_param[np.where(index_param==j)[0]] = arr[0]
 
-        # tie
-        for i in np.unique(index_param)[1:]:
-            newkern.tie_param(np.where(index_param==i)[0])
-
-        return newkern
+        # ties and constrains
+        for i in range(K1.Nparam + K2.Nparam):
+            index = np.where(index_param==i)[0]
+            if index.size > 1:
+                self.tie_param(index)
+        for i in prev_constr_pos:
+            self.constrain_positive(np.where(index_param==i)[0])
+        for i in prev_constr_neg:
+            self.constrain_neg(np.where(index_param==i)[0])
+        for j, i in enumerate(prev_constr_fix):
+            self.constrain_fixed(np.where(index_param==i)[0],prev_constr_fix_values[j])
+        for j, i in enumerate(prev_constr_bou):
+            self.constrain_bounded(np.where(index_param==i)[0],prev_constr_bou_low[j],prev_constr_bou_upp[j])
 
     def _get_params(self):
         return np.hstack([p._get_params() for p in self.parts])
@@ -446,7 +487,7 @@ class kern(parameterised):
             Xnew = np.vstack((xx.flatten(),yy.flatten())).T
             Kx = self.K(Xnew,x,slices2=which_functions)
             Kx = Kx.reshape(resolution,resolution).T
-            pb.contour(xg,yg,Kx,vmin=Kx.min(),vmax=Kx.max(),cmap=pb.cm.jet)
+            pb.contour(xg,yg,Kx,vmin=Kx.min(),vmax=Kx.max(),cmap=pb.cm.jet,*args,**kwargs)
             pb.xlim(xmin[0],xmax[0])
             pb.ylim(xmin[1],xmax[1])
             pb.xlabel("x1")
diff --git a/GPy/kern/periodic_Matern32.py b/GPy/kern/periodic_Matern32.py
index 3e81044c..eeadf1c8 100644
--- a/GPy/kern/periodic_Matern32.py
+++ b/GPy/kern/periodic_Matern32.py
@@ -21,13 +21,14 @@ class periodic_Matern32(kernpart):
     """
 
     def __init__(self,D=1,variance=1.,lengthscale=None,period=2*np.pi,n_freq=10,lower=0.,upper=4*np.pi):
-        assert D==1
+        assert D==1, "Periodic kernels are only defined for D=1"
         self.name = 'periodic_Mat32'
         self.D = D
         if lengthscale is not None:
-            assert lengthscale.shape==(self.D,)
+            lengthscale = np.asarray(lengthscale)
+            assert lengthscale.size == 1, "Wrong size: only one lengthscale needed"
         else:
-            lengthscale = np.ones(self.D)
+            lengthscale = np.ones(1)
         self.lower,self.upper = lower, upper
         self.Nparam = 3
         self.n_freq = n_freq
diff --git a/GPy/kern/periodic_Matern52.py b/GPy/kern/periodic_Matern52.py
index abdd6a7d..2db3d223 100644
--- a/GPy/kern/periodic_Matern52.py
+++ b/GPy/kern/periodic_Matern52.py
@@ -21,13 +21,14 @@ class periodic_Matern52(kernpart):
     """
 
     def __init__(self,D=1,variance=1.,lengthscale=None,period=2*np.pi,n_freq=10,lower=0.,upper=4*np.pi):
-        assert D==1
+        assert D==1, "Periodic kernels are only defined for D=1"
         self.name = 'periodic_Mat52'
         self.D = D
         if lengthscale is not None:
-            assert lengthscale.shape==(self.D,)
+            lengthscale = np.asarray(lengthscale)
+            assert lengthscale.size == 1, "Wrong size: only one lengthscale needed"
         else:
-            lengthscale = np.ones(self.D)
+            lengthscale = np.ones(1)
         self.lower,self.upper = lower, upper
         self.Nparam = 3
         self.n_freq = n_freq
diff --git a/GPy/kern/periodic_exponential.py b/GPy/kern/periodic_exponential.py
index 2637afe5..d99bada8 100644
--- a/GPy/kern/periodic_exponential.py
+++ b/GPy/kern/periodic_exponential.py
@@ -21,13 +21,14 @@ class periodic_exponential(kernpart):
     """
 
     def __init__(self,D=1,variance=1.,lengthscale=None,period=2*np.pi,n_freq=10,lower=0.,upper=4*np.pi):
-        assert D==1
+        assert D==1, "Periodic kernels are only defined for D=1"
         self.name = 'periodic_exp'
         self.D = D
         if lengthscale is not None:
-            assert lengthscale.shape==(self.D,)
+            lengthscale = np.asarray(lengthscale)
+            assert lengthscale.size == 1, "Wrong size: only one lengthscale needed"
         else:
-            lengthscale = np.ones(self.D)
+            lengthscale = np.ones(1)
         self.lower,self.upper = lower, upper
         self.Nparam = 3
         self.n_freq = n_freq
@@ -45,7 +46,7 @@ class periodic_exponential(kernpart):
         r =  np.sqrt(r1**2 + r2**2)
         psi = np.where(r1 != 0, (np.arctan(r2/r1) + (r1<0.)*np.pi),np.arcsin(r2))
         return r,omega[:,0:1], psi
-    
+
     def _int_computation(self,r1,omega1,phi1,r2,omega2,phi2):
         Gint1 = 1./(omega1+omega2.T)*( np.sin((omega1+omega2.T)*self.upper+phi1+phi2.T) - np.sin((omega1+omega2.T)*self.lower+phi1+phi2.T)) + 1./(omega1-omega2.T)*( np.sin((omega1-omega2.T)*self.upper+phi1-phi2.T) - np.sin((omega1-omega2.T)*self.lower+phi1-phi2.T) )
         Gint2 = 1./(omega1+omega2.T)*( np.sin((omega1+omega2.T)*self.upper+phi1+phi2.T) - np.sin((omega1+omega2.T)*self.lower+phi1+phi2.T)) +  np.cos(phi1-phi2.T)*(self.upper-self.lower)
diff --git a/GPy/kern/product.py b/GPy/kern/product.py
new file mode 100644
index 00000000..cefd7458
--- /dev/null
+++ b/GPy/kern/product.py
@@ -0,0 +1,86 @@
+# 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 # This may not be necessary (Nicolas, 20th Feb)
+
+class product(kernpart):
+    """
+    Computes the product of 2 kernels that are defined on the same space
+
+    :param k1, k2: the kernels to multiply
+    :type k1, k2: kernpart
+    :rtype: kernel object
+
+    """
+    def __init__(self,k1,k2):
+        assert k1.D == k2.D, "Error: The input spaces of the kernels to multiply must have the same dimension"
+        self.D = k1.D
+        self.Nparam = k1.Nparam + k2.Nparam
+        self.name = k1.name + '<times>' + k2.name
+        self.k1 = k1
+        self.k2 = k2
+        self._set_params(np.hstack((k1._get_params(),k2._get_params())))
+
+    def _get_params(self):
+        """return the value of the parameters."""
+        return self.params
+
+    def _set_params(self,x):
+        """set the value of the parameters."""
+        self.k1._set_params(x[:self.k1.Nparam])
+        self.k2._set_params(x[self.k1.Nparam:])
+        self.params = x
+
+    def _get_param_names(self):
+        """return parameter names."""
+        return [self.k1.name + '_' + param_name for param_name in self.k1._get_param_names()] + [self.k2.name + '_' + param_name for param_name in self.k2._get_param_names()]
+
+    def K(self,X,X2,target):
+        """Compute the covariance matrix between X and X2."""
+        if X2 is None: X2 = X
+        target1 = np.zeros((X.shape[0],X2.shape[0]))
+        target2 = np.zeros((X.shape[0],X2.shape[0]))
+        self.k1.K(X,X2,target1)
+        self.k2.K(X,X2,target2)
+        target += target1 * target2
+
+    def Kdiag(self,X,target):
+        """Compute the diagonal of the covariance matrix associated to X."""
+        target1 = np.zeros((X.shape[0],))
+        target2 = np.zeros((X.shape[0],))
+        self.k1.Kdiag(X,target1)
+        self.k2.Kdiag(X,target2)
+        target += target1 * target2
+
+    def dK_dtheta(self,partial,X,X2,target):
+        """derivative of the covariance matrix with respect to the parameters."""
+        if X2 is None: X2 = X
+        K1 = np.zeros((X.shape[0],X2.shape[0]))
+        K2 = np.zeros((X.shape[0],X2.shape[0]))
+        self.k1.K(X,X2,K1)
+        self.k2.K(X,X2,K2)
+
+        k1_target = np.zeros(self.k1.Nparam)
+        k2_target = np.zeros(self.k2.Nparam)
+        self.k1.dK_dtheta(partial*K2, X, X2, k1_target)
+        self.k2.dK_dtheta(partial*K1, X, X2, k2_target)
+
+        target[:self.k1.Nparam] += k1_target
+        target[self.k1.Nparam:] += k2_target
+
+    def dK_dX(self,partial,X,X2,target):
+        """derivative of the covariance matrix with respect to X."""
+        if X2 is None: X2 = X
+        K1 = np.zeros((X.shape[0],X2.shape[0]))
+        K2 = np.zeros((X.shape[0],X2.shape[0]))
+        self.k1.K(X,X2,K1)
+        self.k2.K(X,X2,K2)
+
+        self.k1.dK_dX(partial*K2, X, X2, target)
+        self.k2.dK_dX(partial*K1, X, X2, target)
+
+    def dKdiag_dX(self,X,target):
+        pass
diff --git a/GPy/kern/product_orthogonal.py b/GPy/kern/product_orthogonal.py
index 26456d8c..a729c126 100644
--- a/GPy/kern/product_orthogonal.py
+++ b/GPy/kern/product_orthogonal.py
@@ -4,7 +4,7 @@
 from kernpart import kernpart
 import numpy as np
 import hashlib
-from scipy import integrate
+#from scipy import integrate # This may not be necessary (Nicolas, 20th Feb)
 
 class product_orthogonal(kernpart):
     """
@@ -16,13 +16,12 @@ class product_orthogonal(kernpart):
 
     """
     def __init__(self,k1,k2):
-        assert k1._get_param_names()[0] == 'variance' and k2._get_param_names()[0] == 'variance', "Error: The multipication of kernels is only defined when the first parameters of the kernels to multiply is the variance."
         self.D = k1.D + k2.D
-        self.Nparam = k1.Nparam + k2.Nparam - 1
+        self.Nparam = k1.Nparam + k2.Nparam
         self.name = k1.name + '<times>' + k2.name
         self.k1 = k1
         self.k2 = k2
-        self._set_params(np.hstack((k1._get_params()[0]*k2._get_params()[0], k1._get_params()[1:],k2._get_params()[1:])))
+        self._set_params(np.hstack((k1._get_params(),k2._get_params())))
 
     def _get_params(self):
         """return the value of the parameters."""
@@ -30,14 +29,14 @@ class product_orthogonal(kernpart):
 
     def _set_params(self,x):
         """set the value of the parameters."""
-        self.k1._set_params(np.hstack((1.,x[1:self.k1.Nparam])))
-        self.k2._set_params(np.hstack((1.,x[self.k1.Nparam:])))
+        self.k1._set_params(x[:self.k1.Nparam])
+        self.k2._set_params(x[self.k1.Nparam:])
         self.params = x
 
     def _get_param_names(self):
         """return parameter names."""
-        return ['variance']+[self.k1.name + '_' + self.k1._get_param_names()[i+1] for i in range(self.k1.Nparam-1)] +  [self.k2.name + '_' + self.k2._get_param_names()[i+1] for i in range(self.k2.Nparam-1)]
-
+        return [self.k1.name + '_' + param_name for param_name in self.k1._get_param_names()] + [self.k2.name + '_' + param_name for param_name in self.k2._get_param_names()]
+    
     def K(self,X,X2,target):
         """Compute the covariance matrix between X and X2."""
         if X2 is None: X2 = X
@@ -45,7 +44,7 @@ class product_orthogonal(kernpart):
         target2 = np.zeros((X.shape[0],X2.shape[0]))
         self.k1.K(X[:,0:self.k1.D],X2[:,0:self.k1.D],target1)
         self.k2.K(X[:,self.k1.D:],X2[:,self.k1.D:],target2)
-        target += self.params[0]*target1 * target2
+        target += target1 * target2
 
     def Kdiag(self,X,target):
         """Compute the diagonal of the covariance matrix associated to X."""
@@ -53,7 +52,7 @@ class product_orthogonal(kernpart):
         target2 = np.zeros((X.shape[0],))
         self.k1.Kdiag(X[:,0:self.k1.D],target1)
         self.k2.Kdiag(X[:,self.k1.D:],target2)
-        target += self.params[0]*target1 * target2
+        target += target1 * target2
 
     def dK_dtheta(self,partial,X,X2,target):
         """derivative of the covariance matrix with respect to the parameters."""
@@ -68,13 +67,8 @@ class product_orthogonal(kernpart):
         self.k1.dK_dtheta(partial*K2, X[:,:self.k1.D], X2[:,:self.k1.D], k1_target)
         self.k2.dK_dtheta(partial*K1, X[:,self.k1.D:], X2[:,self.k1.D:], k2_target)
 
-        target[0] += np.sum(K1*K2*partial)
-        target[1:self.k1.Nparam] += self.params[0]* k1_target[1:]
-        target[self.k1.Nparam:] += self.params[0]* k2_target[1:]
-        
-    def dKdiag_dtheta(self,partial,X,target):
-        """derivative of the diagonal of the covariance matrix with respect to the parameters."""
-        target[0] += 1
+        target[:self.k1.Nparam] += k1_target
+        target[self.k1.Nparam:] += k2_target
 
     def dK_dX(self,partial,X,X2,target):
         """derivative of the covariance matrix with respect to X."""
diff --git a/doc/Figures/tuto_kern_overview_basicdef.png b/doc/Figures/tuto_kern_overview_basicdef.png
deleted file mode 100644
index bad43b09..00000000
Binary files a/doc/Figures/tuto_kern_overview_basicdef.png and /dev/null differ
diff --git a/doc/Figures/kern-def.png b/doc/Figures/tuto_kern_overview_basicplot.png
similarity index 100%
rename from doc/Figures/kern-def.png
rename to doc/Figures/tuto_kern_overview_basicplot.png
diff --git a/doc/Figures/tuto_kern_overview_mANOVA.png b/doc/Figures/tuto_kern_overview_mANOVA.png
index db49e3bd..f7ff5f50 100644
Binary files a/doc/Figures/tuto_kern_overview_mANOVA.png and b/doc/Figures/tuto_kern_overview_mANOVA.png differ
diff --git a/doc/Figures/tuto_kern_overview_mANOVAdec.png b/doc/Figures/tuto_kern_overview_mANOVAdec.png
index ef154263..6a1b7b73 100644
Binary files a/doc/Figures/tuto_kern_overview_mANOVAdec.png and b/doc/Figures/tuto_kern_overview_mANOVAdec.png differ
diff --git a/doc/Figures/tuto_kern_overview_multadd.png b/doc/Figures/tuto_kern_overview_multadd.png
new file mode 100644
index 00000000..a6aa888a
Binary files /dev/null and b/doc/Figures/tuto_kern_overview_multadd.png differ
diff --git a/doc/Figures/tuto_kern_overview_multperdecay.png b/doc/Figures/tuto_kern_overview_multperdecay.png
new file mode 100644
index 00000000..3970d342
Binary files /dev/null and b/doc/Figures/tuto_kern_overview_multperdecay.png differ
diff --git a/doc/GPy.inference.rst b/doc/GPy.inference.rst
index 357e70c7..f30e7d25 100644
--- a/doc/GPy.inference.rst
+++ b/doc/GPy.inference.rst
@@ -1,6 +1,14 @@
 inference Package
 =================
 
+:mod:`SGD` Module
+-----------------
+
+.. automodule:: GPy.inference.SGD
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
 :mod:`optimization` Module
 --------------------------
 
diff --git a/doc/GPy.kern.rst b/doc/GPy.kern.rst
index d6593939..a3a611b7 100644
--- a/doc/GPy.kern.rst
+++ b/doc/GPy.kern.rst
@@ -113,6 +113,14 @@ kern Package
     :undoc-members:
     :show-inheritance:
 
+:mod:`product` Module
+---------------------
+
+.. automodule:: GPy.kern.product
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
 :mod:`product_orthogonal` Module
 --------------------------------
 
diff --git a/doc/tuto_kernel_overview.rst b/doc/tuto_kernel_overview.rst
index 80e2bee2..a8f5b53d 100644
--- a/doc/tuto_kernel_overview.rst
+++ b/doc/tuto_kernel_overview.rst
@@ -13,20 +13,27 @@ First we import the libraries we will need ::
 For most kernels, the dimension is the only mandatory parameter to define a kernel object. However, it is also possible to specify the values of the parameters. For example, the three following commands are valid for defining a squared exponential kernel (ie rbf or Gaussian) ::
 
     ker1 = GPy.kern.rbf(1)  # Equivalent to ker1 = GPy.kern.rbf(D=1, variance=1., lengthscale=1.)
-    ker2 = GPy.kern.rbf(D=1, variance = 1.5, lengthscale=2.)
+    ker2 = GPy.kern.rbf(D=1, variance = .75, lengthscale=2.)
     ker3 = GPy.kern.rbf(1, .5, .5)
 
-A `plot` and a `print` functions are implemented to represent kernel objects ::
+A ``print`` and a ``plot`` functions are implemented to represent kernel objects. The commands ::
     
-    print ker1
+    print ker2
 
     ker1.plot()
     ker2.plot()
     ker3.plot()
 
-.. figure::  Figures/tuto_kern_overview_basicdef.png
+return::
+
+           Name        |  Value   |  Constraints  |  Ties  
+    -------------------------------------------------------
+       rbf_variance    |  0.7500  |               |        
+      rbf_lengthscale  |  2.0000  |               |        
+
+.. figure::  Figures/tuto_kern_overview_basicplot.png
     :align:   center
-    :height: 350px
+    :height: 300px
 
 Implemented kernels
 ===================
@@ -37,33 +44,131 @@ Many kernels are already implemented in GPy. Here is a summary of most of them:
     :align:  center
     :height: 800px
 
-On the other hand, it is possible to use the `sympy` package to build new kernels. This will be the subject of another tutorial.
+On the other hand, it is possible to use the `sympy` package to build new kernels. This will be the subject of another tutorial.    
 
-Operations to combine kernel
-============================
+Operations to combine kernels
+=============================
 
-In ``GPy``, kernel objects can be combined with the usual ``+`` and ``*`` operators. ::
-    
-    k1 = GPy.kern.rbf(1,variance=1., lengthscale=2)
-    k2 = GPy.kern.Matern32(1,variance=1., lengthscale=2)
+In ``GPy``, kernel objects can be added or multiplied. In both cases, two kinds of operations are possible since one can assume that the kernels to add/multiply are defined on the same space or on different subspaces. In other words, it is possible to use two kernels :math:`k_1,\ k_2` over :math:`\mathbb{R} \times \mathbb{R}` to create 
 
-    ker_add = k1 + k2
-    print ker_add
+    * a kernel over :math:`\mathbb{R} \times \mathbb{R}`:  :math:`k(x,y) = k_1(x,y) \times k_2(x,y)`
+    * a kernel over :math:`\mathbb{R}^2 \times \mathbb{R}^2`:  :math:`k(\mathbf{x},\mathbf{y}) = k_1(x_1,y_1) \times k_2(x_2,y_2)`
 
-    ker_prod = k1 * k2
-    print ker_prod
+These two options are available in GPy under the name ``prod`` and ``prod_orthogonal`` (resp. ``add`` and ``add_orthogonal`` for the addition). Here is a quick example ::
 
-Note that by default, the operator ``+`` adds kernels defined on the same input space whereas ``*`` assumes that the kernels are defined on different input spaces. Here for example ``ker_add.D`` will return ``1`` whereas ``ker_prod.D`` will return ``2``.
+    k1 = GPy.kern.rbf(1,1.,2.)
+    k2 = GPy.kern.Matern32(1, 0.5, 0.2)
 
-In order to add kernels defined on the different input spaces, the required command is::
+    # Product of kernels
+    k_prod = k1.prod(k2)
+    k_prodorth = k1.prod_orthogonal(k2)
 
-    ker_add_orth = k1.add_orthogonal(k2)
+    # Sum of kernels
+    k_add = k1.add(k2)
+    k_addorth = k1.add_orthogonal(k2)    
 
-.. figure::  Figures/tuto_kern_overview_add_orth.png
+..  # plots
+    pb.figure(figsize=(8,8))
+    pb.subplot(2,2,1)
+    k_prod.plot()
+    pb.title('prod')
+    pb.subplot(2,2,2)
+    k_prodorth.plot()
+    pb.title('prod_orthogonal')
+    pb.subplot(2,2,3)
+    k_add.plot()
+    pb.title('add')
+    pb.subplot(2,2,4)
+    k_addorth.plot()
+    pb.title('add_orthogonal')
+    pb.subplots_adjust(wspace=0.3, hspace=0.3)
+
+.. figure::  Figures/tuto_kern_overview_multadd.png
     :align:  center
-    :height: 350px
+    :height: 500px
 
-    Output of ``ker_add_orth.plot(plot_limits=[[-10,-10],[10,10]])``.
+A shortcut for ``add`` and ``prod`` is provided by the usual ``+`` and ``*`` operators. Here is another example where we create a periodic kernel with some decay ::
+    
+    k1 = GPy.kern.rbf(1,1.,2)
+    k2 = GPy.kern.periodic_Matern52(1,variance=1e3, lengthscale=1, period = 1.5, lower=-5., upper = 5)
+
+    k = k1 * k2  # equivalent to k = k1.prod(k2)
+    print k
+
+    # Simulate sample paths
+    X = np.linspace(-5,5,501)[:,None]
+    Y = np.random.multivariate_normal(np.zeros(501),k.K(X),1)
+
+..  # plot
+    pb.figure(figsize=(10,4))
+    pb.subplot(1,2,1)
+    k.plot()
+    pb.subplot(1,2,2)
+    pb.plot(X,Y.T)
+    pb.ylabel("Sample path")
+    pb.subplots_adjust(wspace=0.3)
+
+.. figure::  Figures/tuto_kern_overview_multperdecay.png
+    :align:  center
+    :height: 300px
+
+In general, ``kern`` objects can be seen as a sum of ``kernparts`` objects, where the later are covariance functions denied on the same space. For example, the following code ::
+
+    k = (k1+k2)*(k1+k2)
+    print k.parts[0].name, '\n', k.parts[1].name, '\n', k.parts[2].name, '\n', k.parts[3].name
+
+returns ::
+
+    rbf<times>rbf 
+    rbf<times>periodic_Mat52 
+    periodic_Mat52<times>rbf 
+    periodic_Mat52<times>periodic_Mat52
+
+Constraining the parameters
+===========================
+
+Various constrains can be applied to the parameters of a kernel
+
+    * ``constrain_fixed`` to fix the value of a parameter (the value will not be modified during optimisation)
+    * ``constrain_positive`` to make sure the parameter is greater than 0.
+    * ``constrain_bounded`` to impose the parameter to be in a given range.
+    * ``tie_param`` to impose the value of two (or more) parameters to be equal.
+
+When calling one of these functions, the parameters to constrain can either by specified by a regular expression that matches its name or by a number that corresponds to the rank of the parameter. Here is an example ::
+
+    k1 = GPy.kern.rbf(1)
+    k2 = GPy.kern.Matern32(1)
+    k3 = GPy.kern.white(1)
+
+    k = k1 + k2 + k3
+    print k
+
+    k.constrain_positive('var')
+    k.constrain_fixed(np.array([1]),1.75)
+    k.tie_param('len')
+    k.unconstrain('white')
+    k.constrain_bounded('white',lower=1e-5,upper=.5)
+    print k
+    
+with output::
+
+            Name         |  Value   |  Constraints  |  Ties  
+    ---------------------------------------------------------
+        rbf_variance     |  1.0000  |               |        
+       rbf_lengthscale   |  1.0000  |               |        
+       Mat32_variance    |  1.0000  |               |        
+      Mat32_lengthscale  |  1.0000  |               |        
+       white_variance    |  1.0000  |               |        
+    
+    
+            Name         |  Value   |  Constraints   |  Ties  
+    ----------------------------------------------------------
+        rbf_variance     |  1.0000  |     (+ve)      |        
+       rbf_lengthscale   |  1.7500  |     Fixed      |  (0)   
+       Mat32_variance    |  1.0000  |     (+ve)      |        
+      Mat32_lengthscale  |  1.7500  |                |  (0)   
+       white_variance    |  0.3655  |  (1e-05, 0.5)  |        
+  
 
 Example : Building an ANOVA kernel
 ==================================
@@ -78,23 +183,27 @@ Let us assume that we want to define an ANOVA kernel with a Matern 3/2 kernel fo
 
     k_cst = GPy.kern.bias(1,variance=1.)
     k_mat = GPy.kern.Matern52(1,variance=1., lengthscale=3)
-    Kanova = (k_cst + k_mat) * (k_cst + k_mat)
+    Kanova = (k_cst + k_mat).prod_orthogonal(k_cst + k_mat)
     print Kanova
 
 Printing the resulting kernel outputs the following ::
 
                      Name                  |  Value   |  Constraints  |  Ties  
     ---------------------------------------------------------------------------
-           bias<times>bias_variance        |  1.0000  |               |        
-           bias<times>Mat52_variance       |  1.0000  |               |        
-      bias<times>Mat52_Mat52_lengthscale   |  3.0000  |               |  (1)   
-           Mat52<times>bias_variance       |  1.0000  |               |        
-      Mat52<times>bias_Mat52_lengthscale   |  3.0000  |               |  (0)   
-          Mat52<times>Mat52_variance       |  1.0000  |               |        
-      Mat52<times>Mat52_Mat52_lengthscale  |  3.0000  |               |  (0)   
-      Mat52<times>Mat52_Mat52_lengthscale  |  3.0000  |               |  (1)
+         bias<times>bias_bias_variance     |  1.0000  |               |  (0)   
+         bias<times>bias_bias_variance     |  1.0000  |               |  (3)   
+        bias<times>Mat52_bias_variance     |  1.0000  |               |  (0)   
+        bias<times>Mat52_Mat52_variance    |  1.0000  |               |  (4)   
+      bias<times>Mat52_Mat52_lengthscale   |  3.0000  |               |  (5)   
+        Mat52<times>bias_Mat52_variance    |  1.0000  |               |  (1)   
+      Mat52<times>bias_Mat52_lengthscale   |  3.0000  |               |  (2)   
+        Mat52<times>bias_bias_variance     |  1.0000  |               |  (3)   
+       Mat52<times>Mat52_Mat52_variance    |  1.0000  |               |  (1)   
+      Mat52<times>Mat52_Mat52_lengthscale  |  3.0000  |               |  (2)   
+       Mat52<times>Mat52_Mat52_variance    |  1.0000  |               |  (4)   
+      Mat52<times>Mat52_Mat52_lengthscale  |  3.0000  |               |  (5)   
 
-Note the ties between the lengthscales of ``Kanova`` to keep the number of lengthscales equal to 2. On the other hand, there are four variance terms in the new parameterization: one for each term of the right hand part of the above equation. We can illustrate the use of this kernel on a toy example::
+Note the ties between the parameters of ``Kanova`` that reflect the links between the parameters of the kernparts objects. We can illustrate the use of this kernel on a toy example::
 
     # sample inputs and outputs
     X = np.random.uniform(-3.,3.,(40,2))
@@ -102,12 +211,12 @@ Note the ties between the lengthscales of ``Kanova`` to keep the number of lengt
 
     # Create GP regression model
     m = GPy.models.GP_regression(X,Y,Kanova)
+    pb.figure(figsize=(5,5))
     m.plot()
 
-
 .. figure::  Figures/tuto_kern_overview_mANOVA.png
     :align:  center
-    :height: 350px
+    :height: 300px
 
 As :math:`k_{ANOVA}` corresponds to the sum of 4 kernels, the best predictor can be splited in a sum of 4 functions 
 
@@ -119,7 +228,7 @@ As :math:`k_{ANOVA}` corresponds to the sum of 4 kernels, the best predictor can
 
 The submodels can be represented with the option ``which_function`` of ``plot``: ::
     
-    pb.figure(figsize=(20,5))
+    pb.figure(figsize=(20,3))
     pb.subplots_adjust(wspace=0.5)
     pb.subplot(1,5,1)
     m.plot()
@@ -135,10 +244,11 @@ The submodels can be represented with the option ``which_function`` of ``plot``:
     pb.ylabel("+   ",rotation='horizontal',fontsize='30')
     m.plot(which_functions=[False,False,False,True])
 
+..  pb.savefig('tuto_kern_overview_mANOVAdec.png',bbox_inches='tight')
 
 .. figure::  Figures/tuto_kern_overview_mANOVAdec.png
     :align:  center
-    :height: 200px
+    :height: 250px
 
 
 ..  import pylab as pb