<s> import numpy as np <EOL> import emcee <EOL> import logging <EOL> from scipy . stats import norm <EOL> from robo . acquisition . log_ei import LogEI <EOL> from robo . acquisition . base import AcquisitionFunction <EOL> from robo . initial_design . init_random_uniform import init_random_uniform <EOL> from robo . util import mc_part <EOL> logger = logging . getLogger ( __name__ ) <EOL> class InformationGainMC ( AcquisitionFunction ) : <EOL> def __init__ ( self , model , X_lower , X_upper , <EOL> Nb = <NUM_LIT:50> , Nf = <NUM_LIT> , <EOL> sampling_acquisition = None , <EOL> sampling_acquisition_kw = { "<STR_LIT>" : <NUM_LIT:0.0> } , <EOL> Np = <NUM_LIT:50> , ** kwargs ) : <EOL> """<STR_LIT>""" <EOL> self . Nb = Nb <EOL> super ( InformationGainMC , self ) . __init__ ( model , X_lower , X_upper ) <EOL> self . D = self . X_lower . shape [ <NUM_LIT:0> ] <EOL> self . sn2 = None <EOL> if sampling_acquisition is None : <EOL> sampling_acquisition = LogEI <EOL> self . sampling_acquisition = sampling_acquisition ( <EOL> model , self . X_lower , self . X_upper , ** sampling_acquisition_kw ) <EOL> self . Nf = Nf <EOL> self . Np = Np <EOL> def compute ( self , X , derivative = False , * args ) : <EOL> if derivative : <EOL> raise NotImplementedError <EOL> new_pmin = self . change_pmin_by_innovation ( X ) <EOL> H = - np . sum ( np . multiply ( np . exp ( self . logP ) , ( self . logP + self . lmb ) ) ) <EOL> dHp = - np . sum ( np . multiply ( new_pmin , <EOL> np . add ( np . log ( new_pmin ) , self . lmb ) ) , axis = <NUM_LIT:0> ) - H <EOL> return - np . array ( [ dHp ] ) <EOL> def sample_representer_points ( self ) : <EOL> self . sampling_acquisition . update ( self . model ) <EOL> start_points = init_random_uniform ( self . X_lower , <EOL> self . X_upper , <EOL> self . Nb ) <EOL> sampler = emcee . EnsembleSampler ( self . Nb , <EOL> self . D , <EOL> self . sampling_acquisition_wrapper ) <EOL> self . zb , self . lmb , _ = sampler . run_mcmc ( start_points , <NUM_LIT:200> ) <EOL> if len ( self . zb . shape ) == <NUM_LIT:1> : <EOL> self . zb = self . zb [ : , None ] <EOL> if len ( self . lmb . shape ) == <NUM_LIT:1> : <EOL> self . lmb = self . lmb [ : , None ] <EOL> def sampling_acquisition_wrapper ( self , x ) : <EOL> if np . any ( x < self . X_lower ) or np . any ( x > self . X_upper ) : <EOL> return - np . inf <EOL> return self . sampling_acquisition ( np . array ( [ x ] ) ) [ <NUM_LIT:0> ] <EOL> def update ( self , model ) : <EOL> self . model = model <EOL> self . sn2 = self . model . get_noise ( ) <EOL> self . sampling_acquisition . update ( model ) <EOL> self . sample_representer_points ( ) <EOL> self . W = norm . ppf ( np . linspace ( <NUM_LIT:1.> / ( self . Np + <NUM_LIT:1> ) , <EOL> <NUM_LIT:1> - <NUM_LIT:1.> / ( self . Np + <NUM_LIT:1> ) , <EOL> self . Np ) ) [ np . newaxis , : ] <EOL> self . Mb , self . Vb = self . model . predict ( self . zb , full_cov = True ) <EOL> self . pmin = mc_part . joint_pmin ( self . Mb , self . Vb , self . Nf ) <EOL> self . logP = np . log ( self . pmin ) <EOL> def innovations ( self , x , rep ) : <EOL> _ , v = self . model . predict ( x ) <EOL> v_ = v - self . sn2 <EOL> sigma_x_rep = self . model . predict_variance ( rep , x ) <EOL> norm_cov = np . dot ( sigma_x_rep , np . linalg . inv ( v_ ) ) <EOL> dm_rep = np . dot ( norm_cov , <EOL> np . linalg . cholesky ( v + <NUM_LIT> ) ) <EOL> dm_rep = dm_rep . dot ( self . W ) <EOL> dv_rep = - norm_cov . dot ( sigma_x_rep . T ) <EOL> return dm_rep , dv_rep <EOL> def change_pmin_by_innovation ( self , x ) : <EOL> dmdb , dvdb = self . innovations ( x , self . zb ) <EOL> Mb_new = self . Mb + dmdb <EOL> Vb_new = self . Vb + dvdb <EOL> return mc_part . joint_pmin ( Mb_new , Vb_new , self . Nf ) </s>
