diff --git a/GPy/core/fitc.py b/GPy/core/fitc.py index 515abc29..ef171459 100644 --- a/GPy/core/fitc.py +++ b/GPy/core/fitc.py @@ -3,10 +3,10 @@ import numpy as np import pylab as pb -from ..util.linalg import mdot, jitchol, chol_inv, tdot, symmetrify,pdinv +from ..util.linalg import mdot, jitchol, chol_inv, tdot, symmetrify, pdinv, dtrtrs from ..util.plot import gpplot from .. import kern -from scipy import stats, linalg +from scipy import stats from sparse_gp import SparseGP class FITC(SparseGP): @@ -50,7 +50,7 @@ class FITC(SparseGP): def _computations(self): #factor Kmm self.Lm = jitchol(self.Kmm) - self.Lmi,info = linalg.lapack.flapack.dtrtrs(self.Lm,np.eye(self.num_inducing),lower=1) + self.Lmi,info = dtrtrs(self.Lm,np.eye(self.num_inducing),lower=1) Lmipsi1 = np.dot(self.Lmi,self.psi1) self.Qnn = np.dot(Lmipsi1.T,Lmipsi1).copy() self.Diag0 = self.psi0 - np.diag(self.Qnn) @@ -59,7 +59,7 @@ class FITC(SparseGP): # The rather complex computations of self.A tmp = self.psi1 * (np.sqrt(self.beta_star.flatten().reshape(1, self.num_data))) - tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1) + tmp, _ = dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1) self.A = tdot(tmp) # factor B @@ -68,8 +68,8 @@ class FITC(SparseGP): self.LBi = chol_inv(self.LB) self.psi1V = np.dot(self.psi1, self.V_star) - Lmi_psi1V, info = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(self.psi1V), lower=1, trans=0) - self._LBi_Lmi_psi1V, _ = linalg.lapack.flapack.dtrtrs(self.LB, np.asfortranarray(Lmi_psi1V), lower=1, trans=0) + Lmi_psi1V, info = dtrtrs(self.Lm, np.asfortranarray(self.psi1V), lower=1, trans=0) + self._LBi_Lmi_psi1V, _ = dtrtrs(self.LB, np.asfortranarray(Lmi_psi1V), lower=1, trans=0) Kmmipsi1 = np.dot(self.Lmi.T,Lmipsi1) b_psi1_Ki = self.beta_star * Kmmipsi1.T @@ -188,7 +188,7 @@ class FITC(SparseGP): self.P = Iplus_Dprod_i[:,None] * self.psi1.T self.RPT0 = np.dot(self.Lmi,self.psi1) self.L = np.linalg.cholesky(np.eye(self.num_inducing) + np.dot(self.RPT0,((1. - Iplus_Dprod_i)/self.Diag0)[:,None]*self.RPT0.T)) - self.R,info = linalg.lapack.flapack.dtrtrs(self.L,self.Lmi,lower=1) + self.R,info = dtrtrs(self.L,self.Lmi,lower=1) self.RPT = np.dot(self.R,self.P.T) self.Sigma = np.diag(self.Diag) + np.dot(self.RPT.T,self.RPT) self.w = self.Diag * self.likelihood.v_tilde @@ -210,7 +210,7 @@ class FITC(SparseGP): # = I - [RPT0] * (U*U.T)^-1 * [RPT0].T # = I - V.T * V U = np.linalg.cholesky(np.diag(self.Diag0) + self.Qnn) - V,info = linalg.lapack.flapack.dtrtrs(U,self.RPT0.T,lower=1) + V,info = dtrtrs(U,self.RPT0.T,lower=1) C = np.eye(self.num_inducing) - np.dot(V.T,V) mu_u = np.dot(C,self.RPT0)*(1./self.Diag0[None,:]) #self.C = C diff --git a/GPy/core/gp.py b/GPy/core/gp.py index 7b67e722..5172d9e7 100644 --- a/GPy/core/gp.py +++ b/GPy/core/gp.py @@ -3,10 +3,9 @@ import numpy as np -from scipy import linalg import pylab as pb from .. import kern -from ..util.linalg import pdinv, mdot, tdot +from ..util.linalg import pdinv, mdot, tdot, dpotrs, dtrtrs #from ..util.plot import gpplot, Tango from ..likelihoods import EP from gp_base import GPBase @@ -44,13 +43,13 @@ class GP(GPBase): # the gradient of the likelihood wrt the covariance matrix if self.likelihood.YYT is None: #alpha = np.dot(self.Ki, self.likelihood.Y) - alpha,_ = linalg.lapack.flapack.dpotrs(self.L, self.likelihood.Y,lower=1) + alpha,_ = dpotrs(self.L, self.likelihood.Y,lower=1) self.dL_dK = 0.5 * (tdot(alpha) - self.output_dim * self.Ki) else: #tmp = mdot(self.Ki, self.likelihood.YYT, self.Ki) - tmp, _ = linalg.lapack.flapack.dpotrs(self.L, np.asfortranarray(self.likelihood.YYT), lower=1) - tmp, _ = linalg.lapack.flapack.dpotrs(self.L, np.asfortranarray(tmp.T), lower=1) + tmp, _ = dpotrs(self.L, np.asfortranarray(self.likelihood.YYT), lower=1) + tmp, _ = dpotrs(self.L, np.asfortranarray(tmp.T), lower=1) self.dL_dK = 0.5 * (tmp - self.output_dim * self.Ki) def _get_params(self): @@ -76,7 +75,7 @@ class GP(GPBase): Computes the model fit using YYT if it's available """ if self.likelihood.YYT is None: - tmp, _ = linalg.lapack.flapack.dtrtrs(self.L, np.asfortranarray(self.likelihood.Y), lower=1) + tmp, _ = dtrtrs(self.L, np.asfortranarray(self.likelihood.Y), lower=1) return -0.5 * np.sum(np.square(tmp)) #return -0.5 * np.sum(np.square(np.dot(self.Li, self.likelihood.Y))) else: @@ -108,7 +107,7 @@ class GP(GPBase): """ Kx = self.kern.K(_Xnew,self.X,which_parts=which_parts).T #KiKx = np.dot(self.Ki, Kx) - KiKx, _ = linalg.lapack.flapack.dpotrs(self.L, np.asfortranarray(Kx), lower=1) + KiKx, _ = dpotrs(self.L, np.asfortranarray(Kx), lower=1) mu = np.dot(KiKx.T, self.likelihood.Y) if full_cov: Kxx = self.kern.K(_Xnew, which_parts=which_parts) diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py index df3160f1..22d0b8a2 100644 --- a/GPy/core/sparse_gp.py +++ b/GPy/core/sparse_gp.py @@ -3,7 +3,7 @@ import numpy as np import pylab as pb -from ..util.linalg import mdot, jitchol, tdot, symmetrify, backsub_both_sides, chol_inv +from ..util.linalg import mdot, jitchol, tdot, symmetrify, backsub_both_sides, chol_inv, dtrtrs, dpotrs, dpotri from scipy import linalg from ..likelihoods import Gaussian from gp_base import GPBase @@ -80,7 +80,7 @@ class SparseGP(GPBase): tmp = self.psi1 * (np.sqrt(self.likelihood.precision.flatten().reshape(1, self.num_data))) else: tmp = self.psi1 * (np.sqrt(self.likelihood.precision)) - tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1) + tmp, _ = dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1) self.A = tdot(tmp) @@ -92,10 +92,10 @@ class SparseGP(GPBase): self.psi1V = np.dot(self.psi1, self.likelihood.V) # back substutue C into psi1V - tmp, info1 = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(self.psi1V), lower=1, trans=0) - self._LBi_Lmi_psi1V, _ = linalg.lapack.flapack.dtrtrs(self.LB, np.asfortranarray(tmp), lower=1, trans=0) - tmp, info2 = linalg.lapack.flapack.dpotrs(self.LB, tmp, lower=1) - self.Cpsi1V, info3 = linalg.lapack.flapack.dtrtrs(self.Lm, tmp, lower=1, trans=1) + tmp, info1 = dtrtrs(self.Lm, np.asfortranarray(self.psi1V), lower=1, trans=0) + self._LBi_Lmi_psi1V, _ = dtrtrs(self.LB, np.asfortranarray(tmp), lower=1, trans=0) + tmp, info2 = dpotrs(self.LB, tmp, lower=1) + self.Cpsi1V, info3 = dtrtrs(self.Lm, tmp, lower=1, trans=1) # Compute dL_dKmm tmp = tdot(self._LBi_Lmi_psi1V) @@ -219,7 +219,7 @@ class SparseGP(GPBase): def _raw_predict(self, Xnew, X_variance_new=None, which_parts='all', full_cov=False): """Internal helper function for making predictions, does not account for normalization""" - Bi, _ = linalg.lapack.flapack.dpotri(self.LB, lower=0) # WTH? this lower switch should be 1, but that doesn't work! + Bi, _ = dpotri(self.LB, lower=0) # WTH? this lower switch should be 1, but that doesn't work! symmetrify(Bi) Kmmi_LmiBLmi = backsub_both_sides(self.Lm, np.eye(self.num_inducing) - Bi) diff --git a/GPy/likelihoods/ep.py b/GPy/likelihoods/ep.py index e05e69ea..fb9a55c7 100644 --- a/GPy/likelihoods/ep.py +++ b/GPy/likelihoods/ep.py @@ -1,6 +1,6 @@ import numpy as np -from scipy import stats, linalg -from ..util.linalg import pdinv,mdot,jitchol,chol_inv,DSYR,tdot +from scipy import stats +from ..util.linalg import pdinv,mdot,jitchol,chol_inv,DSYR,tdot,dtrtrs from likelihood import likelihood class EP(likelihood): @@ -124,7 +124,7 @@ class EP(likelihood): Sroot_tilde_K = np.sqrt(self.tau_tilde)[:,None]*K B = np.eye(self.N) + np.sqrt(self.tau_tilde)[None,:]*Sroot_tilde_K L = jitchol(B) - V,info = linalg.lapack.flapack.dtrtrs(L,Sroot_tilde_K,lower=1) + V,info = dtrtrs(L,Sroot_tilde_K,lower=1) Sigma = K - np.dot(V.T,V) mu = np.dot(Sigma,self.v_tilde) epsilon_np1 = sum((self.tau_tilde-self.np1[-1])**2)/self.N @@ -210,7 +210,7 @@ class EP(likelihood): DSYR(LLT,Kmn[:,i].copy(),Delta_tau) #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) + V,info = dtrtrs(L,Kmn,lower=1) Sigma_diag = np.sum(V*V,-2) si = np.sum(V.T*V[:,i],-1) mu += (Delta_v-Delta_tau*mu[i])*si @@ -218,8 +218,8 @@ class EP(likelihood): #Sigma recomputation with Cholesky decompositon LLT = LLT0 + np.dot(Kmn*self.tau_tilde[None,:],Kmn.T) L = jitchol(LLT) - V,info = linalg.lapack.flapack.dtrtrs(L,Kmn,lower=1) - V2,info = linalg.lapack.flapack.dtrtrs(L.T,V,lower=0) + V,info = dtrtrs(L,Kmn,lower=1) + V2,info = dtrtrs(L.T,V,lower=0) Sigma_diag = np.sum(V*V,-2) Knmv_tilde = np.dot(Kmn,self.v_tilde) mu = np.dot(V2.T,Knmv_tilde) @@ -322,7 +322,7 @@ class EP(likelihood): P = Iplus_Dprod_i[:,None] * P0 safe_diag = np.where(Diag0 < self.tau_tilde, self.tau_tilde/(1.+Diag0*self.tau_tilde), (1. - Iplus_Dprod_i)/Diag0) L = jitchol(np.eye(num_inducing) + np.dot(RPT0,safe_diag[:,None]*RPT0.T)) - R,info = linalg.lapack.flapack.dtrtrs(L,R0,lower=1) + R,info = dtrtrs(L,R0,lower=1) RPT = np.dot(R,P.T) Sigma_diag = Diag + np.sum(RPT.T*RPT.T,-1) self.w = Diag * self.v_tilde diff --git a/GPy/testing/examples_tests.py b/GPy/testing/examples_tests.py index a682c547..ec030055 100644 --- a/GPy/testing/examples_tests.py +++ b/GPy/testing/examples_tests.py @@ -39,7 +39,7 @@ def model_instance(model): #assert isinstance(model, GPy.core.model) return isinstance(model, GPy.core.Model) -#@nottest +@nottest def test_models(): examples_path = os.path.dirname(GPy.examples.__file__) # Load modules diff --git a/GPy/util/linalg.py b/GPy/util/linalg.py index c04fc460..d2464128 100644 --- a/GPy/util/linalg.py +++ b/GPy/util/linalg.py @@ -11,6 +11,12 @@ import types import ctypes from ctypes import byref, c_char, c_int, c_double # TODO # import scipy.lib.lapack +import scipy + +if np.all(np.float64((scipy.__version__).split('.')[:2]) >= np.array([0, 10])): + import scipy.linalg.lapack as lapack +else: + import scipy.linalg.lapack.flapack as lapack try: _blaslib = ctypes.cdll.LoadLibrary(np.core._dotblas.__file__) # @UndefinedVariable @@ -18,6 +24,36 @@ try: except: _blas_available = False +def dtrtrs(A, B, lower=0, trans=0, unitdiag=0, overwrite_b=0): + """Wrapper for lapack dtrtrs function + + :param A: Matrix A + :param B: Matrix B + :param lower: is matrix lower (true) or upper (false) + :returns: + """ + return lapack.dtrtrs(A, B, lower=lower, trans=trans, unitdiag=unitdiag, overwrite_b=overwrite_b) + +def dpotrs(A, B, overwrite_b=0, lower=0): + """Wrapper for lapack dpotrs function + + :param A: Matrix A + :param B: Matrix B + :param lower: is matrix lower (true) or upper (false) + :returns: + """ + return lapack.dpotrs(A, B, overwrite_b=overwrite_b, lower=lower) + +def dpotri(A, B, overwrite_b=0, lower=0): + """Wrapper for lapack dpotri function + + :param A: Matrix A + :param B: Matrix B + :param lower: is matrix lower (true) or upper (false) + :returns: + """ + return lapack.dpotri(A, B, overwrite_b=overwrite_b, lower=lower) + def trace_dot(a, b): """ efficiently compute the trace of the matrix product of a and b @@ -56,7 +92,7 @@ def _mdot_r(a, b): def jitchol(A, maxtries=5): A = np.asfortranarray(A) - L, info = linalg.lapack.flapack.dpotrf(A, lower=1) + L, info = lapack.dpotrf(A, lower=1) if info == 0: return L else: @@ -117,7 +153,7 @@ def pdinv(A, *args): L = jitchol(A, *args) logdet = 2.*np.sum(np.log(np.diag(L))) Li = chol_inv(L) - Ai, _ = linalg.lapack.flapack.dpotri(L) + Ai, _ = lapack.dpotri(L) # Ai = np.tril(Ai) + np.tril(Ai,-1).T symmetrify(Ai) @@ -133,7 +169,7 @@ def chol_inv(L): """ - return linalg.lapack.flapack.dtrtri(L, lower=True)[0] + return lapack.dtrtri(L, lower=True)[0] def multiple_pdinv(A): @@ -150,7 +186,7 @@ def multiple_pdinv(A): N = A.shape[-1] chols = [jitchol(A[:, :, i]) for i in range(N)] halflogdets = [np.sum(np.log(np.diag(L[0]))) for L in chols] - invs = [linalg.lapack.flapack.dpotri(L[0], True)[0] for L in chols] + invs = [lapack.dpotri(L[0], True)[0] for L in chols] invs = [np.triu(I) + np.triu(I, 1).T for I in invs] return np.dstack(invs), np.array(halflogdets) @@ -351,9 +387,9 @@ def cholupdate(L, x): def backsub_both_sides(L, X, transpose='left'): """ Return L^-T * X * L^-1, assumuing X is symmetrical and L is lower cholesky""" if transpose == 'left': - 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 + tmp, _ = lapack.dtrtrs(L, np.asfortranarray(X), lower=1, trans=1) + return lapack.dtrtrs(L, np.asfortranarray(tmp.T), lower=1, trans=1)[0].T else: - tmp, _ = linalg.lapack.flapack.dtrtrs(L, np.asfortranarray(X), lower=1, trans=0) - return linalg.lapack.flapack.dtrtrs(L, np.asfortranarray(tmp.T), lower=1, trans=0)[0].T + tmp, _ = lapack.dtrtrs(L, np.asfortranarray(X), lower=1, trans=0) + return lapack.dtrtrs(L, np.asfortranarray(tmp.T), lower=1, trans=0)[0].T