From d444319e59ea3aa445e38f30450b9eea6600aeb3 Mon Sep 17 00:00:00 2001 From: Bruno Abreu Calfa Date: Wed, 10 Jan 2018 13:39:16 -0500 Subject: [PATCH 1/5] Added support for compilers that do not provide __builtin_expect (e.g., Visual Studio) --- pyhsmm/internals/util.h | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/pyhsmm/internals/util.h b/pyhsmm/internals/util.h index 4755b9b3..a000aa9c 100644 --- a/pyhsmm/internals/util.h +++ b/pyhsmm/internals/util.h @@ -1,8 +1,25 @@ #ifndef UTIL_H #define UTIL_H +#ifdef _WIN32 || _WIN64 +#define NO_BUILTIN_EXPECT +#endif + +#ifdef NO_BUILTIN_EXPECT +#if !defined(LIKELY) +#define likely(x) (!!(x),true) +#endif +#if !defined(UNLIKELY) +#define unlikely(x) (!!(x),false) +#endif +#else +#if !defined(LIKELY) #define likely(x) __builtin_expect(!!(x),true) +#endif +#if !defined(UNLIKELY) #define unlikely(x) __builtin_expect(!!(x),false) +#endif +#endif namespace util { using namespace std; From 911320b122d604ca7ca40212684ad9926bd5f6b2 Mon Sep 17 00:00:00 2001 From: Bruno Abreu Calfa Date: Thu, 11 Jan 2018 09:42:46 -0500 Subject: [PATCH 2/5] Added support for compilers that do not provide __builtin_expect --- pyhsmm/internals/util.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyhsmm/internals/util.h b/pyhsmm/internals/util.h index a000aa9c..af9d0605 100644 --- a/pyhsmm/internals/util.h +++ b/pyhsmm/internals/util.h @@ -1,7 +1,7 @@ #ifndef UTIL_H #define UTIL_H -#ifdef _WIN32 || _WIN64 +#ifdef _MSC_VER #define NO_BUILTIN_EXPECT #endif From a5747cdf7ea79c2f011be45b6c88e88ddd57dc10 Mon Sep 17 00:00:00 2001 From: Bruno Abreu Calfa Date: Sat, 13 Jan 2018 20:20:29 -0500 Subject: [PATCH 3/5] Added support for left-right models and example file --- examples/hsmm-left-right.py | 74 +++++++++++++++++++++++++++++++ pyhsmm/internals/initial_state.py | 3 +- pyhsmm/internals/transitions.py | 32 ++++++++----- pyhsmm/models.py | 18 ++++++-- 4 files changed, 110 insertions(+), 17 deletions(-) create mode 100644 examples/hsmm-left-right.py diff --git a/examples/hsmm-left-right.py b/examples/hsmm-left-right.py new file mode 100644 index 00000000..08e6a489 --- /dev/null +++ b/examples/hsmm-left-right.py @@ -0,0 +1,74 @@ + +import numpy as np +np.seterr(divide='ignore') # these warnings are usually harmless for this code +from matplotlib import pyplot as plt +import copy, os + +import pyhsmm +from pyhsmm.util.text import progprint_xrange + +SAVE_FIGURES = False + +print(''' +This demo shows a left-right HDP-HSMM in action. Its iterations are slower than those for +the (Sticky-)HDP-HMM, but explicit duration modeling can be a big advantage for +conditioning the prior or for discovering structure in data. +''') + +############### +# load data # +############### + +T = 1000 +data = np.loadtxt(os.path.join(os.path.dirname(__file__),'example-data.txt'))[:T] + +######################### +# posterior inference # +######################### + +# Set the weak limit truncation level +Nmax = 25 + +# Initialize left-right model parameters +trans_matrix = np.zeros((Nmax, Nmax)) +trans_matrix[-1, -1] = 1.0 +for i in range(Nmax - 1): + trans_matrix[i, i:] = 1.0 / (Nmax - i) # State i -> State i + 1, i + 2, ..., N + +pi_0 = np.zeros(Nmax) +pi_0[0] = 1.0 # Start in the first state + +# and some hyperparameters +obs_dim = data.shape[1] +obs_hypparams = {'mu_0':np.zeros(obs_dim), + 'sigma_0':np.eye(obs_dim), + 'kappa_0':0.25, + 'nu_0':obs_dim+2} +dur_hypparams = {'alpha_0':2*30, + 'beta_0':2} + +obs_distns = [pyhsmm.distributions.Gaussian(**obs_hypparams) for state in range(Nmax)] +dur_distns = [pyhsmm.distributions.PoissonDuration(**dur_hypparams) for state in range(Nmax)] + +posteriormodel = pyhsmm.models.WeakLimitHDPHSMM( + alpha=6.,gamma=6., # these can matter; see concentration-resampling.py + init_state_concentration=6., # pretty inconsequential + obs_distns=obs_distns, + dur_distns=dur_distns, + trans_matrix=trans_matrix, + pi_0=pi_0, + fix_trans_matrix_zeros=True, + fix_init_state_zeros=True) +posteriormodel.add_data(data,trunc=60) # duration truncation speeds things up when it's possible + +for idx in progprint_xrange(150): + posteriormodel.resample_model() + +print(np.array2string(posteriormodel.trans_distn.trans_matrix, + formatter={"float_kind": "{0:.4f}".format})) # Transition matrix +print(np.array2string(posteriormodel.init_state_distn.weights, + formatter={"float_kind": "{0:.4f}".format})) # Initial state probabilities + +posteriormodel.plot() + +plt.show() diff --git a/pyhsmm/internals/initial_state.py b/pyhsmm/internals/initial_state.py index ca0bba83..4a7cfefd 100644 --- a/pyhsmm/internals/initial_state.py +++ b/pyhsmm/internals/initial_state.py @@ -48,8 +48,9 @@ def copy_sample(self, new_model): return new class HMMInitialState(Categorical): - def __init__(self,model,init_state_concentration=None,pi_0=None): + def __init__(self,model,init_state_concentration=None,pi_0=None,fix_init_state_zeros=False): self.model = model + self.fix_init_state_zeros = fix_init_state_zeros if init_state_concentration is not None or pi_0 is not None: self._is_steady_state = False super(HMMInitialState,self).__init__( diff --git a/pyhsmm/internals/transitions.py b/pyhsmm/internals/transitions.py index 5921c66d..a853c58f 100644 --- a/pyhsmm/internals/transitions.py +++ b/pyhsmm/internals/transitions.py @@ -28,8 +28,9 @@ ### HMM class _HMMTransitionsBase(object): - def __init__(self,num_states=None,alpha=None,alphav=None,trans_matrix=None): + def __init__(self,num_states=None,alpha=None,alphav=None,trans_matrix=None, fix_trans_matrix_zeros=False): self.N = num_states + self.fix_trans_matrix_zeros = fix_trans_matrix_zeros if trans_matrix is not None: self._row_distns = [Multinomial(alpha_0=alpha,K=self.N,alphav_0=alphav, @@ -82,7 +83,13 @@ def resample(self,stateseqs=[],trans_counts=None): trans_counts = self._count_transitions(stateseqs) if trans_counts is None \ else trans_counts for distn, counts in zip(self._row_distns,trans_counts): + weights_orig = distn.weights distn.resample(counts) + + # Fix weights at zero (useful for prohibiting transitions by fixing probabilities) + if self.fix_trans_matrix_zeros: + distn.weights = np.where(weights_orig == 0.0, weights_orig, distn.weights) # Set to zero elements that were originally zeros + distn.weights /= np.sum(distn.weights) # Normalize return self class _HMMTransitionsMaxLikelihood(_HMMTransitionsBase): @@ -184,11 +191,12 @@ class HMMTransitionsConc(_ConcentrationResamplingMixin,_HMMTransitionsGibbs): class _HSMMTransitionsBase(_HMMTransitionsBase): def _get_trans_matrix(self): out = self.full_trans_matrix - out.flat[::out.shape[0]+1] = 0 - errs = np.seterr(invalid='ignore') - out /= out.sum(1)[:,na] - out = np.nan_to_num(out) - np.seterr(**errs) + if not self.fix_trans_matrix_zeros: + out.flat[::out.shape[0]+1] = 0 + errs = np.seterr(invalid='ignore') + out /= out.sum(1)[:,na] + out = np.nan_to_num(out) + np.seterr(**errs) return out trans_matrix = property(_get_trans_matrix,_HMMTransitionsBase.trans_matrix.fset) @@ -263,7 +271,7 @@ class HSMMTransitionsConc(_ConcentrationResamplingMixin,_HSMMTransitionsGibbs): ### HDP-HMM class _WeakLimitHDPHMMTransitionsBase(_HMMTransitionsBase): - def __init__(self,gamma,alpha,num_states=None,beta=None,trans_matrix=None): + def __init__(self,gamma,alpha,num_states=None,beta=None,trans_matrix=None,fix_trans_matrix_zeros=False): if num_states is None: assert beta is not None or trans_matrix is not None self.N = len(beta) if beta is not None else trans_matrix.shape[0] @@ -275,7 +283,7 @@ def __init__(self,gamma,alpha,num_states=None,beta=None,trans_matrix=None): super(_WeakLimitHDPHMMTransitionsBase,self).__init__( num_states=self.N,alpha=alpha, - alphav=alpha*self.beta,trans_matrix=trans_matrix) + alphav=alpha*self.beta,trans_matrix=trans_matrix,fix_trans_matrix_zeros=fix_trans_matrix_zeros) @property def beta(self): @@ -339,7 +347,7 @@ class WeakLimitHDPHMMTransitions(_WeakLimitHDPHMMTransitionsGibbs,_HMMTransition class _WeakLimitHDPHMMTransitionsConcBase(_WeakLimitHDPHMMTransitionsBase): def __init__(self,num_states,gamma_a_0,gamma_b_0,alpha_a_0,alpha_b_0, - beta=None,trans_matrix=None,**kwargs): + beta=None,trans_matrix=None,fix_trans_matrix_zeros=False,**kwargs): if num_states is None: assert beta is not None or trans_matrix is not None self.N = len(beta) if beta is not None else trans_matrix.shape[0] @@ -354,7 +362,7 @@ def __init__(self,num_states,gamma_a_0,gamma_b_0,alpha_a_0,alpha_b_0, # because it sets beta_obj in a different way _HMMTransitionsBase.__init__( self, num_states=self.N, alphav=self.alpha*self.beta, - trans_matrix=trans_matrix, **kwargs) + trans_matrix=trans_matrix, fix_trans_matrix_zeros=fix_trans_matrix_zeros, **kwargs) @property def alpha(self): @@ -440,7 +448,7 @@ class WeakLimitStickyHDPHMMTransitionsConc( class _DATruncHDPHMMTransitionsBase(_HMMTransitionsBase): # NOTE: self.beta stores \beta_{1:K}, so \beta_{\text{rest}} is implicit - def __init__(self,gamma,alpha,num_states,beta=None,trans_matrix=None): + def __init__(self,gamma,alpha,num_states,beta=None,trans_matrix=None,fix_trans_matrix_zeros=False): self.N = num_states self.gamma = gamma self._alpha = alpha @@ -452,7 +460,7 @@ def __init__(self,gamma,alpha,num_states,beta=None,trans_matrix=None): betafull = np.concatenate(((beta,(1.-beta.sum(),)))) super(_DATruncHDPHMMTransitionsBase,self).__init__( - num_states=self.N,alphav=alpha*betafull,trans_matrix=trans_matrix) + num_states=self.N,alphav=alpha*betafull,trans_matrix=trans_matrix,fix_trans_matrix_zeros=fix_trans_matrix_zeros) self.beta = beta diff --git a/pyhsmm/models.py b/pyhsmm/models.py index 517e74ec..52314916 100644 --- a/pyhsmm/models.py +++ b/pyhsmm/models.py @@ -37,9 +37,11 @@ def __init__(self, obs_distns, trans_distn=None, alpha=None,alpha_a_0=None,alpha_b_0=None,trans_matrix=None, - init_state_distn=None,init_state_concentration=None,pi_0=None): + init_state_distn=None,init_state_concentration=None,pi_0=None, + fix_init_state_zeros=False): self.obs_distns = obs_distns self.states_list = [] + self.fix_init_state_zeros = fix_init_state_zeros if trans_distn is not None: self.trans_distn = trans_distn @@ -61,7 +63,8 @@ def __init__(self, self.init_state_distn = self._init_state_class( model=self, init_state_concentration=init_state_concentration, - pi_0=pi_0) + pi_0=pi_0, + fix_init_state_zeros=fix_init_state_zeros) self._clear_caches() @@ -458,7 +461,13 @@ def resample_trans_distn(self): self._clear_caches() def resample_init_state_distn(self): + weights_orig = self.init_state_distn.weights self.init_state_distn.resample([s.stateseq[0] for s in self.states_list]) + + # Fix weights at zero (useful for constructing, for example, left-right models) + if self.fix_init_state_zeros: + self.init_state_distn.weights = np.where(weights_orig == 0.0, weights_orig, self.init_state_distn.weights) # Set to zero elements that were originally zeros + self.init_state_distn.weights /= np.sum(self.init_state_distn.weights) # Normalize self._clear_caches() def resample_states(self,num_procs=0): @@ -716,6 +725,7 @@ def __init__(self, obs_distns, trans_distn=None,alpha=None,alpha_a_0=None,alpha_b_0=None, gamma=None,gamma_a_0=None,gamma_b_0=None,trans_matrix=None, + fix_trans_matrix_zeros=False, **kwargs): if trans_distn is not None: @@ -725,11 +735,11 @@ def __init__(self, num_states=len(obs_distns), alpha_a_0=alpha_a_0,alpha_b_0=alpha_b_0, gamma_a_0=gamma_a_0,gamma_b_0=gamma_b_0, - trans_matrix=trans_matrix) + trans_matrix=trans_matrix,fix_trans_matrix_zeros=fix_trans_matrix_zeros) else: trans_distn = self._trans_class( num_states=len(obs_distns),alpha=alpha,gamma=gamma, - trans_matrix=trans_matrix) + trans_matrix=trans_matrix,fix_trans_matrix_zeros=fix_trans_matrix_zeros) super(_WeakLimitHDPMixin,self).__init__( obs_distns=obs_distns,trans_distn=trans_distn,**kwargs) From 4c807445686c156d6c1ebe11cf518a07ccde87b2 Mon Sep 17 00:00:00 2001 From: Bruno Abreu Calfa Date: Sat, 13 Jan 2018 20:25:30 -0500 Subject: [PATCH 4/5] Revert "Added support for left-right models and example file" This reverts commit a5747cdf7ea79c2f011be45b6c88e88ddd57dc10. --- examples/hsmm-left-right.py | 74 ------------------------------- pyhsmm/internals/initial_state.py | 3 +- pyhsmm/internals/transitions.py | 32 +++++-------- pyhsmm/models.py | 18 ++------ 4 files changed, 17 insertions(+), 110 deletions(-) delete mode 100644 examples/hsmm-left-right.py diff --git a/examples/hsmm-left-right.py b/examples/hsmm-left-right.py deleted file mode 100644 index 08e6a489..00000000 --- a/examples/hsmm-left-right.py +++ /dev/null @@ -1,74 +0,0 @@ - -import numpy as np -np.seterr(divide='ignore') # these warnings are usually harmless for this code -from matplotlib import pyplot as plt -import copy, os - -import pyhsmm -from pyhsmm.util.text import progprint_xrange - -SAVE_FIGURES = False - -print(''' -This demo shows a left-right HDP-HSMM in action. Its iterations are slower than those for -the (Sticky-)HDP-HMM, but explicit duration modeling can be a big advantage for -conditioning the prior or for discovering structure in data. -''') - -############### -# load data # -############### - -T = 1000 -data = np.loadtxt(os.path.join(os.path.dirname(__file__),'example-data.txt'))[:T] - -######################### -# posterior inference # -######################### - -# Set the weak limit truncation level -Nmax = 25 - -# Initialize left-right model parameters -trans_matrix = np.zeros((Nmax, Nmax)) -trans_matrix[-1, -1] = 1.0 -for i in range(Nmax - 1): - trans_matrix[i, i:] = 1.0 / (Nmax - i) # State i -> State i + 1, i + 2, ..., N - -pi_0 = np.zeros(Nmax) -pi_0[0] = 1.0 # Start in the first state - -# and some hyperparameters -obs_dim = data.shape[1] -obs_hypparams = {'mu_0':np.zeros(obs_dim), - 'sigma_0':np.eye(obs_dim), - 'kappa_0':0.25, - 'nu_0':obs_dim+2} -dur_hypparams = {'alpha_0':2*30, - 'beta_0':2} - -obs_distns = [pyhsmm.distributions.Gaussian(**obs_hypparams) for state in range(Nmax)] -dur_distns = [pyhsmm.distributions.PoissonDuration(**dur_hypparams) for state in range(Nmax)] - -posteriormodel = pyhsmm.models.WeakLimitHDPHSMM( - alpha=6.,gamma=6., # these can matter; see concentration-resampling.py - init_state_concentration=6., # pretty inconsequential - obs_distns=obs_distns, - dur_distns=dur_distns, - trans_matrix=trans_matrix, - pi_0=pi_0, - fix_trans_matrix_zeros=True, - fix_init_state_zeros=True) -posteriormodel.add_data(data,trunc=60) # duration truncation speeds things up when it's possible - -for idx in progprint_xrange(150): - posteriormodel.resample_model() - -print(np.array2string(posteriormodel.trans_distn.trans_matrix, - formatter={"float_kind": "{0:.4f}".format})) # Transition matrix -print(np.array2string(posteriormodel.init_state_distn.weights, - formatter={"float_kind": "{0:.4f}".format})) # Initial state probabilities - -posteriormodel.plot() - -plt.show() diff --git a/pyhsmm/internals/initial_state.py b/pyhsmm/internals/initial_state.py index 4a7cfefd..ca0bba83 100644 --- a/pyhsmm/internals/initial_state.py +++ b/pyhsmm/internals/initial_state.py @@ -48,9 +48,8 @@ def copy_sample(self, new_model): return new class HMMInitialState(Categorical): - def __init__(self,model,init_state_concentration=None,pi_0=None,fix_init_state_zeros=False): + def __init__(self,model,init_state_concentration=None,pi_0=None): self.model = model - self.fix_init_state_zeros = fix_init_state_zeros if init_state_concentration is not None or pi_0 is not None: self._is_steady_state = False super(HMMInitialState,self).__init__( diff --git a/pyhsmm/internals/transitions.py b/pyhsmm/internals/transitions.py index a853c58f..5921c66d 100644 --- a/pyhsmm/internals/transitions.py +++ b/pyhsmm/internals/transitions.py @@ -28,9 +28,8 @@ ### HMM class _HMMTransitionsBase(object): - def __init__(self,num_states=None,alpha=None,alphav=None,trans_matrix=None, fix_trans_matrix_zeros=False): + def __init__(self,num_states=None,alpha=None,alphav=None,trans_matrix=None): self.N = num_states - self.fix_trans_matrix_zeros = fix_trans_matrix_zeros if trans_matrix is not None: self._row_distns = [Multinomial(alpha_0=alpha,K=self.N,alphav_0=alphav, @@ -83,13 +82,7 @@ def resample(self,stateseqs=[],trans_counts=None): trans_counts = self._count_transitions(stateseqs) if trans_counts is None \ else trans_counts for distn, counts in zip(self._row_distns,trans_counts): - weights_orig = distn.weights distn.resample(counts) - - # Fix weights at zero (useful for prohibiting transitions by fixing probabilities) - if self.fix_trans_matrix_zeros: - distn.weights = np.where(weights_orig == 0.0, weights_orig, distn.weights) # Set to zero elements that were originally zeros - distn.weights /= np.sum(distn.weights) # Normalize return self class _HMMTransitionsMaxLikelihood(_HMMTransitionsBase): @@ -191,12 +184,11 @@ class HMMTransitionsConc(_ConcentrationResamplingMixin,_HMMTransitionsGibbs): class _HSMMTransitionsBase(_HMMTransitionsBase): def _get_trans_matrix(self): out = self.full_trans_matrix - if not self.fix_trans_matrix_zeros: - out.flat[::out.shape[0]+1] = 0 - errs = np.seterr(invalid='ignore') - out /= out.sum(1)[:,na] - out = np.nan_to_num(out) - np.seterr(**errs) + out.flat[::out.shape[0]+1] = 0 + errs = np.seterr(invalid='ignore') + out /= out.sum(1)[:,na] + out = np.nan_to_num(out) + np.seterr(**errs) return out trans_matrix = property(_get_trans_matrix,_HMMTransitionsBase.trans_matrix.fset) @@ -271,7 +263,7 @@ class HSMMTransitionsConc(_ConcentrationResamplingMixin,_HSMMTransitionsGibbs): ### HDP-HMM class _WeakLimitHDPHMMTransitionsBase(_HMMTransitionsBase): - def __init__(self,gamma,alpha,num_states=None,beta=None,trans_matrix=None,fix_trans_matrix_zeros=False): + def __init__(self,gamma,alpha,num_states=None,beta=None,trans_matrix=None): if num_states is None: assert beta is not None or trans_matrix is not None self.N = len(beta) if beta is not None else trans_matrix.shape[0] @@ -283,7 +275,7 @@ def __init__(self,gamma,alpha,num_states=None,beta=None,trans_matrix=None,fix_tr super(_WeakLimitHDPHMMTransitionsBase,self).__init__( num_states=self.N,alpha=alpha, - alphav=alpha*self.beta,trans_matrix=trans_matrix,fix_trans_matrix_zeros=fix_trans_matrix_zeros) + alphav=alpha*self.beta,trans_matrix=trans_matrix) @property def beta(self): @@ -347,7 +339,7 @@ class WeakLimitHDPHMMTransitions(_WeakLimitHDPHMMTransitionsGibbs,_HMMTransition class _WeakLimitHDPHMMTransitionsConcBase(_WeakLimitHDPHMMTransitionsBase): def __init__(self,num_states,gamma_a_0,gamma_b_0,alpha_a_0,alpha_b_0, - beta=None,trans_matrix=None,fix_trans_matrix_zeros=False,**kwargs): + beta=None,trans_matrix=None,**kwargs): if num_states is None: assert beta is not None or trans_matrix is not None self.N = len(beta) if beta is not None else trans_matrix.shape[0] @@ -362,7 +354,7 @@ def __init__(self,num_states,gamma_a_0,gamma_b_0,alpha_a_0,alpha_b_0, # because it sets beta_obj in a different way _HMMTransitionsBase.__init__( self, num_states=self.N, alphav=self.alpha*self.beta, - trans_matrix=trans_matrix, fix_trans_matrix_zeros=fix_trans_matrix_zeros, **kwargs) + trans_matrix=trans_matrix, **kwargs) @property def alpha(self): @@ -448,7 +440,7 @@ class WeakLimitStickyHDPHMMTransitionsConc( class _DATruncHDPHMMTransitionsBase(_HMMTransitionsBase): # NOTE: self.beta stores \beta_{1:K}, so \beta_{\text{rest}} is implicit - def __init__(self,gamma,alpha,num_states,beta=None,trans_matrix=None,fix_trans_matrix_zeros=False): + def __init__(self,gamma,alpha,num_states,beta=None,trans_matrix=None): self.N = num_states self.gamma = gamma self._alpha = alpha @@ -460,7 +452,7 @@ def __init__(self,gamma,alpha,num_states,beta=None,trans_matrix=None,fix_trans_m betafull = np.concatenate(((beta,(1.-beta.sum(),)))) super(_DATruncHDPHMMTransitionsBase,self).__init__( - num_states=self.N,alphav=alpha*betafull,trans_matrix=trans_matrix,fix_trans_matrix_zeros=fix_trans_matrix_zeros) + num_states=self.N,alphav=alpha*betafull,trans_matrix=trans_matrix) self.beta = beta diff --git a/pyhsmm/models.py b/pyhsmm/models.py index 52314916..517e74ec 100644 --- a/pyhsmm/models.py +++ b/pyhsmm/models.py @@ -37,11 +37,9 @@ def __init__(self, obs_distns, trans_distn=None, alpha=None,alpha_a_0=None,alpha_b_0=None,trans_matrix=None, - init_state_distn=None,init_state_concentration=None,pi_0=None, - fix_init_state_zeros=False): + init_state_distn=None,init_state_concentration=None,pi_0=None): self.obs_distns = obs_distns self.states_list = [] - self.fix_init_state_zeros = fix_init_state_zeros if trans_distn is not None: self.trans_distn = trans_distn @@ -63,8 +61,7 @@ def __init__(self, self.init_state_distn = self._init_state_class( model=self, init_state_concentration=init_state_concentration, - pi_0=pi_0, - fix_init_state_zeros=fix_init_state_zeros) + pi_0=pi_0) self._clear_caches() @@ -461,13 +458,7 @@ def resample_trans_distn(self): self._clear_caches() def resample_init_state_distn(self): - weights_orig = self.init_state_distn.weights self.init_state_distn.resample([s.stateseq[0] for s in self.states_list]) - - # Fix weights at zero (useful for constructing, for example, left-right models) - if self.fix_init_state_zeros: - self.init_state_distn.weights = np.where(weights_orig == 0.0, weights_orig, self.init_state_distn.weights) # Set to zero elements that were originally zeros - self.init_state_distn.weights /= np.sum(self.init_state_distn.weights) # Normalize self._clear_caches() def resample_states(self,num_procs=0): @@ -725,7 +716,6 @@ def __init__(self, obs_distns, trans_distn=None,alpha=None,alpha_a_0=None,alpha_b_0=None, gamma=None,gamma_a_0=None,gamma_b_0=None,trans_matrix=None, - fix_trans_matrix_zeros=False, **kwargs): if trans_distn is not None: @@ -735,11 +725,11 @@ def __init__(self, num_states=len(obs_distns), alpha_a_0=alpha_a_0,alpha_b_0=alpha_b_0, gamma_a_0=gamma_a_0,gamma_b_0=gamma_b_0, - trans_matrix=trans_matrix,fix_trans_matrix_zeros=fix_trans_matrix_zeros) + trans_matrix=trans_matrix) else: trans_distn = self._trans_class( num_states=len(obs_distns),alpha=alpha,gamma=gamma, - trans_matrix=trans_matrix,fix_trans_matrix_zeros=fix_trans_matrix_zeros) + trans_matrix=trans_matrix) super(_WeakLimitHDPMixin,self).__init__( obs_distns=obs_distns,trans_distn=trans_distn,**kwargs) From d6e04232989afde6c7c4b059a901f1862bd8b8ff Mon Sep 17 00:00:00 2001 From: Bruno Abreu Calfa Date: Sat, 13 Jan 2018 20:59:00 -0500 Subject: [PATCH 5/5] Added support for left-right models and example file --- examples/hsmm-left-right.py | 74 +++++++++++++++++++++++++++++++ pyhsmm/internals/initial_state.py | 3 +- pyhsmm/internals/transitions.py | 32 ++++++++----- pyhsmm/models.py | 18 ++++++-- 4 files changed, 110 insertions(+), 17 deletions(-) create mode 100644 examples/hsmm-left-right.py diff --git a/examples/hsmm-left-right.py b/examples/hsmm-left-right.py new file mode 100644 index 00000000..08e6a489 --- /dev/null +++ b/examples/hsmm-left-right.py @@ -0,0 +1,74 @@ + +import numpy as np +np.seterr(divide='ignore') # these warnings are usually harmless for this code +from matplotlib import pyplot as plt +import copy, os + +import pyhsmm +from pyhsmm.util.text import progprint_xrange + +SAVE_FIGURES = False + +print(''' +This demo shows a left-right HDP-HSMM in action. Its iterations are slower than those for +the (Sticky-)HDP-HMM, but explicit duration modeling can be a big advantage for +conditioning the prior or for discovering structure in data. +''') + +############### +# load data # +############### + +T = 1000 +data = np.loadtxt(os.path.join(os.path.dirname(__file__),'example-data.txt'))[:T] + +######################### +# posterior inference # +######################### + +# Set the weak limit truncation level +Nmax = 25 + +# Initialize left-right model parameters +trans_matrix = np.zeros((Nmax, Nmax)) +trans_matrix[-1, -1] = 1.0 +for i in range(Nmax - 1): + trans_matrix[i, i:] = 1.0 / (Nmax - i) # State i -> State i + 1, i + 2, ..., N + +pi_0 = np.zeros(Nmax) +pi_0[0] = 1.0 # Start in the first state + +# and some hyperparameters +obs_dim = data.shape[1] +obs_hypparams = {'mu_0':np.zeros(obs_dim), + 'sigma_0':np.eye(obs_dim), + 'kappa_0':0.25, + 'nu_0':obs_dim+2} +dur_hypparams = {'alpha_0':2*30, + 'beta_0':2} + +obs_distns = [pyhsmm.distributions.Gaussian(**obs_hypparams) for state in range(Nmax)] +dur_distns = [pyhsmm.distributions.PoissonDuration(**dur_hypparams) for state in range(Nmax)] + +posteriormodel = pyhsmm.models.WeakLimitHDPHSMM( + alpha=6.,gamma=6., # these can matter; see concentration-resampling.py + init_state_concentration=6., # pretty inconsequential + obs_distns=obs_distns, + dur_distns=dur_distns, + trans_matrix=trans_matrix, + pi_0=pi_0, + fix_trans_matrix_zeros=True, + fix_init_state_zeros=True) +posteriormodel.add_data(data,trunc=60) # duration truncation speeds things up when it's possible + +for idx in progprint_xrange(150): + posteriormodel.resample_model() + +print(np.array2string(posteriormodel.trans_distn.trans_matrix, + formatter={"float_kind": "{0:.4f}".format})) # Transition matrix +print(np.array2string(posteriormodel.init_state_distn.weights, + formatter={"float_kind": "{0:.4f}".format})) # Initial state probabilities + +posteriormodel.plot() + +plt.show() diff --git a/pyhsmm/internals/initial_state.py b/pyhsmm/internals/initial_state.py index ca0bba83..4a7cfefd 100644 --- a/pyhsmm/internals/initial_state.py +++ b/pyhsmm/internals/initial_state.py @@ -48,8 +48,9 @@ def copy_sample(self, new_model): return new class HMMInitialState(Categorical): - def __init__(self,model,init_state_concentration=None,pi_0=None): + def __init__(self,model,init_state_concentration=None,pi_0=None,fix_init_state_zeros=False): self.model = model + self.fix_init_state_zeros = fix_init_state_zeros if init_state_concentration is not None or pi_0 is not None: self._is_steady_state = False super(HMMInitialState,self).__init__( diff --git a/pyhsmm/internals/transitions.py b/pyhsmm/internals/transitions.py index 5921c66d..a853c58f 100644 --- a/pyhsmm/internals/transitions.py +++ b/pyhsmm/internals/transitions.py @@ -28,8 +28,9 @@ ### HMM class _HMMTransitionsBase(object): - def __init__(self,num_states=None,alpha=None,alphav=None,trans_matrix=None): + def __init__(self,num_states=None,alpha=None,alphav=None,trans_matrix=None, fix_trans_matrix_zeros=False): self.N = num_states + self.fix_trans_matrix_zeros = fix_trans_matrix_zeros if trans_matrix is not None: self._row_distns = [Multinomial(alpha_0=alpha,K=self.N,alphav_0=alphav, @@ -82,7 +83,13 @@ def resample(self,stateseqs=[],trans_counts=None): trans_counts = self._count_transitions(stateseqs) if trans_counts is None \ else trans_counts for distn, counts in zip(self._row_distns,trans_counts): + weights_orig = distn.weights distn.resample(counts) + + # Fix weights at zero (useful for prohibiting transitions by fixing probabilities) + if self.fix_trans_matrix_zeros: + distn.weights = np.where(weights_orig == 0.0, weights_orig, distn.weights) # Set to zero elements that were originally zeros + distn.weights /= np.sum(distn.weights) # Normalize return self class _HMMTransitionsMaxLikelihood(_HMMTransitionsBase): @@ -184,11 +191,12 @@ class HMMTransitionsConc(_ConcentrationResamplingMixin,_HMMTransitionsGibbs): class _HSMMTransitionsBase(_HMMTransitionsBase): def _get_trans_matrix(self): out = self.full_trans_matrix - out.flat[::out.shape[0]+1] = 0 - errs = np.seterr(invalid='ignore') - out /= out.sum(1)[:,na] - out = np.nan_to_num(out) - np.seterr(**errs) + if not self.fix_trans_matrix_zeros: + out.flat[::out.shape[0]+1] = 0 + errs = np.seterr(invalid='ignore') + out /= out.sum(1)[:,na] + out = np.nan_to_num(out) + np.seterr(**errs) return out trans_matrix = property(_get_trans_matrix,_HMMTransitionsBase.trans_matrix.fset) @@ -263,7 +271,7 @@ class HSMMTransitionsConc(_ConcentrationResamplingMixin,_HSMMTransitionsGibbs): ### HDP-HMM class _WeakLimitHDPHMMTransitionsBase(_HMMTransitionsBase): - def __init__(self,gamma,alpha,num_states=None,beta=None,trans_matrix=None): + def __init__(self,gamma,alpha,num_states=None,beta=None,trans_matrix=None,fix_trans_matrix_zeros=False): if num_states is None: assert beta is not None or trans_matrix is not None self.N = len(beta) if beta is not None else trans_matrix.shape[0] @@ -275,7 +283,7 @@ def __init__(self,gamma,alpha,num_states=None,beta=None,trans_matrix=None): super(_WeakLimitHDPHMMTransitionsBase,self).__init__( num_states=self.N,alpha=alpha, - alphav=alpha*self.beta,trans_matrix=trans_matrix) + alphav=alpha*self.beta,trans_matrix=trans_matrix,fix_trans_matrix_zeros=fix_trans_matrix_zeros) @property def beta(self): @@ -339,7 +347,7 @@ class WeakLimitHDPHMMTransitions(_WeakLimitHDPHMMTransitionsGibbs,_HMMTransition class _WeakLimitHDPHMMTransitionsConcBase(_WeakLimitHDPHMMTransitionsBase): def __init__(self,num_states,gamma_a_0,gamma_b_0,alpha_a_0,alpha_b_0, - beta=None,trans_matrix=None,**kwargs): + beta=None,trans_matrix=None,fix_trans_matrix_zeros=False,**kwargs): if num_states is None: assert beta is not None or trans_matrix is not None self.N = len(beta) if beta is not None else trans_matrix.shape[0] @@ -354,7 +362,7 @@ def __init__(self,num_states,gamma_a_0,gamma_b_0,alpha_a_0,alpha_b_0, # because it sets beta_obj in a different way _HMMTransitionsBase.__init__( self, num_states=self.N, alphav=self.alpha*self.beta, - trans_matrix=trans_matrix, **kwargs) + trans_matrix=trans_matrix, fix_trans_matrix_zeros=fix_trans_matrix_zeros, **kwargs) @property def alpha(self): @@ -440,7 +448,7 @@ class WeakLimitStickyHDPHMMTransitionsConc( class _DATruncHDPHMMTransitionsBase(_HMMTransitionsBase): # NOTE: self.beta stores \beta_{1:K}, so \beta_{\text{rest}} is implicit - def __init__(self,gamma,alpha,num_states,beta=None,trans_matrix=None): + def __init__(self,gamma,alpha,num_states,beta=None,trans_matrix=None,fix_trans_matrix_zeros=False): self.N = num_states self.gamma = gamma self._alpha = alpha @@ -452,7 +460,7 @@ def __init__(self,gamma,alpha,num_states,beta=None,trans_matrix=None): betafull = np.concatenate(((beta,(1.-beta.sum(),)))) super(_DATruncHDPHMMTransitionsBase,self).__init__( - num_states=self.N,alphav=alpha*betafull,trans_matrix=trans_matrix) + num_states=self.N,alphav=alpha*betafull,trans_matrix=trans_matrix,fix_trans_matrix_zeros=fix_trans_matrix_zeros) self.beta = beta diff --git a/pyhsmm/models.py b/pyhsmm/models.py index 517e74ec..52314916 100644 --- a/pyhsmm/models.py +++ b/pyhsmm/models.py @@ -37,9 +37,11 @@ def __init__(self, obs_distns, trans_distn=None, alpha=None,alpha_a_0=None,alpha_b_0=None,trans_matrix=None, - init_state_distn=None,init_state_concentration=None,pi_0=None): + init_state_distn=None,init_state_concentration=None,pi_0=None, + fix_init_state_zeros=False): self.obs_distns = obs_distns self.states_list = [] + self.fix_init_state_zeros = fix_init_state_zeros if trans_distn is not None: self.trans_distn = trans_distn @@ -61,7 +63,8 @@ def __init__(self, self.init_state_distn = self._init_state_class( model=self, init_state_concentration=init_state_concentration, - pi_0=pi_0) + pi_0=pi_0, + fix_init_state_zeros=fix_init_state_zeros) self._clear_caches() @@ -458,7 +461,13 @@ def resample_trans_distn(self): self._clear_caches() def resample_init_state_distn(self): + weights_orig = self.init_state_distn.weights self.init_state_distn.resample([s.stateseq[0] for s in self.states_list]) + + # Fix weights at zero (useful for constructing, for example, left-right models) + if self.fix_init_state_zeros: + self.init_state_distn.weights = np.where(weights_orig == 0.0, weights_orig, self.init_state_distn.weights) # Set to zero elements that were originally zeros + self.init_state_distn.weights /= np.sum(self.init_state_distn.weights) # Normalize self._clear_caches() def resample_states(self,num_procs=0): @@ -716,6 +725,7 @@ def __init__(self, obs_distns, trans_distn=None,alpha=None,alpha_a_0=None,alpha_b_0=None, gamma=None,gamma_a_0=None,gamma_b_0=None,trans_matrix=None, + fix_trans_matrix_zeros=False, **kwargs): if trans_distn is not None: @@ -725,11 +735,11 @@ def __init__(self, num_states=len(obs_distns), alpha_a_0=alpha_a_0,alpha_b_0=alpha_b_0, gamma_a_0=gamma_a_0,gamma_b_0=gamma_b_0, - trans_matrix=trans_matrix) + trans_matrix=trans_matrix,fix_trans_matrix_zeros=fix_trans_matrix_zeros) else: trans_distn = self._trans_class( num_states=len(obs_distns),alpha=alpha,gamma=gamma, - trans_matrix=trans_matrix) + trans_matrix=trans_matrix,fix_trans_matrix_zeros=fix_trans_matrix_zeros) super(_WeakLimitHDPMixin,self).__init__( obs_distns=obs_distns,trans_distn=trans_distn,**kwargs)