Added support for left-right models

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()

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__(

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
 

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;

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)