Ticket #3726: sage-3726-part2.patch

a/sage/stats/hmm/hmm.pyx

@@ -11 +11 @@
-
-from sage.matrix.all import is_Matrix, matrix
-from sage.rings.all  import RDF
+from sage.misc.randstate import random
-
-from sage.matrix.matrix_real_double_dense cimport Matrix_real_double_dense
-
@@ -57 +58 @@
-    ctypedef struct ghmm_alphabet:
-        pass
-    
+    # Discrete HMM's.
-    ctypedef struct ghmm_dmodel:
-        # Number of states
-        int N   
@@ -133 +135 @@
-
-    int ghmm_dmodel_normalize(ghmm_dmodel *m)
-    int ghmm_dmodel_free(ghmm_dmodel **m)
+    int ghmm_dmodel_logp(ghmm_dmodel *m, int *O, int len, double *log_p)
-
+
+    # Discrete sequences
+    ctypedef struct ghmm_dseq:
+        # sequence array. sequence[i] [j] = j-th symbol of i-th seq
+        int **seq
+        # matrix of state ids, can be used to save the viterbi path during sequence generation.
+        # ATTENTION: is NOT allocated by ghmm_dseq_calloc 
+        int **states
+        # array of sequence length
+        int *seq_len
+        # array of state path lengths
+        int *states_len
+        # array of sequence IDs
+        double *seq_id
+        # positive sequence weights.  default is 1 = no weight
+        double *seq_w
+        # total number of sequences 
+        long seq_number
+        # reserved space for sequences is always >= seq_number
+        long capacity
+        # sum of sequence weights
+        double total_w
+        # matrix of state labels corresponding to seq
+        int **state_labels
+        # number of labels for each sequence
+        int *state_labels_len
+        # internal flags
+        unsigned int flags
+
+    ghmm_dseq *ghmm_dmodel_label_generate_sequences(
+        ghmm_dmodel * mo, int seed, int global_len, long seq_number, int Tmax)
+    ghmm_dseq *ghmm_dmodel_generate_sequences(ghmm_dmodel* mo, int seed, int global_len,
+                                          long seq_number, int Tmax)
-
-cdef class HiddenMarkovModel:
-    pass
@@ -142 +178 @@
-cdef class DiscreteHiddenMarkovModel:
-    cdef ghmm_dmodel* m
-    cdef bint initialized
-
+
+
-    cdef Matrix_real_double_dense A, B
-    
-    def __init__(self, A, B, pi, emission_symbols=None, name=None):
@@ -178 +215 @@
-
-        cdef Py_ssize_t i, j, k
-
-
-
-
-
+
+
+
-
-
-
+
+
+
-        
-        # Set all pointers to NULL
-        self.m.s = NULL; self.m.name = NULL; self.m.silent = NULL
-        self.m.tied_to = NULL; self.m.order = NULL; self.m.background_id = NULL
-        self.m.bp = NULL; self.m.topo_order = NULL; self.m.pow_lookup = NULL;
- = NULL; self.m.label
+
-
-        # Set number of states and number of outputs
-        self.m.N = A.nrows()
-
+_
-        # Set the model type to discrete
-        self.m.model_type = GHMM_kDiscreteHMM
-
-        self.A = A
-        self.B = B
-
-        # Set that no a prior model probabilities are set.
-        self.m.prior = -1
@@ -289 +323 @@
-
-        self.initialized = True
-        
+    def __dealloc__(self):
+        if self.initialized:
+            ghmm_dmodel_free(&self.m)
-
-    def __repr__(self):
-        s = "Discrete Hidden Markov Model%s (%s states, %s outputs)\n"%(
@@ -300 +337 @@
-        return s
-
-    def initial_probabilities(self):
+        """
+        Return the list of initial state probabilities.
+        
+        EXAMPLES:
+            sage: a = hmm.DiscreteHiddenMarkovModel([[0.9,0.1],[0.9,0.1]], [[0.5,0.5,0,0],[0,0,.5,.5]], [0.5,0.5], [1,-1,3,-3])
+            sage: a.initial_probabilities()
+            [0.5, 0.5]
+        """
-        cdef Py_ssize_t i
-        return [self.m.s[i].pi for i in range(self.m.N)]
-
-    def transition_matrix(self, list_only=True):
+        """
+        Return the hidden state transition matrix. 
+        
+        EXAMPLES:
+            sage: a = hmm.DiscreteHiddenMarkovModel([[0.9,0.1],[0.9,0.1]], [[0.5,0.5,0,0],[0,0,.5,.5]], [0.5,0.5], [1,-1,3,-3])
+            sage: a.transition_matrix()
+            [0.9 0.1]
+            [0.9 0.1]
+        """
-        cdef Py_ssize_t i, j
-        for i from 0 <= i < self.m.N:
-            for j from 0 <= j < self.m.s[i].out_states:
@@ -311 +365 @@
-        return self.A
-
-    def emission_matrix(self, list_only=True):
+        """
+        Return the emission probability matrix.
+        
+        EXAMPLES:
+            sage: a = hmm.DiscreteHiddenMarkovModel([[0.9,0.1],[0.9,0.1]], [[0.5,0.5,0,0],[0,0,.5,.5]], [0.5,0.5], [1,-1,3,-3])
+            sage: a.emission_matrix()
+            [0.5 0.5 0.0 0.0]
+            [0.0 0.0 0.5 0.5]
+        """
-        cdef Py_ssize_t i, j
-        for i from 0 <= i < self.m.N:
-            for j from 0 <= j < self.B._ncols:
@@ -320 +383 @@
-    def normalize(self):
-        """
-        Normalize the transition and emission probabilities, if applicable.
+
+        EXAMPLES:
+            sage: a = hmm.DiscreteHiddenMarkovModel([[0.5,1],[1.2,0.9]], [[1,0.5],[0.5,1]], [0.1,1.2])
+            sage: a.normalize()
+            sage: a
+            Discrete Hidden Markov Model (2 states, 2 outputs)
+            Initial probabilities: [0.076923076923076927, 0.92307692307692302]
+            Transition matrix:
+            [0.333333333333 0.666666666667]
+            [0.571428571429 0.428571428571]
+            Emission matrix:
+            [0.666666666667 0.333333333333]
+            [0.333333333333 0.666666666667]
-        """
-        ghmm_dmodel_normalize(self.m)
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-    
-
-
-
-##################################################################################
-# Helper Functions