From 7561c4c2326f9c0a1ece372fdb4ae9cd458a7e93 Mon Sep 17 00:00:00 2001
From: James Hensman <james.hensman@gmail.com>
Date: Thu, 2 May 2013 16:13:39 +0100
Subject: [PATCH 1/2] fixed a bug in all_constrained_indices

---
 GPy/core/parameterised.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/GPy/core/parameterised.py b/GPy/core/parameterised.py
index 2e245c5a..bca242f6 100644
--- a/GPy/core/parameterised.py
+++ b/GPy/core/parameterised.py
@@ -191,8 +191,8 @@ class parameterised(object):
         self.constrain(which, transformations.logistic(lower, upper))
 
     def all_constrained_indices(self):
-        if len(self.constrained_indices):
-            return np.hstack(self.constrained_indices)
+        if len(self.constrained_indices) or len(self.fixed_indices):
+            return np.hstack(self.constrained_indices + self.fixed_indices)
         else:
             return np.empty(shape=(0,))
 

From d9252d0e36a4aba68aeae13b1cde79392d9eb5e0 Mon Sep 17 00:00:00 2001
From: James Hensman <james.hensman@gmail.com>
Date: Fri, 3 May 2013 14:00:22 +0100
Subject: [PATCH 2/2] cholesky update for RA

---
 GPy/likelihoods/EP.py   |  5 +++--
 GPy/models/sparse_GP.py |  5 +++++
 GPy/util/linalg.py      | 27 +++++++++++++++++++++++++++
 3 files changed, 35 insertions(+), 2 deletions(-)

diff --git a/GPy/likelihoods/EP.py b/GPy/likelihoods/EP.py
index 118b226a..8307b6b4 100644
--- a/GPy/likelihoods/EP.py
+++ b/GPy/likelihoods/EP.py
@@ -196,8 +196,9 @@ class EP(likelihood):
                 self.tau_tilde[i] = self.tau_tilde[i] + Delta_tau
                 self.v_tilde[i] = self.v_tilde[i] + Delta_v
                 #Posterior distribution parameters update
-                LLT = LLT + np.outer(Kmn[:,i],Kmn[:,i])*Delta_tau
-                L = jitchol(LLT)
+                #LLT = LLT + np.outer(Kmn[:,i],Kmn[:,i])*Delta_tau
+                #L = jitchol(LLT)
+                cholupdate(L,Kmn[:,i]*np.sqrt(Delta_tau))
                 V,info = linalg.lapack.flapack.dtrtrs(L,Kmn,lower=1)
                 Sigma_diag = np.sum(V*V,-2)
                 si = np.sum(V.T*V[:,i],-1)
diff --git a/GPy/models/sparse_GP.py b/GPy/models/sparse_GP.py
index 14c789b8..c3b9f793 100644
--- a/GPy/models/sparse_GP.py
+++ b/GPy/models/sparse_GP.py
@@ -9,6 +9,11 @@ from .. import kern
 from GP import GP
 from scipy import linalg
 
+def backsub_both_sides(L,X):
+    """ Return L^-T * X * L^-1, assumuing X is symmetrical and L is lower cholesky"""
+    tmp,_ = linalg.lapack.flapack.dtrtrs(L,np.asfortranarray(X),lower=1,trans=1)
+    return linalg.lapack.flapack.dtrtrs(L,np.asfortranarray(tmp.T),lower=1,trans=1)[0].T
+
 class sparse_GP(GP):
     """
     Variational sparse GP model
diff --git a/GPy/util/linalg.py b/GPy/util/linalg.py
index 42a98fea..a62fccb3 100644
--- a/GPy/util/linalg.py
+++ b/GPy/util/linalg.py
@@ -276,3 +276,30 @@ def symmetrify_murray(A):
     nn = A.shape[0]
     A[[range(nn),range(nn)]] /= 2.0
 
+def cholupdate(L,x):
+    """
+    update the LOWER cholesky factor of a pd matrix IN PLACE
+
+    if L is the lower chol. of K, then this function computes L_
+    where L_ is the lower chol of K + x*x^T
+    """
+    support_code = """
+    #include <math.h>
+    """
+    code="""
+    double r,c,s;
+    int j,i;
+    for(j=0; j<N; j++){
+      r = sqrt(L(j,j)*L(j,j) + x(j)*x(j));
+      c = r / L(j,j);
+      s = x(j) / L(j,j);
+      L(j,j) = r;
+      for (i=j+1; i<N; i++){
+        L(i,j) = (L(i,j) + s*x(i))/c;
+        x(i) = c*x(i) - s*L(i,j);
+      }
+    }
+    """
+    x = x.copy()
+    N = x.size
+    weave.inline(code, support_code=support_code, arg_names=['N','L','x'], type_converters=weave.converters.blitz)