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/internals/util.h b/pyhsmm/internals/util.h index 4755b9b3..af9d0605 100644 --- a/pyhsmm/internals/util.h +++ b/pyhsmm/internals/util.h @@ -1,8 +1,25 @@ #ifndef UTIL_H #define UTIL_H +#ifdef _MSC_VER +#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; 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)