diff --git a/GPy/_models/gradient_checker.py b/GPy/_models/gradient_checker.py
index 64b8b2fb..dfd0640f 100644
--- a/GPy/_models/gradient_checker.py
+++ b/GPy/_models/gradient_checker.py
@@ -28,38 +28,37 @@ class GradientChecker(Model):
         :param df: Gradient of function to check
         :param x0:
             Initial guess for inputs x (if it has a shape (a,b) this will be reflected in the parameter names).
-            Can be a list of arrays, if takes a list of arrays. This list will be passed
+            Can be a list of arrays, if f takes a list of arrays. This list will be passed
             to f and df in the same order as given here.
-            If only one argument, make sure not to pass a list!!!
-
+            If f takes only one argument, make sure not to pass a list for x0!!!
         :type x0: [array-like] | array-like | float | int
-        :param names:
+        :param list names:
             Names to print, when performing gradcheck. If a list was passed to x0
             a list of names with the same length is expected.
-        :param args: Arguments passed as f(x, *args, **kwargs) and df(x, *args, **kwargs)
+        :param args kwargs: Arguments passed as f(x, *args, **kwargs) and df(x, *args, **kwargs)
 
         Examples:
         ---------
-            from GPy.models import GradientChecker
-            N, M, Q = 10, 5, 3
+        from GPy.models import GradientChecker
+        N, M, Q = 10, 5, 3
 
-            Sinusoid:
+        Sinusoid:
 
-                X = numpy.random.rand(N, Q)
-                grad = GradientChecker(numpy.sin,numpy.cos,X,'x')
-                grad.checkgrad(verbose=1)
+            X = numpy.random.rand(N, Q)
+            grad = GradientChecker(numpy.sin,numpy.cos,X,'sin_in')
+            grad.checkgrad(verbose=1)
 
-            Using GPy:
+        Using GPy:
 
-                X, Z = numpy.random.randn(N,Q), numpy.random.randn(M,Q)
-                kern = GPy.kern.linear(Q, ARD=True) + GPy.kern.rbf(Q, ARD=True)
-                grad = GradientChecker(kern.K,
-                                       lambda x: 2*kern.dK_dX(numpy.ones((1,1)), x),
-                                       x0 = X.copy(),
-                                       names='X')
-                grad.checkgrad(verbose=1)
-                grad.randomize()
-                grad.checkgrad(verbose=1)
+            X, Z = numpy.random.randn(N,Q), numpy.random.randn(M,Q)
+            kern = GPy.kern.linear(Q, ARD=True) + GPy.kern.rbf(Q, ARD=True)
+            grad = GradientChecker(kern.K,
+                                   lambda x: kern.dK_dX(numpy.ones((1,1)), x),
+                                   x0 = X.copy(),
+                                   names=['X_input'])
+            grad.checkgrad(verbose=1)
+            grad.randomize()
+            grad.checkgrad(verbose=1)
         """
         Model.__init__(self)
         if isinstance(x0, (list, tuple)) and names is None:
diff --git a/GPy/kern/kern.py b/GPy/kern/kern.py
index 46bb01c8..949df5ab 100644
--- a/GPy/kern/kern.py
+++ b/GPy/kern/kern.py
@@ -487,12 +487,11 @@ class kern(Parameterized):
                 p1.psi1(Z, mu, S, psi11)
                 Mu, Sigma = p1._crossterm_mu_S(Z, mu, S)
                 Mu, Sigma = Mu.reshape(NM,self.input_dim), Sigma.reshape(NM,self.input_dim)
-                
+
                 p2.psi1(Z, Mu, Sigma, psi12)
                 eK2 = psi12.reshape(N, M, M)
                 crossterms = eK2 * (psi11[:, :, None] + psi11[:, None, :])
                 target += crossterms
-                #import ipdb;ipdb.set_trace()
             else:
                 raise NotImplementedError, "psi2 cannot be computed for this kernel"
         return target        
@@ -540,15 +539,15 @@ class kern(Parameterized):
                     # turn around to have rbf in front
                     p1, p2 = self.parts[i2], self.parts[i1]
                     ps1, ps2 = self.param_slices[i2], self.param_slices[i1]  
-                
+
                 N, M = mu.shape[0], Z.shape[0]; NM=N*M
 
                 psi11 = np.zeros((N, M))
                 p1.psi1(Z, mu, S, psi11)
-                
+
                 Mu, Sigma = p1._crossterm_mu_S(Z, mu, S)
                 Mu, Sigma = Mu.reshape(NM,self.input_dim), Sigma.reshape(NM,self.input_dim)
-                
+
                 tmp1 = np.zeros_like(target[ps1])
                 tmp2 = np.zeros_like(target[ps2])
 #                 for n in range(N):
@@ -559,7 +558,7 @@ class kern(Parameterized):
 #                             Mu, Sigma= Mu.reshape(N,M,self.input_dim), Sigma.reshape(N,M,self.input_dim)
 #                             p2.dpsi1_dtheta((dL_dpsi2[n:n+1,m:m+1,m_prime:m_prime+1]*(psi11[n:n+1,m_prime:m_prime+1]))[0], Z[m:m+1], Mu[n:n+1,m], Sigma[n:n+1,m], target[ps2])
 #                             p2.dpsi1_dtheta((dL_dpsi2[n:n+1,m:m+1,m_prime:m_prime+1]*(psi11[n:n+1,m:m+1]))[0], Z[m_prime:m_prime+1], Mu[n:n+1, m_prime], Sigma[n:n+1, m_prime], target[ps2])#Z[m_prime:m_prime+1], Mu[n+m:(n+m)+1], Sigma[n+m:(n+m)+1], target[ps2])
-                
+
                 if isinstance(p1, RBF) and isinstance(p2, RBF):
                     psi12 = np.zeros((N, M))
                     p2.psi1(Z, mu, S, psi12)
@@ -571,11 +570,11 @@ class kern(Parameterized):
                 if isinstance(p1, RBF) and isinstance(p2, Linear):
                     #import ipdb;ipdb.set_trace()
                     pass
-                
+
                 p2.dpsi1_dtheta((dL_dpsi2*(psi11[:,:,None] + psi11[:,None,:])).reshape(NM,M), Z, Mu, Sigma, tmp2)
-                
+
                 target[ps1] += tmp1
-                target[ps2] += tmp2                
+                target[ps2] += tmp2
             else:
                 raise NotImplementedError, "psi2 cannot be computed for this kernel"
 
@@ -615,17 +614,17 @@ class kern(Parameterized):
                 psi11 = np.zeros((N, M))
                 psi12 = np.zeros((NM, M))
                 #psi12_t = np.zeros((N,M))
-                
+
                 p1.psi1(Z, mu, S, psi11)
                 Mu, Sigma = p1._crossterm_mu_S(Z, mu, S)
                 Mu, Sigma = Mu.reshape(NM,self.input_dim), Sigma.reshape(NM,self.input_dim)
-                
+
                 p2.psi1(Z, Mu, Sigma, psi12)
                 tmp1 = np.zeros_like(target)
                 p1.dpsi1_dZ((dL_dpsi2*psi12.reshape(N,M,M)).sum(1), Z, mu, S, tmp1)
                 p1.dpsi1_dZ((dL_dpsi2*psi12.reshape(N,M,M)).sum(2), Z, mu, S, tmp1)
                 target += tmp1
-                
+
                 #p2.dpsi1_dtheta((dL_dpsi2*(psi11[:,:,None] + psi11[:,None,:])).reshape(NM,M), Z, Mu, Sigma, target)
                 p2.dpsi1_dZ((dL_dpsi2*(psi11[:,:,None] + psi11[:,None,:])).reshape(NM,M), Z, Mu, Sigma, target)
             else:
@@ -666,21 +665,21 @@ class kern(Parameterized):
                 psi11 = np.zeros((N, M))
                 psi12 = np.zeros((NM, M))
                 #psi12_t = np.zeros((N,M))
-                
+
                 p1.psi1(Z, mu, S, psi11)
                 Mu, Sigma = p1._crossterm_mu_S(Z, mu, S)
                 Mu, Sigma = Mu.reshape(NM,self.input_dim), Sigma.reshape(NM,self.input_dim)
-                
+
                 p2.psi1(Z, Mu, Sigma, psi12)
                 p1.dpsi1_dmuS((dL_dpsi2*psi12.reshape(N,M,M)).sum(1), Z, mu, S, target_mu, target_S)
                 p1.dpsi1_dmuS((dL_dpsi2*psi12.reshape(N,M,M)).sum(2), Z, mu, S, target_mu, target_S)
-                
+
                 #p2.dpsi1_dtheta((dL_dpsi2*(psi11[:,:,None] + psi11[:,None,:])).reshape(NM,M), Z, Mu, Sigma, target)
                 p2.dpsi1_dmuS((dL_dpsi2*(psi11[:,:,None])).sum(1)*2, Z, Mu.reshape(N,M,self.input_dim).sum(1), Sigma.reshape(N,M,self.input_dim).sum(1), target_mu, target_S)
             else:
                 raise NotImplementedError, "psi2 cannot be computed for this kernel"
         return target_mu, target_S
-    
+
     def plot(self, x=None, plot_limits=None, which_parts='all', resolution=None, *args, **kwargs):
         if which_parts == 'all':
             which_parts = [True] * self.num_parts
@@ -737,15 +736,16 @@ class kern(Parameterized):
         else:
             raise NotImplementedError, "Cannot plot a kernel with more than two input dimensions"
 
-from GPy.core.model import Model
-
+from ..core.model import Model
 class Kern_check_model(Model):
     """This is a dummy model class used as a base class for checking that the gradients of a given kernel are implemented correctly. It enables checkgradient() to be called independently on a kernel."""
     def __init__(self, kernel=None, dL_dK=None, X=None, X2=None):
         num_samples = 20
         num_samples2 = 10
         if kernel==None:
+            import GPy
             kernel = GPy.kern.rbf(1)
+            del GPy
         if X==None:
             X = np.random.normal(size=(num_samples, kernel.input_dim))
         if dL_dK==None:
@@ -753,14 +753,14 @@ class Kern_check_model(Model):
                 dL_dK = np.ones((X.shape[0], X.shape[0]))
             else:
                 dL_dK = np.ones((X.shape[0], X2.shape[0]))
-        
+
         self.kernel=kernel
         self.X = X
         self.X2 = X2
         self.dL_dK = dL_dK
         #self.constrained_indices=[]
         #self.constraints=[]
-        Model.__init__(self)
+        super(Kern_check_model, self).__init__()
 
     def is_positive_definite(self):
         v = np.linalg.eig(self.kernel.K(self.X))[0]
@@ -768,7 +768,7 @@ class Kern_check_model(Model):
             return False
         else:
             return True
-        
+
     def _get_params(self):
         return self.kernel._get_params()
 
@@ -783,7 +783,7 @@ class Kern_check_model(Model):
 
     def _log_likelihood_gradients(self):
         raise NotImplementedError, "This needs to be implemented to use the kern_check_model class."
-    
+
 class Kern_check_dK_dtheta(Kern_check_model):
     """This class allows gradient checks for the gradient of a kernel with respect to parameters. """
     def __init__(self, kernel=None, dL_dK=None, X=None, X2=None):
@@ -798,7 +798,7 @@ class Kern_check_dKdiag_dtheta(Kern_check_model):
         Kern_check_model.__init__(self,kernel=kernel,dL_dK=dL_dK, X=X, X2=None)
         if dL_dK==None:
             self.dL_dK = np.ones((self.X.shape[0]))
-        
+
     def log_likelihood(self):
         return (self.dL_dK*self.kernel.Kdiag(self.X)).sum()
 
@@ -815,7 +815,7 @@ class Kern_check_dK_dX(Kern_check_model):
 
     def _get_param_names(self):
         return ['X_'  +str(i) + ','+str(j) for j in range(self.X.shape[1]) for i in range(self.X.shape[0])]
-                
+
     def _get_params(self):
         return self.X.flatten()
 
@@ -837,7 +837,7 @@ class Kern_check_dKdiag_dX(Kern_check_model):
 
     def _get_param_names(self):
         return ['X_'  +str(i) + ','+str(j) for j in range(self.X.shape[1]) for i in range(self.X.shape[0])]
-                
+
     def _get_params(self):
         return self.X.flatten()
 
@@ -861,13 +861,15 @@ def kern_test(kern, X=None, X2=None, output_ind=None, verbose=False, X_positive=
         if X_positive:
             X = abs(X)
         if output_ind is not None:
-            X[:, output_ind] = np.random.randint(kern.parts[0].output_dim, X.shape[0])
+            assert(output_ind<kern.input_dim)
+            X[:, output_ind] = np.random.randint(low=0,high=kern.parts[0].output_dim, size=X.shape[0])
     if X2==None:
         X2 = np.random.randn(20, kern.input_dim)
         if X_positive:
             X2 = abs(X2)
         if output_ind is not None:
-            X2[:, output_ind] = np.random.randint(kern.parts[0].output_dim, X2.shape[0])
+            assert(output_ind<kern.input_dim)
+            X2[:, output_ind] = np.random.randint(low=0, high=kern.parts[0].output_dim, size=X2.shape[0])
 
     if verbose:
         print("Checking covariance function is positive definite.")
@@ -961,3 +963,4 @@ def kern_test(kern, X=None, X2=None, output_ind=None, verbose=False, X_positive=
         return False
 
     return pass_checks
+del Model
diff --git a/GPy/kern/parts/sympy_helpers.cpp b/GPy/kern/parts/sympy_helpers.cpp
index 9f30eea9..56aa6f21 100644
--- a/GPy/kern/parts/sympy_helpers.cpp
+++ b/GPy/kern/parts/sympy_helpers.cpp
@@ -80,7 +80,7 @@ double ln_diff_erf(double x0, double x1){
   else //x0 and x1 non-positive
     return log(erfcx(-x0)-erfcx(-x1)*exp(x0*x0 - x1*x1))-x0*x0;
 }
-
+// TODO: For all these computations of h things are very efficient at the moment. Need to recode sympykern to allow the precomputations to take place and all the gradients to be computed in one function. Not sure of best way forward for that yet. Neil
 double h(double t, double tprime, double d_i, double d_j, double l){
   // Compute the h function for the sim covariance.
   double half_l_di = 0.5*l*d_i;
@@ -170,9 +170,27 @@ double dh_dl(double t, double tprime, double d_i, double d_j, double l){
 }
 
 double dh_dt(double t, double tprime, double d_i, double d_j, double l){
-  return 0.0;
+  // compute gradient of h function with respect to t.
+  double diff_t = t - tprime;
+  double half_l_di = 0.5*l*d_i;
+  double arg_1 = half_l_di + tprime/l;
+  double arg_2 = half_l_di - diff_t/l;
+  double ln_part_1 = ln_diff_erf(arg_1, arg_2);
+  arg_2 = half_l_di - t/l;
+  double ln_part_2 = ln_diff_erf(half_l_di, arg_2);
+  
+  return (d_i*exp(ln_part_2-d_i*t - d_j*tprime) - d_i*exp(ln_part_1-d_i*diff_t) + 2*exp(-d_i*diff_t - pow(half_l_di - diff_t/l, 2))/(sqrt(M_PI)*l) - 2*exp(-d_i*t - d_j*tprime - pow(half_l_di - t/l,2))/(sqrt(M_PI)*l))*exp(half_l_di*half_l_di)/(d_i + d_j);
 }
 
 double dh_dtprime(double t, double tprime, double d_i, double d_j, double l){
-  return 0.0;
+  // compute gradient of h function with respect to tprime.
+  double diff_t = t - tprime;
+  double half_l_di = 0.5*l*d_i;
+  double arg_1 = half_l_di + tprime/l;
+  double arg_2 = half_l_di - diff_t/l;
+  double ln_part_1 = ln_diff_erf(arg_1, arg_2);
+  arg_2 = half_l_di - t/l;
+  double ln_part_2 = ln_diff_erf(half_l_di, arg_2);
+
+  return (d_i*exp(ln_part_1-d_i*diff_t) + d_j*exp(ln_part_2-d_i*t - d_j*tprime) + (-2*exp(-pow(half_l_di - diff_t/l,2)) + 2*exp(-pow(half_l_di + tprime/l,2)))*exp(-d_i*diff_t)/(sqrt(M_PI)*l))*exp(half_l_di*half_l_di)/(d_i + d_j);
 }
diff --git a/GPy/models.py b/GPy/models.py
index 9a847ea0..a56fb305 100644
--- a/GPy/models.py
+++ b/GPy/models.py
@@ -20,3 +20,4 @@ from _models.mrd import MRD#; _mrd = mrd; del mrd
 from _models.gradient_checker import GradientChecker#; _gradient_checker = gradient_checker ; del gradient_checker 
 from _models.gp_multioutput_regression import GPMultioutputRegression#; _gp_multioutput_regression = gp_multioutput_regression ; del gp_multioutput_regression 
 from _models.sparse_gp_multioutput_regression import SparseGPMultioutputRegression#; _sparse_gp_multioutput_regression = sparse_gp_multioutput_regression ; del sparse_gp_multioutput_regression 
+from _models.gradient_checker import GradientChecker
\ No newline at end of file
diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py
index 5d2fbeec..0fceac60 100644
--- a/GPy/testing/kernel_tests.py
+++ b/GPy/testing/kernel_tests.py
@@ -36,7 +36,7 @@ class KernelTests(unittest.TestCase):
     def test_eq_sympykernel(self):
         if SYMPY_AVAILABLE:
             kern = GPy.kern.eq_sympy(5, 3)
-            self.assertTrue(GPy.kern.kern_test(kern, output_ind=3, verbose=verbose))
+            self.assertTrue(GPy.kern.kern_test(kern, output_ind=4, verbose=verbose))
 
     def test_ode1_eqkernel(self):
         if SYMPY_AVAILABLE:
diff --git a/README.md b/README.md
index 0ff3d890..27af0b0d 100644
--- a/README.md
+++ b/README.md
@@ -9,6 +9,44 @@ A Gaussian processes framework in Python.
 
 Continuous integration status: ![CI status](https://travis-ci.org/SheffieldML/GPy.png)
 
+Getting started
+===============
+Installing with pip
+-------------------
+The simplest way to install GPy is using pip. ubuntu users can do:
+
+    sudo apt-get install python-pip
+    pip install gpy
+
+If you'd like to install from source, or want to contribute to the project (e.g. by sending pull requests via github), read on.
+
+Ubuntu
+------
+For the most part, the developers are using ubuntu. To install the required packages:
+
+    sudo apt-get install python-numpy python-scipy python-matplotlib
+
+clone this git repository and add it to your path:
+
+    git clone git@github.com:SheffieldML/GPy.git ~/SheffieldML
+    echo 'PYTHONPATH=$PYTHONPATH:~/SheffieldML' >> ~/.bashrc
+
+
+Windows
+-------
+On windows, we recommend the ![anaconda python distribution](http://continuum.io/downloads). We've also had luck with ![enthought](http://www.enthought.com). git clone or unzip the source to a suitable directory, and add an approptiate PYTHONPATH environment variable. 
+
+On windows 7 (and possibly earlier versions) there's a bug in scipy version 0.13 which tries to write very long filenames. Reverting to scipy 0.12 seems to do the trick:
+
+    conda install scipy=0.12
+
+OSX
+---
+Everything appears to work out-of-the box using ![enthought](http://www.enthought.com) on osx Mavericks. Download/clone GPy, and then add GPy to your PYTHONPATH
+
+    git clone git@github.com:SheffieldML/GPy.git ~/SheffieldML
+    echo 'PYTHONPATH=$PYTHONPATH:~/SheffieldML' >> ~/.profile
+
 
 Compiling documentation:
 ========================