<s> import itertools <EOL> import numpy as np <EOL> import scipy . linalg as linalg <EOL> import scipy . sparse as sp <EOL> from bayespy . utils import misc as utils <EOL> from . import node as EF <EOL> from . import CovarianceFunctions as CF <EOL> class CovarianceMatrix : <EOL> def cholesky ( self ) : <EOL> pass <EOL> def multiply ( A , B ) : <EOL> return np . multiply ( A , B ) <EOL> def gp_posterior_moment_function ( m , k , x , y , k_sparse = None , pseudoinputs = None , noise = None ) : <EOL> mu = m ( x ) [ <NUM_LIT:0> ] <EOL> K_noise = None <EOL> if noise != None : <EOL> if K_noise is None : <EOL> K_noise = noise <EOL> else : <EOL> K_noise += noise <EOL> if k_sparse != None : <EOL> if K_noise is None : <EOL> K_noise = k_sparse ( x , x ) [ <NUM_LIT:0> ] <EOL> else : <EOL> K_noise += k_sparse ( x , x ) [ <NUM_LIT:0> ] <EOL> if pseudoinputs != None : <EOL> p = pseudoinputs <EOL> K_pp = k ( p , p ) [ <NUM_LIT:0> ] <EOL> K_xp = k ( x , p ) [ <NUM_LIT:0> ] <EOL> U = utils . chol ( K_noise ) <EOL> Lambda = K_pp + np . dot ( K_xp . T , utils . chol_solve ( U , K_xp ) ) <EOL> U_lambda = utils . chol ( Lambda ) <EOL> z = utils . chol_solve ( U_lambda , <EOL> np . dot ( K_xp . T , <EOL> utils . chol_solve ( U , <EOL> y - mu ) ) ) <EOL> U = utils . chol ( K_pp ) <EOL> x = p <EOL> else : <EOL> K = K_noise <EOL> if K is None : <EOL> K = k ( x , x ) [ <NUM_LIT:0> ] <EOL> else : <EOL> try : <EOL> K += k ( x , x ) [ <NUM_LIT:0> ] <EOL> except : <EOL> K = K + k ( x , x ) [ <NUM_LIT:0> ] <EOL> N = len ( y ) <EOL> U = None <EOL> z = None <EOL> if N > <NUM_LIT:0> : <EOL> U = utils . chol ( K ) <EOL> z = utils . chol_solve ( U , y - mu ) <EOL> def get_moments ( h , covariance = <NUM_LIT:1> , mean = True ) : <EOL> K_xh = k ( x , h ) [ <NUM_LIT:0> ] <EOL> if k_sparse != None : <EOL> try : <EOL> K_xh += k_sparse ( x , h ) [ <NUM_LIT:0> ] <EOL> except : <EOL> K_xh = K_xh + k_sparse ( x , h ) [ <NUM_LIT:0> ] <EOL> if not sp . issparse ( K_xh ) : <EOL> K_xh = np . asarray ( K_xh ) <EOL> if mean : <EOL> m_h = m ( h ) [ <NUM_LIT:0> ] <EOL> if z != None : <EOL> m_h += K_xh . T . dot ( z ) <EOL> else : <EOL> m_h = None <EOL> if covariance : <EOL> if covariance == <NUM_LIT:1> : <EOL> k_h = k ( h ) [ <NUM_LIT:0> ] <EOL> if k_sparse != None : <EOL> k_h += k_sparse ( h ) [ <NUM_LIT:0> ] <EOL> if U != None : <EOL> if isinstance ( K_xh , np . ndarray ) : <EOL> k_h -= np . einsum ( '<STR_LIT>' , <EOL> K_xh , <EOL> utils . chol_solve ( U , K_xh ) ) <EOL> else : <EOL> k_h -= np . asarray ( K_xh . multiply ( utils . chol_solve ( U , K_xh ) ) ) . sum ( axis = <NUM_LIT:0> ) <EOL> if pseudoinputs != None : <EOL> if isinstance ( K_xh , np . ndarray ) : <EOL> k_h += np . einsum ( '<STR_LIT>' , <EOL> K_xh , <EOL> utils . chol_solve ( U_lambda , K_xh ) ) <EOL> else : <EOL> k_h += np . asarray ( K_xh . multiply ( utils . chol_solve ( U_lambda , K_xh ) ) ) . sum ( axis = <NUM_LIT:0> ) <EOL> k_h [ k_h < <NUM_LIT:0> ] = <NUM_LIT:0> <EOL> return ( m_h , k_h ) <EOL> elif covariance == <NUM_LIT:2> : <EOL> K_hh = k ( h , h ) [ <NUM_LIT:0> ] <EOL> if k_sparse != None : <EOL> K_hh += k_sparse ( h ) [ <NUM_LIT:0> ] <EOL> if U != None : <EOL> K_hh -= K_xh . T . dot ( utils . chol_solve ( U , K_xh ) ) <EOL> if pseudoinputs != None : <EOL> K_hh += K_xh . T . dot ( utils . chol_solve ( U_lambda , K_xh ) ) <EOL> return ( m_h , K_hh ) <EOL> else : <EOL> return ( m_h , None ) <EOL> return get_moments <EOL> class Constant ( EF . Node ) : <EOL> def __init__ ( self , f , ** kwargs ) : <EOL> self . f = f <EOL> EF . Node . __init__ ( self , dims = [ ( np . inf , ) ] , ** kwargs ) <EOL> def message_to_child ( self , gradient = False ) : <EOL> def func ( x , gradient = False ) : <EOL> if gradient : <EOL> return ( [ self . f ( x ) , None ] , [ ] ) <EOL> else : <EOL> return [ self . f ( x ) , None ] <EOL> return func <EOL> class Multiple ( EF . Node ) : <EOL> def __init__ ( self , GPs , ** kwargs ) : <EOL> EF . NodeVariable . __init__ ( self , <EOL> * GPs , <EOL> plates = ( ) , <EOL> dims = [ ( np . inf , ) , ( np . inf , np . inf ) ] , <EOL> ** kwargs ) <EOL> def message_to_parent ( self , index ) : <EOL> raise Exception ( "<STR_LIT>" ) <EOL> def message_to_child ( self , gradient = False ) : <EOL> u = [ parent . message_to_child ( ) for parent in self . parents ] <EOL> def get_moments ( xh , ** kwargs ) : <EOL> mh_all = [ ] <EOL> khh_all = [ ] <EOL> for i in range ( len ( self . parents ) ) : <EOL> xi = np . array ( xh [ i ] ) <EOL> if gradient : <EOL> ( ( mh , khh ) , dm ) = u [ i ] ( xi , ** kwargs ) <EOL> else : <EOL> ( mh , khh ) = u [ i ] ( xi , ** kwargs ) <EOL> mh_all = np . concatenate ( [ mh_all , mh ] ) <EOL> if khh != None : <EOL> print ( khh ) <EOL> raise Exception ( '<STR_LIT>' ) <EOL> if gradient : <EOL> return ( [ mh_all , khh_all ] , [ ] ) <EOL> else : <EOL> return [ mh_all , khh_all ] <EOL> return get_moments <EOL> class GaussianProcess ( EF . Node ) : <EOL> def __init__ ( self , m , k , k_sparse = None , pseudoinputs = None , ** kwargs ) : <EOL> self . x = np . array ( [ ] ) <EOL> self . f = np . array ( [ ] ) <EOL> if pseudoinputs != None : <EOL> pseudoinputs = EF . NodeConstant ( [ pseudoinputs ] , <EOL> dims = [ np . shape ( pseudoinputs ) ] ) <EOL> self . m = None <EOL> self . k = None <EOL> if isinstance ( k , list ) and isinstance ( m , list ) : <EOL> if len ( k ) != len ( m ) : <EOL> raise Exception ( '<STR_LIT>' ) <EOL> k = CF . Multiple ( k ) <EOL> m = Multiple ( m ) <EOL> elif isinstance ( k , list ) : <EOL> D = len ( k ) <EOL> k = CF . Multiple ( k ) <EOL> m = Multiple ( D * [ m ] ) <EOL> elif isinstance ( m , list ) : <EOL> D = len ( m ) <EOL> k = CF . Multiple ( D * [ k ] ) <EOL> m = Multiple ( m ) <EOL> EF . NodeVariable . __init__ ( self , <EOL> m , <EOL> k , <EOL> k_sparse , <EOL> pseudoinputs , <EOL> plates = ( ) , <EOL> dims = [ ( np . inf , ) , ( np . inf , np . inf ) ] , <EOL> ** kwargs ) <EOL> def __call__ ( self , x , covariance = None ) : <EOL> if not covariance : <EOL> return self . u ( x , covariance = False ) [ <NUM_LIT:0> ] <EOL> elif covariance . lower ( ) == '<STR_LIT>' : <EOL> return self . u ( x , covariance = <NUM_LIT:1> ) <EOL> elif covariance . lower ( ) == '<STR_LIT>' : <EOL> return self . u ( x , covariance = <NUM_LIT:2> ) <EOL> else : <EOL> raise Exception ( "<STR_LIT>" ) <EOL> def message_to_parent ( self , index ) : <EOL> if index == <NUM_LIT:0> : <EOL> k = self . parents [ <NUM_LIT:1> ] . message_to_child ( ) [ <NUM_LIT:0> ] <EOL> K = k ( self . x , self . x ) <EOL> return [ self . x , <EOL> self . mu , <EOL> K ] <EOL> if index == <NUM_LIT:1> : <EOL> raise Exception ( "<STR_LIT>" ) <EOL> def message_to_child ( self ) : <EOL> if self . observed : <EOL> raise Exception ( "<STR_LIT>" ) <EOL> return self . u <EOL> def get_parameters ( self ) : <EOL> return self . u <EOL> def observe ( self , x , f ) : <EOL> self . observed = True <EOL> self . x = x <EOL> self . f = f <EOL> def lower_bound_contribution ( self , gradient = False ) : <EOL> m = self . parents [ <NUM_LIT:0> ] . message_to_child ( gradient = gradient ) <EOL> k = self . parents [ <NUM_LIT:1> ] . message_to_child ( gradient = gradient ) <EOL> if self . parents [ <NUM_LIT:2> ] : <EOL> k_sparse = self . parents [ <NUM_LIT:2> ] . message_to_child ( gradient = gradient ) <EOL> else : <EOL> k_sparse = None <EOL> if self . parents [ <NUM_LIT:3> ] : <EOL> pseudoinputs = self . parents [ <NUM_LIT:3> ] . message_to_child ( gradient = gradient ) <EOL> else : <EOL> pseudoinputs = None <EOL> DKs_xx = [ ] <EOL> DKd_xx = [ ] <EOL> DKd_xp = [ ] <EOL> DKd_pp = [ ] <EOL> Dxp = [ ] <EOL> Dmu = [ ] <EOL> if gradient : <EOL> ( ( mu , _ ) , Dmu ) = m ( self . x , gradient = True ) <EOL> if pseudoinputs : <EOL> ( ( Ks_xx , ) , DKs_xx ) = k_sparse ( self . x , self . x , gradient = True ) <EOL> ( ( xp , ) , Dxp ) = pseudoinputs <EOL> ( ( Kd_pp , ) , DKd_pp ) = k ( xp , xp , gradient = True ) <EOL> ( ( Kd_xp , ) , DKd_xp ) = k ( self . x , xp , gradient = True ) <EOL> else : <EOL> ( ( K_xx , ) , DKd_xx ) = k ( self . x , self . x , gradient = True ) <EOL> if k_sparse : <EOL> ( ( Ks_xx , ) , DKs_xx ) = k_sparse ( self . x , self . x , gradient = True ) <EOL> try : <EOL> K_xx += Ks_xx <EOL> except : <EOL> K_xx = K_xx + Ks_xx <EOL> else : <EOL> ( mu , _ ) = m ( self . x ) <EOL> if pseudoinputs : <EOL> ( Ks_xx , ) = k_sparse ( self . x , self . x ) <EOL> ( xp , ) = pseudoinputs <EOL> ( Kd_pp , ) = k ( xp , xp ) <EOL> ( Kd_xp , ) = k ( self . x , xp ) <EOL> else : <EOL> ( K_xx , ) = k ( self . x , self . x ) <EOL> if k_sparse : <EOL> ( Ks_xx , ) = k_sparse ( self . x , self . x ) <EOL> try : <EOL> K_xx += Ks_xx <EOL> except : <EOL> K_xx = K_xx + Ks_xx <EOL> mu = mu [ <NUM_LIT:0> ] <EOL> if self . observed : <EOL> f0 = self . f - mu <EOL> if pseudoinputs : <EOL> try : <EOL> Us_xx = utils . cholesky ( Ks_xx ) <EOL> except linalg . LinAlgError : <EOL> print ( '<STR_LIT>' ) <EOL> return - np . inf <EOL> z = Us_xx . solve ( f0 ) <EOL> Lambda = Kd_pp + np . dot ( Kd_xp . T , <EOL> Us_xx . solve ( Kd_xp ) ) <EOL> try : <EOL> U_Lambda = utils . cholesky ( Lambda ) <EOL> except linalg . LinAlgError : <EOL> print ( '<STR_LIT>' ) <EOL> return - np . inf <EOL> nu = U_Lambda . solve ( np . dot ( Kd_xp . T , z ) ) <EOL> rho = np . dot ( Kd_xp , nu ) <EOL> y2 = np . dot ( f0 , z ) - np . dot ( z , rho ) <EOL> try : <EOL> Ud_pp = utils . cholesky ( Kd_pp ) <EOL> except linalg . LinAlgError : <EOL> print ( '<STR_LIT>' ) <EOL> return - np . inf <EOL> logdet = ( U_Lambda . logdet ( ) <EOL> + Us_xx . logdet ( ) <EOL> - Ud_pp . logdet ( ) ) <EOL> L = gaussian_logpdf ( y2 , <EOL> <NUM_LIT:0> , <EOL> <NUM_LIT:0> , <EOL> logdet , <EOL> np . size ( self . f ) ) <EOL> for ( dmu , func ) in Dmu : <EOL> d = np . nan <EOL> func ( d ) <EOL> for ( dKs_xx , func ) in DKs_xx : <EOL> d = np . nan <EOL> func ( d ) <EOL> for ( dKd_xp , func ) in DKd_xp : <EOL> d = np . nan <EOL> func ( d ) <EOL> V = Ud_pp . solve ( Kd_xp . T ) <EOL> Z = Us_xx . solve ( V . T ) <EOL> for ( dKd_pp , func ) in DKd_pp : <EOL> d = ( <NUM_LIT:0.5> * np . trace ( Ud_pp . solve ( dKd_pp ) ) <EOL> - <NUM_LIT:0.5> * np . trace ( U_Lambda . solve ( dKd_pp ) ) <EOL> + np . dot ( nu , np . dot ( dKd_pp , nu ) ) <EOL> + np . trace ( np . dot ( dKd_pp , <EOL> np . dot ( V , Z ) ) ) ) <EOL> func ( d ) <EOL> for ( dxp , func ) in Dxp : <EOL> d = np . nan <EOL> func ( d ) <EOL> else : <EOL> try : <EOL> U = utils . cholesky ( K_xx ) <EOL> except linalg . LinAlgError : <EOL> print ( '<STR_LIT>' ) <EOL> return - np . inf <EOL> z = U . solve ( f0 ) <EOL> L = utils . gaussian_logpdf ( np . dot ( f0 , z ) , <EOL> <NUM_LIT:0> , <EOL> <NUM_LIT:0> , <EOL> U . logdet ( ) , <EOL> np . size ( self . f ) ) <EOL> for ( dmu , func ) in Dmu : <EOL> d = - np . sum ( z ) <EOL> func ( d ) <EOL> for ( dK , func ) in DKd_xx : <EOL> d = <NUM_LIT:0.5> * ( dK . dot ( z ) . dot ( z ) <EOL> - U . trace_solve_gradient ( dK ) ) <EOL> func ( d ) <EOL> for ( dK , func ) in DKs_xx : <EOL> d = <NUM_LIT:0.5> * ( dK . dot ( z ) . dot ( z ) <EOL> - U . trace_solve_gradient ( dK ) ) <EOL> func ( d ) <EOL> else : <EOL> raise Exception ( '<STR_LIT>' ) <EOL> return L <EOL> def update ( self ) : <EOL> m = self . parents [ <NUM_LIT:0> ] . message_to_child ( ) <EOL> k = self . parents [ <NUM_LIT:1> ] . message_to_child ( ) <EOL> if self . parents [ <NUM_LIT:2> ] : <EOL> k_sparse = self . parents [ <NUM_LIT:2> ] . message_to_child ( ) <EOL> else : <EOL> k_sparse = None <EOL> if self . parents [ <NUM_LIT:3> ] : <EOL> pseudoinputs = self . parents [ <NUM_LIT:3> ] . message_to_child ( ) [ <NUM_LIT:0> ] <EOL> else : <EOL> pseudoinputs = None <EOL> if self . observed : <EOL> self . u = gp_posterior_moment_function ( m , <EOL> k , <EOL> self . x , <EOL> self . f , <EOL> k_sparse = k_sparse , <EOL> pseudoinputs = pseudoinputs ) <EOL> else : <EOL> x = np . array ( [ ] ) <EOL> y = np . array ( [ ] ) <EOL> for ( child , index ) in self . children : <EOL> ( msg , mask ) = child . message_to_parent ( index ) <EOL> x = np . concatenate ( ( x , msg [ <NUM_LIT:0> ] ) , axis = - <NUM_LIT:2> ) <EOL> y = np . concatenate ( ( y , msg [ <NUM_LIT:1> ] ) ) <EOL> V = linalg . block_diag ( V , msg [ <NUM_LIT:2> ] ) <EOL> self . u = gp_posterior_moment_function ( m , k , x , y , covariance = V ) <EOL> self . x = x <EOL> self . f = y </s>
