Ticket #3726: sage-3726-part4.patch

a/sage/stats/hmm/all.py

@@ -1 @@
-
+
+

/dev/null

@@ +1 @@
+
+from hmm cimport HiddenMarkovModel
+
+cdef extern from "ghmm/ghmm.h":
+    cdef int GHMM_kContinuousHMM
+
+cdef extern from "ghmm/smodel.h":
+    #############################################################
+    # Continuous emission density types
+    #############################################################
+    cdef enum ghmm_density_t:
+        normal,
+        normal_right,   # right tail
+        normal_approx,
+        normal_left,    # left tail
+        uniform,
+        binormal,
+        multinormal,
+        density_number
+
+    #############################################################
+    # Continuous emission
+    #############################################################
+    # Mean and variance types
+    cdef union mean_t:
+        double val
+        double *vec
+        
+    cdef union variance_t:
+        double val
+        double *mat
+
+    cdef struct ghmm_c_emission:
+        # Flag for density function for each component of the mixture
+        #   0: normal density
+        #   1: truncated normal (right side) density
+        #   2: approximated normal density
+        #   3: truncated normal (left side)
+        #   4: uniform distribution
+        #   6: multivariate normal
+        int type
+        # dimension > 1 for multivariate normals
+        int dimension
+        # mean for output functions (pointer to mean vector for multivariate)
+        mean_t mean
+        variance_t variance
+        # pointer to inverse of covariance matrix if multivariate normal;
+        # else NULL 
+        double *sigmainv
+        # determinant of covariance matrix for multivariate normal 
+        double det
+        # Cholesky decomposition of covariance matrix A,
+        #   if A = G*G' sigmacd only holds G
+        double *sigmacd
+        # minimum of uniform distribution or left boundary for
+        # right-tail gaussians
+        double min
+        # maximum of uniform distribution or right boundary for
+        # left-tail gaussians
+        double max
+        # if fixed != 0 the parameters of the density are fixed 
+        int fixed
+        
+    #############################################################
+    # Continous Emission States
+    #############################################################
+    ctypedef struct ghmm_cstate:
+        # Number of output densities per state 
+        int M
+        # initial prob.
+        double pi
+        # IDs of successor states
+        int *out_id
+        # IDs of predecessor states
+        int *in_id
+        # transition probs to successor states. It is a
+        # matrix in case of mult. transition matrices (COS > 1)
+        double **out_a
+        # transition probs from predecessor states. It is a
+        # matrix in case of mult. transition matrices (COS > 1)
+        double **in_a
+        # number of  successor states
+        int out_states
+        # number of  predecessor states
+        int in_states
+        # weight vector for output function components
+        double *c
+        # flag for fixation of parameter. If fix = 1 do not change parameters of
+        # output functions, and if fix = 0 do normal training. Default is 0.
+        int fix
+        # vector of ghmm_c_emission
+        # (type and parameters of output function components)
+        ghmm_c_emission *e
+        # contains a description of the state (null terminated utf-8)
+        char *desc
+        # x coordinate position for graph representation plotting
+        int xPosition
+        # y coordinate position for graph representation plotting
+        int yPosition
+
+    #############################################################
+    # Continous Hidden Markov Model
+    #############################################################
+    ctypedef struct ghmm_cmodel
+        
+    ctypedef struct ghmm_cmodel_class_change_context:
+        # Names of class change module/function (for python callback)
+        char *python_module
+        char *python_function
+        # index of current sequence 
+        int k
+        # pointer to class function 
+        int (*get_class) (ghmm_cmodel*, double*, int, int)
+        # space for any data necessary for class switch, USER is RESPONSIBLE
+        void *user_data
+
+    ctypedef struct ghmm_cmodel:
+        # Number of states
+        int N
+        # Maximun number of components in the states
+        int M
+        # Number of dimensions of the emission components.
+        # NOTE: All emissions must have the same number of dimensions.
+        int dim
+        # The ghmm_cmodel includes two continuous models
+        #   cos=1: one transition matrix
+        #   cos>1: (integer) extension for models with several matrices.
+        int cos
+        # Prior for a priori probability of the model.
+        # -1 means no prior specified (all models have equal prob. a priori.)
+        double prior
+        # Contains a arbitrary name for the model (null terminated utf-8)
+        char * name
+        # Contains bit flags for various model extensions such as
+        # kSilentStates (see ghmm.h for a complete list).
+        int model_type
+        # All states of the model. Transition probabilities are part of the states.
+        ghmm_cstate *s
+        # pointer to a ghmm_cmodel_class_change_context struct
+        # necessary for multiple transition classes
+        ghmm_cmodel_class_change_context *class_change
+
+    int ghmm_cmodel_class_change_alloc(ghmm_cmodel * smo)
+
+    #############################################################
+    # Continous Sequence of states
+    #
+    # Sequence structure for double sequences.  Contains an array of
+    # sequences and corresponding data like sequence labels, sequence
+    # weights, etc. Sequences may have different length.
+    # multi-dimension sequences are linearized
+    #############################################################
+    ctypedef struct ghmm_cseq:
+        # sequence array. sequence[i][j] = j-th symbol of i-th seq.
+        # sequence[i][D * j] = first dimension of j-th observation of i-th sequence
+        double **seq
+        # array of sequence length
+        int *seq_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
+        # total number of dimensions 
+        int dim
+        # flags (internal)
+        unsigned int flags
+
+
+    #############################################################
+    # Continous Baum-Welch algorithm context
+    #############################################################
+    ctypedef struct ghmm_cmodel_baum_welch_context:
+        # pointer of continuous model
+        ghmm_cmodel *smo
+        # sequence pointer
+        ghmm_cseq *sqd
+        # calculated log likelihood 
+        double *logp
+        # leave reestimation loop if diff. between successive logp values
+        # is smaller than eps 
+        double eps
+        # max. no of iterations
+        int max_iter
+
+
+    int ghmm_cmodel_free (ghmm_cmodel ** smo)
+    
+        
+cdef class ContinuousHiddenMarkovModel(HiddenMarkovModel):
+    cdef ghmm_cmodel* m
+    cdef bint initialized
+
+cdef class NormalHiddenMarkovModel(ContinuousHiddenMarkovModel):
+    pass
+
+

/dev/null

@@ +1 @@
+
+
+include "../../ext/stdsage.pxi"
+
+include "misc.pxi"
+
+cdef class ContinuousHiddenMarkovModel(HiddenMarkovModel):
+    def __init__(self, A, B, pi, name=None):
+        self.initialized = False
+        HiddenMarkovModel.__init__(self, A, B, pi)
+        self.m = <ghmm_cmodel*> safe_malloc(sizeof(ghmm_cmodel))
+        # Set number of states
+        self.m.N = self.A.nrows()
+
+cdef class NormalHiddenMarkovModel(ContinuousHiddenMarkovModel):
+    """
+    EXAMPLES:
+    The transition matrix:
+        sage: A = [[0,1,0],[0.5,0,0.5],[0.3,0.3,0.4]]
+
+    Parameters of the normal emission distributions in pairs of (mu, sigma):
+        sage: B = [(0,1), (-1,0.5), (1,0.2)]
+
+    The initial probabilities per state:
+        sage: pi = [1,0,0]
+
+    Create the continuous HMM:
+        sage: m = hmm.NormalHiddenMarkovModel(A, B, pi); m
+    """
+    def __init__(self, A, B, pi, name=None):
+        ContinuousHiddenMarkovModel.__init__(self, A, B, pi)
+
+        # Set number of outputs
+        self.m.M = 1
+
+        # Set the model type to continuous
+        self.m.model_type = GHMM_kContinuousHMM
+
+        # Assign model identifier if specified
+        if name is not None: 
+            name = str(name)
+            self.m.name = name
+        else:
+            self.m.name = NULL
+            
+        # 1 transition matrix
+        self.m.cos   =  1
+        # Set that no a prior model probabilities are set.
+        self.m.prior = -1
+        # Dimension is 1
+        self.m.dim   =  1
+
+        # Allocate and initialize states
+        cdef ghmm_cstate* states = <ghmm_cstate*> safe_malloc(sizeof(ghmm_cstate) * self.m.N)
+        cdef ghmm_cstate* state
+        cdef ghmm_c_emission* e
+
+        for i in range(self.m.N):
+            # Parameters of normal distribution
+            mu, sigma = self.B[i]
+            # Get a reference to the i-th state for convenience of the notation below.
+            state = &(states[i])
+            state.desc = NULL
+            state.M = 1
+            e = <ghmm_c_emission*> safe_malloc(sizeof(ghmm_c_emission))
+            e.type = 0  # normal
+            e.dimension = 1
+            e.mean.val = mu
+            e.variance.val = sigma
+            # fixing of emissions is deactivated by default            
+            e.fixed = 0
+            e.sigmacd = NULL
+            e.sigmainv = NULL
+            state.e = e
+            state.c = to_double_array([1.0])
+
+        # Set states
+        self.m.s = states
+
+        self.initialized = True
+
+    def __dealloc__(self):
+        return
+        if self.initialized:
+            ghmm_cmodel_free(&self.m)
+
+    def __repr__(self):
+        """
+        Return string representation of this Continuous HMM.
+
+        OUTPUT:
+            string
+
+        EXAMPLES:
+            sage: m = hmm.NormalHiddenMarkovModel([[0.0,1.0,0],[0.5,0.0,0.5],[0.3,0.3,0.4]], [(0.0,1.0), (-1.0,0.5), (1.0,0.2)], [1,0,0])
+            sage: a.__repr__()
+            "Discrete Hidden Markov Model (2 states, 2 outputs)\nInitial probabilities: [0.5, 0.5]\nTransition matrix:\n[0.1 0.9]\n[0.1 0.9]\nEmission matrix:\n[0.9 0.1]\n[0.1 0.9]\nEmission symbols: [3/4, 'abc']"
+        """
+        s = "Normal Hidden Markov Model%s (%s states, %s outputs)"%(
+            ' ' + self.m.name if self.m.name else '',
+            self.m.N, self.m.M)
+        #s += '\nInitial probabilities: %s'%self.initial_probabilities()
+        #s += '\nTransition matrix:\n%s'%self.transition_matrix()
+        #s += '\nEmission matrix:\n%s'%self.emission_matrix()
+        return s
+

/dev/null

@@ +1 @@
+from sage.matrix.matrix_real_double_dense cimport Matrix_real_double_dense
+
+cdef extern from "ghmm/ghmm.h":
+    cdef int GHMM_kSilentStates
+    cdef int GHMM_kHigherOrderEmissions
+    cdef int GHMM_kDiscreteHMM
+
+cdef extern from "ghmm/model.h":
+    ctypedef struct ghmm_dstate:
+        # Initial probability
+        double pi
+        # Output probability
+        double *b
+
+        # ID's of the following states
+        int *out_id
+        # ID's of the previous states
+        int *in_id
+
+        # Transition probabilities to successor states.
+        double *out_a
+        # Transition probabilities from predecessor states.
+        double *in_a
+
+        # Number of successor states
+        int out_states
+        # Number of precursor states
+        int in_states
+        # if fix == 1 then output probabilities b stay fixed during the training
+        int fix
+        # Contains a description of the state (null terminated utf-8)
+        char * desc
+        # x coordinate position for graph representation plotting
+        int xPosition
+        # y coordinate position for graph representation plotting
+        int yPosition
+
+    ctypedef struct ghmm_dbackground:
+        pass
+
+    ctypedef struct ghmm_alphabet:
+        pass
+    
+    # Discrete HMM's.
+    ctypedef struct ghmm_dmodel:
+        # Number of states
+        int N   
+        # Number of outputs
+        int M   
+        # Vector of states
+        ghmm_dstate *s  
+        # Contains bit flags for varios model extensions such as
+        # kSilentStates, kTiedEmissions (see ghmm.h for a complete list)
+        int model_type
+        # The a priori probability for the model.
+        # A value of -1 indicates that no prior is defined. 
+        # Note: this is not to be confused with priors on emission distributions.
+        double prior
+        # An arbitrary name for the model (null terminated utf-8)
+        char * name
+
+        # Flag variables for each state indicating whether it is emitting or not. 
+        # Note: silent != NULL iff (model_type & kSilentStates) == 1 
+        int *silent
+        
+        # Int variable for the maximum level of higher order emissions.
+        int maxorder
+        
+        # saves the history of emissions as int, 
+        # the nth-last emission is (emission_history * |alphabet|^n+1) % |alphabet|
+        int emission_history
+
+        # Flag variables for each state indicating whether the states emissions
+        # are tied to another state. Groups of tied states are represented
+        # by their tie group leader (the lowest numbered member of the group).
+        # tied_to[s] == kUntied  : s is not a tied state
+        # tied_to[s] == s        : s is a tie group leader
+        # tied_to[t] == s        : t is tied to state s (t>s)
+        # Note: tied_to != NULL iff (model_type & kTiedEmissions) != 0  */
+        int *tied_to
+
+        # Note: State store order information of the emissions.
+        # Classical HMMS have emission order 0; that is, the emission probability
+        # is conditioned only on the state emitting the symbol.
+        # For higher order emissions, the emission are conditioned on the state s
+        # as well as the previous emission_order observed symbols.
+        # The emissions are stored in the state's usual double* b. 
+        # Note: order != NULL iff (model_type & kHigherOrderEmissions) != 0  */
+        int * order
+        
+        # ghmm_dbackground is a pointer to a
+        # ghmm_dbackground structure, which holds (essentially) an
+        # array of background distributions (which are just vectors of floating
+        # point numbers like state.b).
+        # For each state the array background_id indicates which of the background
+        # distributions to use in parameter estimation. A value of kNoBackgroundDistribution
+        # indicates that none should be used.
+        # Note: background_id != NULL iff (model_type & kHasBackgroundDistributions) != 0 
+        int *background_id
+        ghmm_dbackground *bp
+
+        # Topological ordering of silent states 
+        #  Condition: topo_order != NULL iff (model_type & kSilentStates) != 0 
+        int *topo_order
+        int topo_order_length
+
+        # pow_lookup is a array of precomputed powers
+        # It contains in the i-th entry M (alphabet size) to the power of i
+        # The last entry is maxorder+1.
+        int *pow_lookup
+
+        # Store for each state a class label. Limits the possibly state sequence
+        # Note: label != NULL iff (model_type & kLabeledStates) != 0 
+        int* label
+        ghmm_alphabet* label_alphabet
+        ghmm_alphabet* alphabet
+
+
+    # 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)
+    int ghmm_dseq_free (ghmm_dseq ** sq)
+    ghmm_dseq *ghmm_dseq_calloc (long seq_number)
+
+    int ghmm_dmodel_normalize(ghmm_dmodel *m)
+    int ghmm_dmodel_free(ghmm_dmodel **m)
+
+    # Problem 1: Likelihood
+    int ghmm_dmodel_logp(ghmm_dmodel *m, int *O, int len, double *log_p)
+
+    # Problem 2: Decoding
+    int *ghmm_dmodel_viterbi (ghmm_dmodel *m, int *O, int len, int *pathlen, double *log_p)
+
+    # Problem 3: Learning
+    int ghmm_dmodel_baum_welch (ghmm_dmodel *m, ghmm_dseq *sq)
+    
+
+cdef class HiddenMarkovModel:
+    cdef Matrix_real_double_dense A, B
+    cdef list pi
+
+cdef class DiscreteHiddenMarkovModel(HiddenMarkovModel):
+    cdef ghmm_dmodel* m
+    cdef bint initialized
+    cdef object _emission_symbols, _emission_symbols_dict

a/sage/stats/hmm/hmm.pyx

@@ -5 +5 @@
-
-EXAMPLES:
-
+TODO:
+   * make models pickleable (i.e., all parameters should be obtainable
+     using functions to make this easy).
+   * continuous hmm's
-"""
-
-import math
@@ -13 +17 @@
-from sage.rings.all  import RDF
-from sage.misc.randstate import random
-
-from sage.matrix.matrix_real_double_dense cimport Matrix_real_double_dense
-
-include "../../ext/stdsage.pxi"
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
-
-cdef class HiddenMarkovModel:
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
-        if not is_Matrix(A):
-            A = matrix(RDF, len(A), len(A[0]), A)
-        if not is_Matrix(B):
@@ -214 +35 @@
-            A = A.change_ring(RDF)
-        if B.base_ring() != RDF:
-            B = B.change_ring(RDF)
-
+
+
+
-            raise ValueError, "length of pi must equal number of rows of A"
-
-        cdef Py_ssize_t i, j, k
-
-        self.A = A
-        self.B = B
+
+
+cdef class DiscreteHiddenMarkovModel(HiddenMarkovModel):
+    
+    def __init__(self, A, B, pi, emission_symbols=None, name=None):
+        """
+        INPUTS:
+            A  -- square matrix of doubles; the state change probabilities
+            B  -- matrix of doubles; emission probabilities
+            pi -- list of floats; probabilities for each initial state
+            emission_state -- list of B.ncols() symbols (just used for printing)
+            name -- (optional) name of the model
+
+        EXAMPLES:
+        We create a discrete HMM with 2 internal states on an alphabet of
+        size 2.
+        
+            sage: a = hmm.DiscreteHiddenMarkovModel([[0.2,0.8],[0.5,0.5]], [[1,0],[0,1]], [0,1])
+            
+        """
+        self.initialized = False
+        HiddenMarkovModel.__init__(self, A, B, pi)
+
+        cdef Py_ssize_t i, j, k
-        self.set_emission_symbols(emission_symbols)
-
-        self.m = <ghmm_dmodel*> safe_malloc(sizeof(ghmm_dmodel))
@@ -234 +79 @@
-        self.m.label_alphabet = NULL; self.m.alphabet = NULL
-
-        # Set number of states and number of outputs
-
+self.
-        self.m.M = len(self._emission_symbols)
-        # Set the model type to discrete
-        self.m.model_type = GHMM_kDiscreteHMM
-
-        # Set that no a prior model probabilities are set.
-        self.m.prior = -1
+
-        # Assign model identifier if specified
-        if name is not None: 
-            name = str(name)
@@ -256 +102 @@
-        tmp_order     = []
-        
-        for i in range(self.m.N):
-
+self.
-
-            # Get a reference to the i-th state for convenience of the notation below.
-            state = &(states[i])
@@ -274 +120 @@
-                tmp_order.append(order)
-            else:
-                raise ValueError, "number of columns (= %s) of B must be a power of the number of emission symbols (= %s)"%(
-
+self.
-
-            state.b = to_double_array(v)
-
+self.
-
-            silent_states.append( 1 if sum(v) == 0 else 0)
-
-            # Set out probabilities, i.e., the probabilities that each
-            # symbol will be emitted from this state. 
-
+self.
-            nz = v.nonzero_positions()
-            k = len(nz)
-            state.out_states = k
@@ -294 +140 @@
-                state.out_a[j]  = v[nz[j]]
-
-            # Set "in" probabilities
-
+self.
-            nz = v.nonzero_positions()
-            k = len(nz)
-            state.in_states = k
@@ -307 +153 @@
-            state.fix = 0
-                
-        self.m.s = states
-; return
+
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-        
-    def __dealloc__(self):
-        if self.initialized:
@@ -341 +184 @@
-        EXAMPLES:
-            sage: a = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.1,0.9]], [[0.9,0.1],[0.1,0.9]], [0.5,0.5], [3/4, 'abc'])
-            sage: a.__repr__()
-\n
+
-        """
-\n
+
-            ' ' + self.m.name if self.m.name else '',
-            self.m.N, self.m.M)
-        s += '\nInitial probabilities: %s'%self.initial_probabilities()
@@ -439 +282 @@
-        cdef ghmm_dseq *d = ghmm_dmodel_generate_sequences(self.m, seed, length, 1, length)
-        cdef Py_ssize_t i
-        v = [d.seq[0][i] for i in range(length)]
+        ghmm_dseq_free(&d)
-        if self._emission_symbols_dict:
-            w = self._emission_symbols
-            return [w[i] for i in v]
@@ -510 +354 @@
-    ####################################################################
-    def log_likelihood(self, seq):
-        r"""
-Return $\log ( P[seq | model] )$, the log of
-probability of seeing the given sequence given this model,
-using the forward algorithm and assuming independance of the
-
+Likelihood. Return $\log ( P[seq | model] )$,
+log of the probability of seeing the given sequence given
+this model, using the forward algorithm and assuming
+independance of the 
-
-        INPUT:
-            seq -- a list; sequence of observed emissions of the HMM
@@ -566 +410 @@
-    ####################################################################    
-    def viterbi(self, seq):
-        """
-
-
+
+
+
-
-        INPUT:
-            seq -- sequence of emitted symbols
@@ -612 +457 @@
-    # probability of observing O.
-    ####################################################################
-    def baum_welch(self, training_seqs, nsteps=None, log_likelihood_cutoff=None):
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
-
+        INPUT:
+            training_seqs -- a list of lists of emission symbols
+            nsteps -- integer or None (default: None) maximum number
+                      of Baum-Welch steps to take
+            log_likehood_cutoff -- positive float or None (default:
+                      None); the minimal improvement in likelihood
+                      with respect to the last iteration required to
+                      continue. Relative value to log likelihood
+
+        OUTPUT:
+            changes the model in places, or raises a RuntimError
+            exception on error
+
+        EXAMPLES:
+        We make a model that has two states and is equally likely to output
+        0 or 1 in either state and transitions back and forth with equal
+        probability.
+            sage: a = hmm.DiscreteHiddenMarkovModel([[0.5,0.5],[0.5,0.5]], [[0.5,0.5],[0.5,0.5]], [0.5,0.5])
+
+        We give the model some training data this much more likely to
+        be 1 than 0.
+            sage: a.baum_welch([[1,1,1,1,0,1], [1,0,1,1,1,1]])
+
+        Now the model's emission matrix changes since it is much
+        more likely to emit 1 than 0.  
+            sage: a
+            Discrete Hidden Markov Model (2 states, 2 outputs)
+            Initial probabilities: [0.5, 0.5]
+            Transition matrix:
+            [0.5 0.5]
+            [0.5 0.5]
+            Emission matrix:
+            [0.166666666667 0.833333333333]
+            [0.166666666667 0.833333333333]
+
+        Note that 1/6 = 1.666...:
+            sage: 1.0/6
+            0.166666666666667
+
+        TESTS:
+        We test training with non-default string symbols:
+            sage: a = hmm.DiscreteHiddenMarkovModel([[0.5,0.5],[0.5,0.5]], [[0.5,0.5],[0.5,0.5]], [0.5,0.5], ['up','down'])
+            sage: a.baum_welch([['up','up'], ['down','up']])
+            sage: a
+            Discrete Hidden Markov Model (2 states, 2 outputs)
+            Initial probabilities: [0.5, 0.5]
+            Transition matrix:
+            [0.5 0.5]
+            [0.5 0.5]
+            Emission matrix:
+            [0.75 0.25]
+            [0.75 0.25]
+            Emission symbols: ['up', 'down']        
+
+        NOTE: Training for models including silent states is not yet supported.
+
+        REFERENCES:
+            Rabiner, L.R.: "`A Tutorial on Hidden Markov Models and Selected
+            Applications in Speech Recognition"', Proceedings of the IEEE,
+            77, no 2, 1989, pp 257--285. 
+        """
+        if self._emission_symbols_dict:
+            seqs = [[self._emission_symbols_dict[z] for z in x] for x in training_seqs]
+        else:
+            seqs = training_seqs
+            
+        cdef ghmm_dseq* d = malloc_ghmm_dseq(seqs)
-        
+        if ghmm_dmodel_baum_welch(self.m, d):
+            raise RuntimeError, "error running Baum-Welch algorithm"
-        
+        ghmm_dseq_free(&d)
-
-
-
-
+
-##################################################################################
-# Helper Functions
-##################################################################################
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-