<s> import numpy as np <EOL> from scipy . special import gammaln , psi <EOL> from scipy . misc import logsumexp <EOL> from joblib import Parallel , delayed <EOL> from pybasicbayes . abstractions import GibbsSampling , MeanField , MeanFieldSVI <EOL> from pyhawkes . internals . distributions import Bernoulli , Gamma <EOL> from pyhawkes . utils . utils import logistic , logit <EOL> class SpikeAndSlabGammaWeights ( GibbsSampling ) : <EOL> """<STR_LIT>""" <EOL> def __init__ ( self , model , parallel_resampling = True ) : <EOL> """<STR_LIT>""" <EOL> self . model = model <EOL> self . K = model . K <EOL> self . network = model . network <EOL> self . parallel_resampling = parallel_resampling <EOL> self . A = np . ones ( ( self . K , self . K ) ) <EOL> self . W = np . zeros ( ( self . K , self . K ) ) <EOL> self . resample ( ) <EOL> @ property <EOL> def W_effective ( self ) : <EOL> return self . A * self . W <EOL> def log_likelihood ( self , x ) : <EOL> """<STR_LIT>""" <EOL> A , W = x <EOL> assert isinstance ( A , np . ndarray ) and A . shape == ( self . K , self . K ) , "<STR_LIT>" <EOL> assert isinstance ( W , np . ndarray ) and W . shape == ( self . K , self . K ) , "<STR_LIT>" <EOL> rho = np . clip ( self . network . P , <NUM_LIT> , <NUM_LIT:1> - <NUM_LIT> ) <EOL> ll = ( A * np . log ( rho ) + ( <NUM_LIT:1> - A ) * np . log ( <NUM_LIT:1> - rho ) ) . sum ( ) <EOL> ll = np . nan_to_num ( ll ) <EOL> kappa = self . network . kappa <EOL> v = self . network . V <EOL> lp_W = kappa * np . log ( v ) - gammaln ( kappa ) + ( kappa - <NUM_LIT:1> ) * np . log ( W ) - v * W <EOL> ll += ( A * lp_W ) . sum ( ) <EOL> return ll <EOL> def log_probability ( self ) : <EOL> return self . log_likelihood ( ( self . A , self . W ) ) <EOL> def rvs ( self , size = [ ] ) : <EOL> A = np . random . rand ( self . K , self . K ) < self . network . P <EOL> W = np . random . gamma ( self . network . kappa , <NUM_LIT:1.0> / self . network . V , <EOL> size ( self . K , self . K ) ) <EOL> return A , W <EOL> def _joint_resample_A_W ( self ) : <EOL> """<STR_LIT>""" <EOL> raise NotImplementedError ( ) <EOL> def _joblib_resample_A_given_W ( self , data ) : <EOL> """<STR_LIT>""" <EOL> import pyhawkes . internals . parallel_adjacency_resampling as par <EOL> par . model = self . model <EOL> par . data = data <EOL> par . K = self . model . K <EOL> if len ( data ) == <NUM_LIT:0> : <EOL> self . A = np . random . rand ( self . K , self . K ) < self . network . P <EOL> return <EOL> A_cols = Parallel ( n_jobs = - <NUM_LIT:1> , backend = "<STR_LIT>" ) ( <EOL> delayed ( par . _resample_column_of_A ) ( k2 ) for k2 in range ( self . K ) ) <EOL> self . A = np . array ( A_cols ) . T <EOL> def _resample_A_given_W ( self , data ) : <EOL> """<STR_LIT>""" <EOL> p = self . network . P <EOL> for k1 in xrange ( self . K ) : <EOL> for k2 in xrange ( self . K ) : <EOL> if self . model is None : <EOL> ll0 = <NUM_LIT:0> <EOL> ll1 = <NUM_LIT:0> <EOL> else : <EOL> self . A [ k1 , k2 ] = <NUM_LIT:0> <EOL> ll0 = sum ( [ d . log_likelihood_single_process ( k2 ) for d in data ] ) <EOL> self . A [ k1 , k2 ] = <NUM_LIT:1> <EOL> ll1 = sum ( [ d . log_likelihood_single_process ( k2 ) for d in data ] ) <EOL> lp0 = ll0 + np . log ( <NUM_LIT:1.0> - p [ k1 , k2 ] ) <EOL> lp1 = ll1 + np . log ( p [ k1 , k2 ] ) <EOL> Z = logsumexp ( [ lp0 , lp1 ] ) <EOL> self . A [ k1 , k2 ] = np . log ( np . random . rand ( ) ) < lp1 - Z <EOL> def resample_W_given_A_and_z ( self , data = [ ] ) : <EOL> """<STR_LIT>""" <EOL> ss = np . zeros ( ( <NUM_LIT:2> , self . K , self . K ) ) + sum ( [ d . compute_weight_ss ( ) for d in data ] ) <EOL> ss [ <NUM_LIT:1> ] *= self . A <EOL> kappa_post = self . network . kappa + ss [ <NUM_LIT:0> ] <EOL> v_post = self . network . V + ss [ <NUM_LIT:1> ] <EOL> self . W = np . atleast_1d ( np . random . gamma ( kappa_post , <NUM_LIT:1.0> / v_post ) ) . reshape ( ( self . K , self . K ) ) <EOL> def resample ( self , data = [ ] ) : <EOL> """<STR_LIT>""" <EOL> self . resample_W_given_A_and_z ( data ) <EOL> if self . parallel_resampling : <EOL> self . _joblib_resample_A_given_W ( data ) <EOL> else : <EOL> self . _resample_A_given_W ( data ) <EOL> class GammaMixtureWeights ( GibbsSampling , MeanField , MeanFieldSVI ) : <EOL> """<STR_LIT>""" <EOL> def __init__ ( self , model , kappa_0 = <NUM_LIT:0.1> , nu_0 = <NUM_LIT> ) : <EOL> """<STR_LIT>""" <EOL> self . model = model <EOL> self . network = model . network <EOL> self . K = model . K <EOL> self . network = model . network <EOL> assert model . network is not None , "<STR_LIT>" <EOL> self . kappa_0 = kappa_0 <EOL> self . nu_0 = nu_0 <EOL> self . mf_p = np . ones ( ( self . K , self . K ) ) - <NUM_LIT> <EOL> self . mf_kappa_0 = self . kappa_0 * np . ones ( ( self . K , self . K ) ) <EOL> self . mf_v_0 = self . nu_0 * np . ones ( ( self . K , self . K ) ) <EOL> self . mf_kappa_1 = self . network . kappa * np . ones ( ( self . K , self . K ) ) <EOL> self . mf_v_1 = self . network . V * np . ones ( ( self . K , self . K ) ) <EOL> self . A = np . ones ( ( self . K , self . K ) ) <EOL> self . W = np . zeros ( ( self . K , self . K ) ) <EOL> self . resample ( ) <EOL> @ property <EOL> def W_effective ( self ) : <EOL> return self . W <EOL> def log_likelihood ( self , x ) : <EOL> """<STR_LIT>""" <EOL> A , W = x <EOL> assert isinstance ( A , np . ndarray ) and A . shape == ( self . K , self . K ) , "<STR_LIT>" <EOL> assert isinstance ( W , np . ndarray ) and W . shape == ( self . K , self . K ) , "<STR_LIT>" <EOL> rho = self . network . P <EOL> lp_A = ( A * np . log ( rho ) + ( <NUM_LIT:1> - A ) * np . log ( <NUM_LIT:1> - rho ) ) <EOL> kappa = self . network . kappa <EOL> v = self . network . V <EOL> lp_W0 = ( self . kappa_0 * np . log ( self . nu_0 ) - gammaln ( self . kappa_0 ) <EOL> + ( self . kappa_0 - <NUM_LIT:1> ) * np . log ( W ) - self . nu_0 * W ) [ A == <NUM_LIT:0> ] <EOL> lp_W1 = ( kappa * np . log ( v ) - gammaln ( kappa ) <EOL> + ( kappa - <NUM_LIT:1> ) * np . log ( W ) - v * W ) [ A == <NUM_LIT:1> ] <EOL> ll = lp_A . sum ( ) + lp_W0 . sum ( ) + lp_W1 . sum ( ) <EOL> return ll <EOL> def log_probability ( self ) : <EOL> return self . log_likelihood ( ( self . A , self . W ) ) <EOL> def rvs ( self , size = [ ] ) : <EOL> raise NotImplementedError ( ) <EOL> def expected_A ( self ) : <EOL> return self . mf_p <EOL> def expected_W ( self ) : <EOL> """<STR_LIT>""" <EOL> p_A = self . expected_A ( ) <EOL> E_W = p_A * self . expected_W_given_A ( <NUM_LIT:1.0> ) + ( <NUM_LIT:1> - p_A ) * self . expected_W_given_A ( <NUM_LIT:0.0> ) <EOL> if not self . network . allow_self_connections : <EOL> np . fill_diagonal ( E_W , <NUM_LIT:0.0> ) <EOL> return E_W <EOL> def expected_W_given_A ( self , A ) : <EOL> """<STR_LIT>""" <EOL> return A * ( self . mf_kappa_1 / self . mf_v_1 ) + ( <NUM_LIT:1.0> - A ) * ( self . mf_kappa_0 / self . mf_v_0 ) <EOL> def std_A ( self ) : <EOL> """<STR_LIT>""" <EOL> return np . sqrt ( self . mf_p * ( <NUM_LIT:1> - self . mf_p ) ) <EOL> def expected_log_W ( self ) : <EOL> """<STR_LIT>""" <EOL> p_A = self . expected_A ( ) <EOL> E_ln_W = p_A * self . expected_log_W_given_A ( <NUM_LIT:1.0> ) + ( <NUM_LIT:1> - p_A ) * self . expected_log_W_given_A ( <NUM_LIT:0.0> ) <EOL> if not self . network . allow_self_connections : <EOL> np . fill_diagonal ( E_ln_W , - np . inf ) <EOL> return E_ln_W <EOL> def expected_log_W_given_A ( self , A ) : <EOL> """<STR_LIT>""" <EOL> return A * ( psi ( self . mf_kappa_1 ) - np . log ( self . mf_v_1 ) ) + ( <NUM_LIT:1.0> - A ) * ( psi ( self . mf_kappa_0 ) - np . log ( self . mf_v_0 ) ) <EOL> def expected_log_likelihood ( self , x ) : <EOL> raise NotImplementedError ( ) <EOL> def meanfieldupdate_p ( self , stepsize = <NUM_LIT:1.0> ) : <EOL> """<STR_LIT>""" <EOL> logit_p = self . network . expected_log_p ( ) - self . network . expected_log_notp ( ) <EOL> logit_p += self . network . kappa * self . network . expected_log_v ( ) - gammaln ( self . network . kappa ) <EOL> logit_p += gammaln ( self . mf_kappa_1 ) - self . mf_kappa_1 * np . log ( self . mf_v_1 ) <EOL> logit_p += gammaln ( self . kappa_0 ) - self . kappa_0 * np . log ( self . nu_0 ) <EOL> logit_p += self . mf_kappa_0 * np . log ( self . mf_v_0 ) - gammaln ( self . mf_kappa_0 ) <EOL> logit_p_hat = ( <NUM_LIT:1> - stepsize ) * logit ( self . mf_p ) + stepsize * logit_p <EOL> self . mf_p = np . clip ( logistic ( logit_p_hat ) , <NUM_LIT> , <NUM_LIT:1> - <NUM_LIT> ) <EOL> def meanfieldupdate_kappa_v ( self , data = [ ] , minibatchfrac = <NUM_LIT:1.0> , stepsize = <NUM_LIT:1.0> ) : <EOL> """<STR_LIT>""" <EOL> exp_ss = sum ( [ d . compute_exp_weight_ss ( ) for d in data ] ) <EOL> kappa0_hat = self . kappa_0 + exp_ss [ <NUM_LIT:0> ] / minibatchfrac <EOL> kappa1_hat = self . network . kappa + exp_ss [ <NUM_LIT:0> ] / minibatchfrac <EOL> self . mf_kappa_0 = ( <NUM_LIT:1.0> - stepsize ) * self . mf_kappa_0 + stepsize * kappa0_hat <EOL> self . mf_kappa_1 = ( <NUM_LIT:1.0> - stepsize ) * self . mf_kappa_1 + stepsize * kappa1_hat <EOL> v0_hat = self . nu_0 * np . ones ( ( self . K , self . K ) ) + exp_ss [ <NUM_LIT:1> ] / minibatchfrac <EOL> self . mf_v_0 = ( <NUM_LIT:1.0> - stepsize ) * self . mf_v_0 + stepsize * v0_hat <EOL> v1_hat = self . network . expected_v ( ) + exp_ss [ <NUM_LIT:1> ] / minibatchfrac <EOL> self . mf_v_1 = ( <NUM_LIT:1.0> - stepsize ) * self . mf_v_1 + stepsize * v1_hat <EOL> def meanfieldupdate ( self , data = [ ] ) : <EOL> self . meanfieldupdate_kappa_v ( data ) <EOL> self . meanfieldupdate_p ( ) <EOL> def meanfield_sgdstep ( self , data , minibatchfrac , stepsize ) : <EOL> self . meanfieldupdate_kappa_v ( data , minibatchfrac = minibatchfrac , stepsize = stepsize ) <EOL> self . meanfieldupdate_p ( stepsize = stepsize ) <EOL> def get_vlb ( self ) : <EOL> """<STR_LIT>""" <EOL> vlb = <NUM_LIT:0> <EOL> E_A = self . expected_A ( ) <EOL> E_notA = <NUM_LIT:1.0> - E_A <EOL> E_ln_p = self . network . expected_log_p ( ) <EOL> E_ln_notp = self . network . expected_log_notp ( ) <EOL> vlb += Bernoulli ( ) . negentropy ( E_x = E_A , E_notx = E_notA , <EOL> E_ln_p = E_ln_p , E_ln_notp = E_ln_notp ) . sum ( ) <EOL> kappa = self . network . kappa <EOL> E_v = self . network . expected_v ( ) <EOL> E_ln_v = self . network . expected_log_v ( ) <EOL> E_W1 = self . expected_W_given_A ( A = <NUM_LIT:1> ) <EOL> E_ln_W1 = self . expected_log_W_given_A ( A = <NUM_LIT:1> ) <EOL> vlb += ( E_A * Gamma ( kappa ) . negentropy ( E_beta = E_v , E_ln_beta = E_ln_v , <EOL> E_lambda = E_W1 , E_ln_lambda = E_ln_W1 ) ) . sum ( ) <EOL> kappa0 = self . kappa_0 <EOL> v0 = self . nu_0 <EOL> E_W0 = self . expected_W_given_A ( A = <NUM_LIT:0> ) <EOL> E_ln_W0 = self . expected_log_W_given_A ( A = <NUM_LIT:0> ) <EOL> vlb += ( E_notA * Gamma ( kappa0 , v0 ) . negentropy ( E_lambda = E_W0 , E_ln_lambda = E_ln_W0 ) ) . sum ( ) <EOL> vlb -= Bernoulli ( self . mf_p ) . negentropy ( ) . sum ( ) <EOL> vlb -= ( E_A * Gamma ( self . mf_kappa_1 , self . mf_v_1 ) . negentropy ( ) ) . sum ( ) <EOL> vlb -= ( E_notA * Gamma ( self . mf_kappa_0 , self . mf_v_0 ) . negentropy ( ) ) . sum ( ) <EOL> return vlb <EOL> def resample_from_mf ( self ) : <EOL> """<STR_LIT>""" <EOL> self . A = np . random . rand ( self . K , self . K ) < self . mf_p <EOL> self . W = ( <NUM_LIT:1> - self . A ) * np . random . gamma ( self . mf_kappa_0 , <NUM_LIT:1.0> / self . mf_v_0 ) <EOL> self . W += self . A * np . random . gamma ( self . mf_kappa_1 , <NUM_LIT:1.0> / self . mf_v_1 ) <EOL> def resample ( self , data = [ ] ) : <EOL> ss = np . zeros ( ( <NUM_LIT:2> , self . K , self . K ) ) + sum ( [ d . compute_weight_ss ( ) for d in data ] ) <EOL> self . _resample_A ( ss ) <EOL> self . _resample_W_given_A ( ss ) <EOL> def _resample_A ( self , ss ) : <EOL> """<STR_LIT>""" <EOL> p = self . network . P <EOL> v = self . network . V <EOL> kappa0_post = self . kappa_0 + ss [ <NUM_LIT:0> , : , : ] <EOL> v0_post = self . nu_0 + ss [ <NUM_LIT:1> , : , : ] <EOL> kappa1_post = self . network . kappa + ss [ <NUM_LIT:0> , : , : ] <EOL> v1_post = v + ss [ <NUM_LIT:1> , : , : ] <EOL> lp0 = self . kappa_0 * np . log ( self . nu_0 ) - gammaln ( self . kappa_0 ) <EOL> lp0 += gammaln ( kappa0_post ) - kappa0_post * np . log ( v0_post ) <EOL> lp1 = self . network . kappa * np . log ( v ) - gammaln ( self . network . kappa ) <EOL> lp1 += gammaln ( kappa1_post ) - kappa1_post * np . log ( v1_post ) <EOL> lp0 = lp0 + np . log ( <NUM_LIT:1.0> - p ) <EOL> lp1 = lp1 + np . log ( p ) <EOL> Z = logsumexp ( np . concatenate ( ( lp0 [ : , : , None ] , lp1 [ : , : , None ] ) , <EOL> axis = <NUM_LIT:2> ) , <EOL> axis = <NUM_LIT:2> ) <EOL> self . A = np . log ( np . random . rand ( self . K , self . K ) ) < lp1 - Z <EOL> def _resample_W_given_A ( self , ss ) : <EOL> kappa_prior = self . kappa_0 * ( <NUM_LIT:1> - self . A ) + self . network . kappa * self . A <EOL> kappa_cond = kappa_prior + ss [ <NUM_LIT:0> , : , : ] <EOL> v_prior = self . nu_0 * ( <NUM_LIT:1> - self . A ) + self . network . V * self . A <EOL> v_cond = v_prior + ss [ <NUM_LIT:1> , : , : ] <EOL> self . W = np . array ( np . random . gamma ( kappa_cond , <NUM_LIT:1.0> / v_cond ) ) . reshape ( ( self . K , self . K ) ) <EOL> self . W = np . clip ( self . W , <NUM_LIT> , np . inf ) <EOL> def initialize_from_gibbs ( self , A , W , scale = <NUM_LIT:100> ) : <EOL> """<STR_LIT>""" <EOL> conf = <NUM_LIT> <EOL> self . mf_p = conf * A + ( <NUM_LIT:1> - conf ) * ( <NUM_LIT:1> - A ) <EOL> self . mf_kappa_0 = self . kappa_0 <EOL> self . mf_v_0 = self . nu_0 <EOL> self . mf_kappa_1 = scale * W <EOL> self . mf_v_1 = scale <EOL> class SpikeAndSlabContinuousTimeGammaWeights ( GibbsSampling ) : <EOL> """<STR_LIT>""" <EOL> def __init__ ( self , model , parallel_resampling = True ) : <EOL> self . model = model <EOL> self . network = model . network <EOL> self . K = self . model . K <EOL> self . parallel_resampling = parallel_resampling <EOL> self . A = np . ones ( ( self . K , self . K ) ) <EOL> self . W = np . zeros ( ( self . K , self . K ) ) <EOL> self . resample ( ) <EOL> def log_likelihood ( self , x ) : <EOL> raise NotImplementedError <EOL> def log_probability ( self ) : <EOL> return <NUM_LIT:0> <EOL> def rvs ( self , size = [ ] ) : <EOL> raise NotImplementedError <EOL> @ property <EOL> def W_effective ( self ) : <EOL> return self . A * self . W <EOL> def _compute_weighted_impulses_at_events_manual ( self , data ) : <EOL> N , S , C , Z , dt_max = data . N , data . S , data . C , data . Z , self . model . dt_max <EOL> W = self . W <EOL> lmbda = np . zeros ( ( self . K , N ) ) <EOL> for n in xrange ( N ) : <EOL> for par in xrange ( n - <NUM_LIT:1> , - <NUM_LIT:1> , - <NUM_LIT:1> ) : <EOL> dt = S [ n ] - S [ par ] <EOL> if dt == <NUM_LIT:0> : <EOL> continue <EOL> if dt >= dt_max : <EOL> break <EOL> lmbda [ C [ par ] , n ] += W [ C [ par ] , C [ n ] ] * self . model . impulse_model . impulse ( dt , C [ par ] , C [ n ] ) <EOL> return lmbda <EOL> def _compute_weighted_impulses_at_events ( self , data ) : <EOL> from pyhawkes . internals . continuous_time_helpers import compute_weighted_impulses_at_events <EOL> N , S , C , Z , dt_max = data . N , data . S , data . C , data . Z , self . model . dt_max <EOL> W = self . W <EOL> mu , tau = self . model . impulse_model . mu , self . model . impulse_model . tau <EOL> lmbda = np . zeros ( ( N , self . K ) ) <EOL> compute_weighted_impulses_at_events ( S , C , Z , dt_max , W , mu , tau , lmbda ) <EOL> return lmbda <EOL> def _resample_A_given_W ( self , data ) : <EOL> """<STR_LIT>""" <EOL> lmbda_irs = [ self . _compute_weighted_impulses_at_events ( d ) for d in data ] <EOL> lmbda0 = self . model . lambda0 <EOL> def _log_likelihood_single_process ( k ) : <EOL> ll = <NUM_LIT:0> <EOL> for lmbda_ir , d in zip ( lmbda_irs , data ) : <EOL> Ns , C , T = d . Ns , d . C , d . T <EOL> ll -= lmbda0 [ k ] * T <EOL> ll -= self . W_effective [ : , k ] . dot ( Ns ) <EOL> ll += np . log ( lmbda0 [ k ] + np . sum ( self . A [ : , k ] * lmbda_ir [ C == k , : ] , axis = <NUM_LIT:1> ) ) . sum ( ) <EOL> return ll <EOL> p = self . network . P <EOL> for k1 in xrange ( self . K ) : <EOL> for k2 in xrange ( self . K ) : <EOL> if p [ k1 , k2 ] == <NUM_LIT:0.> : <EOL> self . A [ k1 , k2 ] = <NUM_LIT:0> <EOL> continue <EOL> if p [ k1 , k2 ] == <NUM_LIT:1.> : <EOL> self . A [ k1 , k2 ] = <NUM_LIT:1> <EOL> continue <EOL> self . A [ k1 , k2 ] = <NUM_LIT:0> <EOL> ll0 = _log_likelihood_single_process ( k2 ) <EOL> self . A [ k1 , k2 ] = <NUM_LIT:1> <EOL> ll1 = _log_likelihood_single_process ( k2 ) <EOL> lp0 = ll0 + np . log ( <NUM_LIT:1.0> - p [ k1 , k2 ] ) <EOL> lp1 = ll1 + np . log ( p [ k1 , k2 ] ) <EOL> Z = logsumexp ( [ lp0 , lp1 ] ) <EOL> self . A [ k1 , k2 ] = np . log ( np . random . rand ( ) ) < lp1 - Z <EOL> def _joblib_resample_A_given_W ( self , data ) : <EOL> """<STR_LIT>""" <EOL> import pyhawkes . internals . parallel_adjacency_resampling as par <EOL> par . model = self . model <EOL> par . data = data <EOL> par . lambda_irs = [ par . _compute_weighted_impulses_at_events ( d ) for d in data ] <EOL> if len ( data ) == <NUM_LIT:0> : <EOL> self . A = np . random . rand ( self . K , self . K ) < self . network . P <EOL> return <EOL> A_cols = Parallel ( n_jobs = - <NUM_LIT:1> , backend = "<STR_LIT>" ) ( <EOL> delayed ( par . _ct_resample_column_of_A ) ( k2 ) for k2 in range ( self . K ) ) <EOL> self . A = np . array ( A_cols ) . T <EOL> def resample_W_given_A_and_z ( self , N , Zsum ) : <EOL> """<STR_LIT>""" <EOL> kappa_post = self . network . kappa + Zsum <EOL> v_post = self . network . V + N [ : , None ] * self . A <EOL> self . W = np . array ( np . random . gamma ( kappa_post , <NUM_LIT:1.0> / v_post ) ) . reshape ( ( self . K , self . K ) ) <EOL> def resample ( self , data = [ ] ) : <EOL> """<STR_LIT>""" <EOL> assert isinstance ( data , list ) <EOL> N = np . zeros ( ( self . K , ) ) <EOL> Zsum = np . zeros ( ( self . K , self . K ) ) <EOL> for d in data : <EOL> Zsum += d . weight_ss <EOL> N += d . Ns <EOL> self . resample_W_given_A_and_z ( N , Zsum ) <EOL> if self . parallel_resampling : <EOL> self . _joblib_resample_A_given_W ( data ) <EOL> else : <EOL> self . _resample_A_given_W ( data ) </s>
