Ticket #3773: sage-3773-part3.patch

File sage-3773-part3.patch, 14.0 kB (added by was, 5 months ago)

  • a/sage/stats/hmm/chmm.pyx

    old new  
    128128        [(0.0, 1.0), (-1.0, 0.5), (1.0, 0.20000000000000001)] 
    129129        Initial probabilities: [1.0, 0.0, 0.0] 
    130130    """ 
    131     def __init__(self, A, B, pi, name=None): 
     131    def __init__(self, A, B, pi=None, name=None): 
    132132        """ 
    133133        EXAMPLES: 
    134134        We make a very simple model: 
    …… 
    166166        self.initialized = True 
    167167 
    168168    def _initialize_state(self, pi): 
    169         # Allocate and initialize states 
     169        """ 
     170        Allocate and initialize states. 
     171 
     172        INPUT: 
     173            pi -- initial probabilities 
     174 
     175        All other inputs are set as self.A and self.B. 
     176 
     177        EXAMPLES: 
     178        This function is called implicitly during object creation.  It should 
     179        never be called directly by the user, unless they want to LEAKE MEMORY. 
     180 
     181            sage: m = hmm.GaussianHiddenMarkovModel([[1]], [(1,10)], [1]) # indirect test 
     182        """ 
    170183        cdef ghmm_cstate* states = <ghmm_cstate*> safe_malloc(sizeof(ghmm_cstate) * self.m.N) 
    171184        cdef ghmm_cstate* state 
    172185        cdef ghmm_c_emission* e 
    …… 
    338351 
    339352        return 0 
    340353 
    341     def fix_emission_state(self, Py_ssize_t i, bint fixed=True): 
     354    def fix_emissions(self, Py_ssize_t i, bint fixed=True): 
    342355        """ 
    343         Sets the i-th emission state to be either fixed or not fixed. 
    344         If it is fixed, then running the Baum-Welch algorithm will not 
    345         change it. 
     356        Sets the Gaussian emission parameters for the i-th state to be 
     357        either fixed or not fixed.  If it is fixed, then running the 
     358        Baum-Welch algorithm will not change the emission parameters 
     359        for the i-th state. 
    346360         
    347361        INPUT: 
    348362            i -- nonnegative integer < self.m.N 
    …… 
    363377 
    364378        Now we run Baum-Welch with the emission states fixed.  Notice that they don't change.  
    365379            sage: m = hmm.GaussianHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [(0.0,1.0),(1,1)], [1,0]) 
    366             sage: m.fix_emission_state(0); m.fix_emission_state(1) 
     380            sage: m.fix_emissions(0); m.fix_emissions(1) 
    367381            sage: m.baum_welch([0,1]) 
    368382            sage: m 
    369383            Gaussian Hidden Markov Model with 2 States 
    …… 
    377391        if i < 0 or i >= self.m.N: 
    378392            raise IndexError, "index out of range" 
    379393        self.m.s[i].e.fixed = fixed 
    380  
    381     def fix_hidden_state(self, Py_ssize_t i, bint fixed=True): 
    382         """ 
    383         Sets the i-th hidden state to be either fixed or not fixed. 
    384         If it is fixed, then running the Baum-Welch algorithm will not 
    385         change it. 
    386          
    387         INPUT: 
    388             i -- nonnegative integer < self.m.N 
    389             fixed -- bool 
    390  
    391         EXAMPLES: 
    392         """ 
    393         if i < 0 or i >= self.m.N: 
    394             raise IndexError, "index out of range" 
    395         self.m.s[i].fix = fixed 
    396394 
    397395    def __repr__(self): 
    398396        """ 
    …… 
    789787        pi -- list of floats that sums to 1.0; these are 
    790788              the initial probabilities of each hidden state 
    791789        name -- (default: None); a string 
     790 
     791    EXAMPLES: 
     792        sage: A  = [[0.5,0.5],[0.5,0.5]] 
     793        sage: B  = [[(0.5,(0.0,1.0)), (0.1,(1,10000))],[(1,(1,1)), (0,(0,0.1))]] 
     794        sage: pi = [1,0] 
     795        sage: hmm.GaussianMixtureHiddenMarkovModel(A, B, pi) 
     796        Gaussian Hidden Markov Model with 2 States 
     797        Transition matrix: 
     798        [0.5 0.5] 
     799        [0.5 0.5] 
     800        Emission parameters: 
     801        [[(0.5, (0.0, 1.0)), (0.10000000000000001, (1.0, 10000.0))], [(1.0, (1.0, 1.0)), (0.0, (0.0, 0.10000000000000001))]] 
     802        Initial probabilities: [1.0, 0.0] 
     803 
     804    TESTS: 
     805    We test that standard deviations must be positive: 
     806        sage: m = hmm.GaussianMixtureHiddenMarkovModel([[1]], [[(1,(0,0))]], [1]) 
     807        Traceback (most recent call last): 
     808        ... 
     809        ValueError: sigma must be positive (if weight is nonzero) 
     810 
     811    We test that number of mixtures must be positive: 
     812        sage: m = hmm.GaussianMixtureHiddenMarkovModel([[1]], [[]], [1]) 
     813        Traceback (most recent call last): 
     814        ... 
     815        ValueError: number of Gaussian mixtures must be positive 
    792816    """ 
    793     def __init__(self, A, B, pi, name=None): 
     817    def __init__(self, A, B, pi=None, name=None): 
    794818        """ 
    795819        EXAMPLES: 
     820            sage: hmm.GaussianMixtureHiddenMarkovModel([[1]], [[(1,(0,1))]], [1], name='simple') 
     821            Gaussian Hidden Markov Model simple with 1 States 
     822            Transition matrix: 
     823            [1.0] 
     824            Emission parameters: 
     825            [[(1.0, (0.0, 1.0))]] 
     826            Initial probabilities: [1.0] 
    796827        """ 
    797828        # Turn B into a list of lists 
    798829        B = [flatten(x) for x in B] 
    …… 
    800831        if m == 0: 
    801832            raise ValueError, "number of Gaussian mixtures must be positive" 
    802833        B = [x + [0]*(m-len(x)) for x in B] 
    803         GaussianHiddenMarkovModel.__init__(self, A, B, pi) 
    804         print m//3 
     834        GaussianHiddenMarkovModel.__init__(self, A, B, pi, name=name) 
    805835        self.m.M = m//3 
    806836        # Set number of outputs.   
    807837 
    808838    def _initialize_state(self, pi): 
    809         # Allocate and initialize states 
     839        """ 
     840        Allocate and initialize states. 
     841 
     842        INPUT: 
     843            pi -- initial probabilities 
     844 
     845        All other inputs are set as self.A and self.B. 
     846 
     847        EXAMPLES: 
     848        This function is called implicitly during object creation.  It should 
     849        never be called directly by the user, unless they want to LEAKE MEMORY. 
     850 
     851            sage: m = hmm.GaussianMixtureHiddenMarkovModel([[1]], [[(1,(1,10))]], [1]) # indirect test 
     852        """ 
    810853        cdef ghmm_cstate* states = <ghmm_cstate*> safe_malloc(sizeof(ghmm_cstate) * self.m.N) 
    811854        cdef ghmm_cstate* state 
    812855        cdef ghmm_c_emission* e 
    …… 
    831874                e[n].dimension = 1 
    832875                mu             = v[n*3+1] 
    833876                sigma          = v[n*3+2] 
     877                if sigma <= 0 and v[n*3]: 
     878                    raise ValueError, "sigma must be positive (if weight is nonzero)" 
    834879                weights.append(  v[n*3] ) 
    835880                e[n].mean.val     = mu 
    836881                e[n].variance.val = sigma*sigma  # variance! not standard deviation 
    …… 
    890935                 for j in range(self.m.s[i].M)]  for i in range(self.m.N)] 
    891936 
    892937         
     938 
     939    def fix_emissions(self, Py_ssize_t i, Py_ssize_t j, bint fixed=True): 
     940        """ 
     941        Sets the j-th Gaussian of the emission parameters for the i-th 
     942        state to be either fixed or not fixed.  If it is fixed, then 
     943        running the Baum-Welch algorithm will not change the emission 
     944        parameters for the i-th state. 
     945         
     946        INPUT: 
     947            i -- nonnegative integer < self.m.N 
     948            j -- nonnegative integer < self.m.M 
     949            fixed -- bool 
     950 
     951        EXAMPLES: 
     952        We run Baum-Welch once without fixing the emission states: 
     953            sage: m = hmm.GaussianHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [(0.0,1.0),(1,1)], [1,0]) 
     954            sage: m.baum_welch([0,1]) 
     955            sage: m 
     956            Gaussian Hidden Markov Model with 2 States 
     957            Transition matrix: 
     958            [0.0 1.0] 
     959            [0.1 0.9] 
     960            Emission parameters: 
     961            [(0.0, 0.01), (1.0, 0.01)] 
     962            Initial probabilities: [1.0, 0.0] 
     963 
     964        Now we run Baum-Welch with the emission states fixed.  Notice that they don't change.  
     965            sage: m = hmm.GaussianHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [(0.0,1.0),(1,1)], [1,0]) 
     966            sage: m.fix_emissions(0); m.fix_emissions(1) 
     967            sage: m.baum_welch([0,1]) 
     968            sage: m 
     969            Gaussian Hidden Markov Model with 2 States 
     970            Transition matrix: 
     971            [0.000368587006957    0.999631412993] 
     972            [              0.1               0.9] 
     973            Emission parameters: 
     974            [(0.0, 1.0), (1.0, 1.0)] 
     975            Initial probabilities: [1.0, 0.0] 
     976        """ 
     977        if i < 0 or i >= self.m.N: 
     978            raise IndexError, "index i out of range" 
     979        if j < 0 or j >= self.m.M: 
     980            raise IndexError, "index j out of range" 
     981        self.m.s[i].e[j].fixed = fixed 

  • a/sage/stats/hmm/hmm.pyx

    old new  
    22Hidden Markov Models 
    33 
    44AUTHOR: William Stein 
    5  
    6 EXAMPLES: 
    7  
    8 TODO: 
    9    * make models pickleable (i.e., all parameters should be obtainable 
    10      using functions to make this easy). 
    11    * continuous hmm's 
    125""" 
    136 
    147############################################################################# 
    …… 
    4437        sage: sage.stats.hmm.hmm.HiddenMarkovModel([[0.4,0.6],[1,0]], [[1,0],[0.5,0.5]], [0.5,0.5]) 
    4538        <sage.stats.hmm.hmm.HiddenMarkovModel object at ...> 
    4639    """ 
    47     def __init__(self, A, B, pi): 
     40    def __init__(self, A, B, pi=None): 
    4841        """ 
    4942        INPUT: 
    5043            A -- matrix or list 
    …… 
    8477            B = B.change_ring(RDF) 
    8578 
    8679        # Make sure the initial probabilities are all floats.  
    87         self.pi = [float(x) for x in pi] 
    88         if len(self.pi) != A.nrows(): 
    89             raise ValueError, "length of pi must equal number of rows of A" 
     80        if pi is None: 
     81            if A.nrows() == 0: 
     82                self.pi = [] 
     83            else: 
     84                self.pi = [1.0/A.nrows()]*A.nrows() 
     85        else: 
     86            self.pi = [float(x) for x in pi] 
     87            if len(self.pi) != A.nrows(): 
     88                raise ValueError, "length of pi must equal number of rows of A" 
    9089 
    9190        # Record the now validated matrices A and B as attributes. 
    9291        # They get used later as attributes in the constructors for 
    …… 
    9796 
    9897cdef class DiscreteHiddenMarkovModel(HiddenMarkovModel): 
    9998     
    100     def __init__(self, A, B, pi, emission_symbols=None, name=None): 
     99    def __init__(self, A, B, pi=None, emission_symbols=None, name=None): 
    101100        """ 
    102101        INPUTS: 
    103102            A  -- square matrix of doubles; the state change probabilities 
    …… 
    292291        Deallocate memory allocated by the HMM. 
    293292 
    294293        EXAMPLES: 
    295             sage: a = hmm.DiscreteHiddenMarkovModel([[0.2,0.8],[0.5,0.5]], [[1,0],[0,1]], [0,1])  # indirect doctest 
    296             sage: del a 
     294            sage: m = hmm.DiscreteHiddenMarkovModel([[0.2,0.8],[0.5,0.5]], [[1,0],[0,1]], [0,1])  # indirect doctest 
     295            sage: del m 
    297296        """ 
    298297        if self.initialized: 
    299298            ghmm_dmodel_free(&self.m) 
     299 
     300    def fix_emissions(self, Py_ssize_t i, bint fixed=True): 
     301        """ 
     302        Sets the i-th emission parameters to be either fixed or not 
     303        fixed.  If it is fixed, then running the Baum-Welch algorithm 
     304        will not change the emission parameters for the i-th state. 
     305         
     306        INPUT: 
     307            i -- nonnegative integer < self.m.N 
     308            fixed -- bool 
     309 
     310        EXAMPLES: 
     311        First without calling fix_emissions: 
     312            sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]],[[.5,.5],[.5,.5]]) 
     313            sage: m.baum_welch([0,0,0,1,1,1]) 
     314            sage: m.emission_matrix() 
     315            [              1.0               0.0] 
     316            [3.92881039079e-05    0.999960711896] 
     317 
     318        We call fix_emissions on the first state and notice that the first 
     319        row of the emission matrix does not change: 
     320            sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]],[[.5,.5],[.5,.5]]) 
     321            sage: m.fix_emissions(0) 
     322            sage: m.baum_welch([0,0,0,1,1,1]) 
     323            sage: m.emission_matrix() 
     324            [              0.5               0.5] 
     325            [0.000542712675606    0.999457287324] 
     326 
     327        We call fix_emissions on the second state and notice that the second 
     328        row of the emission matrix does not change: 
     329            sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]],[[.5,.5],[.5,.5]]) 
     330            sage: m.fix_emissions(1) 
     331            sage: m.baum_welch([0,0,0,1,1,1]) 
     332            sage: m.emission_matrix() 
     333            [   0.999999904763 9.52366620142e-08] 
     334            [              0.5               0.5] 
     335 
     336        TESTS: 
     337        Make sure that out of range indices are handled correctly with an IndexError.  
     338            sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]],[[.5,.5],[.5,.5]]) 
     339            sage: m.fix_emissions(2) 
     340            Traceback (most recent call last): 
     341            ... 
     342            IndexError: index out of range 
     343            sage: m.fix_emissions(-1) 
     344            Traceback (most recent call last): 
     345            ... 
     346            IndexError: index out of range 
     347        """ 
     348        if i < 0 or i >= self.m.N: 
     349            raise IndexError, "index out of range" 
     350        self.m.s[i].fix = fixed 
    300351 
    301352    def __repr__(self): 
    302353        """ 
    …… 
    574625                    Viterbi path. 
    575626            float -- log of the probability that the sequence of hidden 
    576627                     states actually produced the given sequence seq. 
    577                      [[TODO: I do not understand precisely what this means.]] 
    578628 
    579629        EXAMPLES: 
    580630            sage: a = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.1,0.9]], [[0.9,0.1],[0.1,0.9]], [0.5,0.5])