Merge branch 'devel' into mrd

This commit is contained in:
Max Zwiessele 2013-05-03 14:39:48 +01:00
commit 02b3f41680
4 changed files with 37 additions and 4 deletions

View file

@ -191,8 +191,8 @@ class parameterised(object):
self.constrain(which, transformations.logistic(lower, upper)) self.constrain(which, transformations.logistic(lower, upper))
def all_constrained_indices(self): def all_constrained_indices(self):
if len(self.constrained_indices): if len(self.constrained_indices) or len(self.fixed_indices):
return np.hstack(self.constrained_indices) return np.hstack(self.constrained_indices + self.fixed_indices)
else: else:
return np.empty(shape=(0,)) return np.empty(shape=(0,))

View file

@ -196,8 +196,9 @@ class EP(likelihood):
self.tau_tilde[i] = self.tau_tilde[i] + Delta_tau self.tau_tilde[i] = self.tau_tilde[i] + Delta_tau
self.v_tilde[i] = self.v_tilde[i] + Delta_v self.v_tilde[i] = self.v_tilde[i] + Delta_v
#Posterior distribution parameters update #Posterior distribution parameters update
LLT = LLT + np.outer(Kmn[:,i],Kmn[:,i])*Delta_tau #LLT = LLT + np.outer(Kmn[:,i],Kmn[:,i])*Delta_tau
L = jitchol(LLT) #L = jitchol(LLT)
cholupdate(L,Kmn[:,i]*np.sqrt(Delta_tau))
V,info = linalg.lapack.flapack.dtrtrs(L,Kmn,lower=1) V,info = linalg.lapack.flapack.dtrtrs(L,Kmn,lower=1)
Sigma_diag = np.sum(V*V,-2) Sigma_diag = np.sum(V*V,-2)
si = np.sum(V.T*V[:,i],-1) si = np.sum(V.T*V[:,i],-1)

View file

@ -9,6 +9,11 @@ from .. import kern
from GP import GP from GP import GP
from scipy import linalg 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): class sparse_GP(GP):
""" """
Variational sparse GP model Variational sparse GP model

View file

@ -276,3 +276,30 @@ def symmetrify_murray(A):
nn = A.shape[0] nn = A.shape[0]
A[[range(nn),range(nn)]] /= 2.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)