<s> """<STR_LIT>""" <EOL> from __future__ import print_function <EOL> import os <EOL> import re <EOL> import gzip <EOL> import time <EOL> import cPickle <EOL> import operator <EOL> import collections <EOL> import numpy as np <EOL> from scipy . misc import logsumexp <EOL> from sklearn . feature_extraction . text import CountVectorizer <EOL> import matplotlib . pyplot as plt <EOL> from hips . plotting . layout import create_axis_at_location , create_figure <EOL> import brewer2mpl <EOL> from pgmult . lds import MultinomialLDS <EOL> from pgmult . particle_lds import LogisticNormalMultinomialLDS , ParticleSBMultinomialLDS <EOL> from pgmult . hmm import MultinomialHMM <EOL> from pgmult . utils import pi_to_psi <EOL> from pylds . models import DefaultLDS , NonstationaryLDS <EOL> from pybasicbayes . distributions import GaussianFixed , Multinomial , Regression <EOL> from pybasicbayes . util . text import progprint_xrange <EOL> from autoregressive . distributions import AutoRegression <EOL> colors = brewer2mpl . get_map ( "<STR_LIT>" , "<STR_LIT>" , <NUM_LIT:9> ) . mpl_colors <EOL> goodcolors = np . array ( [ <NUM_LIT:0> , <NUM_LIT:1> , <NUM_LIT:2> , <NUM_LIT:4> , <NUM_LIT:6> , <NUM_LIT:7> , <NUM_LIT:8> ] ) <EOL> colors = np . array ( colors ) [ goodcolors ] <EOL> np . seterr ( invalid = "<STR_LIT>" ) <EOL> np . random . seed ( <NUM_LIT:0> ) <EOL> np . seterr ( invalid = "<STR_LIT>" ) <EOL> np . random . seed ( <NUM_LIT:0> ) <EOL> K = <NUM_LIT:200> <EOL> def load ( filename = None ) : <EOL> if filename is None : <EOL> bigstr = download_ap ( ) <EOL> else : <EOL> with open ( filename , '<STR_LIT:r>' ) as infile : <EOL> bigstr = infile . read ( ) <EOL> docs = re . findall ( r'<STR_LIT>' , bigstr . translate ( None , '<STR_LIT:\n>' ) ) <EOL> vectorizer = CountVectorizer ( stop_words = '<STR_LIT>' , max_features = K ) . fit ( docs ) <EOL> docs = [ make_onehot_seq ( doc , vectorizer ) for doc in docs ] <EOL> words = vectorizer . get_feature_names ( ) <EOL> return docs , words <EOL> def download_ap ( ) : <EOL> from cStringIO import StringIO <EOL> from urllib2 import urlopen <EOL> import tarfile <EOL> print ( "<STR_LIT>" ) <EOL> response = urlopen ( '<STR_LIT>' ) <EOL> tar = tarfile . open ( fileobj = StringIO ( response . read ( ) ) ) <EOL> return tar . extractfile ( '<STR_LIT>' ) . read ( ) <EOL> def filter_wordseq ( doc , vectorizer ) : <EOL> return [ w for w in doc if w in vectorizer . vocabulary_ ] <EOL> def make_onehot_seq ( doc , vectorizer ) : <EOL> lst = filter_wordseq ( vectorizer . build_analyzer ( ) ( doc ) , vectorizer ) <EOL> indices = { word : idx for idx , word in enumerate ( vectorizer . vocabulary_ . keys ( ) ) } <EOL> out = np . zeros ( ( len ( lst ) , len ( indices ) ) ) <EOL> for wordidx , word in enumerate ( lst ) : <EOL> out [ wordidx , indices [ word ] ] = <NUM_LIT:1> <EOL> return out <EOL> Results = collections . namedtuple ( "<STR_LIT>" , [ "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" ] ) <EOL> def fit_lds_model ( Xs , Xtest , D , N_samples = <NUM_LIT:100> ) : <EOL> Nx = len ( Xs ) <EOL> assert len ( Xtest ) == Nx <EOL> model = MultinomialLDS ( K , D , <EOL> init_dynamics_distn = GaussianFixed ( mu = np . zeros ( D ) , sigma = <NUM_LIT:1> * np . eye ( D ) ) , <EOL> dynamics_distn = AutoRegression ( nu_0 = D + <NUM_LIT:1> , S_0 = <NUM_LIT:1> * np . eye ( D ) , M_0 = np . zeros ( ( D , D ) ) , K_0 = <NUM_LIT:1> * np . eye ( D ) ) , <EOL> sigma_C = <NUM_LIT:1.> <EOL> ) <EOL> for X in Xs : <EOL> model . add_data ( X ) <EOL> model . resample_parameters ( ) <EOL> compute_pred_ll = lambda : sum ( [ model . predictive_log_likelihood ( Xt , data_index = i , M = <NUM_LIT:10> ) [ <NUM_LIT:0> ] <EOL> for i , Xt in enumerate ( Xtest ) ] ) <EOL> init_results = ( <NUM_LIT:0> , None , model . log_likelihood ( ) , <EOL> np . nan , <EOL> compute_pred_ll ( ) ) <EOL> def resample ( ) : <EOL> tic = time . time ( ) <EOL> model . resample_model ( ) <EOL> toc = time . time ( ) - tic <EOL> return toc , None , model . log_likelihood ( ) , np . nan , compute_pred_ll ( ) <EOL> times , samples , lls , test_lls , pred_lls = map ( np . array , zip ( * ( [ init_results ] + <EOL> [ resample ( ) for _ in progprint_xrange ( N_samples , perline = <NUM_LIT:5> ) ] ) ) ) <EOL> timestamps = np . cumsum ( times ) <EOL> return Results ( lls , test_lls , pred_lls , samples , timestamps ) <EOL> def fit_hmm ( Xs , Xtest , D_hmm , N_samples = <NUM_LIT:100> ) : <EOL> print ( "<STR_LIT>" % D_hmm ) <EOL> model = MultinomialHMM ( K , D_hmm , alpha_0 = <NUM_LIT> ) <EOL> for X in Xs : <EOL> model . add_data ( X ) <EOL> compute_pred_ll = lambda : sum ( [ model . log_likelihood ( np . vstack ( ( Xs [ i ] , Xtest [ i ] ) ) ) <EOL> - model . log_likelihood ( Xs [ i ] ) <EOL> for i , Xt in enumerate ( Xtest ) ] ) <EOL> init_results = ( <NUM_LIT:0> , None , model . log_likelihood ( ) , <EOL> model . log_likelihood ( Xtest ) , <EOL> compute_pred_ll ( ) ) <EOL> def resample ( ) : <EOL> tic = time . time ( ) <EOL> model . resample_model ( ) <EOL> toc = time . time ( ) - tic <EOL> return toc , None , model . log_likelihood ( ) , np . nan , compute_pred_ll ( ) <EOL> times , samples , lls , test_lls , pred_lls = map ( np . array , zip ( * ( [ init_results ] + <EOL> [ resample ( ) for _ in progprint_xrange ( N_samples , perline = <NUM_LIT:5> ) ] ) ) ) <EOL> timestamps = np . cumsum ( times ) <EOL> return Results ( lls , test_lls , pred_lls , samples , timestamps ) <EOL> def fit_gaussian_lds_model ( Xs , Xtest , D_gauss_lds , N_samples = <NUM_LIT:100> ) : <EOL> print ( "<STR_LIT>" % D_gauss_lds ) <EOL> from pylds . models import NonstationaryLDS <EOL> model = NonstationaryLDS ( <EOL> init_dynamics_distn = GaussianFixed ( mu = np . zeros ( D ) , sigma = <NUM_LIT:1> * np . eye ( D ) ) , <EOL> dynamics_distn = AutoRegression ( nu_0 = D + <NUM_LIT:1> , S_0 = <NUM_LIT:1> * np . eye ( D ) , M_0 = np . zeros ( ( D , D ) ) , K_0 = <NUM_LIT:1> * np . eye ( D ) ) , <EOL> emission_distn = Regression ( nu_0 = K + <NUM_LIT:1> , S_0 = K * np . eye ( K ) , M_0 = np . zeros ( ( K , D ) ) , K_0 = K * np . eye ( D ) ) ) <EOL> Xs_centered = [ X - np . mean ( X , axis = <NUM_LIT:0> ) [ None , : ] + <NUM_LIT> * np . random . randn ( * X . shape ) for X in Xs ] <EOL> for X in Xs_centered : <EOL> model . add_data ( X ) <EOL> def compute_pred_ll ( ) : <EOL> pred_ll = <NUM_LIT:0> <EOL> for i , Xt in enumerate ( Xtest ) : <EOL> Npred = <NUM_LIT:10> <EOL> Tpred = Xt . shape [ <NUM_LIT:0> ] <EOL> preds = model . sample_predictions ( Xs_centered [ i ] , Tpred , Npred = Npred ) <EOL> inds = np . argmax ( preds , axis = <NUM_LIT:1> ) <EOL> pi = np . array ( [ np . bincount ( inds [ t ] , minlength = K ) for t in xrange ( Tpred ) ] ) / float ( Npred ) <EOL> assert np . allclose ( pi . sum ( axis = <NUM_LIT:1> ) , <NUM_LIT:1.0> ) <EOL> pi = np . clip ( pi , <NUM_LIT> , <NUM_LIT:1.0> ) <EOL> pi /= pi . sum ( axis = <NUM_LIT:1> ) [ : , None ] <EOL> pred_ll += np . sum ( [ Multinomial ( weights = pi [ t ] , K = K ) . log_likelihood ( Xt [ t ] [ None , : ] ) <EOL> for t in xrange ( Tpred ) ] ) <EOL> return pred_ll <EOL> init_results = ( <NUM_LIT:0> , None , np . nan , np . nan , compute_pred_ll ( ) ) <EOL> def resample ( ) : <EOL> tic = time . time ( ) <EOL> model . resample_model ( ) <EOL> toc = time . time ( ) - tic <EOL> return toc , None , np . nan , np . nan , compute_pred_ll ( ) <EOL> times , samples , lls , test_lls , pred_lls = map ( np . array , zip ( * ( [ init_results ] + <EOL> [ resample ( ) for _ in progprint_xrange ( N_samples , perline = <NUM_LIT:5> ) ] ) ) ) <EOL> timestamps = np . cumsum ( times ) <EOL> return Results ( lls , test_lls , pred_lls , samples , timestamps ) <EOL> def fit_ln_lds_model ( Xs , Xtest , D , N_samples = <NUM_LIT:100> ) : <EOL> """<STR_LIT>""" <EOL> print ( "<STR_LIT>" % D ) <EOL> model = LogisticNormalMultinomialLDS ( <EOL> init_dynamics_distn = GaussianFixed ( mu = np . zeros ( D ) , sigma = <NUM_LIT:1> * np . eye ( D ) ) , <EOL> dynamics_distn = AutoRegression ( nu_0 = D + <NUM_LIT:1> , S_0 = D * np . eye ( D ) , M_0 = np . zeros ( ( D , D ) ) , K_0 = D * np . eye ( D ) ) , <EOL> emission_distn = Regression ( nu_0 = K + <NUM_LIT:1> , S_0 = K * np . eye ( K ) , M_0 = np . zeros ( ( K , D ) ) , K_0 = K * np . eye ( D ) ) ) <EOL> model . A = <NUM_LIT:0.5> * np . eye ( D ) <EOL> model . sigma_states = np . eye ( D ) <EOL> model . C = <NUM_LIT> * np . random . randn ( K , D ) <EOL> model . sigma_obs = <NUM_LIT:0.1> * np . eye ( K ) <EOL> for X in Xs : <EOL> model . add_data ( X ) <EOL> compute_pred_ll = lambda : sum ( [ model . predictive_log_likelihood ( Xt , data_index = i , Npred = <NUM_LIT:10> ) [ <NUM_LIT:0> ] <EOL> for i , Xt in enumerate ( Xtest ) ] ) <EOL> init_results = ( <NUM_LIT:0> , None , model . log_likelihood ( ) , <EOL> np . nan , compute_pred_ll ( ) ) <EOL> def resample ( ) : <EOL> tic = time . time ( ) <EOL> model . resample_model ( ) <EOL> toc = time . time ( ) - tic <EOL> return toc , None , model . log_likelihood ( ) , np . nan , compute_pred_ll ( ) <EOL> times , samples , lls , test_lls , pred_lls = map ( np . array , zip ( * ( [ init_results ] + <EOL> [ resample ( ) for _ in progprint_xrange ( N_samples , perline = <NUM_LIT:5> ) ] ) ) ) <EOL> timestamps = np . cumsum ( times ) <EOL> return Results ( lls , test_lls , pred_lls , samples , timestamps ) <EOL> def fit_lds_model_with_pmcmc ( Xs , Xtest , D , N_samples = <NUM_LIT:100> ) : <EOL> """<STR_LIT>""" <EOL> print ( "<STR_LIT>" % D ) <EOL> model = ParticleSBMultinomialLDS ( <EOL> init_dynamics_distn = GaussianFixed ( mu = np . zeros ( D ) , sigma = <NUM_LIT:1> * np . eye ( D ) ) , <EOL> dynamics_distn = AutoRegression ( nu_0 = D + <NUM_LIT:1> , S_0 = D * np . eye ( D ) , M_0 = np . zeros ( ( D , D ) ) , K_0 = D * np . eye ( D ) ) , <EOL> emission_distn = Regression ( nu_0 = K + <NUM_LIT:1> , S_0 = K * np . eye ( K ) , M_0 = np . zeros ( ( K , D ) ) , K_0 = K * np . eye ( D ) ) , <EOL> mu = pi_to_psi ( np . ones ( K ) / K ) ) <EOL> model . A = <NUM_LIT:0.5> * np . eye ( D ) <EOL> model . sigma_states = np . eye ( D ) <EOL> model . C = np . random . randn ( K - <NUM_LIT:1> , D ) <EOL> model . sigma_obs = <NUM_LIT:0.1> * np . eye ( K ) <EOL> for X in Xs : <EOL> model . add_data ( X ) <EOL> compute_pred_ll = lambda : sum ( [ model . predictive_log_likelihood ( Xt , data_index = i , Npred = <NUM_LIT:10> ) [ <NUM_LIT:0> ] <EOL> for i , Xt in enumerate ( Xtest ) ] ) <EOL> init_results = ( <NUM_LIT:0> , None , model . log_likelihood ( ) , <EOL> np . nan , compute_pred_ll ( ) ) <EOL> def resample ( ) : <EOL> tic = time . time ( ) <EOL> model . resample_model ( ) <EOL> toc = time . time ( ) - tic <EOL> return toc , None , model . log_likelihood ( ) , np . nan , compute_pred_ll ( ) <EOL> times , samples , lls , test_lls , pred_lls = map ( np . array , zip ( * ( [ init_results ] + [ resample ( ) for _ in progprint_xrange ( N_samples ) ] ) ) ) <EOL> timestamps = np . cumsum ( times ) <EOL> return Results ( lls , test_lls , pred_lls , samples , timestamps ) <EOL> def plot_log_likelihood ( results , names , run , outname = "<STR_LIT>" ) : <EOL> plt . figure ( figsize = ( <NUM_LIT:3> , <NUM_LIT> ) ) <EOL> for i , ( result , name ) in enumerate ( zip ( results , names ) ) : <EOL> plt . plot ( result . timestamps , result . lls , lw = <NUM_LIT:2> , color = colors [ i ] , label = name ) <EOL> plt . legend ( loc = "<STR_LIT>" ) <EOL> plt . xlabel ( '<STR_LIT>' ) <EOL> plt . ylabel ( "<STR_LIT>" ) <EOL> plt . savefig ( os . path . join ( "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" % run , outname ) ) <EOL> plt . tight_layout ( ) <EOL> def plot_pred_log_likelihood ( results , names , run , outname = "<STR_LIT>" , smooth = True ) : <EOL> plt . figure ( figsize = ( <NUM_LIT:3> , <NUM_LIT> ) ) <EOL> for i , ( result , name ) in enumerate ( zip ( results , names ) ) : <EOL> if result . pred_lls . ndim == <NUM_LIT:2> : <EOL> pred_ll = result . pred_lls [ : , <NUM_LIT:0> ] <EOL> else : <EOL> pred_ll = result . pred_lls <EOL> if smooth : <EOL> win = <NUM_LIT:10> <EOL> pad_pred_ll = np . concatenate ( ( pred_ll [ <NUM_LIT:0> ] * np . ones ( win ) , pred_ll ) ) <EOL> smooth_pred_ll = np . array ( [ logsumexp ( pad_pred_ll [ j - win : j + <NUM_LIT:1> ] ) - np . log ( win ) <EOL> for j in xrange ( win , pad_pred_ll . size ) ] ) <EOL> plt . plot ( np . clip ( result . timestamps , <NUM_LIT> , np . inf ) , smooth_pred_ll , <EOL> lw = <NUM_LIT:2> , color = colors [ i ] , label = name ) <EOL> else : <EOL> plt . plot ( np . clip ( result . timestamps , <NUM_LIT> , np . inf ) , result . pred_lls , <EOL> lw = <NUM_LIT:2> , color = colors [ i ] , label = name ) <EOL> plt . xlabel ( '<STR_LIT>' ) <EOL> plt . xscale ( "<STR_LIT>" ) <EOL> plt . ylabel ( "<STR_LIT>" ) <EOL> plt . ylim ( - <NUM_LIT> , - <NUM_LIT> ) <EOL> plt . savefig ( os . path . join ( "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" % run , outname ) ) <EOL> plt . tight_layout ( ) <EOL> def plot_pred_ll_vs_D ( all_results , Ds , Xtrain , Xtest , run ) : <EOL> N = len ( Ds ) <EOL> M = len ( all_results [ <NUM_LIT:0> ] ) <EOL> T = len ( all_results [ <NUM_LIT:0> ] [ <NUM_LIT:0> ] . pred_lls ) <EOL> pred_lls = np . zeros ( ( N , M , T ) ) <EOL> for n in xrange ( N ) : <EOL> for m in xrange ( M ) : <EOL> if all_results [ n ] [ m ] . pred_lls . ndim == <NUM_LIT:2> : <EOL> pred_lls [ n , m ] = all_results [ n ] [ m ] . pred_lls [ : , <NUM_LIT:0> ] <EOL> else : <EOL> pred_lls [ n , m ] = all_results [ n ] [ m ] . pred_lls <EOL> burnin = T // <NUM_LIT:2> <EOL> pred_ll_mean = logsumexp ( pred_lls [ : , : , burnin : ] , axis = - <NUM_LIT:1> ) - np . log ( T - burnin ) <EOL> pred_ll_std = np . zeros_like ( pred_ll_mean ) <EOL> for n in xrange ( N ) : <EOL> for m in xrange ( M ) : <EOL> samples = np . random . choice ( pred_lls [ n , m , burnin : ] , size = ( <NUM_LIT:100> , ( T - burnin ) ) , replace = True ) <EOL> pll_samples = logsumexp ( samples , axis = <NUM_LIT:1> ) - np . log ( T - burnin ) <EOL> pred_ll_std [ n , m ] = pll_samples . std ( ) <EOL> baseline = <NUM_LIT:0> <EOL> for Xtr , Xte in zip ( Xtrain , Xtest ) : <EOL> pi_emp = Xtr . sum ( <NUM_LIT:0> ) / float ( Xtr . sum ( ) ) <EOL> pi_emp = np . clip ( pi_emp , <NUM_LIT> , np . inf ) <EOL> pi_emp /= pi_emp . sum ( ) <EOL> baseline += Multinomial ( weights = pi_emp , K = Xtr . shape [ <NUM_LIT:1> ] ) . log_likelihood ( Xte ) . sum ( ) <EOL> fig = plt . figure ( figsize = ( <NUM_LIT:3> , <NUM_LIT:3> ) ) <EOL> ax = fig . add_subplot ( <NUM_LIT> ) <EOL> width = np . min ( np . diff ( Ds ) ) / ( M + <NUM_LIT:1.0> ) if len ( Ds ) > <NUM_LIT:1> else <NUM_LIT:1.> <EOL> for m in xrange ( M ) : <EOL> ax . bar ( Ds + m * width , <EOL> pred_ll_mean [ : , m ] - baseline , <EOL> yerr = pred_ll_std [ : , m ] , <EOL> width = <NUM_LIT> * width , color = colors [ m ] , ecolor = '<STR_LIT:k>' ) <EOL> ax . plot ( [ Ds . min ( ) - width , Ds . max ( ) + ( M + <NUM_LIT:1> ) * width ] , np . zeros ( <NUM_LIT:2> ) , '<STR_LIT>' ) <EOL> ax . set_xlim ( Ds . min ( ) - width , Ds . max ( ) + ( M + <NUM_LIT:1> ) * width ) <EOL> ax . set_xticks ( Ds + ( M * width ) / <NUM_LIT> ) <EOL> ax . set_xticklabels ( Ds ) <EOL> ax . set_xlabel ( "<STR_LIT:D>" ) <EOL> ax . set_ylabel ( "<STR_LIT>" ) <EOL> ax . set_title ( "<STR_LIT>" ) <EOL> plt . savefig ( os . path . join ( "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" % run , "<STR_LIT>" ) ) <EOL> def fit_joint_corpus ( ) : <EOL> run = <NUM_LIT> <EOL> results_dir = os . path . join ( "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" % run ) <EOL> from pgmult . internals . utils import mkdir <EOL> if not os . path . exists ( results_dir ) : <EOL> print ( "<STR_LIT>" , results_dir ) <EOL> mkdir ( results_dir ) <EOL> Xs , words = load ( ) <EOL> N_docs = <NUM_LIT:10> <EOL> T_split = <NUM_LIT:10> <EOL> Xfilt = filter ( lambda X : X . shape [ <NUM_LIT:0> ] > <NUM_LIT:2> * T_split , Xs ) <EOL> Xtrain = [ X [ : - T_split ] for X in Xfilt [ : N_docs ] ] <EOL> Xtest = [ X [ - T_split : ] for X in Xfilt [ : N_docs ] ] <EOL> N_samples = <NUM_LIT> <EOL> all_results = [ ] <EOL> Ds = np . array ( [ <NUM_LIT:10> ] ) <EOL> models = [ "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" ] <EOL> methods = [ fit_lds_model , fit_hmm , fit_gaussian_lds_model , fit_ln_lds_model , fit_lds_model_with_pmcmc ] <EOL> for D in Ds : <EOL> D_results = [ ] <EOL> for model , method in zip ( models , methods ) : <EOL> results_file = os . path . join ( results_dir , "<STR_LIT>" % ( model , D ) ) <EOL> if os . path . exists ( results_file ) : <EOL> print ( "<STR_LIT>" , results_file ) <EOL> with gzip . open ( results_file , "<STR_LIT:r>" ) as f : <EOL> D_model_results = cPickle . load ( f ) <EOL> else : <EOL> print ( "<STR_LIT>" , model , "<STR_LIT>" , D ) <EOL> D_model_results = method ( Xtrain , Xtest , D , N_samples ) <EOL> with gzip . open ( results_file , "<STR_LIT:w>" ) as f : <EOL> print ( "<STR_LIT>" , results_file ) <EOL> cPickle . dump ( D_model_results , f , protocol = - <NUM_LIT:1> ) <EOL> D_results . append ( D_model_results ) <EOL> all_results . append ( D_results ) <EOL> res_index = <NUM_LIT:0> <EOL> plot_log_likelihood ( all_results [ res_index ] , <EOL> models , <EOL> run , <EOL> outname = "<STR_LIT>" % Ds [ res_index ] ) <EOL> plot_pred_log_likelihood ( all_results [ res_index ] , <EOL> models , <EOL> run , <EOL> outname = "<STR_LIT>" % Ds [ res_index ] ) <EOL> plot_pred_ll_vs_D ( all_results , Ds , Xtrain , Xtest , run ) <EOL> plt . show ( ) <EOL> if __name__ == "<STR_LIT:__main__>" : <EOL> fit_joint_corpus ( ) </s>
