diff --git a/AUTHORS.txt b/AUTHORS.txt
index fc9e0fb7..f81db5ec 100644
--- a/AUTHORS.txt
+++ b/AUTHORS.txt
@@ -3,4 +3,5 @@ Nicolo Fusi
 Ricardo Andrade
 Nicolas Durrande
 Alan Saul
+Max Zwiessele
 Neil D. Lawrence
diff --git a/GPy/core/model.py b/GPy/core/model.py
index 093a56a3..574bdcb1 100644
--- a/GPy/core/model.py
+++ b/GPy/core/model.py
@@ -447,7 +447,7 @@ class model(parameterised):
         assert isinstance(self.likelihood, likelihoods.EP), "EPEM is only available for EP likelihoods"
         ll_change = epsilon + 1.
         iteration = 0
-        last_ll = -np.exp(1000)
+        last_ll = -np.inf
 
         convergence = False
         alpha = 0
diff --git a/GPy/kern/rbf.py b/GPy/kern/rbf.py
index c413b469..8b644241 100644
--- a/GPy/kern/rbf.py
+++ b/GPy/kern/rbf.py
@@ -56,6 +56,13 @@ class rbf(kernpart):
         self._Z, self._mu, self._S = np.empty(shape=(3,1))
         self._X, self._X2, self._params = np.empty(shape=(3,1))
 
+        #a set of optional args to pass to weave
+        self.weave_options = {'headers'           : ['<omp.h>'],
+                         'extra_compile_args': ['-fopenmp -O3'],  #-march=native'],
+                         'extra_link_args'   : ['-lgomp']}
+
+
+
     def _get_params(self):
         return np.hstack((self.variance,self.lengthscale))
 
@@ -85,8 +92,43 @@ class rbf(kernpart):
         self._K_computations(X,X2)
         target[0] += np.sum(self._K_dvar*dL_dK)
         if self.ARD:
-            if X2 is None: X2 = X
-            [np.add(target[1+q:2+q],(self.variance/self.lengthscale[q]**3)*np.sum(self._K_dvar*dL_dK*np.square(X[:,q][:,None]-X2[:,q][None,:])),target[1+q:2+q]) for q in range(self.D)]
+            dvardLdK = self._K_dvar*dL_dK
+            var_len3 = self.variance/np.power(self.lengthscale,3)
+            if X2 is None:
+                #save computation for the symmetrical case
+                dvardLdK += dvardLdK.T
+                code = """
+                int q,i,j;
+                double tmp;
+                for(q=0; q<D; q++){
+                  tmp = 0;
+                  for(i=0; i<N; i++){
+                    for(j=0; j<i; j++){
+                      tmp += (X(i,q)-X(j,q))*(X(i,q)-X(j,q))*dvardLdK(i,j);
+                    }
+                  }
+                  target(q+1) += var_len3(q)*tmp;
+                }
+                """
+                N,M,D = X.shape[0], X.shape[0], self.D
+            else:
+                code = """
+                int q,i,j;
+                double tmp;
+                for(q=0; q<D; q++){
+                  tmp = 0;
+                  for(i=0; i<N; i++){
+                    for(j=0; j<M; j++){
+                      tmp += (X(i,q)-X2(j,q))*(X(i,q)-X2(j,q))*dvardLdK(i,j);
+                    }
+                  }
+                  target(q+1) += var_len3(q)*tmp;
+                }
+                """
+                N,M,D = X.shape[0], X2.shape[0], self.D
+            #[np.add(target[1+q:2+q],var_len3[q]*np.sum(dvardLdK*np.square(X[:,q][:,None]-X2[:,q][None,:])),target[1+q:2+q]) for q in range(self.D)]
+            weave.inline(code, arg_names=['N','M','D','X','X2','target','dvardLdK','var_len3'],
+                 type_converters=weave.converters.blitz,**self.weave_options)
         else:
             target[1] += (self.variance/self.lengthscale)*np.sum(self._K_dvar*self._K_dist2*dL_dK)
 
@@ -227,10 +269,6 @@ class rbf(kernpart):
             self._Z, self._mu, self._S = Z, mu,S
 
     def weave_psi2(self,mu,Zhat):
-        weave_options = {'headers'           : ['<omp.h>'],
-                         'extra_compile_args': ['-fopenmp -O3'],  #-march=native'],
-                         'extra_link_args'   : ['-lgomp']}
-
         N,Q = mu.shape
         M = Zhat.shape[0]
 
@@ -288,6 +326,6 @@ class rbf(kernpart):
         """
         weave.inline(code, support_code=support_code, libraries=['gomp'],
                      arg_names=['N','M','Q','mu','Zhat','mudist_sq','mudist','lengthscale2','_psi2_denom','psi2_Zdist_sq','psi2_exponent','half_log_psi2_denom','psi2','variance_sq'],
-                     type_converters=weave.converters.blitz,**weave_options)
+                     type_converters=weave.converters.blitz,**self.weave_options)
 
         return mudist,mudist_sq, psi2_exponent, psi2
diff --git a/GPy/models/sparse_GP.py b/GPy/models/sparse_GP.py
index 7ce767ce..3da7e096 100644
--- a/GPy/models/sparse_GP.py
+++ b/GPy/models/sparse_GP.py
@@ -245,3 +245,40 @@ class sparse_GP(GP):
                 var = Kxx - np.sum(np.sum(psi2*Kmmi_LmiBLmi[None,:,:],1),1)
 
         return mu, var[:, None]
+
+    def predict(self, Xnew, X_variance_new=None, which_parts='all', full_cov=False):
+        """
+        Predict the function(s) at the new point(s) Xnew.
+
+        Arguments
+        ---------
+        :param Xnew: The points at which to make a prediction
+        :type Xnew: np.ndarray, Nnew x self.Q
+        :param X_variance_new: The uncertainty in the prediction points
+        :type X_variance_new: np.ndarray, Nnew x self.Q
+        :param which_parts:  specifies which outputs kernel(s) to use in prediction
+        :type which_parts: ('all', list of bools)
+        :param full_cov: whether to return the folll covariance matrix, or just the diagonal
+        :type full_cov: bool
+        :rtype: posterior mean,  a Numpy array, Nnew x self.D
+        :rtype: posterior variance, a Numpy array, Nnew x 1 if full_cov=False, Nnew x Nnew otherwise
+        :rtype: lower and upper boundaries of the 95% confidence intervals, Numpy arrays,  Nnew x self.D
+
+
+           If full_cov and self.D > 1, the return shape of var is Nnew x Nnew x self.D. If self.D == 1, the return shape is Nnew x Nnew.
+           This is to allow for different normalizations of the output dimensions.
+
+        """
+        # normalize X values
+        Xnew = (Xnew.copy() - self._Xmean) / self._Xstd
+        if X_variance_new is not None:
+            X_variance_new = X_variance_new / self._Xstd**2
+
+        #here's the actual prediction by the GP model
+        mu, var = self._raw_predict(Xnew, X_variance_new, full_cov=full_cov, which_parts=which_parts)
+
+        # now push through likelihood
+        mean, var, _025pm, _975pm = self.likelihood.predictive_values(mu, var, full_cov)
+
+        return mean, var, _025pm, _975pm
+
diff --git a/GPy/util/datasets/mocap/cmu/README.txt b/GPy/util/datasets/mocap/cmu/README.txt
new file mode 100644
index 00000000..83bec0b7
--- /dev/null
+++ b/GPy/util/datasets/mocap/cmu/README.txt
@@ -0,0 +1 @@
+this otherwise empty directory is for storing mnocap data, which we don;t distribute
diff --git a/GPy/util/linalg.py b/GPy/util/linalg.py
index 994cadbf..79cd63d7 100644
--- a/GPy/util/linalg.py
+++ b/GPy/util/linalg.py
@@ -313,7 +313,10 @@ def symmetrify(A,upper=False):
     elif A.flags['F_CONTIGUOUS'] and not upper:
         weave.inline(f_contig_code,['A','N'], extra_compile_args=['-O3'])
     else:
-        tmp = np.tril(A)
+        if upper:
+            tmp = np.tril(A.T)
+        else:
+            tmp = np.tril(A)
         A[:] = 0.0
         A += tmp
         A += np.tril(tmp,-1).T