Ticket #8547: trac-8547-take2.patch

File trac-8547-take2.patch, 243.4 KB (added by was, 12 years ago)

this completely replaces the previous patch

  • module_list.py 

    # HG changeset patch
# User William Stein <wstein@gmail.com>
# Date 1269812759 25200
# Node ID e964f241c1d92393c7c5bc9b24687e72cdb9a218
# Parent  a0a95c2cf60608c1bdedef3475e8a515ba62fead
trac 8547 -- take 2 -- implement hidden markov models in Cython from scratch

diff --git a/module_list.py b/module_list.py
    a b  
    13341334    ##
    13351335    ################################
    13361336
    1337     Extension('sage.stats.hmm.chmm',
    1338               sources = ['sage/stats/hmm/chmm.pyx'],
    1339               libraries = ['ghmm'],
    1340               include_dirs=numpy_include_dirs),
     1337    Extension('sage.stats.hmm.util',
     1338              sources = ['sage/stats/hmm/util.pyx']),
     1339
     1340    Extension('sage.stats.hmm.distributions',
     1341              sources = ['sage/stats/hmm/distributions.pyx']),
    13411342   
    13421343    Extension('sage.stats.hmm.hmm',
    1343               sources = ['sage/stats/hmm/hmm.pyx'],
    1344               libraries = ['ghmm'],
    1345               include_dirs=numpy_include_dirs),
     1344              sources = ['sage/stats/hmm/hmm.pyx']),
     1345
     1346    Extension('sage.stats.hmm.chmm',
     1347              sources = ['sage/stats/hmm/chmm.pyx']),
     1348   
     1349    Extension('sage.stats.intlist',
     1350              sources = ['sage/stats/intlist.pyx']),
    13461351
    13471352    ################################
    13481353    ##

  • sage/finance/time_series.pxd 

    diff --git a/sage/finance/time_series.pxd b/sage/finance/time_series.pxd
    a b  
    11cdef class TimeSeries:
    22    cdef double* _values
    33    cdef Py_ssize_t _length
     4    cpdef rescale(self, double s)
     5    cpdef double sum(self)

  • sage/finance/time_series.pyx 

    diff --git a/sage/finance/time_series.pyx b/sage/finance/time_series.pyx
    a b  
    6868        """
    6969        self._values = NULL
    7070       
    71     def __init__(self, values):
     71    def __init__(self, values, bint initialize=True):
    7272        """
    7373        Initialize new time series.
    7474       
    7575        INPUT:
    76             values -- integer (number of values) or an iterable of floats
     76
     77            - values -- integer (number of values) or an iterable of
     78              floats
     79           
     80            - initialize -- bool (default: True); if False, do not
     81              bother to zero out the entries of the new time series.
     82              For large series that you are going to just fill in,
     83              this can be way faster.
    7784
    7885        EXAMPLES:
    79         This implicitly calls init.
     86
     87        This implicitly calls init.::
     88       
    8089            sage: finance.TimeSeries([pi, 3, 18.2])
    8190            [3.1416, 3.0000, 18.2000]
    8291
    83         Conversion from a numpy 1d array, which is very fast.
     92        Conversion from a numpy 1d array, which is very fast.::
     93       
    8494            sage: v = finance.TimeSeries([1..5])
    8595            sage: w = v.numpy()
    8696            sage: finance.TimeSeries(w)
    8797            [1.0000, 2.0000, 3.0000, 4.0000, 5.0000]
    8898
    89         Conversion from an n-dimensional numpy array also works:
     99        Conversion from an n-dimensional numpy array also works::
     100       
    90101            sage: import numpy
    91102            sage: v = numpy.array([[1,2], [3,4]], dtype=float); v
    92103            array([[ 1.,  2.],
    … …  
    98109            sage: u = numpy.array([[1,2],[3,4]])
    99110            sage: finance.TimeSeries(u)
    100111            [1.0000, 2.0000, 3.0000, 4.0000]
     112
     113        For speed purposes we don't initialize (so value is garbage)::
     114
     115            sage: t = finance.TimeSeries(10, initialize=False)       
    101116        """
    102117        cdef Vector_real_double_dense z
    103118        cdef cnumpy.ndarray np
    104119        cdef double *np_data
    105120        cdef unsigned int j
    106121        if isinstance(values, (int, long, Integer)):
    107             values = [0.0]*values
     122            self._length = values
     123            values = None
    108124        elif PY_TYPE_CHECK(values, Vector_real_double_dense) or PY_TYPE_CHECK(values, cnumpy.ndarray):
    109125            if PY_TYPE_CHECK(values, Vector_real_double_dense):
    110126                np  = values._vector_numpy
    … …  
    130146            return
    131147        else:
    132148            values = [float(x) for x in values]
    133         self._length = len(values)
     149            self._length = len(values)
     150           
    134151        self._values = <double*> sage_malloc(sizeof(double) * self._length)
    135152        if self._values == NULL:
    136153            raise MemoryError
     154        if not initialize: return
    137155        cdef Py_ssize_t i
    138         for i from 0 <= i < self._length:
    139             self._values[i] = values[i]
     156        if values is not None:
     157            for i from 0 <= i < self._length:
     158                self._values[i] = values[i]
     159        else:
     160            for i from 0 <= i < self._length:
     161                self._values[i] = 0
    140162
    141163    def __reduce__(self):
    142164        """
    … …  
    265287        """
    266288        if prec is None: prec = digits
    267289        format = '%.' + str(prec) + 'f'
    268         v = self.list()
    269         if len(v) > max_print:
    270             v0 = v[:max_print//2]
    271             v1 = v[-max_print//2:]
     290        if len(self) > max_print:
     291            v0 = self[:max_print//2]
     292            v1 = self[-max_print//2:]
    272293            return '[' + ', '.join([format%x for x in v0]) + ' ... ' + \
    273294                         ', '.join([format%x for x in v1]) + ']'
    274295        else:
    275             return '[' + ', '.join([format%x for x in v]) + ']'
     296            return '[' + ', '.join([format%x for x in self]) + ']'
    276297
    277298    def __len__(self):
    278299        """
    … …  
    363384                memcpy(t._values, self._values + start, sizeof(double)*t._length)
    364385            return t
    365386        else:
    366             if i < 0:
    367                 i += self._length
    368                 if i < 0:
     387            j = i
     388            if j < 0:
     389                j += self._length
     390                if j < 0:
    369391                    raise IndexError, "TimeSeries index out of range"
    370             elif i >= self._length:
     392            elif j >= self._length:
    371393                raise IndexError, "TimeSeries index out of range"
    372             return self._values[i]
     394            return self._values[j]
    373395
    374396    def __setitem__(self, Py_ssize_t i, double x):
    375397        """
    … …  
    643665            sage: v.list()
    644666            [1.0, -4.0, 3.0, -2.5, -4.0, 3.0]
    645667        """
    646         v = [0.0]*self._length
    647668        cdef Py_ssize_t i
    648         for i from 0 <= i < self._length:
    649             v[i] = self._values[i]
    650         return v
     669        return [self._values[i] for i in range(self._length)]
    651670
    652671    def log(self):
    653672        """
    … …  
    666685            [2.7183, 0.0183, 20.0855, 0.0821, 0.0183, 20.0855]
    667686            sage: v.exp().log()
    668687            [1.0000, -4.0000, 3.0000, -2.5000, -4.0000, 3.0000]
     688
     689        Log of 0 gives -inf:
     690
     691            sage: finance.TimeSeries([1,0,3]).log()
     692            [0.0000, -inf, 1.0986]
    669693        """
    670694        cdef Py_ssize_t i
    671         for i from 0 <= i < self._length:
    672             if self._values[i] <= 0:
    673                 raise ValueError, "every entry of self must be positive."
    674 
    675695        cdef TimeSeries t = new_time_series(self._length)
    676696        for i from 0 <= i < self._length:
    677697            t._values[i] = log(self._values[i])
    … …  
    801821            t._values[i] = self._values[i*k]
    802822        return t
    803823
     824    cpdef rescale(self, double s):
     825        """
     826        Change self by multiplying every value in the series by s.
     827
     828        INPUT:
     829            s -- float
     830
     831        EXAMPLES:
     832            sage: v = finance.TimeSeries([5,4,1.3,2,8,10,3,-5]); v
     833            [5.0000, 4.0000, 1.3000, 2.0000, 8.0000, 10.0000, 3.0000, -5.0000]
     834            sage: v.rescale(0.5)
     835            sage: v
     836            [2.5000, 2.0000, 0.6500, 1.0000, 4.0000, 5.0000, 1.5000, -2.5000]
     837        """
     838        for i from 0 <= i < self._length:
     839            self._values[i] = self._values[i] * s
     840
    804841    def scale(self, double s):
    805842        """
    806843        Return new time series obtained by multiplying every value in the series by s.
    … …  
    10671104            t._values[i] = s
    10681105        return t
    10691106
    1070     def sum(self):
     1107    cpdef double sum(self):
    10711108        """
    10721109        Return the sum of all the entries of self.  If self has
    10731110        length 0, returns 0.
    … …  
    15551592        """
    15561593        if self._length == 0:
    15571594            raise ValueError, "min() arg is an empty sequence"
    1558             return 0.0, -1
    15591595        cdef Py_ssize_t i, j
    15601596        cdef double s = self._values[0]
    15611597        j = 0
    … …  
    15891625        """
    15901626        if self._length == 0:
    15911627            raise ValueError, "max() arg is an empty sequence"
    1592         cdef Py_ssize_t i, j
     1628        cdef Py_ssize_t i, j = 0
    15931629        cdef double s = self._values[0]
    15941630        for i from 1 <= i < self._length:
    15951631            if self._values[i] > s:
    … …  
    17431779            bins -- positive integer (default: 50)
    17441780            normalize -- (default: True) whether to normalize so the total
    17451781                         area in the bars of the histogram is 1.
    1746             **kwds -- passed to the bar_chart function
    17471782        OUTPUT:
    17481783            a histogram plot
    17491784
    … …  
    17511786            sage: v = finance.TimeSeries([1..50])
    17521787            sage: v.plot_histogram(bins=10)
    17531788        """
    1754         from sage.plot.all import bar_chart, polygon
     1789        from sage.plot.all import polygon
    17551790        counts, intervals = self.histogram(bins, normalize=normalize)
    17561791        s = 0
    17571792        for i, (x0,x1) in enumerate(intervals):
    … …  
    22392274
    22402275def unpickle_time_series_v1(v, Py_ssize_t n):
    22412276    """
    2242     Version 0 unpickle method.
     2277    Version 1 unpickle method.
    22432278
    22442279    INPUT:
    22452280        v -- a raw char buffer

  • sage/matrix/matrix_complex_double_dense.pxd 

    diff --git a/sage/matrix/matrix_complex_double_dense.pxd b/sage/matrix/matrix_complex_double_dense.pxd
    a b  
     1import matrix_double_dense
     2
    13cimport matrix_double_dense
    24
    35cdef class Matrix_complex_double_dense(matrix_double_dense.Matrix_double_dense):

  • sage/matrix/matrix_complex_double_dense.pyx 

    diff --git a/sage/matrix/matrix_complex_double_dense.pyx b/sage/matrix/matrix_complex_double_dense.pyx
    a b  
    3636#  The full text of the GPL is available at:
    3737#                  http://www.gnu.org/licenses/
    3838##############################################################################
     39
     40import matrix_double_dense
     41
    3942from sage.rings.complex_double import CDF
    4043
    4144cimport numpy as cnumpy

  • sage/matrix/matrix_double_dense.pyx 

    diff --git a/sage/matrix/matrix_double_dense.pyx b/sage/matrix/matrix_double_dense.pyx
    a b  
    5252cnumpy.import_array()
    5353
    5454
    55 from matrix_complex_double_dense import Matrix_complex_double_dense
    56 
    5755cdef class Matrix_double_dense(matrix_dense.Matrix_dense):
    5856    """
    5957    Base class for matrices over the Real Double Field and the Complex

  • sage/stats/all.py 

    diff --git a/sage/stats/all.py b/sage/stats/all.py
    a b  
    33
    44import hmm.all as hmm
    55
     6from sage.finance.time_series import TimeSeries
     7from intlist import IntList
     8

  • sage/stats/hmm/__init__.py 

    diff --git a/sage/stats/hmm/__init__.py b/sage/stats/hmm/__init__.py
    a b  
    1 # Hidden Markov Models
     1# hmm

  • sage/stats/hmm/all.py 

    diff --git a/sage/stats/hmm/all.py b/sage/stats/hmm/all.py
    a b  
    11#############################################################################
    2 #       Copyright (C) 2008 William Stein <wstein@gmail.com>
    3 #  Distributed under the terms of the GNU General Public License (GPL),
    4 #  version 2 or any later version at your option.
     2#       Copyright (C) 2010 William Stein <wstein@gmail.com>
     3#  Distributed under the terms of the GNU General Public License (GPL), v2+.
    54#  The full text of the GPL is available at:
    65#                  http://www.gnu.org/licenses/
    76#############################################################################
    87
    98from hmm  import DiscreteHiddenMarkovModel
    10 from chmm import GaussianHiddenMarkovModel, GaussianMixtureHiddenMarkovModel
     9
     10from chmm import (GaussianHiddenMarkovModel,
     11                  GaussianMixtureHiddenMarkovModel)
     12
     13from distributions import GaussianMixtureDistribution
     14

  • deleted file sage/stats/hmm/chmm.pxd 

    diff --git a/sage/stats/hmm/chmm.pxd b/sage/stats/hmm/chmm.pxd
deleted file mode 100644
    + -  
    1 #############################################################################
    2 #       Copyright (C) 2008 William Stein <wstein@gmail.com>
    3 #  Distributed under the terms of the GNU General Public License (GPL),
    4 #  version 2 or any later version at your option.
    5 #  The full text of the GPL is available at:
    6 #                  http://www.gnu.org/licenses/
    7 #############################################################################
    8 
    9 from hmm cimport HiddenMarkovModel
    10 
    11 cdef extern from "ghmm/ghmm.h":
    12     cdef int GHMM_kContinuousHMM
    13 
    14 cdef extern from "ghmm/smodel.h":
    15     #############################################################
    16     # Continuous emission density types
    17     #############################################################
    18     cdef enum ghmm_density_t:
    19         normal,
    20         normal_right,   # right tail
    21         normal_approx,
    22         normal_left,    # left tail
    23         uniform,
    24         binormal,
    25         multinormal,
    26         density_number
    27 
    28     #############################################################
    29     # Continuous emission
    30     #############################################################
    31     # Mean and variance types
    32     cdef union mean_t:
    33         double val
    34         double *vec
    35        
    36     cdef union variance_t:
    37         double val
    38         double *mat
    39 
    40     cdef struct ghmm_c_emission:
    41         # Flag for density function for each component of the mixture
    42         #   0: normal density
    43         #   1: truncated normal (right side) density
    44         #   2: approximated normal density
    45         #   3: truncated normal (left side)
    46         #   4: uniform distribution
    47         #   6: multivariate normal
    48         int type
    49         # dimension > 1 for multivariate normals
    50         int dimension
    51         # mean for output functions (pointer to mean vector for multivariate)
    52         mean_t mean
    53         variance_t variance
    54         # pointer to inverse of covariance matrix if multivariate normal;
    55         # else NULL
    56         double *sigmainv
    57         # determinant of covariance matrix for multivariate normal
    58         double det
    59         # Cholesky decomposition of covariance matrix A,
    60         #   if A = G*G' sigmacd only holds G
    61         double *sigmacd
    62         # minimum of uniform distribution or left boundary for
    63         # right-tail gaussians
    64         double min
    65         # maximum of uniform distribution or right boundary for
    66         # left-tail gaussians
    67         double max
    68         # if fixed != 0 the parameters of the density are fixed
    69         int fixed
    70        
    71     #############################################################
    72     # Continous Emission States
    73     #############################################################
    74     ctypedef struct ghmm_cstate:
    75         # Number of output densities per state
    76         int M
    77         # initial prob.
    78         double pi
    79         # IDs of successor states
    80         int *out_id
    81         # IDs of predecessor states
    82         int *in_id
    83         # transition probs to successor states. It is a
    84         # matrix in case of mult. transition matrices (cos > 1)
    85         double **out_a
    86         # transition probs from predecessor states. It is a
    87         # matrix in case of mult. transition matrices (cos > 1)
    88         double **in_a
    89         # number of  successor states
    90         int out_states
    91         # number of  predecessor states
    92         int in_states
    93         # weight vector for output function components
    94         double *c
    95         # flag for fixation of parameter. If fix = 1 do not change parameters of
    96         # output functions, and if fix = 0 do normal training. Default is 0.
    97         int fix
    98         # vector of ghmm_c_emission
    99         # (type and parameters of output function components)
    100         ghmm_c_emission *e
    101         # contains a description of the state (null terminated utf-8)
    102         char *desc
    103         # x coordinate position for graph representation plotting
    104         int xPosition
    105         # y coordinate position for graph representation plotting
    106         int yPosition
    107 
    108 
    109     double **ighmm_cmatrix_alloc (int nrows, int ncols)
    110 
    111     #############################################################
    112     # Continous Hidden Markov Model
    113     #############################################################
    114     ctypedef struct ghmm_cmodel
    115     ctypedef struct ghmm_cmodel_class_change_context
    116        
    117     ctypedef struct ghmm_cmodel:
    118         # Number of states
    119         int N
    120         # Maximun number of components in the states
    121         int M
    122         # Number of dimensions of the emission components.
    123         # NOTE: All emissions must have the same number of dimensions.
    124         int dim
    125         # The ghmm_cmodel includes two continuous models
    126         #   cos=1: one transition matrix
    127         #   cos>1: (integer) extension for models with several matrices.
    128         int cos
    129         # Prior for a priori probability of the model.
    130         # -1 means no prior specified (all models have equal prob. a priori.)
    131         double prior
    132         # Contains a arbitrary name for the model (null terminated utf-8)
    133         char * name
    134         # Contains bit flags for various model extensions such as
    135         # kSilentStates (see ghmm.h for a complete list).
    136         int model_type
    137         # All states of the model. Transition probabilities are part of the states.
    138         ghmm_cstate *s
    139         # pointer to a ghmm_cmodel_class_change_context struct
    140         # necessary for multiple transition classes
    141         ghmm_cmodel_class_change_context *class_change
    142 
    143     int ghmm_cmodel_class_change_alloc(ghmm_cmodel * smo)
    144 
    145     #############################################################
    146     # Continous Sequence of states
    147     #
    148     # Sequence structure for double sequences.  Contains an array of
    149     # sequences and corresponding data like sequence labels, sequence
    150     # weights, etc. Sequences may have different length.
    151     # multi-dimension sequences are linearized
    152     #############################################################
    153     ctypedef struct ghmm_cseq:
    154         # sequence array. sequence[i][j] = j-th symbol of i-th seq.
    155         # sequence[i][D * j] = first dimension of j-th observation of i-th sequence
    156         double **seq
    157         # array of sequence lengths
    158         int *seq_len
    159         # array of sequence IDs
    160         double *seq_id
    161         # positive! sequence weights.  default is 1 = no weight
    162         double *seq_w
    163         # total number of sequences
    164         long seq_number
    165         # reserved space for sequences is always >= seq_number
    166         long capacity
    167         # sum of sequence weights
    168         double total_w
    169         # total number of dimensions
    170         int dim
    171         # flags (internal)
    172         unsigned int flags
    173         # obsolete
    174         int* seq_label
    175 
    176     int ghmm_cseq_free (ghmm_cseq ** sq)
    177 
    178     #############################################################
    179     # Continous Baum-Welch algorithm context
    180     #############################################################
    181     ctypedef struct ghmm_cmodel_baum_welch_context:
    182         # pointer of continuous model
    183         ghmm_cmodel *smo
    184         # sequence pointer
    185         ghmm_cseq *sqd
    186         # calculated log likelihood
    187         double *logp
    188         # leave reestimation loop if diff. between successive logp values
    189         # is smaller than eps
    190         double eps
    191         # max. no of iterations
    192         int max_iter
    193 
    194 
    195     int ghmm_cmodel_free (ghmm_cmodel ** smo)
    196 
    197     # Normalizes initial and transition probabilities and mixture weights.
    198     # 0 on success / -1 on error
    199     int ghmm_cmodel_normalize(ghmm_cmodel *smo)   
    200 
    201    
    202     # Produces sequences computed using a given model. All memory that
    203     # is needed for the sequences is allocated inside the function. It
    204     # is possible to define the length of the sequences global
    205     # (global_len > 0) or it can be set inside the function, when a
    206     # final state in the model is reached (a state with no output). If
    207     # the model has no final state, the sequences will have length
    208     # MAX_SEQ_LEN.
    209     #
    210     #  INPUT:
    211     #      smo       -- model
    212     #      seed      -- initial parameter for the random number generator;
    213     #                   if seed == 0, don't reinitialize random number generator
    214     #      global_en -- length of sequence.  If 0, sequence ended automatically
    215     #                   when the final state is reached
    216     #      seq_number-- number of sequences
    217     #      Tmax      -- maximal sequence length (or -1)
    218     #
    219     #  OUTPUT:
    220     #      pointer to an array of sequences
    221     #
    222 
    223     ghmm_cseq *ghmm_cmodel_generate_sequences (ghmm_cmodel * smo, int seed,
    224                                              int global_len, long seq_number,
    225                                              int Tmax)
    226 
    227     # Computes sum over all sequence of seq_w * log( P ( O|lambda )).
    228     # If a sequence can't be generated by smo error cost of seq_w *
    229     # PENALTY_LOGP are imposed.
    230     #   OUTPUT: number of evaluated sequences; -1: error
    231     int ghmm_cmodel_likelihood (ghmm_cmodel * smo, ghmm_cseq * sqd, double *log_p)
    232 
    233    
    234     # Viterbi algorithm: calculation of the viterbi path (best possible
    235     # state sequence for a given sequence for a given model (smo)). Also
    236     # calculates log(p) according to this path, the matrices in the local_store
    237     # struct are allocated using stat_matrix_d_alloc.
    238     int *ghmm_cmodel_viterbi (ghmm_cmodel * smo, double *o, int n, double *log_p)
    239 
    240 cdef extern from "ghmm/sreestimate.h":
    241     ctypedef struct ghmm_cmodel_class_change_context:
    242         # Names of class change module/function (for python callback)
    243         char *python_module
    244         char *python_function
    245         # index of current sequence
    246         int k
    247         # pointer to class function
    248         int (*get_class) (ghmm_cmodel*, double*, int, int)
    249         # space for any data necessary for class switch, USER is RESPONSIBLE
    250         void *user_data
    251 
    252     # Baum-Welch Algorithm for continuous HMMs
    253     # Training of model parameter with multiple double sequences (incl. scaling).
    254     # New parameters set directly in hmm (no storage of previous values!). Matrices
    255     # are allocated with stat_matrix_d_alloc.
    256     int ghmm_cmodel_baum_welch (ghmm_cmodel_baum_welch_context * cs)
    257 
    258 
    259 cdef class ContinuousHiddenMarkovModel(HiddenMarkovModel):
    260     cdef ghmm_cmodel* m
    261     cdef bint initialized
    262 
    263 cdef class GaussianHiddenMarkovModel(ContinuousHiddenMarkovModel):
    264     pass
    265 
    266 

  • sage/stats/hmm/chmm.pyx 

    diff --git a/sage/stats/hmm/chmm.pyx b/sage/stats/hmm/chmm.pyx
    a b  
    1 """nodoctest
    2 Continuous Hidden Markov Models
     1"""
     2Continuous Emission Hidden Markov Models
    33
    4 AUTHORS:
    5     -- William Stein (2008)
    6     -- The authors of GHMM http://ghmm.sourceforge.net/
     4AUTHOR:
     5
     6   - William Stein, 2010-03
    77"""
    88
    99#############################################################################
    10 #       Copyright (C) 2008 William Stein <wstein@gmail.com>
    11 #  Distributed under the terms of the GNU General Public License (GPL),
    12 #  version 2 or any later version at your option.
     10#       Copyright (C) 2010 William Stein <wstein@gmail.com>
     11#  Distributed under the terms of the GNU General Public License (GPL)
    1312#  The full text of the GPL is available at:
    1413#                  http://www.gnu.org/licenses/
    1514#############################################################################
    1615
    1716include "../../ext/stdsage.pxi"
     17include "../../ext/cdefs.pxi"
     18include "../../ext/interrupt.pxi"
    1819
    19 include "misc.pxi"
     20cdef extern from "math.h":
     21    double log(double)
     22    double sqrt(double)
     23    double exp(double)
     24    int isnormal(double)
     25    int isfinite(double)
    2026
    21 from sage.misc.randstate import random
    22 from sage.misc.flatten   import flatten
     27import math
     28
     29from sage.misc.flatten  import flatten
     30from sage.matrix.all import is_Matrix
     31
     32cdef double sqrt2pi = sqrt(2*math.pi)
     33cdef double NaN = 1.0/0.0
    2334
    2435from sage.finance.time_series cimport TimeSeries
     36from sage.stats.intlist cimport IntList
    2537
    26 cdef extern from "math.h":
    27     double sqrt(double)
    28    
    29 cdef class ContinuousHiddenMarkovModel(HiddenMarkovModel):
     38from hmm cimport HiddenMarkovModel
     39from util cimport HMM_Util
     40from distributions cimport GaussianMixtureDistribution
     41
     42cdef HMM_Util util = HMM_Util()
     43
     44from sage.misc.randstate cimport current_randstate, randstate
     45
     46
     47# DELETE THIS FUNCTION WHEN MOVE Gaussian stuff to distributions.pyx!!!
     48cdef double random_normal(double mean, double std, randstate rstate):
     49    # Ported from http://users.tkk.fi/~nbeijar/soft/terrain/source_o2/boxmuller.c
     50    # This the box muller algorithm.
     51    # Client code can get the current random state from:
     52    #         cdef randstate rstate = current_randstate()
     53
     54    cdef double x1, x2, w, y1, y2
     55    while True:
     56        x1 = 2*rstate.c_rand_double() - 1
     57        x2 = 2*rstate.c_rand_double() - 1
     58        w = x1*x1 + x2*x2
     59        if w < 1: break
     60    w = sqrt( (-2*log(w))/w )
     61    y1 = x1 * w
     62    return mean + y1*std
     63
     64cdef class GaussianHiddenMarkovModel(HiddenMarkovModel):
    3065    """
    31     Abstract base class for continuous hidden Markov models.
     66    GaussianHiddenMarkovModel(A, B, pi)
    3267
    33        
     68    Gaussian emissions Hidden Markov Model.
     69
    3470    INPUT:
    35         A -- square matrix (or list of lists)
    36         B -- list or matrix with numbers of rows equal to number of rows of A;
    37              each row defines the emissions
    38         pi -- list
    39         name -- (default: None); a string
     71        - A  -- matrix; the N x N transition matrix
     72        - B -- list of pairs (mu,sigma) that define the distributions
     73        - pi -- initial state probabilities
    4074
    41     EXAMPLES:
    42         sage: sage.stats.hmm.chmm.ContinuousHiddenMarkovModel([[1.0]], [(-0.0,10.0)], [1], "model")
    43         <sage.stats.hmm.chmm.ContinuousHiddenMarkovModel object at ...>
     75    EXAMPLES::
     76
     77    We illustrate the primary functions with an example 2-state Gaussian HMM::
     78
     79        sage: m = hmm.GaussianHiddenMarkovModel([[.1,.9],[.5,.5]], [(1,1), (-1,1)], [.5,.5]); m
     80        Gaussian Hidden Markov Model with 2 States
     81        Transition matrix:
     82        [0.1 0.9]
     83        [0.5 0.5]
     84        Emission parameters:
     85        [(1.0, 1.0), (-1.0, 1.0)]
     86        Initial probabilities: [0.5000, 0.5000]
     87
     88    We query the defining transition matrix, emission parameters, and
     89    initial state probabilities::
     90
     91        sage: m.transition_matrix()
     92        [0.1 0.9]
     93        [0.5 0.5]
     94        sage: m.emission_parameters()
     95        [(1.0, 1.0), (-1.0, 1.0)]
     96        sage: m.initial_probabilities()
     97        [0.5000, 0.5000]
     98
     99    We obtain a sample sequence with 10 entries in it, and compute the
     100    logarithm of the probability of obtaining his sequence, given the
     101    model::
     102
     103        sage: obs = m.sample(10); obs
     104        [-1.6835, 0.0635, -2.1688, 0.3043, -0.3188, -0.7835, 1.0398, -1.3558, 1.0882, 0.4050]
     105        sage: m.log_likelihood(obs)
     106        -15.2262338077988...
     107
     108    We compute the Viterbi path, and probability that the given path
     109    of states produced obs::
     110
     111        sage: m.viterbi(obs)
     112        ([1, 0, 1, 0, 1, 1, 0, 1, 0, 1], -16.677382701707881)
     113
     114    We use the Baum-Welch iterative algorithm to find another model
     115    for which our observation sequence is more likely::
     116
     117        sage: m.baum_welch(obs)
     118        (-10.610333495739708, 14)
     119        sage: m.log_likelihood(obs)
     120        -10.610333495739708
     121
     122    Notice that running Baum-Welch changed our model::
     123
     124        sage: m
     125        Gaussian Hidden Markov Model with 2 States                                                   
     126        Transition matrix:                                                                           
     127        [   0.415498136619    0.584501863381]                                                         
     128        [   0.999999317425 6.82574625899e-07]                                                         
     129        Emission parameters:                                                                         
     130        [(0.417888242712, 0.517310966436), (-1.50252086313, 0.508551283606)]                         
     131        Initial probabilities: [0.0000, 1.0000]   
    44132    """
    45     def __init__(self, A, B, pi, name):
     133    cdef TimeSeries B, prob
     134    cdef int n_out
     135
     136    def __init__(self, A, B, pi, bint normalize=True):
    46137        """
    47         Constructor for continuous Hidden Markov model abstract base class.
     138        Create a Gaussian emissions HMM with transition probability
     139        matrix A, normal emissions given by B, and initial state
     140        probability distribution pi.
    48141
    49         EXAMPLES:
    50         This class is an abstract base class, so shouldn't ever be constructed by users.
    51             sage: sage.stats.hmm.chmm.ContinuousHiddenMarkovModel([[1.0]], [(0.0,1.0)], [1], None)
    52             <sage.stats.hmm.chmm.ContinuousHiddenMarkovModel object at ...>
     142        INPUT::
     143
     144           - A -- a list of lists or a square N x N matrix, whose
     145             (i,j) entry gives the probability of transitioning from
     146             state i to state j.
     147
     148           - B -- a list of N pairs (mu,std), where if B[i]=(mu,std),
     149             then the probability distribution associated with state i
     150             normal with mean mu and standard deviation std.
     151
     152           - pi -- the probabilities of starting in each initial
     153             state, i.e,. pi[i] is the probability of starting in
     154             state i.
     155
     156           - normalize --bool (default: True); if given, input is
     157             normalized to define valid probability distributions,
     158             e.g., the entries of A are made nonnegative and the rows
     159             sum to 1.
     160
     161        EXAMPLES::
     162
     163
     164            sage: hmm.GaussianHiddenMarkovModel([[.1,.9],[.5,.5]], [(1,1), (-1,1)], [.5,.5])
     165            Gaussian Hidden Markov Model with 2 States
     166            Transition matrix:
     167            [0.1 0.9]
     168            [0.5 0.5]
     169            Emission parameters:
     170            [(1.0, 1.0), (-1.0, 1.0)]
     171            Initial probabilities: [0.5000, 0.5000]
     172
     173        We input a model in which both A and pi have to be
     174        renormalized to define valid probability distributions::
     175
     176            sage: hmm.GaussianHiddenMarkovModel([[-1,.7],[.3,.4]], [(1,1), (-1,1)], [-1,.3])
     177            Gaussian Hidden Markov Model with 2 States
     178            Transition matrix:
     179            [           0.0            1.0]
     180            [0.428571428571 0.571428571429]
     181            Emission parameters:
     182            [(1.0, 1.0), (-1.0, 1.0)]
     183            Initial probabilities: [0.0000, 1.0000]
     184
     185        Bad things can happen::
     186
     187            sage: hmm.GaussianHiddenMarkovModel([[-1,.7],[.3,.4]], [(1,1), (-1,1)], [-1,.3], normalize=False)
     188            Gaussian Hidden Markov Model with 2 States
     189            Transition matrix:
     190            [-1.0  0.7]
     191            [ 0.3  0.4]
     192            ...
    53193        """
    54         self.initialized = False
    55         HiddenMarkovModel.__init__(self, A, B, pi)
    56         self.m = <ghmm_cmodel*> safe_malloc(sizeof(ghmm_cmodel))
    57         # Set number of states
    58         self.m.N = self.A.nrows()
    59         # Assign model identifier (the name) if specified
    60         if name is not None:
    61             name = str(name)
    62             self.m.name = <char*> safe_malloc(len(name)+1)
    63             strcpy(self.m.name, name)
     194        self.pi = util.initial_probs_to_TimeSeries(pi, normalize)
     195        self.N = len(self.pi)
     196        self.A = util.state_matrix_to_TimeSeries(A, self.N, normalize)
     197
     198        # B should be a matrix of N rows, with column 0 the mean and 1
     199        # the standard deviation.
     200        if is_Matrix(B):
     201            B = B.list()
    64202        else:
    65             self.m.name = NULL
    66            
     203            B = flatten(B)
     204        self.B = TimeSeries(B)
     205        self.probability_init()
    67206
    68     def name(self):
     207    def __cmp__(self, other):
    69208        """
    70         Return the name of this model.
     209        Compare self and other, which must both be GaussianHiddenMarkovModel's.
     210
     211        EXAMPLES::
     212
     213            sage: m = hmm.GaussianHiddenMarkovModel([[1]], [(0,1)], [1])
     214            sage: n = hmm.GaussianHiddenMarkovModel([[1]], [(1,1)], [1])
     215            sage: m < n
     216            True
     217            sage: m == m
     218            True
     219            sage: n > m
     220            True
     221            sage: n < m
     222            False
     223        """
     224        if not isinstance(other, GaussianHiddenMarkovModel):
     225            raise ValueError
     226        return cmp(self.__reduce__()[1], other.__reduce__()[1])
     227
     228    def __getitem__(self, Py_ssize_t i):
     229        """
     230        Return the mean and standard distribution for the i-th state.
     231
     232        INPUT:
     233
     234            - i -- integer
    71235
    72236        OUTPUT:
    73             string or None
    74237
    75         EXAMPLES:
    76             sage: m = hmm.GaussianHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [(0.0,1.0),(1,1)], [1,0], 'Test Model')
    77             sage: m.name()
    78             'Test Model'
     238            - 2 floats
    79239
    80         If the model is not explicitly named then this function returns None:
    81             sage: m = hmm.GaussianHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [(0.0,1.0),(1,1)], [1,0])
    82             sage: m.name()
    83             sage: m.name() is None
     240        EXAMPLES::
     241
     242            sage: m = hmm.GaussianHiddenMarkovModel([[.1,.9],[.5,.5]], [(1,.5), (-2,.3)], [.5,.5])
     243            sage: m[0]
     244            (1.0, 0.5)
     245            sage: m[1]
     246            (-2.0, 0.29999999999999999)
     247            sage: m[-1]
     248            (-2.0, 0.29999999999999999)
     249            sage: m[3]
     250            Traceback (most recent call last):
     251            ...
     252            IndexError: index out of range
     253            sage: m[-3]
     254            Traceback (most recent call last):
     255            ...
     256            IndexError: index out of range
     257        """
     258        if i < 0:
     259            i += self.N
     260        if i < 0 or i >= self.N:
     261            raise IndexError, 'index out of range'
     262
     263        # TODO: change to be a normal distribution class
     264        return self.B[2*i], self.B[2*i+1]
     265
     266    def __reduce__(self):
     267        """
     268        Used in pickling.
     269
     270        EXAMPLES::
     271
     272            sage: G = hmm.GaussianHiddenMarkovModel([[1]], [(0,1)], [1])
     273            sage: loads(dumps(G)) == G
    84274            True
    85275        """
    86         if self.m.name:
    87             s = str(self.m.name)
    88             return s
     276        return unpickle_gaussian_hmm_v1, \
     277               (self.A, self.B, self.pi, self.prob, self.n_out)
     278
     279    def emission_parameters(self):
     280        """
     281        Return the parameters that define the normal distributions
     282        associated to all of the states.
     283
     284        OUTPUT:
     285
     286            - a list B of pairs B[i] = (mu, std), such that the
     287              distribution associated to state i is normal with mean
     288              mu and standard deviation std.
     289
     290        EXAMPLES::
     291
     292            sage: hmm.GaussianHiddenMarkovModel([[.1,.9],[.5,.5]], [(1,.5), (-1,3)], [.1,.9]).emission_parameters()
     293            [(1.0, 0.5), (-1.0, 3.0)]
     294        """
     295        cdef Py_ssize_t i
     296        from sage.rings.all import RDF
     297        return [(RDF(self.B[2*i]),RDF(self.B[2*i+1])) for i in range(self.N)]
     298
     299    def __repr__(self):
     300        """
     301        Return string representation.
     302
     303        EXAMPLES::
     304
     305            sage: hmm.GaussianHiddenMarkovModel([[.1,.9],[.5,.5]], [(1,.5), (-1,3)], [.1,.9]).__repr__()
     306            'Gaussian Hidden Markov Model with 2 States\nTransition matrix:\n[0.1 0.9]\n[0.5 0.5]\nEmission parameters:\n[(1.0, 0.5), (-1.0, 3.0)]\nInitial probabilities: [0.1000, 0.9000]'
     307        """
     308        s = "Gaussian Hidden Markov Model with %s States"%self.N
     309        s += '\nTransition matrix:\n%s'%self.transition_matrix()
     310        s += '\nEmission parameters:\n%s'%self.emission_parameters()
     311        s += '\nInitial probabilities: %s'%self.initial_probabilities()
     312        return s
     313
     314
     315    def generate_sequence(self, Py_ssize_t length):
     316        """
     317        Return a sample of the given length from this HMM.
     318
     319        INPUT:
     320
     321            - length -- positive integer
     322
     323        OUTPUT:
     324
     325            - an IntList or list of emission symbols
     326            - TimeSeries of emissions
     327
     328        EXAMPLES::
     329
     330            sage: m =hmm.GaussianHiddenMarkovModel([[.1,.9],[.5,.5]], [(1,.5), (-1,3)], [.1,.9])
     331            sage: m.generate_sequence(5)
     332            ([-3.0505, 0.5317, -4.5065, 0.6521, 1.0435], [1, 0, 1, 0, 1])
     333            sage: m.generate_sequence(0)
     334            ([], [])
     335            sage: m.generate_sequence(-1)
     336            Traceback (most recent call last):
     337            ...
     338            ValueError: length must be nonnegative
     339        """
     340        if length < 0:
     341            raise ValueError, "length must be nonnegative"
     342
     343        # Create Integer lists for states and observations
     344        cdef IntList states = IntList(length)
     345        cdef TimeSeries obs = TimeSeries(length)
     346        if length == 0:
     347            return states, obs
     348
     349        # Setup variables, including random state.
     350        cdef Py_ssize_t i, j
     351        cdef randstate rstate = current_randstate()
     352        cdef int q = 0
     353        cdef double r, accum
     354
     355        # Choose the starting state.
     356        # See the remark in hmm.pyx about how this should get
     357        # replaced by some general fast discrete distribution code.
     358        r = rstate.c_rand_double()
     359        accum = 0
     360        for i in range(self.N):
     361            if r < self.pi._values[i] + accum:
     362                q = i
     363            else:
     364                accum += self.pi._values[i]
     365
     366        states._values[0] = q
     367        obs._values[0] = self.random_sample(q, rstate)
     368
     369        cdef double* row
     370        cdef int O
     371        _sig_on
     372        for i in range(1, length):
     373            accum = 0
     374            row = self.A._values + q*self.N
     375            r = rstate.c_rand_double()
     376            for j in range(self.N):
     377                if r < row[j] + accum:
     378                    q = j
     379                    break
     380                else:
     381                    accum += row[j]
     382            states._values[i] = q
     383            obs._values[i] = self.random_sample(q, rstate)
     384        _sig_off
     385
     386        return obs, states
     387
     388    cdef probability_init(self):
     389        """
     390        Used internally to compute caching information that makes
     391        certain computations in the Baum-Welch algorithm faster.
     392        """
     393        self.prob = TimeSeries(2*self.N)
     394        cdef int i
     395        for i in range(self.N):
     396            self.prob[2*i] = 1.0/(sqrt2pi*self.B[2*i+1])
     397            self.prob[2*i+1] = -1.0/(2*self.B[2*i+1]*self.B[2*i+1])
     398
     399    cdef double random_sample(self, int state, randstate rstate):
     400        """
     401        Return a random sample from the normal distribution associated
     402        to the given state.
     403
     404        This is only used internally, and no bounds or other error
     405        checking is done, so calling this improperly can lead to seg
     406        faults.
     407
     408        INPUT:
     409
     410            - state -- integer
     411            - rstate -- randstate instance
     412
     413        OUTPUT:
     414
     415            - double
     416        """
     417        return random_normal(self.B._values[state*2], self.B._values[state*2+1], rstate)
     418
     419    cdef double probability_of(self, int state, double observation):
     420        """
     421        Return the probability b_j(o) of seeing the given observation
     422        given that we're in the given state j (=state).
     423
     424        This is a continuous probability, so this really returns a
     425        number p such that the probability of a value in the interval
     426        [o,o+d] is p*d.
     427
     428        INPUT:
     429
     430            - state -- integer
     431            - observation -- double
     432
     433        OUTPUT:
     434
     435            - double
     436        """
     437        # The code below is an optimized version of the following code:
     438        #     cdef double mean = self.B._values[2*state], \
     439        #                 std = self.B._values[2*state+1]
     440        #     return 1/(sqrt2pi*std) * \
     441        #             exp(-(observation-mean)*(observation-mean)/(2*std*std))
     442
     443        cdef double x = observation - self.B._values[2*state] # observation - mean
     444        return self.prob._values[2*state] * exp(x*x*self.prob._values[2*state+1])
     445
     446    def log_likelihood(self, obs):
     447        """
     448        Return the logarithm of the probability that this model
     449        produced the given observation sequence.
     450
     451        INPUT:
     452
     453            - obs -- sequence of observations
     454
     455        OUTPUT:
     456
     457            - float
     458
     459        EXAMPLES::
     460
     461            sage: m = hmm.GaussianHiddenMarkovModel([[.1,.9],[.5,.5]], [(1,.5), (-1,3)], [.1,.9])
     462            sage: m.log_likelihood([1,1,1])
     463            -4.2978807660724856
     464            sage: set_random_seed(0); s = m.sample(20)
     465            sage: m.log_likelihood(s)
     466            -40.115714129484...
     467        """
     468        if len(obs) == 0:
     469            return 1.0
     470        if not isinstance(obs, TimeSeries):
     471            obs = TimeSeries(obs)
     472        return self._forward_scale(obs)
     473
     474    def _forward_scale(self, TimeSeries obs):
     475        """
     476        Memory-efficient implementation of the forward algorithm (with scaling).
     477
     478        INPUT:
     479
     480            - obs -- an integer list of observation states.
     481
     482        OUTPUT:
     483
     484            - float -- the log of the probability that the model
     485              produced this sequence
     486
     487        EXAMPLES::
     488
     489            sage: m = hmm.GaussianHiddenMarkovModel([[.1,.9],[.5,.5]], [(1,.5), (-1,3)], [.1,.9])
     490            sage: m._forward_scale(stats.TimeSeries([1,-1,-1,1]))
     491            -7.641988207069133
     492        """
     493        cdef Py_ssize_t i, j, t, T = len(obs)
     494
     495        # The running some of the log probabilities
     496        cdef double log_probability = 0, sum, a
     497
     498        cdef TimeSeries alpha  = TimeSeries(self.N), \
     499                        alpha2 = TimeSeries(self.N)
     500
     501        # Initialization
     502        sum = 0
     503        for i in range(self.N):
     504            a = self.pi[i] * self.probability_of(i, obs._values[0])
     505            alpha[i] = a
     506            sum += a
     507
     508        log_probability = log(sum)
     509        alpha.rescale(1/sum)
     510
     511        # Induction
     512        cdef double s
     513        for t in range(1, T):
     514            sum = 0
     515            for j in range(self.N):
     516                s = 0
     517                for i in range(self.N):
     518                    s += alpha._values[i] * self.A._values[i*self.N + j]
     519                a = s * self.probability_of(j, obs._values[t])
     520                alpha2._values[j] = a
     521                sum += a
     522
     523            log_probability += log(sum)
     524            for j in range(self.N):
     525                alpha._values[j] = alpha2._values[j] / sum
     526
     527        # Termination
     528        return log_probability
     529
     530    def viterbi(self, obs):
     531        """
     532        Determine "the" hidden sequence of states that is most likely
     533        to produce the given sequence seq of observations, along with
     534        the probability that this hidden sequence actually produced
     535        the observation.
     536
     537        INPUT:
     538
     539            - seq -- sequence of emitted ints or symbols
     540
     541        OUTPUT:
     542
     543            - list -- "the" most probable sequence of hidden states, i.e.,
     544              the Viterbi path.
     545
     546            - float -- log of probability that the observed sequence
     547              was produced by the Viterbi sequence of states.
     548
     549        EXAMPLES::
     550
     551        We find the optimal state sequence for a given model::
     552
     553            sage: m = hmm.GaussianHiddenMarkovModel([[0.5,0.5],[0.5,0.5]], [(0,1),(10,1)], [0.5,0.5])
     554            sage: m.viterbi([0,1,10,10,1])
     555            ([0, 0, 1, 1, 0], -9.0604285688230899)
     556
     557        Another example in which the most likely states change based
     558        on the last observation::
     559
     560            sage: m = hmm.GaussianHiddenMarkovModel([[.1,.9],[.5,.5]], [(1,.5), (-1,3)], [.1,.9])
     561            sage: m.viterbi([-2,-1,.1,0.1])
     562            ([1, 1, 0, 1], -9.61823698847639...)
     563            sage: m.viterbi([-2,-1,.1,0.3])
     564            ([1, 1, 1, 0], -9.5660236533785135)
     565        """
     566        cdef TimeSeries _obs
     567        if not isinstance(obs, TimeSeries):
     568            _obs = TimeSeries(obs)
    89569        else:
    90             return None
    91            
     570            _obs = obs
    92571
    93 cdef class GaussianHiddenMarkovModel(ContinuousHiddenMarkovModel):
    94     r"""
    95     Create a Gaussian Hidden Markov Model.  The probability
    96     distribution associated with each state is a Gaussian
    97     distribution.
     572        # The algorithm is the same as _viterbi above, except
     573        # we take the logarithms of everything first, and add
     574        # instead of multiply.
     575        cdef Py_ssize_t t, T = _obs._length
     576        cdef IntList state_sequence = IntList(T)
     577        if T == 0:
     578            return state_sequence, 0.0
    98579
    99     GaussianHiddenMarkovModel(A, B, pi, name)
     580        cdef int i, j, N = self.N
     581
     582        # delta[i] is the maximum of the probabilities over all
     583        # paths ending in state i.
     584        cdef TimeSeries delta = TimeSeries(N), delta_prev = TimeSeries(N)
     585
     586        # We view psi as an N x T matrix of ints. The quantity
     587        #          psi[N*t + j]
     588        # is a most probable hidden state at time t-1, given
     589        # that j is a most probable state at time j.
     590        cdef IntList psi = IntList(N * T)  # initialized to 0 by default
     591
     592        # Log Preprocessing
     593        cdef TimeSeries A = self.A.log()
     594        cdef TimeSeries pi = self.pi.log()
     595
     596        # Initialization:
     597        for i in range(N):
     598            delta._values[i] = pi._values[i] + log(self.probability_of(i, _obs._values[0]))
     599
     600        # Recursion:
     601        cdef double mx, tmp, minus_inf = float('-inf')
     602        cdef int index
     603
     604        for t in range(1, T):
     605            delta_prev, delta = delta, delta_prev
     606            for j in range(N):
     607                # Compute delta_t[j] = max_i(delta_{t-1}(i) a_{i,j}) * b_j(_obs[t])
     608                mx = minus_inf
     609                index = -1
     610                for i in range(N):
     611                    tmp = delta_prev._values[i] + A._values[i*N+j]
     612                    if tmp > mx:
     613                        mx = tmp
     614                        index = i
     615                delta._values[j] = mx + log(self.probability_of(j, _obs._values[t]))
     616                psi._values[N*t + j] = index
     617
     618        # Termination:
     619        mx, index = delta.max(index=True)
     620
     621        # Path (state sequence) backtracking:
     622        state_sequence._values[T-1] = index
     623        t = T-2
     624        while t >= 0:
     625            state_sequence._values[t] = psi._values[N*(t+1) + state_sequence._values[t+1]]
     626            t -= 1
     627
     628        return state_sequence, mx
     629
     630    cdef TimeSeries _backward_scale_all(self, TimeSeries obs, TimeSeries scale):
     631        """
     632        This function returns the matrix beta_t(i), and is used
     633        internally as part of the Baum-Welch algorithm.
     634
     635        The quantity beta_t(i) is the probability of observing the
     636        sequence obs[t+1:] assuming that the model is in state i at
     637        time t.
     638
     639        INPUT:
     640
     641            - obs -- TimeSeries
     642            - scale -- TimeSeries
     643
     644        OUTPUT:
     645
     646            - TimeSeries beta such that beta_t(i) = beta[t*N + i]
     647            - scale is also changed by this function
     648        """
     649        cdef Py_ssize_t t, T = obs._length
     650        cdef int N = self.N, i, j
     651        cdef double s
     652        cdef TimeSeries beta = TimeSeries(N*T, initialize=False)
     653
     654        # 1. Initialization
     655        for i in range(N):
     656            beta._values[(T-1)*N + i] = 1 / scale._values[T-1]
     657
     658        # 2. Induction
     659        t = T-2
     660        while t >= 0:
     661            for i in range(N):
     662                s = 0
     663                for j in range(N):
     664                    s += self.A._values[i*N+j] * \
     665                         self.probability_of(j, obs._values[t+1]) * beta._values[(t+1)*N+j]
     666                beta._values[t*N + i] = s/scale._values[t]
     667            t -= 1
     668        return beta
     669
     670    cdef _forward_scale_all(self, TimeSeries obs):
     671        """
     672        Return scaled values alpha_t(i), the sequence of scalings, and
     673        the log probability.
     674
     675        The quantity alpha_t(i) is the probability of observing the
     676        sequence obs[:t+1] assuming that the model is in state i at
     677        time t.
     678
     679        INPUT:
     680
     681            - obs -- TimeSeries
     682
     683        OUTPUT:
     684
     685            - TimeSeries alpha with alpha_t(i) = alpha[t*N + i]
     686            - TimeSeries scale with scale[t] the scaling at step t
     687            - float -- log_probability of the observation sequence
     688              being produced by the model.
     689        """
     690        cdef Py_ssize_t i, j, t, T = len(obs)
     691        cdef int N = self.N
     692
     693        # The running some of the log probabilities
     694        cdef double log_probability = 0, sum, a
     695
     696        cdef TimeSeries alpha = TimeSeries(N*T, initialize=False)
     697        cdef TimeSeries scale = TimeSeries(T, initialize=False)
     698
     699        # Initialization
     700        sum = 0
     701        for i in range(self.N):
     702            a = self.pi._values[i] * self.probability_of(i, obs._values[0])
     703            alpha._values[0*N + i] = a
     704            sum += a
     705
     706        scale._values[0] = sum
     707        log_probability = log(sum)
     708        for i in range(self.N):
     709            alpha._values[0*N + i] /= sum
     710
     711        # Induction
     712        # The code below is just an optimized version of:
     713        #     alpha = (alpha * A).pairwise_product(B[O[t+1]])
     714        #     alpha = alpha.scale(1/alpha.sum())
     715        # along with keeping track of the log of the scaling factor.
     716        cdef double s
     717        for t in range(1,T):
     718            sum = 0
     719            for j in range(N):
     720                s = 0
     721                for i in range(N):
     722                    s += alpha._values[(t-1)*N + i] * self.A._values[i*N + j]
     723                a = s * self.probability_of(j, obs._values[t])
     724                alpha._values[t*N + j] = a
     725                sum += a
     726
     727            log_probability += log(sum)
     728            scale._values[t] = sum
     729            for j in range(self.N):
     730                alpha._values[t*N + j] /= sum
     731
     732        # Termination
     733        return alpha, scale, log_probability
     734
     735    cdef TimeSeries _baum_welch_gamma(self, TimeSeries alpha, TimeSeries beta):
     736        # TODO: factor out to hmm.pyx
     737        cdef int j, N = self.N
     738        cdef Py_ssize_t t, T = alpha._length//N
     739        cdef TimeSeries gamma = TimeSeries(alpha._length, initialize=False)
     740        cdef double denominator
     741        for t in range(T):
     742            denominator = 0
     743            for j in range(N):
     744                gamma._values[t*N + j] = alpha._values[t*N + j] * beta._values[t*N + j]
     745                denominator += gamma._values[t*N + j]
     746            for j in range(N):
     747                gamma._values[t*N + j] /= denominator
     748        return gamma
     749
     750    cdef TimeSeries _baum_welch_xi(self, TimeSeries alpha, TimeSeries beta, TimeSeries obs):
     751        # TODO: factor out to hmm.pyx
     752        """
     753        Return 3-dimensional array of xi values.
     754        We have x[t,i,j] = x[t*N*N + i*N + j].
     755        """
     756        cdef int i, j, N = self.N
     757        cdef double sum
     758        cdef Py_ssize_t t, T = alpha._length//N
     759        cdef TimeSeries xi = TimeSeries(T*N*N, initialize=False)
     760        for t in range(T-1):
     761            sum = 0.0
     762            for i in range(N):
     763                for j in range(N):
     764                    xi._values[t*N*N+i*N+j] = alpha._values[t*N+i]*beta._values[(t+1)*N+j]*\
     765                       self.A._values[i*N+j] * self.probability_of(j, obs._values[t+1])
     766                    sum += xi._values[t*N*N+i*N+j]
     767            for i in range(N):
     768                for j in range(N):
     769                    xi._values[t*N*N+i*N+j] /= sum
     770        return xi
     771
     772    def baum_welch(self, obs, int max_iter=500, double log_likelihood_cutoff=1e-4, double eps=1e-12, bint fix_emissions=False):
     773        """
     774        Given an observation sequence obs, improve this HMM using the
     775        Baum-Welch algorithm to increase the probability of observing obs.
     776
     777        INPUT:
     778
     779            - obs -- a time series of emissions
     780
     781            - max_iter -- integer (default: 500) maximum number
     782              of Baum-Welch steps to take
     783
     784            - log_likehood_cutoff -- positive float (default: 1e-4);
     785              the minimal improvement in likelihood with respect to
     786              the last iteration required to continue. Relative value
     787              to log likelihood.
     788
     789            - eps -- very small positive float (default: 1e-12); used
     790              during the iteration to replace numbers, e.g., standard
     791              deviations that are 0 but are not allowed to be 0.
     792
     793            - fix_emissions -- bool (default: False); if True, do not
     794              change emissions when updating
     795
     796        OUTPUT:
     797
     798            - changes the model in places, and returns the log
     799              likelihood and number of iterations.
     800
     801        EXAMPLES::
     802
     803            sage: m = hmm.GaussianHiddenMarkovModel([[.1,.9],[.5,.5]], [(1,.5), (-1,3)], [.1,.9])
     804            sage: m.log_likelihood([-2,-1,.1,0.1])
     805            -8.8582822159862751
     806            sage: m.baum_welch([-2,-1,.1,0.1])
     807            (22.164539478647512, 8)
     808            sage: m.log_likelihood([-2,-1,.1,0.1])
     809            22.164539478647512
     810            sage: m
     811            Gaussian Hidden Markov Model with 2 States
     812            Transition matrix:
     813            [              1.0 1.99514384486e-21]
     814            [   0.499999997782    0.500000002218]
     815            Emission parameters:
     816            [(0.09999999999..., 1.48492424379e-06), (-1.4999999929, 0.500000010249)]
     817            Initial probabilities: [0.0000, 1.0000]
     818
     819        We watch the log likelihoods of the model converge, step by step::
     820
     821            sage: m = hmm.GaussianHiddenMarkovModel([[.1,.9],[.5,.5]], [(1,.5), (-1,3)], [.1,.9])
     822            sage: v = m.sample(10)
     823            sage: stats.TimeSeries([m.baum_welch(v,max_iter=1)[0] for _ in range(len(v))])
     824            [-20.1167, -17.7611, -16.9814, -16.9364, -16.9314, -16.9309, -16.9309, -16.9309, -16.9309, -16.9309]
     825
     826        We illustrate fixing emissions::
     827
     828            sage: m = hmm.GaussianHiddenMarkovModel([[.1,.9],[.9,.1]], [(1,2),(-1,.5)], [.3,.7])
     829            sage: set_random_seed(0); v = m.sample(100)
     830            sage: m.baum_welch(v,fix_emissions=True)
     831            (-164.72944548204..., 23)
     832            sage: m.emission_parameters()
     833            [(1.0, 2.0), (-1.0, 0.5)]
     834            sage: m = hmm.GaussianHiddenMarkovModel([[.1,.9],[.9,.1]], [(1,2),(-1,.5)], [.3,.7])
     835            sage: m.baum_welch(v)
     836            (-162.85437039799817, 49)
     837            sage: m.emission_parameters()
     838            [(1.27224191726, 2.37136875176), (-0.948617467518, 0.576236038512)]
     839        """
     840        if not isinstance(obs, TimeSeries):
     841            obs = TimeSeries(obs)
     842        cdef TimeSeries _obs = obs
     843        cdef TimeSeries alpha, beta, scale, gamma, xi
     844        cdef double log_probability, log_probability0, log_probability_prev, delta
     845        cdef int i, j, k, N, n_iterations
     846        cdef Py_ssize_t t, T
     847        cdef double denominator_A, numerator_A, denominator_B, numerator_mean, numerator_std
     848
     849        # Initialization
     850        alpha, scale, log_probability0 = self._forward_scale_all(_obs)
     851        if not isfinite(log_probability0):
     852            return (0.0, 0)
     853        log_probability = log_probability0
     854        beta = self._backward_scale_all(_obs, scale)
     855        gamma = self._baum_welch_gamma(alpha, beta)
     856        xi = self._baum_welch_xi(alpha, beta, _obs)
     857        log_probability_prev = log_probability
     858        N = self.N
     859        n_iterations = 0
     860        T = len(_obs)
     861
     862        # Re-estimation
     863        while True:
     864
     865            # Reestimate
     866            for i in range(N):
     867                if not gamma._values[0*N+i]:
     868                    gamma._values[0*N+i] = eps
     869                elif not isfinite(gamma._values[0*N+i]):
     870                    util.normalize_probability_TimeSeries(self.pi, 0, self.pi._length)
     871                    raise RuntimeError, "impossible to compute gamma during reestimation"
     872                self.pi._values[i] = gamma._values[0*N+i]
     873
     874            # Update the probabilities pi to define a valid discrete distribution
     875            util.normalize_probability_TimeSeries(self.pi, 0, self.pi._length)
     876
     877            # Reestimate transition matrix and emission probabilities in
     878            # each state.
     879            for i in range(N):
     880                # Compute the updated transition matrix
     881                denominator_A = 0.0
     882                for t in range(T-1):
     883                    denominator_A += gamma._values[t*N+i]
     884                if not isnormal(denominator_A):
     885                    raise RuntimeError, "unable to re-estimate transition matrix"
     886                for j in range(N):
     887                    numerator_A = 0.0
     888                    for t in range(T-1):
     889                        numerator_A += xi._values[t*N*N+i*N+j]
     890                    self.A._values[i*N+j] = numerator_A / denominator_A
     891
     892                # Rescale the i-th row of the transition matrix to be
     893                # a valid stochastic matrix:
     894                util.normalize_probability_TimeSeries(self.A, i*N, (i+1)*N)
     895
     896                if not fix_emissions:
     897                    denominator_B = denominator_A + gamma._values[(T-1)*N + i]
     898                    if not isnormal(denominator_B):
     899                        raise RuntimeError, "unable to re-estimate emission probabilities"
     900
     901                    numerator_mean = 0.0
     902                    numerator_std = 0.0
     903                    for t in range(T):
     904                        numerator_mean += gamma._values[t*N + i] * _obs._values[t]
     905                        numerator_std  += gamma._values[t*N + i] * \
     906                               (_obs._values[t] - self.B._values[2*i])*(_obs._values[t] - self.B._values[2*i])
     907                    # restimated mean
     908                    self.B._values[2*i] = numerator_mean / denominator_B
     909                    # restimated standard deviation
     910                    self.B._values[2*i+1] = sqrt(numerator_std / denominator_B)
     911                    if not self.B._values[2*i+1]:
     912                        self.B._values[2*i+1] = eps
     913                    self.probability_init()
     914
     915            n_iterations += 1
     916            if n_iterations >= max_iter: break
     917
     918            # Initialization for next iteration
     919            alpha, scale, log_probability0 = self._forward_scale_all(_obs)
     920            if not isfinite(log_probability0): break
     921            log_probability = log_probability0
     922            beta = self._backward_scale_all(_obs, scale)
     923            gamma = self._baum_welch_gamma(alpha, beta)
     924            xi = self._baum_welch_xi(alpha, beta, _obs)
     925
     926            # Compute the difference between the log probability of
     927            # two iterations.
     928            delta = log_probability - log_probability_prev
     929            log_probability_prev = log_probability
     930
     931            # If the log probability does not change by more than delta,
     932            # then terminate
     933            if delta >= 0 and delta <= log_likelihood_cutoff:
     934                break
     935
     936        return log_probability, n_iterations
     937
     938
     939cdef class GaussianMixtureHiddenMarkovModel(GaussianHiddenMarkovModel):
     940    """
     941    GaussianMixtureHiddenMarkovModel(A, B, pi)
     942
     943    Gaussian mixture Hidden Markov Model.
    100944
    101945    INPUT:
    102         A  -- matrix; the transition matrix (n x n)
    103         B  -- list of n pairs (mu, sigma) that define the
    104               Gaussian distributions associated to each state,
    105               where mu is the mean and sigma the standard deviation.
    106         pi -- list of floats that sums to 1.0; these are
    107               the initial probabilities of each hidden state
    108         name -- (default: None); a string
    109         normalize -- (optional; default=True) whether or not to
    110                      normalize the model so the probabilities add to 1
     946        - A  -- matrix; the N x N transition matrix
     947        - B -- list of pairs (mu,sigma) that define the distributions
     948        - pi -- initial state probabilities
    111949
    112     EXAMPLES:
    113     Define the transition matrix:
    114         sage: A = [[0,1,0],[0.5,0,0.5],[0.3,0.3,0.4]]
    115950
    116     Parameters of the normal emission distributions in pairs of (mu, sigma):
    117         sage: B = [(0,1), (-1,0.5), (1,0.2)]
     951    EXAMPLES::
    118952
    119     The initial probabilities per state:
    120         sage: pi = [1,0,0]
     953        sage: A  = [[0.5,0.5],[0.5,0.5]]
     954        sage: B  = [[(0.9,(0.0,1.0)), (0.1,(1,10000))],[(1,(1,1)), (0,(0,0.1))]]
     955        sage: hmm.GaussianMixtureHiddenMarkovModel(A, B, [1,0])
     956        Gaussian Mixture Hidden Markov Model with 2 States
     957        Transition matrix:
     958        [0.5 0.5]
     959        [0.5 0.5]
     960        Emission parameters:
     961        [0.9*N(0.0,1.0) + 0.1*N(1.0,10000.0), 1.0*N(1.0,1.0) + 0.0*N(0.0,0.1)]
     962        Initial probabilities: [1.0000, 0.0000]
    121963
    122     We create the continuous Gaussian hidden Markov model defined by $A,B,\pi$:
    123         sage: hmm.GaussianHiddenMarkovModel(A, B, pi)
    124         Gaussian Hidden Markov Model with 3 States
     964    TESTS:
     965
     966    If a standard deviation is 0, it is normalized to be slightly bigger than 0.::
     967
     968        sage: hmm.GaussianMixtureHiddenMarkovModel([[1]], [[(1,(0,0))]], [1])
     969        Gaussian Mixture Hidden Markov Model with 1 States
    125970        Transition matrix:
    126         [0.0 1.0 0.0]
    127         [0.5 0.0 0.5]
    128         [0.3 0.3 0.4]
     971        [1.0]
    129972        Emission parameters:
    130         [(0.0, 1.0), (-1.0, 0.5), (1.0, 0.20000000000000001)]
    131         Initial probabilities: [1.0, 0.0, 0.0]
     973        [1.0*N(0.0,1e-08)]
     974        Initial probabilities: [1.0000]
     975
     976    We test that number of emission distributions must be the same as the number of states::
     977
     978        sage: hmm.GaussianMixtureHiddenMarkovModel([[1]], [], [1])
     979        Traceback (most recent call last):
     980        ...
     981        ValueError: number of GaussianMixtures must be the same as number of entries of pi
     982
     983        sage: hmm.GaussianMixtureHiddenMarkovModel([[1]], [[]], [1])
     984        Traceback (most recent call last):
     985        ...
     986        ValueError: must specify at least one component of the mixture model
     987
     988    We test that the number of initial probabilities must equal the number of states::
     989
     990        sage: hmm.GaussianMixtureHiddenMarkovModel([[1]], [[]], [1,2])
     991        Traceback (most recent call last):
     992        ...
     993        ValueError: number of entries of transition matrix A must be the square of the number of entries of pi
    132994    """
    133     def __init__(self, A, B, pi=None, name=None, normalize=True):
     995
     996    cdef object mixture # mixture
     997
     998    def __init__(self, A, B, pi=None, bint normalize=True):
    134999        """
    135         EXAMPLES:
    136         We make a very simple model:
    137             sage: hmm.GaussianHiddenMarkovModel([[1]], [(0,1)], [1], 'simple')
    138             Gaussian Hidden Markov Model simple with 1 States
     1000        EXAMPLES::
     1001
     1002            sage: hmm.GaussianMixtureHiddenMarkovModel([[.9,.1],[.4,.6]], [[(.4,(0,1)), (.6,(1,0.1))],[(1,(0,1))]], [.7,.3])
     1003            Gaussian Mixture Hidden Markov Model with 2 States
    1391004            Transition matrix:
    140             [1.0]
     1005            [0.9 0.1]
     1006            [0.4 0.6]
    1411007            Emission parameters:
    142             [(0.0, 1.0)]
    143             Initial probabilities: [1.0]
     1008            [0.4*N(0.0,1.0) + 0.6*N(1.0,0.1), 1.0*N(0.0,1.0)]
     1009            Initial probabilities: [0.7000, 0.3000]
     1010        """
     1011        self.pi = util.initial_probs_to_TimeSeries(pi, normalize)
     1012        self.N = len(self.pi)
     1013        self.A = util.state_matrix_to_TimeSeries(A, self.N, normalize)
     1014        if self.N*self.N != len(self.A):
     1015            raise ValueError, "number of entries of transition matrix A must be the square of the number of entries of pi"
    1441016
    145         We test a degenerate case:
    146             sage: hmm.GaussianHiddenMarkovModel([], [], [], 'simple')
    147             Gaussian Hidden Markov Model simple with 0 States
    148             Transition matrix:
    149             []
    150             Emission parameters:
    151             []
    152             Initial probabilities: []
     1017        self.mixture = [b if isinstance(b, GaussianMixtureDistribution) else \
     1018                            GaussianMixtureDistribution([flatten(x) for x in b]) for b in B]
     1019        if len(self.mixture) != self.N:
     1020            raise ValueError, "number of GaussianMixtures must be the same as number of entries of pi"
    1531021
    154         TESTS:
    155         We test normalize=False:
    156             sage: hmm.GaussianHiddenMarkovModel([[1,100],[0,1]], [(0,1),(0,2)], [1/2,1/2], normalize=False).transition_matrix()
    157             [  1.0 100.0]
    158             [  0.0   1.0]           
     1022    def __repr__(self):
    1591023        """
    160         ContinuousHiddenMarkovModel.__init__(self, A, B, pi, name=name)
     1024        Return string representation.
    1611025
    162         # Set number of outputs. 
    163         self.m.M = 1
    164         # Set the model type to continuous
    165         self.m.model_type = GHMM_kContinuousHMM
    166         # 1 transition matrix
    167         self.m.cos   =  1
    168         # Set that no a prior model probabilities are set.
    169         self.m.prior = -1
    170         # Dimension is 1
    171         self.m.dim   =  1
    172         self._initialize_state(pi)
    173         self.m.class_change = NULL
    174         self.initialized = True
    175         if normalize:
    176             self.normalize()
     1026        EXAMPLES::
    1771027
    178     def _initialize_state(self, pi):
     1028            sage: hmm.GaussianMixtureHiddenMarkovModel([[.9,.1],[.4,.6]], [[(.4,(0,1)), (.6,(1,0.1))],[(1,(0,1))]], [.7,.3]).__repr__()
     1029            'Gaussian Mixture Hidden Markov Model with 2 States\nTransition matrix:\n[0.9 0.1]\n[0.4 0.6]\nEmission parameters:\n[0.4*N(0.0,1.0) + 0.6*N(1.0,0.1), 1.0*N(0.0,1.0)]\nInitial probabilities: [0.7000, 0.3000]'
    1791030        """
    180         Allocate and initialize states.
    181 
    182         INPUT:
    183             pi -- initial probabilities
    184 
    185         All other inputs are set as self.A and self.B.
    186 
    187         EXAMPLES:
    188         This function is called implicitly during object creation.  It should
    189         never be called directly by the user, unless they want to LEAK MEMORY.
    190 
    191             sage: m = hmm.GaussianHiddenMarkovModel([[1]], [(1,10)], [1]) # indirect test
    192         """
    193         cdef ghmm_cstate* states = <ghmm_cstate*> safe_malloc(sizeof(ghmm_cstate) * self.m.N)
    194         cdef ghmm_cstate* state
    195         cdef ghmm_c_emission* e
    196         cdef Py_ssize_t i, j, k
    197 
    198         for i in range(self.m.N):
    199             # Parameters of normal distribution
    200             mu, sigma   = self.B[i]
    201             # Get a reference to the i-th state for convenience of the notation below.
    202             state = &(states[i])
    203             state.M     = 1
    204             state.pi    = pi[i]
    205             state.desc  = NULL
    206             state.fix   = 0
    207             e = <ghmm_c_emission*> safe_malloc(sizeof(ghmm_c_emission))
    208             e.type      = 0  # normal
    209             e.dimension = 1
    210             e.mean.val  = mu
    211             e.variance.val = sigma*sigma  # variance! not standard deviation
    212             # fixing of emissions is deactivated by default           
    213             e.fixed     = 0
    214             e.sigmacd   = NULL
    215             e.sigmainv  = NULL
    216             state.e     = e
    217             state.c     = to_double_array([1.0])
    218 
    219             #########################################################
    220             # Initialize state transition data.
    221             # NOTE: This code is similar to a block of code in hmm.pyx.
    222            
    223             # Set "out" probabilities, i.e., the probabilities to
    224             # transition to another hidden state from this state.
    225             v = self.A[i]
    226             k = self.m.N
    227             state.out_states = k
    228             state.out_id = <int*> safe_malloc(sizeof(int)*k)
    229             state.out_a  = ighmm_cmatrix_alloc(1, k)
    230             for j in range(k):
    231                 state.out_id[j] = j
    232                 state.out_a[0][j]  = v[j]
    233 
    234             # Set "in" probabilities
    235             v = self.A.column(i)
    236             state.in_states = k
    237             state.in_id = <int*> safe_malloc(sizeof(int)*k)
    238             state.in_a  = ighmm_cmatrix_alloc(1, k)
    239             for j in range(k):
    240                 state.in_id[j] = j
    241                 state.in_a[0][j]  = v[j]
    242 
    243             #########################################################               
    244         # Set states
    245         self.m.s = states
     1031        s = "Gaussian Mixture Hidden Markov Model with %s States"%self.N
     1032        s += '\nTransition matrix:\n%s'%self.transition_matrix()
     1033        s += '\nEmission parameters:\n%s'%self.emission_parameters()
     1034        s += '\nInitial probabilities: %s'%self.initial_probabilities()
     1035        return s
    2461036
    2471037    def __reduce__(self):
    2481038        """
    2491039        Used in pickling.
    2501040
    2511041        EXAMPLES:
    252             sage: m = hmm.GaussianHiddenMarkovModel([[1]], [(0,1)], [1], 'test')
    253             sage: f,g = m.__reduce__()
    254             sage: f(*g) == m
     1042            sage: m = hmm.GaussianMixtureHiddenMarkovModel([[1]], [[(.4,(0,1)), (.6,(1,0.1))]], [1])
     1043            sage: loads(dumps(m)) == m
    2551044            True
    2561045        """
    257         return unpickle_gaussian_hmm_v0, (self.transition_matrix(), self.emission_parameters(),
    258                              self.initial_probabilities(), self.name())
     1046        return unpickle_gaussian_mixture_hmm_v1, \
     1047               (self.A, self.B, self.pi, self.mixture)
    2591048
    260     def __dealloc__(self):
    261         """
    262         Dealloc the memory used by this Gaussian HMM, but only if the
    263         HMM was successfully initialized.
    264 
    265         EXAMPLES: 
    266             sage: m = hmm.GaussianHiddenMarkovModel([[1.0]], [(0.0,1.0)], [1])   # implicit doctest
    267             sage: del m
    268         """
    269         if self.initialized:
    270             ghmm_cmodel_free(&self.m)
    271 
    272     def __copy__(self):
    273         """
    274         Return a copy of this Gaussian HMM.
    275        
    276         EXAMPLES:
    277             sage: m = hmm.GaussianHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [(0.0,1.0),(1,1)], [1,0], 'NAME')
    278             sage: copy(m)
    279             Gaussian Hidden Markov Model NAME with 2 States
    280             Transition matrix:
    281             [0.4 0.6]
    282             [0.1 0.9]
    283             Emission parameters:
    284             [(0.0, 1.0), (1.0, 1.0)]
    285             Initial probabilities: [1.0, 0.0]
    286         """
    287        
    288         return GaussianHiddenMarkovModel(self.transition_matrix(), self.emission_parameters(),
    289                                          self.initial_probabilities(), self.name())
    2901049
    2911050    def __cmp__(self, other):
    2921051        """
    293         Compare two Gaussian HMM's.
     1052        Compare self and other, which must both be GaussianMixtureHiddenMarkovModel's.
     1053
     1054        EXAMPLES::
     1055
     1056            sage: m = hmm.GaussianMixtureHiddenMarkovModel([[1]], [[(.4,(0,1)), (.6,(1,0.1))]], [1])
     1057            sage: n = hmm.GaussianMixtureHiddenMarkovModel([[1]], [[(.5,(0,1)), (.5,(1,0.1))]], [1])
     1058            sage: m < n
     1059            True
     1060            sage: m == m
     1061            True
     1062            sage: n > m
     1063            True
     1064            sage: n < m
     1065            False
     1066        """
     1067        if not isinstance(other, GaussianMixtureHiddenMarkovModel):
     1068            raise ValueError
     1069        return cmp(self.__reduce__()[1], other.__reduce__()[1])
     1070
     1071    def __getitem__(self, Py_ssize_t i):
     1072        """
     1073        Return the Gaussian mixture distribution associated to the
     1074        i-th state.
    2941075
    2951076        INPUT:
    296             self, other -- objects; if other is not a Gaussian HMM compare types.
    297         OUTPUT:
    298             -1,0,1
    2991077
    300         The transition matrices are compared, then the emission
    301         parameters, and the initial probabilities.  The names are not
    302         compared, so two GHMM's with the same parameters but different
    303         names compare to be equal.
    304 
    305         EXAMPLES:
    306             sage: m = hmm.GaussianHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [(0.0,1.0),(1,1)], [1,2], "Test 1", normalize=False)
    307             sage: m.__cmp__(m)
    308             0
    309 
    310         Note that the name doesn't matter:
    311             sage: n = hmm.GaussianHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [(0.0,1.0),(1,1)], [1,2], "Test 2", normalize=False)
    312             sage: m.__cmp__(n)
    313             0
    314 
    315         Normalizing fixes the initial probabilities, hence m and n are no longer equal.
    316             sage: n.normalize()
    317             sage: m.__cmp__(n)
    318             1
    319         """
    320         if not isinstance(other, GaussianHiddenMarkovModel):
    321             return cmp(type(self), type(other))
    322 
    323         if self is other: return 0  # easy special case
    324        
    325         cdef GaussianHiddenMarkovModel o = other
    326         if self.m.N < o.m.N:
    327             return -1
    328         elif self.m.N > o.m.N:
    329             return 1
    330         cdef Py_ssize_t i, j
    331 
    332         # The code below is somewhat long and tedious because I want it to be
    333         # very fast.  All it does is explicitly loop through the transition
    334         # matrix, emission parameters and initial state probabilities checking
    335         # that they agree, and if not returning -1 or 1.
    336         # Compare transition matrices
    337         for i from 0 <= i < self.m.N:
    338             for j from 0 <= j < self.m.s[i].out_states:
    339                 if self.m.s[i].out_a[0][j] < o.m.s[i].out_a[0][j]:
    340                     return -1
    341                 elif self.m.s[i].out_a[0][j] > o.m.s[i].out_a[0][j]:
    342                     return 1
    343 
    344         # Compare emissions parameters
    345         for i from 0 <= i < self.m.N:
    346             if self.m.s[i].e.mean.val < o.m.s[i].e.mean.val:
    347                 return -1
    348             elif self.m.s[i].e.mean.val > o.m.s[i].e.mean.val:
    349                 return 1
    350             if self.m.s[i].e.variance.val < o.m.s[i].e.variance.val:
    351                 return -1
    352             elif self.m.s[i].e.variance.val > o.m.s[i].e.variance.val:
    353                 return 1
    354            
    355         # Compare initial state probabilities
    356         for 0 <= i < self.m.N:
    357             if self.m.s[i].pi < o.m.s[i].pi:
    358                 return -1
    359             elif self.m.s[i].pi > o.m.s[i].pi:
    360                 return 1
    361 
    362         return 0
    363 
    364     def fix_emissions(self, Py_ssize_t i, bint fixed=True):
    365         """
    366         Sets the Gaussian emission parameters for the i-th state to be
    367         either fixed or not fixed.  If it is fixed, then running the
    368         Baum-Welch algorithm will not change the emission parameters
    369         for the i-th state.
    370        
    371         INPUT:
    372             i -- nonnegative integer < self.m.N
    373             fixed -- bool
    374 
    375         EXAMPLES:
    376         We run Baum-Welch once without fixing the emission states:
    377             sage: m = hmm.GaussianHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [(0.0,1.0),(1,1)], [1,0])
    378             sage: m.baum_welch([0,1])
    379             sage: m
    380             Gaussian Hidden Markov Model with 2 States
    381             Transition matrix:
    382             [0.0 1.0]
    383             [0.1 0.9]
    384             Emission parameters:
    385             [(0.0, 0.01), (1.0, 0.01)]
    386             Initial probabilities: [1.0, 0.0]
    387 
    388         Now we run Baum-Welch with the emission states fixed.  Notice that they don't change.
    389             sage: m = hmm.GaussianHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [(0.0,1.0),(1,1)], [1,0])
    390             sage: m.fix_emissions(0); m.fix_emissions(1)
    391             sage: m.baum_welch([0,1])
    392             sage: m
    393             Gaussian Hidden Markov Model with 2 States
    394             Transition matrix:
    395             [0.000368587006957    0.999631412993]
    396             [              0.1               0.9]
    397             Emission parameters:
    398             [(0.0, 1.0), (1.0, 1.0)]
    399             Initial probabilities: [..., 0.0]
    400         """
    401         if i < 0 or i >= self.m.N:
    402             raise IndexError, "index out of range"
    403         self.m.s[i].e.fixed = fixed
    404 
    405     def __repr__(self):
    406         """
    407         Return string representation of this Continuous HMM.
     1078            - i -- integer
    4081079
    4091080        OUTPUT:
    410             string
    4111081
    412         EXAMPLES:
    413             sage: m = hmm.GaussianHiddenMarkovModel([[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])
    414             sage: m.__repr__()
    415             'Gaussian Hidden Markov Model with 3 States\nTransition matrix:\n[0.0 1.0 0.0]\n[0.5 0.0 0.5]\n[0.3 0.3 0.4]\nEmission parameters:\n[(0.0, 1.0), (-1.0, 0.5), (1.0, 0.20000000000000001)]\nInitial probabilities: [1.0, 0.0, 0.0]'
     1082            - a Gaussian mixture distribution object
     1083
     1084        EXAMPLES::
     1085
     1086            sage: m = hmm.GaussianMixtureHiddenMarkovModel([[.9,.1],[.4,.6]], [[(.4,(0,1)), (.6,(1,0.1))],[(1,(0,1))]], [.7,.3])
     1087            sage: m[0]
     1088            0.4*N(0.0,1.0) + 0.6*N(1.0,0.1)
     1089            sage: m[1]
     1090            1.0*N(0.0,1.0)
     1091
     1092        Negative indexing works::
     1093
     1094            sage: m[-1]
     1095            1.0*N(0.0,1.0)
     1096
     1097        Bounds are checked::
     1098
     1099            sage: m[2]
     1100            Traceback (most recent call last):
     1101            ...
     1102            IndexError: index out of range
     1103            sage: m[-3]
     1104            Traceback (most recent call last):
     1105            ...
     1106            IndexError: index out of range
    4161107        """
    417         s = "Gaussian Hidden Markov Model%s with %s States"%(
    418             ' ' + self.m.name if self.m.name else '',
    419             self.m.N)
    420         s += '\nTransition matrix:\n%s'%self.transition_matrix()
    421         s += '\nEmission parameters:\n%s'%self.emission_parameters()
    422         s += '\nInitial probabilities: %s'%self.initial_probabilities()
    423         return s
    424 
    425     def initial_probabilities(self):
    426         """
    427         Return the list of initial state probabilities.
    428 
    429         OUTPUT:
    430             list of floats
    431        
    432         EXAMPLES:
    433             sage: m = hmm.GaussianHiddenMarkovModel([[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)], [0.4,0.3,0.3])
    434             sage: m.initial_probabilities()
    435             [0.40000000000000002, 0.29999999999999999, 0.29999999999999999]
    436         """
    437         cdef Py_ssize_t i
    438         return [self.m.s[i].pi for i in range(self.m.N)]
    439    
    440     def transition_matrix(self):
    441         """
    442         Return the hidden state transition matrix.
    443 
    444         OUTPUT:
    445             matrix whose rows give the transition probabilities of the
    446             Hidden Markov Model states.
    447        
    448         EXAMPLES:
    449             sage: m = hmm.GaussianHiddenMarkovModel([[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])
    450             sage: m.transition_matrix()
    451             [0.0 1.0 0.0]
    452             [0.5 0.0 0.5]
    453             [0.3 0.3 0.4]
    454         """
    455         cdef Py_ssize_t i, j
    456         # Update the state of the "immutable" matrix A, then return a reference to it.
    457         for i from 0 <= i < self.m.N:
    458             for j from 0 <= j < self.m.s[i].out_states:
    459                 self.A.set_unsafe_double(i,j,self.m.s[i].out_a[0][j])
    460         return self.A
     1108        if i < 0:
     1109            i += self.N
     1110        if i < 0 or i >= self.N:
     1111            raise IndexError, 'index out of range'
     1112        return self.mixture[i]
    4611113
    4621114    def emission_parameters(self):
    4631115        """
    464         Return the emission parameters list.
     1116        Returns a list of all the emission distributions.
    4651117
    4661118        OUTPUT:
    467             list of tuples (mu, sigma) that define Gaussian distributions associated to each state.
    468             Here mu is the mean and sigma the standard deviation.
    469        
    470         EXAMPLES:
    471             sage: m = hmm.GaussianHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [(1.5,2),(-1,3)], [1,0], 'NAME')
     1119
     1120            - list of Gaussian mixtures
     1121
     1122        EXAMPLES::
     1123
     1124            sage: m = hmm.GaussianMixtureHiddenMarkovModel([[.9,.1],[.4,.6]], [[(.4,(0,1)), (.6,(1,0.1))],[(1,(0,1))]], [.7,.3])
    4721125            sage: m.emission_parameters()
    473             [(1.5, 2.0), (-1.0, 3.0)]
     1126            [0.4*N(0.0,1.0) + 0.6*N(1.0,0.1), 1.0*N(0.0,1.0)]
    4741127        """
    475         cdef Py_ssize_t i       
    476         return [(self.m.s[i].e.mean.val, sqrt(self.m.s[i].e.variance.val)) for i in range(self.m.N)]
     1128        return list(self.mixture)
    4771129
    478     def normalize(self):
     1130    cdef double random_sample(self, int state, randstate rstate):
    4791131        """
    480         Normalize the transition and emission probabilities, if
    481         applicable.  This changes self in place.
     1132        Return a random sample from the normal distribution associated
     1133        to the given state.
    4821134
    483         EXAMPLES:
    484             sage: m = hmm.GaussianHiddenMarkovModel([[1.0,1.2],[0,0.1]], [(0.0,1.0),(1,1)], [1,2])
    485             sage: m.normalize()
    486             sage: m
    487             Gaussian Hidden Markov Model with 2 States
    488             Transition matrix:
    489             [0.454545454545 0.545454545455]
    490             [           0.0            1.0]
    491             Emission parameters:
    492             [(0.0, 1.0), (1.0, 1.0)]
    493             Initial probabilities: [0.33333333333333331, 0.66666666666666663]
    494         """
    495         if ghmm_cmodel_normalize(self.m):
    496             raise RuntimeError, "error normalizing model (note that model may still have been partly changed)"
    497 
    498     def sample(self, long length, number=None):
    499         """
    500         Return number samples from this HMM of given length.
     1135        This is only used internally, and no bounds or other error
     1136        checking is done, so calling this improperly can lead to seg
     1137        faults.
    5011138
    5021139        INPUT:
    503             length -- positive integer
    504            
     1140
     1141            - state -- integer
     1142            - rstate -- randstate instance
     1143
    5051144        OUTPUT:
    506             if number is not given, return a single TimeSeries.
    507             if number is given, return a list of TimeSeries.
    5081145
    509         EXAMPLES:
    510             sage: m = hmm.GaussianHiddenMarkovModel([[1]], [(0,1)], [1])
    511             sage: set_random_seed(0)
    512             sage: m.sample(5,2)
    513             [[-2.2808, -0.0699, 0.1858, 1.3624, 1.8252],
    514              [0.0080, 0.1244, 0.5098, 0.9961, 0.4710]]
     1146            - double
     1147        """
     1148        cdef GaussianMixtureDistribution G = self.mixture[state]
     1149        return G._sample(rstate)
    5151150
    516             sage: m = hmm.GaussianHiddenMarkovModel([[1]], [(0,1)], [1])
    517             sage: set_random_seed(0)
    518             sage: m.sample(5)
    519             [-2.2808, -0.0699, 0.1858, 1.3624, 1.8252]
     1151    cdef double probability_of(self, int state, double observation):
    5201152        """
    521         if number is None:
    522             single = True
    523             number = 1
    524         else:
    525             single = False
    526         seed = random()
    527         cdef ghmm_cseq *d = ghmm_cmodel_generate_sequences(self.m, seed, length, number, length)
    528         cdef Py_ssize_t i, j
    529         cdef TimeSeries T
    530         ans = []
    531         for j from 0 <= j < number:
    532             T = TimeSeries(length)
    533             for i from 0 <= i < length:
    534                 T._values[i] = d.seq[j][i]
    535             ans.append(T)
    536         if single:
    537             return ans[0]
    538         return ans
    539         # The line below would replace the above by something that returns a list of lists.
    540         #return [[d.seq[j][i] for i in range(length)] for j in range(number)]
     1153        Return the probability b_j(o) of see the given observation o
     1154        (=observation) given that we're in the given state j (=state).
    5411155
    542 
    543            
    544 
    545     ####################################################################
    546     # HMM Problem 1 -- Computing likelihood: Given the parameter set
    547     # lambda of an HMM model and an observation sequence O, determine
    548     # the likelihood P(O | lambda).
    549     ####################################################################
    550     def log_likelihood(self, seq):
    551         r"""
    552         HMM Problem 1: Likelihood.
    553         Computes sum over all sequence of seq_w * log( P ( O|lambda )).
     1156        This is a continuous probability, so this really returns a
     1157        number p such that the probability of a value in the interval
     1158        [o,o+d] is p*d.
    5541159
    5551160        INPUT:
    556             seq -- a single TimeSeries, or
    557                    a list of object z where z is either a TimeSeries
    558                    or a pair (TimeSeries, float), where float is a positive
    559                    weight.    When the weight isn't given it defaults to 1.
     1161
     1162            - state -- integer
     1163            - observation -- double
    5601164
    5611165        OUTPUT:
    562             float -- the weight sum of logs of the probability of
    563                      observing these sequences given the model
    5641166
    565         EXAMPLES:
    566         We create a very simple GHMM that generates random numbers with mean 0
    567         and standard deviation 1.
    568             sage: m = hmm.GaussianHiddenMarkovModel([[1]], [(0,1)], [1])
     1167            - double
     1168        """
     1169        cdef GaussianMixtureDistribution G = self.mixture[state]
     1170        return G.prob(observation)
    5691171
    570         We compute the log probability of a certain series:
    571             sage: m.log_likelihood(finance.TimeSeries([1,0,1,1]))
    572             -5.1757541328186907
     1172    # TODO: factor out to hmm.pyx
     1173    cdef TimeSeries _baum_welch_gamma(self, TimeSeries alpha, TimeSeries beta):
     1174        """
     1175        Returns gamma and gamma_all, where gamma_all contains the
     1176        (double-precision float) numbers gamma_t(j,m), and gamma
     1177        contains just the gamma_t(j).
     1178        """
     1179        cdef int j, N = self.N
     1180        cdef Py_ssize_t t, T = alpha._length//N
     1181        cdef TimeSeries gamma = TimeSeries(alpha._length, initialize=False)
     1182        cdef double denominator
     1183        for t in range(T):
     1184            denominator = 0
     1185            for j in range(N):
     1186                gamma._values[t*N + j] = alpha._values[t*N + j] * beta._values[t*N + j]
     1187                denominator += gamma._values[t*N + j]
     1188            for j in range(N):
     1189                gamma._values[t*N + j] /= denominator
    5731190
    574         We compute the log probability of another much less likely sequence:
    575             sage: m.log_likelihood(finance.TimeSeries([1,0,1,20]))
    576             -204.6757541328187
     1191        return gamma
    5771192
    578         We compute weighted sum of log probabilities of two sequences:
    579             sage: m.log_likelihood([ ([1,0,1,1], 10),  ([1,0,1,20], 0.1)  ])
    580             -72.2251167414687...
    581 
    582         We verify that the above weight computation is right.
    583             sage: 10 * m.log_likelihood([1,0,1,1]) + 0.1 * m.log_likelihood([1,0,1,20])
    584             -72.2251167414688
    585 
    586         We generate two normally distributed sequences and see that the one
    587         with mean 0 and standard deviation 1 is vastly more likely given
    588         the model than the one with mean 10 and standard deviation 1.
    589             sage: set_random_seed(0)
    590             sage: m.log_likelihood(finance.TimeSeries(100).randomize('normal',0,1))
    591             -129.87209711900121
    592             sage: m.log_likelihood(finance.TimeSeries(100).randomize('normal',10,1))
    593             -5010.151947016132
    594 
    595         However, if we change the model then the situation reverses:
    596             sage: m = hmm.GaussianHiddenMarkovModel([[1]], [(10,1)], [1])
    597             sage: m.log_likelihood(finance.TimeSeries(100).randomize('normal',10,1))
    598             -149.68737352182211
    599             sage: m.log_likelihood(finance.TimeSeries(100).randomize('normal',0,1))
    600             -5275.30829407876...
     1193    cdef TimeSeries _baum_welch_mixed_gamma(self, TimeSeries alpha, TimeSeries beta,
     1194                                            TimeSeries obs, int j):
    6011195        """
    602 
    603         cdef double log_p
    604         cdef ghmm_cseq* sqd
    605         try:
    606             sqd = to_cseq(seq)
    607         except RuntimeError:
    608             # no sequences
    609             return float(0)
    610         cdef int ret = ghmm_cmodel_likelihood(self.m, sqd, &log_p)
    611         ghmm_cseq_free(&sqd)
    612         if ret == -1:
    613             raise RuntimeError, "unable to compute log likelihood"
    614         return log_p
    615        
    616 
    617     ####################################################################
    618     # HMM Problem 2 -- Decoding: Given the complete parameter set that
    619     # defines the model and an observation sequence seq, determine the
    620     # best hidden sequence Q.  Use the Viterbi algorithm.
    621     ####################################################################   
    622     def viterbi(self, seq):
    623         """
    624         HMM Problem 2: Decoding.  Determine a hidden sequence of
    625         states that is most likely to produce the given sequence seq
    626         of observations.
     1196        Let gamma_t(j,m) be the m-component (in the mixture) of the
     1197        probability of being in state j at time t, given the
     1198        observation sequence.  This function outputs a TimeSeries v
     1199        such that v[m*T + t] gives gamma_t(j, m) where T is the number
     1200        of time steps.
    6271201
    6281202        INPUT:
    629             seq -- TimeSeries
    630            
     1203
     1204            - alpha -- TimeSeries
     1205            - beta -- TimeSeries
     1206            - obs -- TimeSeries
     1207            - j -- int
     1208
    6311209        OUTPUT:
    632             list -- a most probable sequence of hidden states, i.e., the
    633                     Viterbi path.
    634             float -- log of the probability that the sequence of hidden
    635                      states actually produced the given sequence seq.
    6361210
    637         EXAMPLES:
    638         We find the optimal state sequence for a given model.
    639             sage: m = hmm.GaussianHiddenMarkovModel([[0.5,0.5],[0.5,0.5]], [(0,1),(10,1)], [0.5,0.5])
    640             sage: m.viterbi([0,1,10,10,1])
    641             ([0, 0, 1, 1, 0], -9.0604285688230899)
     1211            - TimeSeries
    6421212        """
    643         cdef double log_p
    644         if not isinstance(seq, TimeSeries):
    645             seq = TimeSeries(seq)
    646         cdef TimeSeries T = seq
    647         cdef int* path = ghmm_cmodel_viterbi(self.m, T._values, T._length, &log_p)
    648         if not path:
    649             raise RuntimeError, "sequence can't be built from model"
    650         cdef Py_ssize_t i
    651         v = [path[i] for i in range(T._length)]
    652         sage_free(path)
    653         return v, log_p
     1213        cdef int i, k, m, N = self.N
     1214        cdef Py_ssize_t t, T = alpha._length//N
    6541215
    655     ####################################################################
    656     # HMM Problem 3 -- Learning: Given an observation sequence O and
    657     # the set of states in the HMM, improve the HMM to increase the
    658     # probability of observing O.
    659     ####################################################################
    660     def baum_welch(self, training_seqs, max_iter=10000, log_likelihood_cutoff=0.00001):
     1216        cdef double numer, alpha_minus, P, s, prob
     1217        cdef GaussianMixtureDistribution G = self.mixture[j]
     1218        cdef int M = len(G)
     1219        cdef TimeSeries mixed_gamma = TimeSeries(T*M)
     1220
     1221        for t in range(T):
     1222            prob = self.probability_of(j, obs._values[t])
     1223            if prob == 0:
     1224                # If the probability of observing obs[t] in state j is 0, then
     1225                # all of the m-mixture components *have* to automatically be 0,
     1226                # since prob is the sum of those and they are all nonnegative.
     1227                for m in range(M):
     1228                    mixed_gamma._values[m*T + t] = 0
     1229            else:
     1230                # Compute the denominator we used when scaling gamma.
     1231                # The key thing is that this is consistent between
     1232                # gamma and mixed_gamma.
     1233                P = 0
     1234                for k in range(N):
     1235                    P += alpha._values[t*N+k]*beta._values[t*N+k]
     1236
     1237                # Divide out the total probability, so we can multiply back in
     1238                # the m-components of the probability.
     1239                alpha_minus = alpha._values[t*N + j] / prob
     1240                for m in range(M):
     1241                    numer =  alpha_minus * G.prob_m(obs._values[t], m) * beta._values[t*N + j]
     1242                    mixed_gamma._values[m*T + t] = numer / P
     1243
     1244        return mixed_gamma
     1245
     1246    def baum_welch(self, obs, int max_iter=1000, double log_likelihood_cutoff=1e-12, double eps=1e-12, bint fix_emissions=False):
    6611247        """
    662         HMM Problem 3: Learning.  Given an observation sequence O and
    663         the set of states in the HMM, improve the HMM using the
    664         Baum-Welch algorithm to increase the probability of observing O.
     1248        Given an observation sequence obs, improve this HMM using the
     1249        Baum-Welch algorithm to increase the probability of observing obs.
    6651250
    6661251        INPUT:
    667             training_seqs -- a list of lists of emission symbols, where all sequences of
    668                       length 0 are ignored; or, a list of pairs
    669                             (sample_sequence, weight),
    670                       where sample_sequence is a list or TimeSeries, and weight is
    671                       a positive real number.
    672             max_iter -- integer or None (default: 10000) maximum number
    673                       of Baum-Welch steps to take
    674             log_likehood_cutoff -- positive float or None (default: 0.00001);
    675                       the minimal improvement in likelihood
    676                       with respect to the last iteration required to
    677                       continue. Relative value to log likelihood
     1252
     1253            - obs -- a time series of emissions
     1254            - max_iter -- integer (default: 1000) maximum number
     1255              of Baum-Welch steps to take
     1256            - log_likehood_cutoff -- positive float (default: 1e-12);
     1257              the minimal improvement in likelihood with respect to
     1258              the last iteration required to continue. Relative value
     1259              to log likelihood.
     1260            - eps -- very small positive float (default: 1e-12); used
     1261              during the iteration to replace numbers, e.g., standard
     1262              deviations that are 0 but are not allowed to be 0.
     1263            - fix_emissions -- bool (default: False); if True, do not
     1264              change emissions when updating
    6781265
    6791266        OUTPUT:
    680             changes the model in places, or raises a RuntimError
    681             exception on error
    6821267
    683         EXAMPLES:
    684         We train a very simple model:
    685             sage: m = hmm.GaussianHiddenMarkovModel([[1]], [(0,1)], [1])
    686             sage: m.baum_welch([1,1,1,1])
     1268            - changes the model in places, and returns the log
     1269              likelihood and number of iterations.
    6871270
    688         Notice that after training the mean in the emission parameters changes
    689         form 0 to 1.
     1271        EXAMPLES::
     1272
     1273            sage: m = hmm.GaussianMixtureHiddenMarkovModel([[.9,.1],[.4,.6]], [[(.4,(0,1)), (.6,(1,0.1))],[(1,(0,1))]], [.7,.3])
     1274            sage: set_random_seed(0); v = m.sample(10); v
     1275            [0.3576, -0.9365, 0.9449, -0.6957, 1.0217, 0.9644, 0.9987, -0.5950, -1.0219, 0.6477]
     1276            sage: m.log_likelihood(v)
     1277            -8.31408655939536...
     1278            sage: m.baum_welch(v)
     1279            (2.18905068682301..., 15)
     1280            sage: m.log_likelihood(v)
     1281            2.18905068682301...
    6901282            sage: m
    691             Gaussian Hidden Markov Model with 1 States
     1283            Gaussian Mixture Hidden Markov Model with 2 States
    6921284            Transition matrix:
    693             [1.0]
     1285            [   0.874636333977    0.125363666023]
     1286            [              1.0 1.45168520229e-40]
    6941287            Emission parameters:
    695             [(1.0, 0.01)]
    696             Initial probabilities: [1.0]
     1288            [0.500161629343*N(-0.812298726239,0.173329026744) + 0.499838370657*N(0.982433690378,0.029719932009), 1.0*N(0.503260056832,0.145881515324)]
     1289            Initial probabilities: [0.0000, 1.0000]
    6971290
    698         We train using a list of lists:
    699             sage: m = hmm.GaussianHiddenMarkovModel([[1,0],[0,1]], [(0,1),(0,2)], [1/2,1/2])
    700             sage: m.baum_welch([[1,2], [3,2]])
    701             sage: m
    702             Gaussian Hidden Markov Model with 2 States
    703             Transition matrix:
    704             [1.0 0.0]
    705             [0.0 1.0]
    706             Emission parameters:
    707             [(1.94653953598434..., 0.705082962992410...), (2.02081569132933..., 0.70680033099099593)]
    708             Initial probabilities: [0.28024729110782..., 0.71975270889217...]
     1291        We illustrate fixing all emissions::
    7091292
    710         We train the same model, but waiting one of the lists more than the other.
    711             sage: m = hmm.GaussianHiddenMarkovModel([[1,0],[0,1]], [(0,1),(0,2)], [1/2,1/2])
    712             sage: m.baum_welch([([1,2,],10), ([3,2],1)])
    713             sage: m
    714             Gaussian Hidden Markov Model with 2 States
    715             Transition matrix:
    716             [1.0 0.0]
    717             [0.0 1.0]
    718             Emission parameters:
    719             [(1.58517861517798..., 0.572645807401051...), (1.594503566606473..., 0.57928632238916...)]
    720             Initial probabilities: [0.385468570528119..., 0.61453142947188...]
    721        
     1293            sage: m = hmm.GaussianMixtureHiddenMarkovModel([[.9,.1],[.4,.6]], [[(.4,(0,1)), (.6,(1,0.1))],[(1,(0,1))]], [.7,.3])
     1294            sage: set_random_seed(0); v = m.sample(10)
     1295            sage: m.baum_welch(v, fix_emissions=True)
     1296            (-7.5865685899788922, 36)
     1297            sage: m.emission_parameters()
     1298            [0.4*N(0.0,1.0) + 0.6*N(1.0,0.1), 1.0*N(0.0,1.0)]
     1299        """
     1300        if not isinstance(obs, TimeSeries):
     1301            obs = TimeSeries(obs)
     1302        cdef TimeSeries _obs = obs
     1303        cdef TimeSeries alpha, beta, scale, gamma, mixed_gamma, mixed_gamma_m, xi
     1304        cdef double log_probability, log_probability0, log_probability_prev, delta
     1305        cdef int i, j, k, m, N, n_iterations
     1306        cdef Py_ssize_t t, T
     1307        cdef double denominator_A, numerator_A, denominator_B, numerator_mean, numerator_std, \
     1308             numerator_c, c, mu, std, numer, denom, new_mu, new_std, new_c, s
     1309        cdef GaussianMixtureDistribution G
    7221310
    723         Training sequences of length 0 are gracefully ignored:
    724             sage: m.baum_welch([])
    725             sage: m.baum_welch([([],1)])
    726         """
    727         cdef ghmm_cmodel_baum_welch_context cs
     1311        # Initialization
     1312        alpha, scale, log_probability0 = self._forward_scale_all(_obs)
     1313        if not isfinite(log_probability0):
     1314            return (0.0, 0)
     1315        log_probability = log_probability0
     1316        beta = self._backward_scale_all(_obs, scale)
     1317        gamma = self._baum_welch_gamma(alpha, beta)
     1318        xi = self._baum_welch_xi(alpha, beta, _obs)
     1319        log_probability_prev = log_probability
     1320        N = self.N
     1321        n_iterations = 0
     1322        T = len(_obs)
    7281323
    729         cs.smo      = self.m
    730         try:
    731             cs.sqd      = to_cseq(training_seqs)
    732         except RuntimeError:
    733             # No nonempty sequences
    734             return
    735         cs.logp     = <double*> safe_malloc(sizeof(double))
    736         cs.eps      = log_likelihood_cutoff
    737         cs.max_iter = max_iter
    738        
    739         if ghmm_cmodel_baum_welch(&cs):
    740             raise RuntimeError, "error running Baum-Welch algorithm"
    741         ghmm_cseq_free(&cs.sqd)
    742        
    743        
    744 cdef ghmm_cseq* to_cseq(seq) except NULL:
    745     """
    746     Return a pointer to a ghmm_cseq C struct.
     1324        # Re-estimation
     1325        while True:
    7471326
    748     All empty sequences are ignored.  If there are no nonempty
    749     sequences a RuntimeError is raised, since GHMM doesn't treat
    750     this degenerate case well.   If there are any nonpositive
    751     weights, then a ValueError is raised.
    752     """
    753     if isinstance(seq, list) and len(seq) > 0 and not isinstance(seq[0], (list, tuple, TimeSeries)):
    754         seq = TimeSeries(seq)
    755     if isinstance(seq, TimeSeries):
    756         seq = [(seq,float(1))]
    757     cdef Py_ssize_t i, n
    758     for i in range(len(seq)):
    759         z = seq[i]
    760         if isinstance(z, tuple) and len(z) == 2:
    761             if isinstance(z[0],TimeSeries):
    762                 z = (z[0], float(z[1]))
    763             else:
    764                 z = (TimeSeries(z[0]), float(z[1]))
    765         else:
    766             if isinstance(z, TimeSeries):
    767                 z = (z, float(1))
    768             else:
    769                 z = (TimeSeries(z), float(1))
    770         seq[i] = z
    771     seq = [x for x in seq if len(x[0]) > 0]
     1327            # Reestimate frequency of state i in time t=0
     1328            for i in range(N):
     1329                if not gamma._values[0*N+i]:
     1330                    gamma._values[0*N+i] = eps
     1331                elif not isnormal(gamma._values[0*N+i]):
     1332                    util.normalize_probability_TimeSeries(self.pi, 0, self.pi._length)
     1333                    raise RuntimeError, "impossible to compute gamma during reestimation"
     1334                self.pi._values[i] = gamma._values[0*N+i]
    7721335
    773     n = len(seq)
    774     if n == 0:
    775         raise RuntimeError, "there must be at least one nonempty sequence"
     1336            # Update the probabilities pi to define a valid discrete distribution
     1337            util.normalize_probability_TimeSeries(self.pi, 0, self.pi._length)
    7761338
    777     for _, w in seq:
    778         if w <= 0:
    779             raise ValueError, "each weight must be positive"
    780    
    781     cdef ghmm_cseq* sqd = <ghmm_cseq*>safe_malloc(sizeof(ghmm_cseq))
    782     sqd.seq        = <double**>safe_malloc(sizeof(double*) * n)
    783     sqd.seq_len    = to_int_array([len(v) for v,_ in seq])
    784     sqd.seq_id     = to_double_array([0]*n)
    785     sqd.seq_label  = to_int_array([])  # obsolete but no choice
    786     weights        = [w for _, w in seq]
    787     sqd.seq_w      = to_double_array(weights)
    788     sqd.seq_number = n
    789     sqd.total_w    = sum(weights)
    790     sqd.dim        = 1
    791     sqd.flags      = 0
    792     sqd.capacity   = n
     1339            # Reestimate transition matrix and emission probabilities in
     1340            # each state.
     1341            for i in range(N):
     1342                # Reestimate the state transition matrix
     1343                denominator_A = 0.0
     1344                for t in range(T-1):
     1345                    denominator_A += gamma._values[t*N+i]
     1346                if not isnormal(denominator_A):
     1347                    raise RuntimeError, "unable to re-estimate pi (1)"
     1348                for j in range(N):
     1349                    numerator_A = 0.0
     1350                    for t in range(T-1):
     1351                        numerator_A += xi._values[t*N*N+i*N+j]
     1352                    self.A._values[i*N+j] = numerator_A / denominator_A
    7931353
    794     cdef TimeSeries T
    795     for i from 0 <= i < len(seq):
    796         T = seq[i][0]
    797         sqd.seq[i] = <double*>safe_malloc(sizeof(double) * T._length)
    798         memcpy(sqd.seq[i], T._values , sizeof(double)*T._length)
     1354                # Rescale the i-th row of the transition matrix to be
     1355                # a valid stochastic matrix:
     1356                util.normalize_probability_TimeSeries(self.A, i*N, (i+1)*N)
    7991357
    800     return sqd
     1358                ########################################################################
     1359                # Re-estimate the emission probabilities
     1360                ########################################################################
     1361                G = self.mixture[i]
     1362                if not fix_emissions and not G.is_fixed():
     1363                    mixed_gamma = self._baum_welch_mixed_gamma(alpha, beta, _obs, i)
     1364                    new_G = []
     1365                    for m in range(len(G)):
     1366                        if G.fixed._values[m]:
     1367                            new_G.append(G[m])
     1368                            continue
    8011369
     1370                        # Compute re-estimated mu_{j,m}
     1371                        numer = 0
     1372                        denom = 0
     1373                        for t in range(T):
     1374                            numer += mixed_gamma._values[m*T + t] * _obs._values[t]
     1375                            denom += mixed_gamma._values[m*T + t]
     1376                        new_mu = numer / denom
     1377
     1378                        # Compute re-estimated standard deviation
     1379                        numer = 0
     1380                        mu = G[m][1]
     1381                        for t in range(T):
     1382                            numer += mixed_gamma._values[m*T + t] * \
     1383                                     (_obs._values[t] - mu)*(_obs._values[t] - mu)
     1384
     1385                        new_std = sqrt(numer / denom)
     1386
     1387                        # Compute re-estimated weighting coefficient
     1388                        new_c = denom
     1389                        s = 0
     1390                        for t in range(T):
     1391                            s += gamma._values[t*N + i]
     1392                        new_c /= s
     1393
     1394                        new_G.append((new_c,new_mu,new_std))
     1395
     1396                    self.mixture[i] = GaussianMixtureDistribution(new_G)
     1397
     1398            n_iterations += 1
     1399            if n_iterations >= max_iter: break
     1400
     1401            ########################################################################
     1402            # Initialization for next iteration
     1403            ########################################################################
     1404            alpha, scale, log_probability0 = self._forward_scale_all(_obs)
     1405            if not isfinite(log_probability0): break
     1406            log_probability = log_probability0
     1407            beta = self._backward_scale_all(_obs, scale)
     1408            gamma = self._baum_welch_gamma(alpha, beta)
     1409            xi = self._baum_welch_xi(alpha, beta, _obs)
     1410
     1411            # Compute the difference between the log probability of
     1412            # two iterations.
     1413            delta = log_probability - log_probability_prev
     1414            log_probability_prev = log_probability
     1415
     1416            # If the log probability does not improve by more than
     1417            # delta, then terminate
     1418            if delta >= 0 and delta <= log_likelihood_cutoff:
     1419                break
     1420
     1421        return log_probability, n_iterations
     1422
     1423
     1424##################################################
     1425# For Unpickling
     1426##################################################
     1427
     1428# We keep the _v0 function for backwards compatible.
    8021429def unpickle_gaussian_hmm_v0(A, B, pi, name):
    8031430    """
    804     EXAMPLES:
    805         sage: m = hmm.GaussianHiddenMarkovModel([[1]], [(0,1)], [1], 'test')
    806         sage: loads(dumps(m)) == m
    807         True
     1431    EXAMPLES::
     1432
     1433        sage: m = hmm.GaussianHiddenMarkovModel([[1]], [(0,1)], [1])
    8081434        sage: sage.stats.hmm.chmm.unpickle_gaussian_hmm_v0(m.transition_matrix(), m.emission_parameters(), m.initial_probabilities(), 'test')
    809         Gaussian Hidden Markov Model test with 1 States
     1435        Gaussian Hidden Markov Model with 1 States
    8101436        Transition matrix:
    8111437        [1.0]
    8121438        Emission parameters:
    8131439        [(0.0, 1.0)]
    814         Initial probabilities: [1.0]
     1440        Initial probabilities: [1.0000]
    8151441    """
    816     return GaussianHiddenMarkovModel(A,B,pi,name)
     1442    return GaussianHiddenMarkovModel(A,B,pi)
    8171443
    8181444
    819 cdef class GaussianMixtureHiddenMarkovModel(GaussianHiddenMarkovModel):
     1445def unpickle_gaussian_hmm_v1(A, B, pi, prob, n_out):
    8201446    """
    821     GaussianMixtureHiddenMarkovModel(A, B, pi, name)
     1447    EXAMPLES::
    8221448
    823     INPUT:
    824         A  -- matrix; the transition matrix (n x n)
    825         B  -- list of lists of pairs (w, (mu, sigma)) that define the
    826               Gaussian mixture associated to each state, where w is
    827               the weight, mu is the mean and sigma the standard
    828               deviation.
    829         pi -- list of floats that sums to 1.0; these are
    830               the initial probabilities of each hidden state
    831         name -- (default: None); a string
    832         normalize -- (optional; default=True) whether or not to normalize
    833                      the model so the probabilities add to 1
     1449        sage: m = hmm.GaussianHiddenMarkovModel([[1]], [(0,1)], [1])
     1450        sage: loads(dumps(m)) == m   # indirect test
     1451        True
     1452    """
     1453    cdef GaussianHiddenMarkovModel m = PY_NEW(GaussianHiddenMarkovModel)
     1454    m.A = A
     1455    m.B = B
     1456    m.pi = pi
     1457    m.prob = prob
     1458    m.n_out = n_out
     1459    return m
    8341460
    835     EXAMPLES:
    836         sage: A  = [[0.5,0.5],[0.5,0.5]]
    837         sage: B  = [[(0.5,(0.0,1.0)), (0.1,(1,10000))],[(1,(1,1)), (0,(0,0.1))]]
    838         sage: pi = [1,0]
    839         sage: hmm.GaussianMixtureHiddenMarkovModel(A, B, pi)
    840         Gaussian Hidden Markov Model with 2 States
    841         Transition matrix:
    842         [0.5 0.5]
    843         [0.5 0.5]
    844         Emission parameters:
    845         [[(1.0, (0.0, 1.0)), (0.10000000000000001, (1.0, 10000.0))], [(1.0, (1.0, 1.0)), (0.0, (0.0, 0.10000000000000001))]]
    846         Initial probabilities: [1.0, 0.0]
    847        
     1461def unpickle_gaussian_mixture_hmm_v1(A, B, pi, mixture):
     1462    """
     1463    EXAMPLES::
    8481464
    849     TESTS:
    850     We test that standard deviations must be positive:
    851         sage: m = hmm.GaussianMixtureHiddenMarkovModel([[1]], [[(1,(0,0))]], [1])
    852         Traceback (most recent call last):
    853         ...
    854         ValueError: sigma must be positive (if weight is nonzero)
     1465        sage: m = hmm.GaussianMixtureHiddenMarkovModel([[1]], [[(.4,(0,1)), (.6,(1,0.1))]], [1])
     1466        sage: loads(dumps(m)) == m   # indirect test
     1467        True
     1468    """
     1469    cdef GaussianMixtureHiddenMarkovModel m = PY_NEW(GaussianMixtureHiddenMarkovModel)
     1470    m.A = A
     1471    m.B = B
     1472    m.pi = pi
     1473    m.mixture = mixture
     1474    return m
    8551475
    856     We test that number of mixtures must be positive:
    857         sage: m = hmm.GaussianMixtureHiddenMarkovModel([[1]], [[]], [1])
    858         Traceback (most recent call last):
    859         ...
    860         ValueError: number of Gaussian mixtures must be positive
    861     """
    862     def __init__(self, A, B, pi=None, name=None, normalize=True):
    863         """
    864         EXAMPLES:
    865             sage: hmm.GaussianMixtureHiddenMarkovModel([[1]], [[(1,(0,1))]], [1], name='simple')
    866             Gaussian Hidden Markov Model simple with 1 States
    867             Transition matrix:
    868             [1.0]
    869             Emission parameters:
    870             [[(1.0, (0.0, 1.0))]]
    871             Initial probabilities: [1.0]
    872         """
    873         # Turn B into a list of lists
    874         B = [flatten(x) for x in B]
    875         m = max([len(x) for x in B])
    876         if m == 0:
    877             raise ValueError, "number of Gaussian mixtures must be positive"
    878         B = [x + [0]*(m-len(x)) for x in B]
    879         GaussianHiddenMarkovModel.__init__(self, A, B, pi, name=name)
    880         self.m.M = m//3
    881         # Set number of outputs. 
    882 
    883     def _initialize_state(self, pi):
    884         """
    885         Allocate and initialize states.
    886 
    887         INPUT:
    888             pi -- initial probabilities
    889 
    890         All other inputs are set as self.A and self.B.
    891 
    892         EXAMPLES:
    893         This function is called implicitly during object creation.  It should
    894         never be called directly by the user, unless they want to LEAK MEMORY.
    895 
    896             sage: m = hmm.GaussianMixtureHiddenMarkovModel([[1]], [[(1,(1,10))]], [1]) # indirect test
    897         """
    898         cdef ghmm_cstate* states = <ghmm_cstate*> safe_malloc(sizeof(ghmm_cstate) * self.m.N)
    899         cdef ghmm_cstate* state
    900         cdef ghmm_c_emission* e
    901         cdef Py_ssize_t i, j, k, M, n
    902 
    903         for i in range(self.m.N):
    904             # Parameters of Gaussian distributions
    905             v = self.B[i]
    906             M = len(v)//3   # number of distinct Gaussians
    907 
    908             # Get a reference to the i-th state for convenience of the notation below.
    909             state = &(states[i])
    910             state.M     = M
    911             state.pi    = pi[i]
    912             state.desc  = NULL
    913             state.fix   = 0
    914             e           = <ghmm_c_emission*> safe_malloc(sizeof(ghmm_c_emission)*M)
    915             weights     = []
    916            
    917             for n in range(M):
    918                 e[n].type      = 0  # Gaussian
    919                 e[n].dimension = 1
    920                 mu             = v[n*3+1]
    921                 sigma          = v[n*3+2]
    922                 if sigma <= 0 and v[n*3]:
    923                     raise ValueError, "sigma must be positive (if weight is nonzero)"
    924                 weights.append(  v[n*3] )
    925                 e[n].mean.val     = mu
    926                 e[n].variance.val = sigma*sigma  # variance! not standard deviation
    927                
    928                 # fixing of emissions is deactivated by default           
    929                 e[n].fixed     = 0
    930                 e[n].sigmacd   = NULL
    931                 e[n].sigmainv  = NULL
    932                
    933             state.e     = e
    934             state.c     = to_double_array(weights)
    935 
    936             #########################################################
    937             # Initialize state transition data.
    938             # NOTE: This code is similar to a block of code in hmm.pyx.
    939            
    940             # Set "out" probabilities, i.e., the probabilities to
    941             # transition to another hidden state from this state.
    942             v = self.A[i]
    943             k = self.m.N
    944             state.out_states = k
    945             state.out_id = <int*> safe_malloc(sizeof(int)*k)
    946             state.out_a  = ighmm_cmatrix_alloc(1, k)
    947             for j in range(k):
    948                 state.out_id[j] = j
    949                 state.out_a[0][j]  = v[j]
    950 
    951             # Set "in" probabilities
    952             v = self.A.column(i)
    953             state.in_states = k
    954             state.in_id = <int*> safe_malloc(sizeof(int)*k)
    955             state.in_a  = ighmm_cmatrix_alloc(1, k)
    956             for j in range(k):
    957                 state.in_id[j] = j
    958                 state.in_a[0][j]  = v[j]
    959 
    960             #########################################################               
    961         # Set states
    962         self.m.s = states
    963 
    964     def emission_parameters(self):
    965         """
    966         Return the emission parameters list.
    967 
    968         OUTPUT:
    969            list of lists of tuples (weight, (mu, sigma))
    970 
    971         EXAMPLES:
    972             sage: m = hmm.GaussianMixtureHiddenMarkovModel([[0.5,0.5],[0.5,0.5]], [[(0.5,(0.0,1.0)), (0.1,(1,10000))],[(1,(1,1))]], [1,0])
    973             sage: m.emission_parameters()
    974             [[(1.0, (0.0, 1.0)), (0.10000000000000001, (1.0, 10000.0))], [(1.0, (1.0, 1.0)), (0.0, (0.0, 0.0))]]
    975         """
    976         cdef Py_ssize_t i,j
    977        
    978         return [[(self.m.s[i].c[j], (self.m.s[i].e[j].mean.val, sqrt(self.m.s[i].e[j].variance.val)))
    979                  for j in range(self.m.s[i].M)]  for i in range(self.m.N)]
    980 
    981        
    982 
    983     def fix_emissions(self, Py_ssize_t i, Py_ssize_t j, bint fixed=True):
    984         """
    985         Sets the j-th Gaussian of the emission parameters for the i-th
    986         state to be either fixed or not fixed.  If it is fixed, then
    987         running the Baum-Welch algorithm will not change the emission
    988         parameters for the i-th state.
    989        
    990         INPUT:
    991             i -- nonnegative integer < self.m.N
    992             j -- nonnegative integer < self.m.M
    993             fixed -- bool
    994 
    995         EXAMPLES:
    996         We run Baum-Welch once without fixing the emission states:
    997             sage: m = hmm.GaussianHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [(0.0,1.0),(1,1)], [1,0])
    998             sage: m.baum_welch([0,1])
    999             sage: m
    1000             Gaussian Hidden Markov Model with 2 States
    1001             Transition matrix:
    1002             [0.0 1.0]
    1003             [0.1 0.9]
    1004             Emission parameters:
    1005             [(0.0, 0.01), (1.0, 0.01)]
    1006             Initial probabilities: [1.0, 0.0]
    1007 
    1008         Now we run Baum-Welch with the emission states fixed.  Notice that they don't change.
    1009             sage: m = hmm.GaussianHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [(0.0,1.0),(1,1)], [1,0])
    1010             sage: m.fix_emissions(0); m.fix_emissions(1)
    1011             sage: m.baum_welch([0,1])
    1012             sage: m
    1013             Gaussian Hidden Markov Model with 2 States
    1014             Transition matrix:
    1015             [0.000368587006957    0.999631412993]
    1016             [              0.1               0.9]
    1017             Emission parameters:
    1018             [(0.0, 1.0), (1.0, 1.0)]
    1019             Initial probabilities: [..., 0.0]
    1020         """
    1021         if i < 0 or i >= self.m.N:
    1022             raise IndexError, "index i out of range"
    1023         if j < 0 or j >= self.m.M:
    1024             raise IndexError, "index j out of range"
    1025         self.m.s[i].e[j].fixed = fixed

  • new file sage/stats/hmm/distributions.pxd 

    diff --git a/sage/stats/hmm/distributions.pxd b/sage/stats/hmm/distributions.pxd
new file mode 100644
    - +  
     1#############################################################################
     2#       Copyright (C) 2010 William Stein <wstein@gmail.com>
     3#  Distributed under the terms of the GNU General Public License (GPL)
     4#  The full text of the GPL is available at:
     5#                  http://www.gnu.org/licenses/
     6#############################################################################
     7
     8from sage.finance.time_series cimport TimeSeries
     9from sage.stats.intlist cimport IntList
     10from sage.misc.randstate cimport randstate
     11
     12cdef class Distribution:
     13    pass
     14
     15cdef class DiscreteDistribution(Distribution):
     16    cdef object v
     17
     18cdef class GaussianDistribution(Distribution):
     19    pass
     20
     21cdef class GaussianMixtureDistribution(Distribution):
     22    cdef TimeSeries c0, c1, param
     23    cdef IntList fixed
     24   
     25    cdef double _sample(self, randstate rstate)
     26    cpdef double prob(self, double x)
     27    cpdef double prob_m(self, double x, int m)
     28    cpdef is_fixed(self, i=?)
     29
     30
     31

  • new file sage/stats/hmm/distributions.pyx 

    diff --git a/sage/stats/hmm/distributions.pyx b/sage/stats/hmm/distributions.pyx
new file mode 100644
    - +  
     1"""
     2Distributions used in implementing Hidden Markov Models
     3
     4These distribution classes are designed specifically for HMM's and not
     5for general use in statistics. For example, they have fixed or
     6non-fixed status, which only make sense relative to being used in a
     7hidden Markov model.
     8
     9AUTHOR:
     10
     11   - William Stein, 2010-03
     12"""
     13
     14#############################################################################
     15#       Copyright (C) 2010 William Stein <wstein@gmail.com>
     16#  Distributed under the terms of the GNU General Public License (GPL)
     17#  The full text of the GPL is available at:
     18#                  http://www.gnu.org/licenses/
     19#############################################################################
     20
     21include "../../ext/stdsage.pxi"
     22
     23cdef extern from "math.h":
     24    double exp(double)
     25    double log(double)
     26    double sqrt(double)
     27
     28import math
     29cdef double sqrt2pi = sqrt(2*math.pi)
     30
     31from sage.misc.randstate cimport current_randstate, randstate
     32from sage.finance.time_series cimport TimeSeries
     33
     34
     35
     36cdef double random_normal(double mean, double std, randstate rstate):
     37    """
     38    Return a floating point number chosen from the normal distribution
     39    with given mean and standard deviation, using the given randstate.
     40    The computation uses the box muller algorithm.
     41
     42    INPUT:
     43
     44        - mean -- float; the mean
     45        - std -- float; the standard deviation
     46        - rstate -- randstate; the random number generator state
     47
     48    OUTPUT:
     49
     50        - double
     51    """
     52    # Ported from http://users.tkk.fi/~nbeijar/soft/terrain/source_o2/boxmuller.c
     53    # This the box muller algorithm.
     54    # Client code can get the current random state from:
     55    #         cdef randstate rstate = current_randstate()
     56    cdef double x1, x2, w, y1, y2
     57    while True:
     58        x1 = 2*rstate.c_rand_double() - 1
     59        x2 = 2*rstate.c_rand_double() - 1
     60        w = x1*x1 + x2*x2
     61        if w < 1: break
     62    w = sqrt( (-2*log(w))/w )
     63    y1 = x1 * w
     64    return mean + y1*std
     65
     66# Abstract base class for distributions used for hidden Markov models.
     67
     68cdef class Distribution:
     69    """
     70    A distribution.
     71    """
     72    def sample(self, n=None):
     73        """
     74        Return either a single sample (the default) or n samples from
     75        this probability distribution.
     76
     77        INPUT:
     78
     79           - n -- None or a positive integer
     80
     81        OUTPUT:
     82
     83           - a single sample if n is 1; otherwise many samples
     84
     85        EXAMPLES::
     86
     87        This method must be defined in a derived class::
     88       
     89            sage: sage.stats.hmm.distributions.Distribution().sample()
     90            Traceback (most recent call last):
     91            ...
     92            NotImplementedError
     93        """
     94        raise NotImplementedError
     95
     96    def prob(self, x):
     97        """
     98        The probability density function evaluated at x.
     99
     100        INPUT:
     101
     102           - x -- object
     103
     104        OUTPUT:
     105
     106           - float
     107
     108        EXAMPLES::
     109
     110        This method must be defined in a derived class::
     111
     112            sage: sage.stats.hmm.distributions.Distribution().prob(0)
     113            Traceback (most recent call last):
     114            ...
     115            NotImplementedError
     116        """
     117        raise NotImplementedError
     118
     119    def plot(self, *args, **kwds):
     120        """
     121        Return a plot of the probability density function.
     122
     123        INPUT:
     124
     125            - args and kwds, passed to the Sage plot function
     126
     127        OUTPUT:
     128
     129            - a Graphics object
     130
     131        EXAMPLES::
     132       
     133            sage: P = hmm.GaussianMixtureDistribution([(.2,-10,.5),(.6,1,1),(.2,20,.5)])
     134            sage: P.plot(-10,30)
     135        """
     136        from sage.plot.all import plot
     137        return plot(self.prob, *args, **kwds)
     138
     139cdef class GaussianMixtureDistribution(Distribution):
     140    """
     141    A probability distribution defined by taking a weighted linear
     142    combination of Gaussian distributions.
     143
     144    EXAMPLES::
     145
     146        sage: P = hmm.GaussianMixtureDistribution([(.3,1,2),(.7,-1,1)]); P
     147        0.3*N(1.0,2.0) + 0.7*N(-1.0,1.0)
     148        sage: P[0]
     149        (0.29999999999999999, 1.0, 2.0)
     150        sage: P.is_fixed()
     151        False
     152        sage: P.fix(1)
     153        sage: P.is_fixed(0)
     154        False
     155        sage: P.is_fixed(1)
     156        True
     157        sage: P.unfix(1)
     158        sage: P.is_fixed(1)
     159        False
     160    """
     161    def __init__(self, B, eps=1e-8, bint normalize=True):
     162        """
     163        INPUT:
     164
     165            - `B` -- a list of triples `(c_i, mean_i, std_i)`, where
     166              the `c_i` and `std_i` are positive and the sum of the
     167              `c_i` is `1`.
     168
     169            - eps -- positive real number; any standard deviation in B
     170              less than eps is replaced by eps.
     171
     172            - normalize -- if True, ensure that the c_i are nonnegative
     173
     174        EXAMPLES::
     175
     176            sage: hmm.GaussianMixtureDistribution([(.3,1,2),(.7,-1,1)])
     177            0.3*N(1.0,2.0) + 0.7*N(-1.0,1.0)
     178            sage: hmm.GaussianMixtureDistribution([(1,-1,0)], eps=1e-3)
     179            1.0*N(-1.0,0.001)
     180        """
     181        B = [[c if c>=0 else 0,  mu,  std if std>0 else eps] for c,mu,std in B]
     182        if len(B) == 0:
     183            raise ValueError, "must specify at least one component of the mixture model"
     184        cdef double s
     185        if normalize:
     186            s = sum([a[0] for a in B])
     187            if s != 1:
     188                if s == 0:
     189                    s = 1.0/len(B)
     190                    for a in B:
     191                        a[0] = s
     192                else:
     193                    for a in B:
     194                        a[0] /= s
     195        self.c0 = TimeSeries([c/(sqrt2pi*std) for c,_,std in B])
     196        self.c1 = TimeSeries([-1.0/(2*std*std) for _,_,std in B])
     197        self.param = TimeSeries(sum([list(x) for x in B],[]))
     198        self.fixed = IntList(self.c0._length)
     199
     200    def __getitem__(self, Py_ssize_t i):
     201        """
     202        Returns triple (coefficient, mu, std).
     203
     204        INPUT:
     205
     206            - i -- integer
     207
     208        OUTPUT:
     209
     210            - triple of floats
     211
     212        EXAMPLES::
     213
     214            sage: P = hmm.GaussianMixtureDistribution([(.2,-10,.5),(.6,1,1),(.2,20,.5)])
     215            sage: P[0]
     216            (0.20000000000000001, -10.0, 0.5)
     217            sage: P[2]
     218            (0.20000000000000001, 20.0, 0.5)
     219            sage: [-1]
     220            [-1]
     221            sage: P[-1]
     222            (0.20000000000000001, 20.0, 0.5)
     223            sage: P[3]
     224            Traceback (most recent call last):
     225            ...
     226            IndexError: index out of range
     227            sage: P[-4]
     228            Traceback (most recent call last):
     229            ...
     230            IndexError: index out of range
     231        """
     232        if i < 0: i += self.param._length//3
     233        if i < 0 or i >= self.param._length//3: raise IndexError, "index out of range"
     234        return self.param._values[3*i], self.param._values[3*i+1], self.param._values[3*i+2]
     235
     236    def __reduce__(self):
     237        """
     238        Used in pickling.
     239
     240        EXAMPLES::
     241
     242            sage: G = sage.stats.hmm.distributions.GaussianMixtureDistribution([(.1,1,2), (.9,0,1)])
     243            sage: loads(dumps(G)) == G
     244            True
     245        """
     246        return unpickle_gaussian_mixture_distribution_v1, (
     247            self.c0, self.c1, self.param, self.fixed)
     248
     249    def __cmp__(self, other):
     250        """
     251        EXAMPLES::
     252
     253            sage: G = sage.stats.hmm.distributions.GaussianMixtureDistribution([(.1,1,2), (.9,0,1)])
     254            sage: H = sage.stats.hmm.distributions.GaussianMixtureDistribution([(.3,1,2), (.7,1,5)])
     255            sage: G < H
     256            True
     257            sage: H > G
     258            True
     259            sage: G == H
     260            False
     261            sage: G == G
     262            True
     263        """
     264        if not isinstance(other, GaussianMixtureDistribution):
     265            raise ValueError
     266        return cmp(self.__reduce__()[1], other.__reduce__()[1])
     267
     268    def __len__(self):
     269        """
     270        Return the number of components of this GaussianMixtureDistribution.
     271
     272        EXAMPLES::
     273       
     274            sage: len(hmm.GaussianMixtureDistribution([(.2,-10,.5),(.6,1,1),(.2,20,.5)]))
     275            3
     276        """
     277        return self.c0._length
     278
     279    cpdef is_fixed(self, i=None):
     280        """
     281        Return whether or not this GaussianMixtureDistribution is
     282        fixed when using Baum-Welch to update the corresponding HMM.
     283
     284        INPUT:
     285
     286            - i - None (default) or integer; if given, only return
     287              whether the i-th component is fixed
     288
     289        EXAMPLES::
     290
     291            sage: P = hmm.GaussianMixtureDistribution([(.2,-10,.5),(.6,1,1),(.2,20,.5)])
     292            sage: P.is_fixed()
     293            False
     294            sage: P.is_fixed(0)
     295            False
     296            sage: P.fix(0); P.is_fixed()
     297            False
     298            sage: P.is_fixed(0)
     299            True
     300            sage: P.fix(); P.is_fixed()
     301            True
     302        """
     303        if i is None:
     304            return bool(self.fixed.prod())
     305        else:
     306            return bool(self.fixed[i])
     307
     308    def fix(self, i=None):
     309        """
     310        Set that this GaussianMixtureDistribution (or its ith
     311        component) is fixed when using Baum-Welch to update
     312        the corresponding HMM.
     313
     314        INPUT:
     315
     316            - i - None (default) or integer; if given, only fix the
     317              i-th component
     318
     319        EXAMPLES::
     320       
     321            sage: P = hmm.GaussianMixtureDistribution([(.2,-10,.5),(.6,1,1),(.2,20,.5)])
     322            sage: P.fix(1); P.is_fixed()
     323            False
     324            sage: P.is_fixed(1)
     325            True
     326            sage: P.fix(); P.is_fixed()
     327            True
     328        """
     329        cdef int j
     330        if i is None:
     331            for j in range(self.c0._length):
     332                self.fixed[j] = 1
     333        else:
     334            self.fixed[i] = 1
     335
     336    def unfix(self, i=None):
     337        """
     338        Set that this GaussianMixtureDistribution (or its ith
     339        component) is not fixed when using Baum-Welch to update the
     340        corresponding HMM.
     341
     342        INPUT:
     343
     344            - i - None (default) or integer; if given, only fix the
     345              i-th component
     346
     347        EXAMPLES::
     348       
     349            sage: P = hmm.GaussianMixtureDistribution([(.2,-10,.5),(.6,1,1),(.2,20,.5)])
     350            sage: P.fix(1); P.is_fixed(1)
     351            True
     352            sage: P.unfix(1); P.is_fixed(1)
     353            False
     354            sage: P.fix(); P.is_fixed()
     355            True
     356            sage: P.unfix(); P.is_fixed()
     357            False
     358
     359        """
     360        cdef int j
     361        if i is None:
     362            for j in range(self.c0._length):
     363                self.fixed[j] = 0
     364        else:
     365            self.fixed[i] = 0
     366
     367
     368    def __repr__(self):
     369        """
     370        Return string representation of this mixed Gaussian distribution.
     371
     372        EXAMPLES::
     373
     374            sage: hmm.GaussianMixtureDistribution([(.2,-10,.5),(.6,1,1),(.2,20,.5)]).__repr__()
     375            '0.2*N(-10.0,0.5) + 0.6*N(1.0,1.0) + 0.2*N(20.0,0.5)'       
     376        """
     377        return ' + '.join(["%s*N(%s,%s)"%x for x in self])
     378
     379    def sample(self, n=None):
     380        """
     381        Return a single sample from this distribution (by default), or
     382        if n>1, return a TimeSeries of samples.
     383       
     384        INPUT:
     385
     386            - n -- integer or None (default: None)
     387
     388        OUTPUT:
     389
     390            - float if n is None (default); otherwise a TimeSeries
     391           
     392        EXAMPLES::
     393
     394            sage: P = hmm.GaussianMixtureDistribution([(.2,-10,.5),(.6,1,1),(.2,20,.5)])
     395            sage: P.sample()
     396            19.658243610875129
     397            sage: P.sample(1)
     398            [-10.4683]
     399            sage: P.sample(5)
     400            [-0.1688, -10.3479, 1.6812, 20.1083, -9.9801]
     401            sage: P.sample(0)
     402            []
     403            sage: P.sample(-3)
     404            Traceback (most recent call last):
     405            ...
     406            ValueError: n must be nonnegative
     407        """
     408        cdef randstate rstate = current_randstate()
     409        cdef Py_ssize_t i
     410        cdef TimeSeries T
     411        if n is None:
     412            return self._sample(rstate)
     413        else:
     414            _n = n
     415            if _n < 0:
     416                raise ValueError, "n must be nonnegative"
     417            T = TimeSeries(_n)
     418            for i in range(_n):
     419                T._values[i] = self._sample(rstate)
     420            return T
     421
     422    cdef double _sample(self, randstate rstate):
     423        """
     424        Used internally to compute a sample from this distribution quickly.
     425
     426        INPUT:
     427
     428            - rstate -- a randstate object
     429
     430        OUTPUT:
     431
     432            - double
     433        """
     434        cdef double accum, r
     435        cdef int n
     436        accum = 0
     437        r = rstate.c_rand_double()
     438       
     439        # See the remark in hmm.pyx about using GSL to remove this
     440        # silly way of sampling from a discrete distribution.
     441        for n in range(self.c0._length):
     442            accum += self.param._values[3*n]
     443            if r <= accum:
     444                return random_normal(self.param._values[3*n+1], self.param._values[3*n+2], rstate)
     445        raise RuntimeError, "invalid probability distribution"
     446
     447    cpdef double prob(self, double x):
     448        """
     449        Return the probability of x.
     450
     451        Since this is a continuous distribution, this is defined to be
     452        the limit of the p's such that the probability of [x,x+h] is p*h.
     453
     454        INPUT:
     455
     456            - x -- float
     457
     458        OUTPUT:
     459
     460            - float
     461
     462        EXAMPLES::
     463
     464            sage: P = hmm.GaussianMixtureDistribution([(.2,-10,.5),(.6,1,1),(.2,20,.5)])
     465            sage: P.prob(.5)
     466            0.21123919605857971
     467            sage: P.prob(-100)
     468            0.0
     469            sage: P.prob(20)
     470            0.1595769121605731
     471        """
     472        # The tricky-looking code below is a fast version of this:
     473        #       return sum([c/(sqrt(2*math.pi)*std) * \
     474        #                  exp(-(x-mean)*(x-mean)/(2*std*std)) for
     475        #                  c, mean, std in self.B])
     476        cdef double s=0, mu
     477        cdef int n
     478        for n in range(self.c0._length):
     479            mu = self.param._values[3*n+1]
     480            s += self.c0._values[n]*exp((x-mu)*(x-mu)*self.c1._values[n])
     481        return s
     482
     483    cpdef double prob_m(self, double x, int m):
     484        """
     485        Return the probability of x using just the m-th summand.
     486
     487        INPUT:
     488
     489            - x -- float
     490            - m -- integer
     491
     492        OUTPUT:
     493
     494            - float
     495
     496        EXAMPLES::
     497
     498            sage: P = hmm.GaussianMixtureDistribution([(.2,-10,.5),(.6,1,1),(.2,20,.5)])
     499            sage: P.prob_m(.5, 0)
     500            2.760811768050888...e-97
     501            sage: P.prob_m(.5, 1)
     502            0.21123919605857971
     503            sage: P.prob_m(.5, 2)
     504            0.0
     505        """
     506        cdef double s, mu
     507        if m < 0 or m >= self.param._length//3:
     508            raise IndexError, "index out of range"
     509        mu = self.param._values[3*m+1]
     510        return self.c0._values[m]*exp((x-mu)*(x-mu)*self.c1._values[m])
     511
     512def unpickle_gaussian_mixture_distribution_v1(TimeSeries c0, TimeSeries c1,
     513                                              TimeSeries param, IntList fixed):
     514    """
     515    Used in unpickling GaussianMixtureDistribution's.
     516
     517    EXAMPLES::
     518
     519        sage: P = hmm.GaussianMixtureDistribution([(.2,-10,.5),(.6,1,1),(.2,20,.5)])
     520        sage: loads(dumps(P)) == P          # indirect doctest
     521        True
     522    """
     523    cdef GaussianMixtureDistribution G = PY_NEW(GaussianMixtureDistribution)
     524    G.c0 = c0
     525    G.c1 = c1
     526    G.param = param
     527    G.fixed = fixed
     528    return G

  • sage/stats/hmm/hmm.pxd 

    diff --git a/sage/stats/hmm/hmm.pxd b/sage/stats/hmm/hmm.pxd
    a b  
    11#############################################################################
    2 #       Copyright (C) 2008 William Stein <wstein@gmail.com>
    3 #  Distributed under the terms of the GNU General Public License (GPL),
    4 #  version 2 or any later version at your option.
     2#       Copyright (C) 2010 William Stein <wstein@gmail.com>
     3#  Distributed under the terms of the GNU General Public License (GPL) v2+.
    54#  The full text of the GPL is available at:
    65#                  http://www.gnu.org/licenses/
    76#############################################################################
    87
    9 from sage.matrix.matrix_real_double_dense cimport Matrix_real_double_dense
    108
    11 cdef extern from "ghmm/ghmm.h":
    12     cdef int GHMM_kSilentStates
    13     cdef int GHMM_kHigherOrderEmissions
    14     cdef int GHMM_kDiscreteHMM
    15     cdef int GHMM_EPS_ITER_BW
    16     cdef int GHMM_MAX_ITER_BW
    17 
    18 cdef extern from "ghmm/model.h":
    19     ctypedef struct ghmm_dstate:
    20         # Initial probability
    21         double pi
    22         # Output probability
    23         double *b
    24 
    25         # ID's of the following states
    26         int *out_id
    27         # ID's of the previous states
    28         int *in_id
    29 
    30         # Transition probabilities to successor states.
    31         double *out_a
    32         # Transition probabilities from predecessor states.
    33         double *in_a
    34 
    35         # Number of successor states
    36         int out_states
    37         # Number of precursor states
    38         int in_states
    39         # if fix == 1 then output probabilities b stay fixed during the training
    40         int fix
    41         # Contains a description of the state (null terminated utf-8)
    42         char * desc
    43         # x coordinate position for graph representation plotting
    44         int xPosition
    45         # y coordinate position for graph representation plotting
    46         int yPosition
    47 
    48     ctypedef struct ghmm_dbackground:
    49         pass
    50 
    51     ctypedef struct ghmm_alphabet:
    52         pass
    53    
    54     # Discrete HMM's.
    55     ctypedef struct ghmm_dmodel:
    56         # Number of states
    57         int N   
    58         # Number of outputs
    59         int M   
    60         # Vector of states
    61         ghmm_dstate *s 
    62         # Contains bit flags for varios model extensions such as
    63         # kSilentStates, kTiedEmissions (see ghmm.h for a complete list)
    64         int model_type
    65         # The a priori probability for the model.
    66         # A value of -1 indicates that no prior is defined.
    67         # Note: this is not to be confused with priors on emission distributions.
    68         double prior
    69         # An arbitrary name for the model (null terminated utf-8)
    70         char * name
    71 
    72         # Flag variables for each state indicating whether it is emitting or not.
    73         # Note: silent != NULL iff (model_type & kSilentStates) == 1
    74         int *silent
    75        
    76         # Int variable for the maximum level of higher order emissions.
    77         int maxorder
    78        
    79         # saves the history of emissions as int,
    80         # the nth-last emission is (emission_history * |alphabet|^n+1) % |alphabet|
    81         int emission_history
    82 
    83         # Flag variables for each state indicating whether the states emissions
    84         # are tied to another state. Groups of tied states are represented
    85         # by their tie group leader (the lowest numbered member of the group).
    86         # tied_to[s] == kUntied  : s is not a tied state
    87         # tied_to[s] == s        : s is a tie group leader
    88         # tied_to[t] == s        : t is tied to state s (t>s)
    89         # Note: tied_to != NULL iff (model_type & kTiedEmissions) != 0  */
    90         int *tied_to
    91 
    92         # Note: State store order information of the emissions.
    93         # Classical HMMS have emission order 0; that is, the emission probability
    94         # is conditioned only on the state emitting the symbol.
    95         # For higher order emissions, the emission are conditioned on the state s
    96         # as well as the previous emission_order observed symbols.
    97         # The emissions are stored in the state's usual double* b.
    98         # Note: order != NULL iff (model_type & kHigherOrderEmissions) != 0  */
    99         int * order
    100        
    101         # ghmm_dbackground is a pointer to a
    102         # ghmm_dbackground structure, which holds (essentially) an
    103         # array of background distributions (which are just vectors of floating
    104         # point numbers like state.b).
    105         # For each state the array background_id indicates which of the background
    106         # distributions to use in parameter estimation. A value of kNoBackgroundDistribution
    107         # indicates that none should be used.
    108         # Note: background_id != NULL iff (model_type & kHasBackgroundDistributions) != 0
    109         int *background_id
    110         ghmm_dbackground *bp
    111 
    112         # Topological ordering of silent states
    113         #  Condition: topo_order != NULL iff (model_type & kSilentStates) != 0
    114         int *topo_order
    115         int topo_order_length
    116 
    117         # pow_lookup is a array of precomputed powers
    118         # It contains in the i-th entry M (alphabet size) to the power of i
    119         # The last entry is maxorder+1.
    120         int *pow_lookup
    121 
    122         # Store for each state a class label. Limits the possibly state sequence
    123         # Note: label != NULL iff (model_type & kLabeledStates) != 0
    124         int* label
    125         ghmm_alphabet* label_alphabet
    126         ghmm_alphabet* alphabet
    127 
    128 
    129     # Discrete sequences
    130     ctypedef struct ghmm_dseq:
    131         # sequence array. sequence[i] [j] = j-th symbol of i-th seq
    132         int **seq
    133         # matrix of state ids, can be used to save the viterbi path during sequence generation.
    134         # ATTENTION: is NOT allocated by ghmm_dseq_calloc
    135         int **states
    136         # array of sequence length
    137         int *seq_len
    138         # array of state path lengths
    139         int *states_len
    140         # array of sequence IDs
    141         double *seq_id
    142         # positive sequence weights.  default is 1 = no weight
    143         double *seq_w
    144         # total number of sequences
    145         long seq_number
    146         # reserved space for sequences is always >= seq_number
    147         long capacity
    148         # sum of sequence weights
    149         double total_w
    150         # matrix of state labels corresponding to seq
    151         int **state_labels
    152         # number of labels for each sequence
    153         int *state_labels_len
    154         # internal flags
    155         unsigned int flags
    156 
    157     ghmm_dseq *ghmm_dmodel_label_generate_sequences(
    158         ghmm_dmodel * mo, int seed, int global_len, long seq_number, int Tmax)
    159     ghmm_dseq *ghmm_dmodel_generate_sequences(ghmm_dmodel* mo, int seed, int global_len,
    160                                           long seq_number, int Tmax)
    161     int ghmm_dseq_free (ghmm_dseq ** sq)
    162     ghmm_dseq *ghmm_dseq_calloc (long seq_number)
    163 
    164     int ghmm_dmodel_normalize(ghmm_dmodel *m)
    165     int ghmm_dmodel_free(ghmm_dmodel **m)
    166 
    167     # Problem 1: Likelihood
    168     int ghmm_dmodel_logp(ghmm_dmodel *m, int *O, int len, double *log_p)
    169 
    170     # Problem 2: Decoding
    171     int *ghmm_dmodel_viterbi (ghmm_dmodel *m, int *O, int len, int *pathlen, double *log_p)
    172 
    173     # Problem 3: Learning
    174     int ghmm_dmodel_baum_welch (ghmm_dmodel *m, ghmm_dseq *sq)
    175     int ghmm_dmodel_baum_welch_nstep (ghmm_dmodel *m, ghmm_dseq *sq,
    176                                 int nsteps, double log_likelihood_cutoff)
     9from sage.finance.time_series cimport TimeSeries
    17710
    17811cdef class HiddenMarkovModel:
    179     cdef Matrix_real_double_dense A, B
    180     cdef list pi
     12    cdef int N
     13    cdef TimeSeries A, pi
    18114
    182 cdef class DiscreteHiddenMarkovModel(HiddenMarkovModel):
    183     cdef ghmm_dmodel* m
    184     cdef bint initialized
    185     cdef object _emission_symbols, _emission_symbols_dict

  • sage/stats/hmm/hmm.pyx 

    diff --git a/sage/stats/hmm/hmm.pyx b/sage/stats/hmm/hmm.pyx
    a b  
    1 r"""nodoctest
     1"""
    22Hidden Markov Models
    33
    4 AUTHOR: William Stein
     4This is a complete pure-Cython optimized implementation of Hidden
     5Markov Models.  It fully supports Discrete, Gaussian, and Mixed
     6Gaussian emissions.
     7
     8The best references for the basic HMM algorithms implemented here are:
     9
     10   -  Tapas Kanungo's "Hidden Markov Models"
     11
     12   -  Jackson's HMM tutorial:
     13        http://personal.ee.surrey.ac.uk/Personal/P.Jackson/tutorial/
     14
     15LICENSE: Some of the code in this file is based on reading Kanungo's
     16GPLv2+ implementation of discrete HMM's, hence the present code must
     17be licensed with a GPLv2+ compatible license.
     18
     19AUTHOR:
     20
     21   - William Stein, 2010-03
    522"""
    623
    724#############################################################################
    8 #       Copyright (C) 2008 William Stein <wstein@gmail.com>
    9 #  Distributed under the terms of the GNU General Public License (GPL),
    10 #  version 2 or any later version at your option.
     25#       Copyright (C) 2010 William Stein <wstein@gmail.com>
     26#  Distributed under the terms of the GNU General Public License (GPL) v2+.
    1127#  The full text of the GPL is available at:
    1228#                  http://www.gnu.org/licenses/
    1329#############################################################################
    1430
    15 import math
     31include "../../ext/stdsage.pxi"
     32include "../../ext/cdefs.pxi"
     33include "../../ext/interrupt.pxi"
    1634
     35cdef extern from "math.h":
     36    double log(double)
     37
     38from sage.finance.time_series cimport TimeSeries
     39from sage.stats.intlist cimport IntList
    1740from sage.matrix.all import is_Matrix, matrix
    18 from sage.rings.all  import RDF
    19 from sage.misc.randstate import random
     41from sage.misc.randstate cimport current_randstate, randstate
    2042
    21 include "../../ext/stdsage.pxi"
     43from util cimport HMM_Util
    2244
    23 include "misc.pxi"
     45cdef HMM_Util util = HMM_Util()
     46
     47###########################################
    2448
    2549cdef class HiddenMarkovModel:
    2650    """
    2751    Abstract base class for all Hidden Markov Models.
     52    """
     53    def initial_probabilities(self):
     54        """
     55        Return the initial probabilities, which as a TimeSeries of
     56        length N, where N is the number of states of the Markov model.
    2857
    29     INPUT:
    30         A -- matrix or list
    31         B -- matrix or list
    32         pi -- list of floats
     58        EXAMPLES::
    3359
    34     EXAMPLES:
    35     One shouldn't directly call this constructor since this class is an abstract
    36     base class. 
    37         sage: sage.stats.hmm.hmm.HiddenMarkovModel([[0.4,0.6],[1,0]], [[1,0],[0.5,0.5]], [0.5,0.5])
    38         <sage.stats.hmm.hmm.HiddenMarkovModel object at ...>
    39     """
    40     def __init__(self, A, B, pi=None):
     60            sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.1,0.9],[0.5,0.5]], [.2,.8])
     61            sage: pi = m.initial_probabilities(); pi
     62            [0.2000, 0.8000]
     63            sage: type(pi)
     64            <type 'sage.finance.time_series.TimeSeries'>
     65
     66        The returned time series is a copy, so changing it does not
     67        change the model.
     68
     69            sage: pi[0] = .1; pi[1] = .9
     70            sage: m.initial_probabilities()
     71            [0.2000, 0.8000]
     72
     73        Some other models::
     74       
     75            sage: hmm.GaussianHiddenMarkovModel([[.1,.9],[.5,.5]], [(1,1), (-1,1)], [.1,.9]).initial_probabilities()
     76            [0.1000, 0.9000]
     77            sage: hmm.GaussianMixtureHiddenMarkovModel([[.9,.1],[.4,.6]], [[(.4,(0,1)), (.6,(1,0.1))],[(1,(0,1))]], [.7,.3]).initial_probabilities()
     78            [0.7000, 0.3000]
    4179        """
    42         INPUT:
    43             A -- matrix or list
    44             B -- matrix or list
    45             pi -- list of floats
     80        return TimeSeries(self.pi)
    4681
    47         EXAMPLES:
    48             sage: hmm.DiscreteHiddenMarkovModel([[1]], [[0.0,1.0]], [1])
    49             Discrete Hidden Markov Model with 1 States and 2 Emissions
    50             Transition matrix:
    51             [1.0]
    52             Emission matrix:
    53             [0.0 1.0]
    54             Initial probabilities: [1.0]
     82    def transition_matrix(self):
    5583        """
    56         cdef Py_ssize_t n
    57        
    58         # Convert A to a matrix
    59         if not is_Matrix(A):
    60             n = len(A)
    61             A = matrix(RDF, n, len(A[0]) if n > 0 else 0, A)
     84        Return the state transition matrix.
    6285
    63         # Convert B to a matrix           
    64         if not is_Matrix(B):
    65             n = len(B)
    66             B = matrix(RDF, n, len(B[0]) if n > 0 else 0, B)
     86        OUTPUT:
    6787
    68         # Some consistency checks
    69         if not A.is_square():
    70             print A.parent()
    71             raise ValueError, "A must be square"
    72         if A.nrows() != B.nrows():
    73             raise ValueError, "number of rows of A and B must be the same"
     88            - a Sage matrix with real double precision (RDF) entries.
    7489
    75         # Make sure A and B are over RDF.
    76         if A.base_ring() != RDF:
    77             A = A.change_ring(RDF)
    78         if B.base_ring() != RDF:
    79             B = B.change_ring(RDF)
     90        EXAMPLES::
    8091
    81         # Make sure the initial probabilities are all floats.
    82         if pi is None:
    83             if A.nrows() == 0:
    84                 self.pi = []
    85             else:
    86                 self.pi = [1.0/A.nrows()]*A.nrows()
    87         else:
    88             self.pi = [float(x) for x in pi]
    89             if len(self.pi) != A.nrows():
    90                 raise ValueError, "length of pi must equal number of rows of A"
     92            sage: M = hmm.DiscreteHiddenMarkovModel([[0.7,0.3],[0.9,0.1]], [[0.5,.5],[.1,.9]], [0.3,0.7])
     93            sage: T = M.transition_matrix(); T
     94            [0.7 0.3]
     95            [0.9 0.1]
    9196
    92         # Record the now validated matrices A and B as attributes.
    93         # They get used later as attributes in the constructors for
    94         # derived classes.
    95         self.A = A
    96         self.B = B
     97        The returned matrix is mutable, but changing it does not
     98        change the transition matrix for the model::
    9799
    98         # Check that all entries of A are nonnegative and raise a
    99         # ValueError otherwise. Note that we don't check B since it
    100         # has entries that are potentially negative in the continuous
    101         # case.  But GHMM prints clear warnings when the emission
    102         # probabilities are negative, i.e., it does not silently give
    103         # wrong results like it does for negative transition
    104         # probabilities.
    105         cdef Py_ssize_t i, j
    106         for i from 0 <= i < self.A._nrows:
    107             for j from 0 <= j < self.A._ncols:
    108                 if self.A.get_unsafe_double(i,j) < 0:
    109                     raise ValueError, "each transition probability must be nonnegative"
    110                
     100            sage: T[0,0] = .1; T[0,1] = .9
     101            sage: M.transition_matrix()
     102            [0.7 0.3]
     103            [0.9 0.1]
    111104
     105        Transition matrices for other types of models::
    112106
    113 cdef class DiscreteHiddenMarkovModel(HiddenMarkovModel):
    114     """
    115     Create a discrete hidden Markov model.
     107            sage: hmm.GaussianHiddenMarkovModel([[.1,.9],[.5,.5]], [(1,1), (-1,1)], [.5,.5]).transition_matrix()
     108            [0.1 0.9]
     109            [0.5 0.5]
     110            sage: hmm.GaussianMixtureHiddenMarkovModel([[.9,.1],[.4,.6]], [[(.4,(0,1)), (.6,(1,0.1))],[(1,(0,1))]], [.7,.3]).transition_matrix()
     111            [0.9 0.1]
     112            [0.4 0.6]
     113        """
     114        from sage.matrix.constructor import matrix
     115        from sage.rings.all import RDF
     116        return matrix(RDF, self.N, self.A.list())
    116117
    117     hmm.DiscreteHiddenMarkovModel(A, B, pi=None, emission_symbols=None, name=None, normalize=True)
    118 n   
    119     INPUTS:
    120         A  -- square matrix of doubles; the state change probabilities
    121         B  -- matrix of doubles; emission probabilities
    122         pi -- list of floats; probabilities for each initial state
    123         emission_symbols -- list of B.ncols() symbols (just used for printing)
    124         name -- (optional) name of the model
    125         normalize -- (optional; default=True) whether or not to normalize
    126                      the model so the probabilities add to 1
    127 
    128     EXAMPLES:
    129     We create a discrete HMM with 2 internal states on an alphabet of size 2.
    130         sage: hmm.DiscreteHiddenMarkovModel([[0.2,0.8],[0.5,0.5]], [[1,0],[0,1]], [0,1])
    131         Discrete Hidden Markov Model with 2 States and 2 Emissions
    132         Transition matrix:
    133         [0.2 0.8]
    134         [0.5 0.5]
    135         Emission matrix:
    136         [1.0 0.0]
    137         [0.0 1.0]
    138         Initial probabilities: [0.0, 1.0]
    139 
    140     The transition probabilities must be nonnegative:
    141         sage: hmm.DiscreteHiddenMarkovModel([[-0.2,0.8],[0.5,0.5]], [[1,0],[0,1]], [0,1])
    142         Traceback (most recent call last):
    143         ...
    144         ValueError: each transition probability must be nonnegative
    145 
    146     The transition probabilities are by default automatically normalized:
    147         sage: a = hmm.DiscreteHiddenMarkovModel([[0.2,0.3],[0.5,0.5]], [[1,0],[0,1]], [0,1])
    148         sage: a.transition_matrix()
    149         [0.4 0.6]
    150         [0.5 0.5]
    151     """
    152     def __init__(self, A, B, pi=None, emission_symbols=None, name=None, normalize=True):
     118    def graph(self, eps=1e-3):
    153119        """
    154         EXAMPLES:
    155             sage: hmm.DiscreteHiddenMarkovModel([[1]], [[0.0,1.0]], [1])
    156             Discrete Hidden Markov Model with 1 States and 2 Emissions
    157             Transition matrix:
    158             [1.0]
    159             Emission matrix:
    160             [0.0 1.0]
    161             Initial probabilities: [1.0]
    162         """
    163         # memory has not all been setup yet.
    164         self.initialized = False 
    165 
    166         # This sets self.A, self.B and pi after doing appropriate coercions, etc.
    167         HiddenMarkovModel.__init__(self, A, B, pi)
    168 
    169         self.set_emission_symbols(emission_symbols)
    170 
    171         self.m = <ghmm_dmodel*> safe_malloc(sizeof(ghmm_dmodel))
    172 
    173         self.m.label = to_int_array(range(len(self._emission_symbols)))
    174        
    175         # Set all pointers to NULL
    176         self.m.s = NULL; self.m.name = NULL; self.m.silent = NULL
    177         self.m.tied_to = NULL; self.m.order = NULL; self.m.background_id = NULL
    178         self.m.bp = NULL; self.m.topo_order = NULL; self.m.pow_lookup = NULL;
    179         self.m.label_alphabet = NULL; self.m.alphabet = NULL
    180 
    181         # Set number of states and number of outputs
    182         self.m.N = self.A.nrows()
    183         self.m.M = len(self._emission_symbols)
    184         # Set the model type to discrete
    185         self.m.model_type = GHMM_kDiscreteHMM
    186 
    187         # Set that no a prior model probabilities are set.
    188         self.m.prior = -1
    189 
    190         # Assign model identifier if specified
    191         if name is not None:
    192             name = str(name)
    193             self.m.name = <char*> safe_malloc(len(name)+1)
    194             strcpy(self.m.name, name)
    195         else:
    196             self.m.name = NULL
    197 
    198         # Allocate and initialize states
    199         cdef ghmm_dstate* states = <ghmm_dstate*> safe_malloc(sizeof(ghmm_dstate) * self.m.N)
    200         cdef ghmm_dstate* state
    201 
    202         silent_states = []
    203         tmp_order     = []
    204        
    205         cdef Py_ssize_t i, j, k
    206        
    207         for i from 0 <= i < self.m.N:
    208             v = self.B[i]
    209 
    210             # Get a reference to the i-th state for convenience of the notation below.
    211             state = &(states[i])
    212             state.desc = NULL
    213 
    214             # Compute state order
    215             if self.m.M > 1:
    216                 order = math.log( len(v), self.m.M ) - 1
    217             else:
    218                 order = len(v) - 1
    219 
    220             # Check for valid number of emission parameters
    221             order = int(order)
    222             if self.m.M**(order+1) == len(v):
    223                 tmp_order.append(order)
    224             else:
    225                 raise ValueError, "number of columns (= %s) of B must be a power of the number of emission symbols (= %s)"%(
    226                     self.B.ncols(), len(emission_symbols))
    227 
    228             state.b = to_double_array(v)
    229             state.pi = self.pi[i]
    230 
    231             silent_states.append( 1 if sum(v) == 0 else 0)
    232 
    233             # Set "out" probabilities, i.e., the probabilities to
    234             # transition to another hidden state from this state.
    235             v = self.A[i]
    236             k = self.m.N
    237             state.out_states = k
    238             state.out_id = <int*> safe_malloc(sizeof(int)*k)
    239             state.out_a  = <double*> safe_malloc(sizeof(double)*k)
    240             for j from 0 <= j < k:
    241                 state.out_id[j] = j
    242                 state.out_a[j]  = v[j]
    243 
    244             # Set "in" probabilities
    245             v = self.A.column(i)
    246             state.in_states = k
    247             state.in_id = <int*> safe_malloc(sizeof(int)*k)
    248             state.in_a  = <double*> safe_malloc(sizeof(double)*k)
    249             for j from 0 <= j < k:
    250                 state.in_id[j] = j
    251                 state.in_a[j]  = v[j]
    252                
    253             state.fix = 0
    254                
    255         self.m.s = states
    256         self.initialized = True
    257         if normalize:
    258             self.normalize()
    259 
    260     def __cmp__(self, other):
    261         """
    262         Compare two Discrete HMM's.
     120        Create a weighted directed graph from the transition matrix,
     121        not including any edge with a probability less than eps.
    263122
    264123        INPUT:
    265             self, other -- objects; if other is not a Discrete HMM compare types.
    266         OUTPUT:
    267             -1,0,1
    268124
    269         The transition matrices are compared, then the emission
    270         parameters, and the initial probabilities.  The names are not
    271         compared, so two GHMM's with the same parameters but different
    272         names compare to be equal.
    273 
    274         EXAMPLES:
    275             sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.0,1.0],[1,1]], [1,2], normalize=False)
    276             sage: m.__cmp__(m)
    277             0
    278 
    279         Note that the name doesn't matter:
    280             sage: n = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.0,1.0],[1,1]], [1,2], normalize=False)
    281             sage: m.__cmp__(n)
    282             0
    283 
    284         Normalizing fixes the initial probabilities, hence m and n are no longer equal.
    285             sage: n.normalize()
    286             sage: m.__cmp__(n)
    287             1
    288         """
    289         if not isinstance(other, DiscreteHiddenMarkovModel):
    290             return cmp(type(self), type(other))
    291 
    292         if self is other: return 0  # easy special case
    293        
    294         cdef DiscreteHiddenMarkovModel o = other
    295         if self.m.N < o.m.N:
    296             return -1
    297         elif self.m.N > o.m.N:
    298             return 1
    299         cdef Py_ssize_t i, j
    300 
    301         # This code is similar to code in chmm.pyx, but with several small differences.
    302        
    303         # Compare transition matrices
    304         for i from 0 <= i < self.m.N:
    305             for j from 0 <= j < self.m.s[i].out_states:
    306                 if self.m.s[i].out_a[j] < o.m.s[i].out_a[j]:
    307                     return -1
    308                 elif self.m.s[i].out_a[j] > o.m.s[i].out_a[j]:
    309                     return 1
    310 
    311         # Compare emission matrices
    312         for i from 0 <= i < self.m.N:
    313             for j from 0 <= j < self.B._ncols:
    314                 if self.m.s[i].b[j] < o.m.s[i].b[j]:
    315                     return -1
    316                 elif self.m.s[i].b[j] > o.m.s[i].b[j]:
    317                     return 1
    318                
    319         # Compare initial state probabilities
    320         for 0 <= i < self.m.N:
    321             if self.m.s[i].pi < o.m.s[i].pi:
    322                 return -1
    323             elif self.m.s[i].pi > o.m.s[i].pi:
    324                 return 1
    325 
    326         # Compare emission symbols
    327         return cmp(self._emission_symbols, o._emission_symbols)
    328 
    329     def __reduce__(self):
    330         """
    331         Used in pickling.
    332 
    333         EXAMPLES:
    334             sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.0,1.0],[1,1]], [0,1], ['a','b'], name='test model')
    335             sage: f,g = m.__reduce__()
    336             sage: f(*g) == m
    337             True
    338         """
    339         return unpickle_discrete_hmm_v0, (self.transition_matrix(), self.emission_matrix(),
    340                       self.initial_probabilities(), self._emission_symbols, self.name())
    341 
    342     def __dealloc__(self):
    343         """
    344         Deallocate memory allocated by the HMM.
    345 
    346         EXAMPLES:
    347             sage: m = hmm.DiscreteHiddenMarkovModel([[0.2,0.8],[0.5,0.5]], [[1,0],[0,1]], [0,1])  # indirect doctest
    348             sage: del m
    349         """
    350         if self.initialized:
    351             ghmm_dmodel_free(&self.m)
    352 
    353     def fix_emissions(self, Py_ssize_t i, bint fixed=True):
    354         """
    355         Sets the i-th emission parameters to be either fixed or not
    356         fixed.  If it is fixed, then running the Baum-Welch algorithm
    357         will not change the emission parameters for the i-th state.
    358        
    359         INPUT:
    360             i -- nonnegative integer < self.m.N
    361             fixed -- bool
    362 
    363         EXAMPLES:
    364         First without calling fix_emissions:
    365             sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]],[[.5,.5],[.5,.5]])
    366             sage: m.baum_welch([0,0,0,1,1,1])
    367             sage: m.emission_matrix()
    368             [              1.0               0.0]
    369             [3.92881039079e-05    0.999960711896]
    370 
    371         We call fix_emissions on the first state and notice that the first
    372         row of the emission matrix does not change:
    373             sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]],[[.5,.5],[.5,.5]])
    374             sage: m.fix_emissions(0)
    375             sage: m.baum_welch([0,0,0,1,1,1])
    376             sage: m.emission_matrix()
    377             [              0.5               0.5]
    378             [0.000542712675606    0.999457287324]
    379 
    380         We call fix_emissions on the second state and notice that the second
    381         row of the emission matrix does not change:
    382             sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]],[[.5,.5],[.5,.5]])
    383             sage: m.fix_emissions(1)
    384             sage: m.baum_welch([0,0,0,1,1,1])
    385             sage: m.emission_matrix()
    386             [   0.999999904763 9.52366620142e-08]
    387             [              0.5               0.5]
    388 
    389         TESTS:
    390         Make sure that out of range indices are handled correctly with an IndexError.
    391             sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]],[[.5,.5],[.5,.5]])
    392             sage: m.fix_emissions(2)
    393             Traceback (most recent call last):
    394             ...
    395             IndexError: index out of range
    396             sage: m.fix_emissions(-1)
    397             Traceback (most recent call last):
    398             ...
    399             IndexError: index out of range
    400         """
    401         if i < 0 or i >= self.m.N:
    402             raise IndexError, "index out of range"
    403         self.m.s[i].fix = fixed
    404 
    405     def __repr__(self):
    406         """
    407         Return string representation of this HMM.
     125            - eps -- nonnegative real number
    408126
    409127        OUTPUT:
    410             string
    411128
    412         EXAMPLES:
    413             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'])
    414             sage: a.__repr__()
    415             "Discrete Hidden Markov Model with 2 States and 2 Emissions\nTransition matrix:\n[0.1 0.9]\n[0.1 0.9]\nEmission matrix:\n[0.9 0.1]\n[0.1 0.9]\nInitial probabilities: [0.5, 0.5]\nEmission symbols: [3/4, 'abc']"
     129            - a digraph
     130
     131        EXAMPLES::
     132
     133            sage: m = hmm.DiscreteHiddenMarkovModel([[.3,0,.7],[0,0,1],[.5,.5,0]], [[.5,.5,.2]]*3, [1/3]*3)
     134            sage: G = m.graph(); G
     135            Looped multi-digraph on 3 vertices
     136            sage: G.edges()
     137            [(0, 0, 0.3), (0, 2, 0.7), (1, 2, 1.0), (2, 0, 0.5), (2, 1, 0.5)]
     138            sage: G.plot()
    416139        """
    417         s = "Discrete Hidden Markov Model%s with %s States and %s Emissions"%(
    418             ' ' + self.m.name if self.m.name else '',
    419             self.m.N, self.m.M)
    420         s += '\nTransition matrix:\n%s'%self.transition_matrix()
    421         s += '\nEmission matrix:\n%s'%self.emission_matrix()
    422         s += '\nInitial probabilities: %s'%self.initial_probabilities()
    423         if self._emission_symbols_dict:
    424             s += '\nEmission symbols: %s'%self._emission_symbols
    425         return s
     140        cdef int i, j
     141        m = self.transition_matrix()
     142        for i in range(self.N):
     143            for j in range(self.N):
     144                if m[i,j] < eps:
     145                    m[i,j] = 0
     146        from sage.graphs.all import DiGraph
     147        return DiGraph(m, weighted=True)
    426148
    427     def name(self):
    428         """
    429         Return the name of this model.
    430 
    431         OUTPUT:
    432             string or None
    433 
    434         EXAMPLES:
    435             sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.0,1.0],[1,1]], [1,2], name='test model')
    436             sage: m.name()
    437             'test model'
    438 
    439         If the model is not explicitly named then this function returns None:
    440             sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.0,1.0],[1,1]], [1,2])
    441             sage: m.name() is None
    442             True
    443         """
    444         if self.m.name:
    445             s = str(self.m.name)
    446             return s
    447         else:
    448             return None
    449 
    450     def initial_probabilities(self):
    451         """
    452         Return the list of initial state probabilities.
    453        
    454         OUTPUT:
    455             list of floats
    456 
    457         EXAMPLES:
    458             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])
    459             sage: a.initial_probabilities()
    460             [0.5, 0.5]
    461         """
    462         cdef Py_ssize_t i
    463         return [self.m.s[i].pi for i in range(self.m.N)]
    464 
    465     def transition_matrix(self, list_only=True):
    466         """
    467         Return the hidden state transition matrix.
    468        
    469         EXAMPLES:
    470             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])
    471             sage: a.transition_matrix()
    472             [0.9 0.1]
    473             [0.9 0.1]
    474         """
    475         cdef Py_ssize_t i, j
    476         for i from 0 <= i < self.m.N:
    477             for j from 0 <= j < self.m.s[i].out_states:
    478                 self.A.set_unsafe_double(i,j,self.m.s[i].out_a[j])
    479         return self.A
    480 
    481     def emission_matrix(self, list_only=True):
    482         """
    483         Return the emission probability matrix.
    484        
    485         EXAMPLES:
    486             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])
    487             sage: a.emission_matrix()
    488             [0.5 0.5 0.0 0.0]
    489             [0.0 0.0 0.5 0.5]
    490         """
    491         cdef Py_ssize_t i, j
    492         for i from 0 <= i < self.m.N:
    493             for j from 0 <= j < self.B._ncols:
    494                 self.B.set_unsafe_double(i,j,self.m.s[i].b[j])
    495         return self.B
    496 
    497     def normalize(self):
    498         """
    499         Normalize the transition and emission probabilities, if applicable.
    500 
    501         EXAMPLES:
    502             sage: a = hmm.DiscreteHiddenMarkovModel([[0.5,1],[1.2,0.9]], [[1,0.5],[0.5,1]], [0.1,1.2])
    503             sage: a.normalize()
    504             sage: a
    505             Discrete Hidden Markov Model with 2 States and 2 Emissions
    506             Transition matrix:
    507             [0.333333333333 0.666666666667]
    508             [0.571428571429 0.428571428571]
    509             Emission matrix:
    510             [0.666666666667 0.333333333333]
    511             [0.333333333333 0.666666666667]
    512             Initial probabilities: [0.076923076923076927, 0.92307692307692302]
    513         """
    514         ghmm_dmodel_normalize(self.m)
    515 
    516     def sample(self, long length, number=None):
     149    def sample(self, Py_ssize_t length, number=None):
    517150        """
    518151        Return number samples from this HMM of given length.
    519152
    … …  
    529162            sage: set_random_seed(0)
    530163            sage: a = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.1,0.9]], [[1,0],[0,1]], [0,1])
    531164            sage: print a.sample(10, 3)
    532             [[1, 0, 1, 1, 0, 1, 1, 0, 1, 0], [1, 1, 1, 1, 1, 1, 1, 1, 1, 0], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1]]
     165            [[1, 0, 1, 1, 1, 1, 0, 1, 1, 1], [1, 1, 0, 0, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 0, 1, 0, 1, 1, 1]]
    533166            sage: a.sample(15)
    534             [1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 1]
     167            [1, 1, 1, 1, 0 ... 1, 1, 1, 1, 1]
     168            sage: a.sample(3, 1)
     169            [[1, 1, 1]]
    535170            sage: list(a.sample(1000)).count(0)
    536             95
     171            88
    537172
    538173        If the emission symbols are set
    539174            sage: set_random_seed(0)
    540175            sage: a = hmm.DiscreteHiddenMarkovModel([[0.5,0.5],[0.1,0.9]], [[1,0],[0,1]], [0,1], ['up', 'down'])
    541176            sage: a.sample(10)
    542             ['down', 'up', 'down', 'down', 'up', 'down', 'down', 'up', 'down', 'up']
     177            ['down', 'up', 'down', 'down', 'down', 'down', 'up', 'up', 'up', 'up']
    543178        """
    544         seed = random()
    545179        if number is None:
    546             number = 1
    547             single = True
    548         else:
    549             single = False
    550         cdef ghmm_dseq *d = ghmm_dmodel_generate_sequences(self.m, seed, length, number, length)
    551         cdef Py_ssize_t i, j
    552         v = [[d.seq[j][i] for i in range(length)] for j in range(number)]
    553         if self._emission_symbols_dict:
    554             w = self._emission_symbols
    555             v = [[w[i] for i in z] for z in v]
    556         if single:
    557             return v[0]
    558         return v
     180            return self.generate_sequence(length)[0]
    559181
    560     def emission_symbols(self):
    561         """
    562         Return a copy of the list of emission symbols of self.
     182        cdef Py_ssize_t i
     183        return [self.generate_sequence(length)[0] for i in range(number)]
    563184
    564         Use set_emission_symbols to set the list of emission symbols.
     185cdef class DiscreteHiddenMarkovModel(HiddenMarkovModel):
     186    """
     187    A discrete Hidden Markov model implemented using double precision
     188    floating point arithmetic.
    565189
    566         EXAMPLES:
    567             sage: a = hmm.DiscreteHiddenMarkovModel([[0.5,0.5],[0.1,0.9]], [[1,0],[0,1]], [0,1], ['up', -3/179])
    568             sage: a.emission_symbols()
    569             ['up', -3/179]
    570         """
    571         return list(self._emission_symbols)
     190    EXAMPLES::
    572191
    573     def set_emission_symbols(self, emission_symbols):
    574         """
    575         Set the list of emission symbols.
    576 
    577         EXAMPLES:
    578             sage: set_random_seed(0)
    579             sage: a = hmm.DiscreteHiddenMarkovModel([[0.5,0.5],[0.1,0.9]], [[1,0],[0,1]], [0,1], ['up', 'down'])
    580             sage: a.set_emission_symbols([3,5])
    581             sage: a.emission_symbols()
    582             [3, 5]
    583             sage: a.sample(10)
    584             [5, 3, 5, 5, 3, 5, 5, 3, 5, 3]
    585             sage: a.set_emission_symbols([pi,5/9+e])
    586             sage: a.sample(10)
    587             [e + 5/9, e + 5/9, e + 5/9, e + 5/9, e + 5/9, e + 5/9, pi, pi, e + 5/9, pi]
    588         """
    589         if emission_symbols is None:
    590             self._emission_symbols = range(self.B.ncols())
    591             self._emission_symbols_dict = None
    592         else:
    593             self._emission_symbols = list(emission_symbols)
    594             if self._emission_symbols != range(self.B.ncols()):
    595                 self._emission_symbols_dict = dict([(x,i) for i, x in enumerate(emission_symbols)])
    596 
    597 
    598     ####################################################################
    599     # HMM Problem 1 -- Computing likelihood: Given the parameter set
    600     # lambda of an HMM model and an observation sequence O, determine
    601     # the likelihood P(O | lambda).
    602     ####################################################################
    603     def log_likelihood(self, seq):
    604         r"""
    605         HMM Problem 1: Likelihood. Return $\log ( P[seq | model] )$,
    606         the log of the probability of seeing the given sequence given
    607         this model, using the forward algorithm and assuming
    608         independence of the sequence seq.
    609 
    610         INPUT:
    611             seq -- a list; sequence of observed emissions of the HMM
    612 
    613         OUTPUT:
    614             float -- the log of the probability of seeing this sequence
    615                      given the model
    616 
    617         WARNING: By convention we return -inf for 0 probability events.
    618        
    619         EXAMPLES:
    620             sage: a = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.1,0.9]], [[1,0],[0,1]], [0,1])
    621             sage: a.log_likelihood([1,1])
    622             -0.10536051565782635
    623             sage: a.log_likelihood([1,0])
    624             -2.3025850929940459
    625 
    626         Notice that any sequence starting with 0 can't occur, since
    627         the above model always starts in a state that produces 1 with
    628         probability 1.  By convention log(probability 0) is -inf.
    629             sage: a.log_likelihood([0,0])
    630             -inf
    631 
    632         Here's a special case where each sequence is equally probable.
    633             sage: a = hmm.DiscreteHiddenMarkovModel([[0.5,0.5],[0.5,0.5]], [[1,0],[0,1]], [0.5,0.5])
    634             sage: a.log_likelihood([0,0])
    635             -1.3862943611198906
    636             sage: log(0.25)
    637             -1.38629436111989
    638             sage: a.log_likelihood([0,1])
    639             -1.3862943611198906
    640             sage: a.log_likelihood([1,0])
    641             -1.3862943611198906
    642             sage: a.log_likelihood([1,1])
    643             -1.3862943611198906
    644         """
    645         if self._emission_symbols_dict:
    646             seq = [self._emission_symbols_dict[z] for z in seq]
    647         cdef double log_p
    648         cdef int* O = to_int_array(seq)
    649         cdef int ret = ghmm_dmodel_logp(self.m, O, len(seq), &log_p)
    650         sage_free(O)
    651         if ret == -1:
    652             # forward returned -1: sequence can't be built
    653             return -float('Inf')
    654         return log_p
    655 
    656     ####################################################################
    657     # HMM Problem 2 -- Decoding: Given the complete parameter set that
    658     # defines the model and an observation sequence seq, determine the
    659     # best hidden sequence Q.  Use the Viterbi algorithm.
    660     ####################################################################   
    661     def viterbi(self, seq):
    662         """
    663         HMM Problem 2: Decoding.  Determine a hidden sequence of
    664         states that is most likely to produce the given sequence seq
    665         of observations.
    666 
    667         INPUT:
    668             seq -- sequence of emitted symbols
    669            
    670         OUTPUT:
    671             list -- a most probable sequence of hidden states, i.e., the
    672                     Viterbi path.
    673             float -- log of the probability that the sequence of hidden
    674                      states actually produced the given sequence seq.
    675 
    676         EXAMPLES:
    677             sage: a = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.1,0.9]], [[0.9,0.1],[0.1,0.9]], [0.5,0.5])
    678             sage: a.viterbi([1,0,0,1,0,0,1,1])
    679             ([1, 0, 0, 1, ..., 0, 1, 1], -11.06245322477221...)
    680            
    681         We predict the state sequence when the emissions are 3/4 and 'abc'.
    682             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'])
    683 
    684         Note that state 0 is common below, despite the model trying hard to
    685         switch to state 1:
    686             sage: a.viterbi([3/4, 'abc', 'abc'] + [3/4]*10)
    687             ([0, 1, 1, 0, 0, ..., 0, 0, 0, 0, 0, 0, 0], -25.299405845367794)
    688         """
    689         if len(seq) == 0:
    690             return [], 0.0
    691         if self._emission_symbols_dict:
    692             seq = [self._emission_symbols_dict[z] for z in seq]
    693         cdef int* path
    694         cdef int* O = to_int_array(seq)
    695         cdef int pathlen
    696         cdef double log_p
    697 
    698         path = ghmm_dmodel_viterbi(self.m, O, len(seq), &pathlen, &log_p)
    699         sage_free(O)
    700         if not path:
    701             raise RuntimeError, "error computing viterbi path"
    702         p = [path[i] for i in range(pathlen)]
    703         sage_free(path)
    704 
    705         return p, log_p
    706 
    707     ####################################################################
    708     # HMM Problem 3 -- Learning: Given an observation sequence O and
    709     # the set of states in the HMM, improve the HMM to increase the
    710     # probability of observing O.
    711     ####################################################################
    712     def baum_welch(self, training_seqs, nsteps=None, log_likelihood_cutoff=None):
    713         """
    714         HMM Problem 3: Learning.  Given an observation sequence O and
    715         the set of states in the HMM, improve the HMM using the
    716         Baum-Welch algorithm to increase the probability of observing O.
    717 
    718         INPUT:
    719             training_seqs -- a list of lists of emission symbols (or a single list)
    720             nsteps -- integer or None (default: None) maximum number
    721                       of Baum-Welch steps to take
    722             log_likehood_cutoff -- positive float or None (default:
    723                       None); the minimal improvement in likelihood
    724                       with respect to the last iteration required to
    725                       continue. Relative value to log likelihood
    726 
    727         OUTPUT:
    728             changes the model in places, or raises a RuntimError
    729             exception on error
    730 
    731         EXAMPLES:
    732         We make a model that has two states and is equally likely to output
    733         0 or 1 in either state and transitions back and forth with equal
    734         probability.
    735             sage: a = hmm.DiscreteHiddenMarkovModel([[0.5,0.5],[0.5,0.5]], [[0.5,0.5],[0.5,0.5]], [0.5,0.5])
    736 
    737         We give the model some training data this much more likely to
    738         be 1 than 0.
    739             sage: a.baum_welch([[1,1,1,1,0,1], [1,0,1,1,1,1]])
    740 
    741         Now the model's emission matrix changes since it is much
    742         more likely to emit 1 than 0. 
    743             sage: a
    744             Discrete Hidden Markov Model with 2 States and 2 Emissions
    745             Transition matrix:
    746             [0.5 0.5]
    747             [0.5 0.5]
    748             Emission matrix:
    749             [0.166666666667 0.833333333333]
    750             [0.166666666667 0.833333333333]
    751             Initial probabilities: [0.5, 0.5]
    752 
    753         Note that 1/6 = 1.666...:
    754             sage: 1.0/6
    755             0.166666666666667
    756 
    757         We run Baum-Welch on a single sequence:
    758             sage: a = hmm.DiscreteHiddenMarkovModel([[0.5,0.5],[0.5,0.5]], [[0.5,0.5],[0.5,0.5]], [0.5,0.5])
    759             sage: a.baum_welch([1,0,1]*10)
    760             sage: a
    761             Discrete Hidden Markov Model with 2 States and 2 Emissions
    762             Transition matrix:
    763             [0.5 0.5]
    764             [0.5 0.5]
    765             Emission matrix:
    766             [0.333333333333 0.666666666667]
    767             [0.333333333333 0.666666666667]
    768             Initial probabilities: [0.5, 0.5]
    769 
    770         We compare using a non-default number of steps and non-default log likelihood cutoff:
    771             sage: h = hmm.DiscreteHiddenMarkovModel([[.1,.9],[.4,.6]], [[1/2]*2]*2, [1/2]*2)
    772             sage: h.baum_welch([1,1,1,1,1,0,1,1,1,1,1,0,0])
    773             sage: h
    774             Discrete Hidden Markov Model with 2 States and 2 Emissions
    775             Transition matrix:
    776             [  0.643888431046   0.356111568954]
    777             [0.00232031442167   0.997679685578]
    778             Emission matrix:
    779             [           0.0            1.0]
    780             [0.296080269308 0.703919730692]
    781             Initial probabilities: [1.0, 3.662034729231...e-20]
    782             sage: h = hmm.DiscreteHiddenMarkovModel([[.1,.9],[.4,.6]], [[1/2]*2]*2, [1/2]*2)
    783             sage: h.baum_welch([1,1,1,1,1,0,1,1,1,1,1,0,0], nsteps=1)
    784             sage: h
    785             Discrete Hidden Markov Model with 2 States and 2 Emissions
    786             Transition matrix:
    787             [0.1 0.9]
    788             [0.4 0.6]
    789             Emission matrix:
    790             [0.222426510432 0.777573489568]
    791             [0.234678486789 0.765321513211]
    792             Initial probabilities: [0.5, 0.5]
    793             sage: h = hmm.DiscreteHiddenMarkovModel([[.1,.9],[.4,.6]], [[1/2]*2]*2, [1/2]*2)
    794             sage: h.baum_welch([1,1,1,1,1,0,1,1,1,1,1,0,0], log_likelihood_cutoff=0.01)
    795             sage: h
    796             Discrete Hidden Markov Model with 2 States and 2 Emissions
    797             Transition matrix:
    798             [0.0999471960754  0.900052803925]
    799             [ 0.399248483887  0.600751516113]
    800             Emission matrix:
    801             [0.209328925617 0.790671074383]
    802             [ 0.24081056723  0.75918943277]
    803             Initial probabilities: [0.508779821929586..., 0.491220178070413...]       
    804 
    805         TESTS:
    806         We test training with non-default string symbols:
    807             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'])
    808             sage: a.baum_welch([['up','up'], ['down','up']])
    809             sage: a
    810             Discrete Hidden Markov Model with 2 States and 2 Emissions
    811             Transition matrix:
    812             [0.5 0.5]
    813             [0.5 0.5]
    814             Emission matrix:
    815             [0.75 0.25]
    816             [0.75 0.25]
    817             Initial probabilities: [0.5, 0.5]
    818             Emission symbols: ['up', 'down']
    819 
    820         NOTE: Training for models including silent states is not yet supported.
    821 
    822         REFERENCES:
    823             Rabiner, L.R.: "`A Tutorial on Hidden Markov Models and Selected
    824             Applications in Speech Recognition"', Proceedings of the IEEE,
    825             77, no 2, 1989, pp 257--285.
    826         """
    827         if len(training_seqs) > 0 and not isinstance(training_seqs[0], (list, tuple)):
    828             training_seqs = [training_seqs]
    829            
    830         if self._emission_symbols_dict:
    831             seqs = [[self._emission_symbols_dict[z] for z in x] for x in training_seqs]
    832         else:
    833             seqs = training_seqs
    834            
    835         cdef ghmm_dseq* d = malloc_ghmm_dseq(seqs)
    836        
    837         if nsteps or log_likelihood_cutoff:
    838             if nsteps is None:
    839                 nsteps = GHMM_MAX_ITER_BW
    840             if log_likelihood_cutoff is None:
    841                 log_likelihood_cutoff = GHMM_EPS_ITER_BW
    842             if ghmm_dmodel_baum_welch_nstep(self.m, d, nsteps, log_likelihood_cutoff):
    843                 raise RuntimeError, "error running Baum-Welch algorithm"
    844         else:
    845             if ghmm_dmodel_baum_welch(self.m, d):
    846                 raise RuntimeError, "error running Baum-Welch algorithm"
    847        
    848         ghmm_dseq_free(&d)
    849 
    850            
    851 ##################################################################################
    852 # Helper Functions
    853 ##################################################################################
    854 
    855 cdef ghmm_dseq* malloc_ghmm_dseq(seqs) except NULL:
    856     """
    857     Allocate a discrete sequence of samples.
    858 
    859     INPUT:
    860         seqs -- a list of sequences
    861 
    862     OUTPUT:
    863         C pointer to ghmm_dseq
    864     """
    865     cdef ghmm_dseq* d = ghmm_dseq_calloc(len(seqs))
    866     if d == NULL:
    867         raise MemoryError
    868     cdef int i, j, m, n
    869     m = len(seqs)
    870     d.seq_number = m
    871     d.capacity = m
    872     d.total_w = m
    873     for i from 0 <= i < m:
    874         v = seqs[i]
    875         n = len(v)
    876         d.seq[i] = <int*> safe_malloc(sizeof(int) * n)
    877         for j from 0 <= j < n:
    878             d.seq[i][j] = v[j]
    879         d.seq_len[i] = n
    880         d.seq_id[i] = i
    881         d.seq_w[i] = 1
    882     d.flags = 0
    883     return d
    884 
    885 
    886 def unpickle_discrete_hmm_v0(A, B, pi, emission_symbols,name):
    887     """
    888     TESTS:
    889         sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.0,1.0],[0.5,0.5]], [1,0], name='test model')
    890         sage: loads(dumps(m)) == m
    891         True
    892         sage: loads(dumps(m)).name()
    893         'test model'
    894         sage: sage.stats.hmm.hmm.unpickle_discrete_hmm_v0(m.transition_matrix(), m.emission_matrix(), m.initial_probabilities(), ['a','b'], m.name())
    895         Discrete Hidden Markov Model test model with 2 States and 2 Emissions
     192        sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.1,0.9],[0.5,0.5]], [.5,.5]); m
     193        Discrete Hidden Markov Model with 2 States and 2 Emissions
    896194        Transition matrix:
    897195        [0.4 0.6]
    898196        [0.1 0.9]
    899197        Emission matrix:
    900         [0.0 1.0]
     198        [0.1 0.9]
    901199        [0.5 0.5]
    902         Initial probabilities: [1.0, 0.0]
    903         Emission symbols: ['a', 'b']
     200        Initial probabilities: [0.5000, 0.5000]
     201        sage: m.log_likelihood([0,1,0,1,0,1])
     202        -4.66693474691329...
     203        sage: m.viterbi([0,1,0,1,0,1])
     204        ([1, 1, 1, 1, 1, 1], -5.3788328422087481)
     205        sage: m.baum_welch([0,1,0,1,0,1])
     206        (0.0, 22)
     207        sage: m
     208        Discrete Hidden Markov Model with 2 States and 2 Emissions
     209        Transition matrix:
     210        [1.0134345614...e-70               1.0]
     211        [              1.0 3.99743527136e-19]
     212        Emission matrix:
     213        [7.3802215662...e-54               1.0]
     214        [              1.0  3.9974352626e-19]
     215        Initial probabilities: [0.0000, 1.0000]
     216        sage: m.sample(10)
     217        [0, 1, 0, 1, 0, 1, 0, 1, 0, 1]
     218        sage: m.graph().plot()
    904219    """
    905     return DiscreteHiddenMarkovModel(A,B,pi,emission_symbols,name)
     220    cdef TimeSeries B
     221    cdef int n_out
     222    cdef object _emission_symbols, _emission_symbols_dict
     223
     224    def __init__(self, A, B, pi, emission_symbols=None, bint normalize=True):
     225        """
     226        Create a discrete emissions HMM with transition probability
     227        matrix A, emission probabilities given by B, initial state
     228        probabilities pi, and given emission symbols (which default
     229        to the first few nonnegative integers).
     230
     231        INPUT::
     232
     233           - A -- a list of lists or a square N x N matrix, whose
     234             (i,j) entry gives the probability of transitioning from
     235             state i to state j.
     236
     237           - B -- a list of N lists or a matrix with N rows, such that
     238             B[i,k] gives the probability of emitting symbol k while
     239             in state i.
     240
     241           - pi -- the probabilities of starting in each initial
     242             state, i.e,. pi[i] is the probability of starting in
     243             state i.
     244
     245           - emission_symbols -- None or list (default: None); if
     246             None, the emission_symbols are the ints [0..N-1], where N
     247             is the number of states.  Otherwise, they are the entries
     248             of the list emissions_symbols, which must all be hashable.
     249
     250           - normalize --bool (default: True); if given, input is
     251             normalized to define valid probability distributions,
     252             e.g., the entries of A are made nonnegative and the rows
     253             sum to 1.
     254
     255        EXAMPLES::
     256
     257            sage: hmm.DiscreteHiddenMarkovModel([.5,0,-1,.5], [[1],[1]],[.5,.5]).transition_matrix()
     258            [1.0 0.0]
     259            [0.0 1.0]
     260            sage: hmm.DiscreteHiddenMarkovModel([0,0,.1,.9], [[1],[1]],[.5,.5]).transition_matrix()
     261            [0.5 0.5]
     262            [0.1 0.9]
     263            sage: hmm.DiscreteHiddenMarkovModel([-1,-2,.1,.9], [[1],[1]],[.5,.5]).transition_matrix()
     264            [0.5 0.5]
     265            [0.1 0.9]
     266            sage: hmm.DiscreteHiddenMarkovModel([1,2,.1,1.2], [[1],[1]],[.5,.5]).transition_matrix()
     267            [ 0.333333333333  0.666666666667]
     268            [0.0769230769231  0.923076923077]
     269        """
     270        self.pi = util.initial_probs_to_TimeSeries(pi, normalize)
     271        self.N = len(self.pi)
     272        self.A = util.state_matrix_to_TimeSeries(A, self.N, normalize)
     273        self._emission_symbols = emission_symbols
     274        if self._emission_symbols is not None:
     275            self._emission_symbols_dict = dict([(y,x) for x,y in enumerate(emission_symbols)])
     276
     277        # TODO : normalization
     278        if not is_Matrix(B):
     279            B = matrix(B)
     280        if B.nrows() != self.N:
     281            raise ValueError, "number of rows of B must equal number of states"
     282        self.B = TimeSeries(B.list())
     283        self.n_out = B.ncols()
     284        if emission_symbols is not None and len(emission_symbols) != self.n_out:
     285            raise ValueError, "number of emission symbols must equal number of output states"
     286        cdef Py_ssize_t i
     287        if normalize:
     288            for i in range(self.N):
     289                util.normalize_probability_TimeSeries(self.B, i*self.n_out, (i+1)*self.n_out)
     290
     291    def __reduce__(self):
     292        """
     293        Used in pickling.
     294
     295        EXAMPLES:
     296            sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.0,1.0],[1,1]], [0,1], ['a','b'])
     297            sage: loads(dumps(m)) == m
     298            True
     299        """
     300        return unpickle_discrete_hmm_v1, \
     301               (self.A, self.B, self.pi, self.n_out, self._emission_symbols, self._emission_symbols_dict)
     302
     303    def __cmp__(self, other):
     304        """
     305        EXAMPLES::
     306
     307            sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.0,1.0],[0.5,0.5]], [.5,.5])
     308            sage: n = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.0,1.0],[0.5,0.5]], [1,0])
     309            sage: m == n
     310            False
     311            sage: m == m
     312            True
     313            sage: m < n
     314            True
     315            sage: n < m
     316            False
     317        """
     318        if not isinstance(other, DiscreteHiddenMarkovModel):
     319            raise ValueError
     320        return cmp(self.__reduce__()[1], other.__reduce__()[1])
     321
     322
     323    def emission_matrix(self):
     324        """
     325        Return the matrix whose i-th row specifies the emission
     326        probability distribution for the i-th state.  More precisely,
     327        the i,j entry of the matrix is the probability of the Markov
     328        model outputing the j-th symbol when it is in the i-th state.
     329
     330        OUTPUT:
     331
     332            - a Sage matrix with real double precision (RDF) entries.
     333
     334        EXAMPLES::
     335
     336            sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.1,0.9],[0.5,0.5]], [.5,.5])
     337            sage: E = m.emission_matrix(); E
     338            [0.1 0.9]
     339            [0.5 0.5]
     340
     341        The returned matrix is mutable, but changing it does not
     342        change the transition matrix for the model::
     343
     344            sage: E[0,0] = 0; E[0,1] = 1
     345            sage: m.emission_matrix()
     346            [0.1 0.9]
     347            [0.5 0.5]
     348        """
     349        from sage.matrix.constructor import matrix
     350        from sage.rings.all import RDF
     351        return matrix(RDF, self.N, self.n_out, self.B.list())
     352
     353
     354    def __repr__(self):
     355        """
     356        Return string representation of this discrete hidden Markov model.
     357
     358        EXAMPLES::
     359
     360            sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.1,0.9],[0.5,0.5]], [.2,.8])
     361            sage: m.__repr__()
     362            'Discrete Hidden Markov Model with 2 States and 2 Emissions\nTransition matrix:\n[0.4 0.6]\n[0.1 0.9]\nEmission matrix:\n[0.1 0.9]\n[0.5 0.5]\nInitial probabilities: [0.2000, 0.8000]'
     363        """
     364        s = "Discrete Hidden Markov Model with %s States and %s Emissions"%(
     365            self.N, self.n_out)
     366        s += '\nTransition matrix:\n%s'%self.transition_matrix()
     367        s += '\nEmission matrix:\n%s'%self.emission_matrix()
     368        s += '\nInitial probabilities: %s'%self.initial_probabilities()
     369        if self._emission_symbols is not None:
     370            s += '\nEmission symbols: %s'%self._emission_symbols
     371        return s
     372
     373    def _emission_symbols_to_IntList(self, obs):
     374        """
     375        Internal function used to convert a list of emission symbols to an IntList.
     376
     377        INPUT:
     378
     379            - obs -- a list of objects
     380
     381        OUTPUT:
     382
     383            - an IntList
     384
     385        EXAMPLES::
     386
     387            sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.1,0.9],[0.5,0.5]], [.2,.8], ['smile', 'frown'])
     388            sage: m._emission_symbols_to_IntList(['frown','smile'])
     389            [1, 0]
     390        """
     391        d = self._emission_symbols_dict
     392        return IntList([d[x] for x in obs])
     393
     394    def _IntList_to_emission_symbols(self, obs):
     395        """
     396        Internal function used to convert a list of emission symbols to an IntList.
     397
     398        INPUT:
     399
     400            - obs -- a list of objects
     401
     402        OUTPUT:
     403
     404            - an IntList
     405
     406        EXAMPLES::
     407
     408            sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.1,0.9],[0.5,0.5]], [.2,.8], ['smile', 'frown'])
     409            sage: m._IntList_to_emission_symbols(stats.IntList([0,0,1,0]))
     410            ['smile', 'smile', 'frown', 'smile']
     411        """
     412        d = self._emission_symbols
     413        return [d[x] for x in obs]
     414
     415    def log_likelihood(self, obs, bint scale=True):
     416        """
     417        Return the logarithm of the probability that this model produced the given
     418        observation sequence.  Thus the output is a non-positive number.
     419
     420        INPUT:
     421
     422            - obs -- sequence of observations
     423
     424            - scale -- boolean (default: True); if True, use rescaling
     425              to overoid loss of precision due to the very limited
     426              dynamic range of floats.  You should leave this as True
     427              unless the obs sequence is very small.
     428
     429        EXAMPLES::
     430
     431            sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.1,0.9],[0.5,0.5]], [.2,.8])
     432            sage: m.log_likelihood([0, 1, 0, 1, 1, 0, 1, 0, 0, 0])
     433            -7.3301308009370825
     434            sage: m.log_likelihood([0, 1, 0, 1, 1, 0, 1, 0, 0, 0], scale=False)
     435            -7.3301308009370816
     436            sage: m.log_likelihood([])
     437            0.0
     438
     439            sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.1,0.9],[0.5,0.5]], [.2,.8], ['happy','sad'])
     440            sage: m.log_likelihood(['happy','happy'])
     441            -1.6565295199679506
     442            sage: m.log_likelihood(['happy','sad'])
     443            -1.4731602941415523
     444
     445        Overflow from not using the scale option::
     446
     447            sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.1,0.9],[0.5,0.5]], [.2,.8])
     448            sage: m.log_likelihood([0,1]*1000, scale=True)
     449            -1433.8206666527281
     450            sage: m.log_likelihood([0,1]*1000, scale=False)
     451            -inf
     452        """
     453        if len(obs) == 0:
     454            return 0.0
     455        if self._emission_symbols is not None:
     456            obs = self._emission_symbols_to_IntList(obs)
     457        if not isinstance(obs, IntList):
     458            obs = IntList(obs)
     459        if scale:
     460            return self._forward_scale(obs)
     461        else:
     462            return self._forward(obs)
     463
     464    def _forward(self, IntList obs):
     465        """
     466        Memory-efficient implementation of the forward algorithm, without scaling.
     467
     468        INPUT:
     469
     470            - obs -- an integer list of observation states.
     471
     472        OUTPUT:
     473
     474            - float -- the log of the probability that the model produced this sequence
     475
     476        EXAMPLES::
     477
     478            sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.1,0.9],[0.5,0.5]], [.2,.8])
     479            sage: m._forward(stats.IntList([0,1]*10))
     480            -14.378663512219456
     481
     482        The forward algorithm computes the log likelihood::
     483
     484            sage: m.log_likelihood(stats.IntList([0,1]*10), scale=False)
     485            -14.378663512219456
     486
     487        But numerical overflow will happen (without scaling) for long sequences::
     488
     489            sage: m._forward(stats.IntList([0,1]*1000))
     490            -inf
     491        """
     492        if obs.max() > self.N or obs.min() < 0:
     493            raise ValueError, "invalid observation sequence, since it includes unknown states"
     494
     495        cdef Py_ssize_t i, j, t, T = len(obs)
     496
     497        cdef TimeSeries alpha  = TimeSeries(self.N), \
     498                        alpha2 = TimeSeries(self.N)
     499
     500        # Initialization
     501        for i in range(self.N):
     502            alpha[i] = self.pi[i] * self.B[i*self.n_out + obs._values[0]]
     503        alpha[i] = self.pi[i] * self.B[i*self.n_out + obs._values[0]]
     504
     505        # Induction
     506        cdef double s
     507        for t in range(1, T):
     508            for j in range(self.N):
     509                s = 0
     510                for i in range(self.N):
     511                    s += alpha._values[i] * self.A._values[i*self.N + j]
     512                alpha2._values[j] = s * self.B._values[j*self.n_out+obs._values[t]]
     513            for j in range(self.N):
     514                alpha._values[j] = alpha2._values[j]
     515
     516        # Termination
     517        return log(alpha.sum())
     518
     519    def _forward_scale(self, IntList obs):
     520        """
     521        Memory-efficient implementation of the forward algorithm, with scaling.
     522
     523        INPUT:
     524
     525            - obs -- an integer list of observation states.
     526
     527        OUTPUT:
     528
     529            - float -- the log of the probability that the model produced this sequence
     530
     531        EXAMPLES::
     532
     533            sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.1,0.9],[0.5,0.5]], [.2,.8])
     534            sage: m._forward_scale(stats.IntList([0,1]*10))
     535            -14.378663512219454
     536
     537        The forward algorithm computes the log likelihood::
     538
     539            sage: m.log_likelihood(stats.IntList([0,1]*10))
     540            -14.378663512219454
     541
     542        Note that the scale algorithm ensures that numerical overflow
     543        won't happen for long sequences like it does with the forward
     544        non-scaled algorithm::
     545
     546            sage: m._forward_scale(stats.IntList([0,1]*1000))
     547            -1433.8206666527281
     548            sage: m._forward(stats.IntList([0,1]*1000))
     549            -inf
     550
     551        A random sequence produced by the model is more likely::
     552
     553            sage: set_random_seed(0); v = m.sample(1000)
     554            sage: m._forward_scale(v)
     555            -686.87531893650555
     556        """
     557        # This is just like self._forward(obs) above, except at every step of the
     558        # algorithm, we rescale the vector alpha so that the sum of
     559        # the entries in alpha is 1.  Then, at the end of the
     560        # algorithm, the sum of probabilities (alpha.sum()) is of
     561        # course just 1.  However, the true probability that we want
     562        # is the product of the numbers that we divided by when
     563        # rescaling, since at each step of the iteration the next term
     564        # depends linearly on the previous one.  Instead of returning
     565        # the product, we return the sum of the logs, which avoid
     566        # numerical error.
     567        cdef Py_ssize_t i, j, t, T = len(obs)
     568
     569        # The running some of the log probabilities
     570        cdef double log_probability = 0, sum, a
     571
     572        cdef TimeSeries alpha  = TimeSeries(self.N), \
     573                        alpha2 = TimeSeries(self.N)
     574
     575        # Initialization
     576        sum = 0
     577        for i in range(self.N):
     578            a = self.pi[i] * self.B[i*self.n_out + obs._values[0]]
     579            alpha[i] = a
     580            sum += a
     581
     582        log_probability = log(sum)
     583        alpha.rescale(1/sum)
     584
     585        # Induction
     586        # The code below is just an optimized version of:
     587        #     alpha = (alpha * A).pairwise_product(B[O[t+1]])
     588        #     alpha = alpha.scale(1/alpha.sum())
     589        # along with keeping track of the log of the scaling factor.
     590        cdef double s
     591        for t in range(1, T):
     592            sum = 0
     593            for j in range(self.N):
     594                s = 0
     595                for i in range(self.N):
     596                    s += alpha._values[i] * self.A._values[i*self.N + j]
     597                a = s * self.B._values[j*self.n_out+obs._values[t]]
     598                alpha2._values[j] = a
     599                sum += a
     600
     601            log_probability += log(sum)
     602            for j in range(self.N):
     603                alpha._values[j] = alpha2._values[j] / sum
     604
     605        # Termination
     606        return log_probability
     607
     608    def generate_sequence(self, Py_ssize_t length):
     609        """
     610        Return a sample of the given length from this HMM.
     611
     612        INPUT:
     613
     614            - length -- positive integer
     615
     616        OUTPUT:
     617
     618            - an IntList or list of emission symbols
     619            - IntList of the actual states the model was in when
     620              emitting the corresponding symbols
     621
     622        EXAMPLES::
     623
     624        In this example, the emission symbols are not set::
     625
     626            sage: set_random_seed(0)
     627            sage: a = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.1,0.9]], [[1,0],[0,1]], [0,1])
     628            sage: a.generate_sequence(5)
     629            ([1, 0, 1, 1, 1], [1, 0, 1, 1, 1])
     630            sage: list(a.generate_sequence(1000)[0]).count(0)
     631            90
     632
     633        Here the emission symbols are set::
     634
     635            sage: set_random_seed(0)
     636            sage: a = hmm.DiscreteHiddenMarkovModel([[0.5,0.5],[0.1,0.9]], [[1,0],[0,1]], [0,1], ['up', 'down'])
     637            sage: a.generate_sequence(5)
     638            (['down', 'up', 'down', 'down', 'down'], [1, 0, 1, 1, 1])
     639        """
     640        if length < 0:
     641            raise ValueError, "length must be nonnegative"
     642
     643        # Create Integer lists for states and observations
     644        cdef IntList states = IntList(length)
     645        cdef IntList obs = IntList(length)
     646        if length == 0:
     647            # A special case
     648            if self._emission_symbols is None:
     649                return states, obs
     650            else:
     651                return states, []
     652
     653        # Setup variables, including random state.
     654        cdef Py_ssize_t i, j
     655        cdef randstate rstate = current_randstate()
     656        cdef int q = 0
     657        cdef double r, accum
     658        r = rstate.c_rand_double()
     659
     660        # Now choose random variables from our discrete distribution.
     661
     662        # This standard naive algorithm has complexity that is linear
     663        # in the number of states.  It might be possible to replace it
     664        # by something that is more efficient.  However, make sure to
     665        # refactor this code into distribution.pyx before doing so.
     666        # Note that state switching involves the same algorithm
     667        # (below).  Use GSL as described here to make this O(1):
     668        #    http://www.gnu.org/software/gsl/manual/html_node/General-Discrete-Distributions.html
     669
     670        # Choose initial state:
     671        accum = 0
     672        for i in range(self.N):
     673            if r < self.pi._values[i] + accum:
     674                q = i
     675            else:
     676                accum += self.pi._values[i]
     677
     678        states._values[0] = q
     679        # Generate a symbol from state q
     680        obs._values[0] = self._gen_symbol(q, rstate.c_rand_double())
     681
     682        cdef double* row
     683        cdef int O
     684        _sig_on
     685        for i in range(1, length):
     686            # Choose next state
     687            accum = 0
     688            row = self.A._values + q*self.N
     689            r = rstate.c_rand_double()
     690            for j in range(self.N):
     691                if r < row[j] + accum:
     692                    q = j
     693                    break
     694                else:
     695                    accum += row[j]
     696            states._values[i] = q
     697            # Generate symbol from this new state q
     698            obs._values[i] = self._gen_symbol(q, rstate.c_rand_double())
     699        _sig_off
     700
     701        if self._emission_symbols is None:
     702            # No emission symbol mapping
     703            return obs, states
     704        else:
     705            # Emission symbol mapping, so change our intlist into a list of symbols
     706            return self._IntList_to_emission_symbols(obs), states
     707
     708    cdef int _gen_symbol(self, int q, double r):
     709        """
     710        Generate a symbol in state q using the randomly chosen
     711        floating point number r, which should be between 0 and 1.
     712
     713        INPUT:
     714
     715            - q -- a nonnegative integer, which specifies a state
     716
     717            - r -- a real number between 0 and 1
     718
     719        OUTPUT:
     720
     721            - a nonnegative int
     722
     723        EXAMPLES:
     724
     725            sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.1,0.9],[0.9,0.1]], [.2,.8])
     726            sage: set_random_seed(0)
     727            sage: m.generate_sequence(10)  # indirect test
     728            ([0, 1, 0, 0, 0, 0, 1, 0, 1, 1], [1, 0, 1, 1, 1, 1, 0, 0, 0, 0])
     729
     730        """
     731        cdef Py_ssize_t j
     732        cdef double a, accum = 0
     733        # See the comments above about switching to GSL for this; also note
     734        # that this should get factored out into a call to something
     735        # defined in distributions.pyx.
     736        for j in range(self.n_out):
     737            a = self.B._values[q*self.N+j]
     738            if r < a + accum:
     739                return j
     740            else:
     741                accum += a
     742        # None of the values was selected: shouldn't this only happen if the
     743        # distribution is invalid?   Anyway, this will get factored out.
     744        return self.n_out - 1
     745
     746    def viterbi(self, obs, log_scale=True):
     747        """
     748        Determine "the" hidden sequence of states that is most likely
     749        to produce the given sequence seq of observations, along with
     750        the probability that this hidden sequence actually produced
     751        the observation.
     752
     753        INPUT:
     754
     755            - seq -- sequence of emitted ints or symbols
     756
     757            - log_scale -- bool (default: True) whether to scale the
     758              sequence in order to avoid numerical overflow.
     759
     760        OUTPUT:
     761
     762            - list -- "the" most probable sequence of hidden states, i.e.,
     763              the Viterbi path.
     764
     765            - float -- log of probability that the observed sequence
     766              was produced by the Viterbi sequence of states.
     767
     768        EXAMPLES::
     769
     770            sage: a = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.1,0.9]], [[0.9,0.1],[0.1,0.9]], [0.5,0.5])
     771            sage: a.viterbi([1,0,0,1,0,0,1,1])
     772            ([1, 0, 0, 1, ..., 0, 1, 1], -11.06245322477221...)
     773
     774        We predict the state sequence when the emissions are 3/4 and 'abc'.::
     775
     776            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'])
     777
     778        Note that state 0 is common below, despite the model trying hard to
     779        switch to state 1::
     780
     781            sage: a.viterbi([3/4, 'abc', 'abc'] + [3/4]*10)
     782            ([0, 1, 1, 0, 0 ... 0, 0, 0, 0, 0], -25.299405845367794)
     783        """
     784        if self._emission_symbols is not None:
     785            obs = self._emission_symbols_to_IntList(obs)
     786        elif not isinstance(obs, IntList):
     787            obs = IntList(obs)
     788        if log_scale:
     789            return self._viterbi_scale(obs)
     790        else:
     791            return self._viterbi(obs)
     792
     793    cpdef _viterbi(self, IntList obs):
     794        """
     795        Used internally to compute the viterbi path, without
     796        rescaling.  This can be useful for short sequences.
     797
     798        INPUT:
     799
     800            - obs -- IntList
     801
     802        OUTPUT:
     803
     804            - IntList (most likely state sequence)
     805
     806            - log of probability that the observed sequence was
     807              produced by the Viterbi sequence of states.
     808
     809        EXAMPLES::
     810
     811            sage: m = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.9,0.1]], [[1,0],[0,1]], [.2,.8])
     812            sage: m._viterbi(stats.IntList([1]*5))
     813            ([1, 1, 1, 1, 1], -9.4334839232903924)
     814            sage: m._viterbi(stats.IntList([0]*5))
     815            ([0, 0, 0, 0, 0], -10.819778284410283)
     816
     817        Long sequences will overflow::
     818
     819            sage: m._viterbi(stats.IntList([0]*1000))
     820            ([0, 0, 0, 0, 0 ... 0, 0, 0, 0, 0], -inf)
     821        """
     822        cdef Py_ssize_t t, T = obs._length
     823        cdef IntList state_sequence = IntList(T)
     824        if T == 0:
     825            return state_sequence, 0.0
     826
     827        cdef int i, j, N = self.N
     828
     829        # delta[i] is the maximum of the probabilities over all
     830        # paths ending in state i.
     831        cdef TimeSeries delta = TimeSeries(N), delta_prev = TimeSeries(N)
     832
     833        # We view psi as an N x T matrix of ints. The quantity
     834        #          psi[N*t + j]
     835        # is a most probable hidden state at time t-1, given
     836        # that j is a most probable state at time j.
     837        cdef IntList psi = IntList(N * T)  # initialized to 0 by default
     838
     839        # Initialization:
     840        for i in range(N):
     841            delta._values[i] = self.pi._values[i] * self.B._values[self.n_out*i + obs._values[0]]
     842
     843        # Recursion:
     844        cdef double mx, tmp
     845        cdef int index
     846        for t in range(1, T):
     847            delta_prev, delta = delta, delta_prev
     848            for j in range(N):
     849                # delta_t[j] = max_i(delta_{t-1}(i) a_{i,j}) * b_j(obs[t])
     850                mx = -1
     851                index = -1
     852                for i in range(N):
     853                    tmp = delta_prev._values[i]*self.A._values[i*N+j]
     854                    if tmp > mx:
     855                        mx = tmp
     856                        index = i
     857                delta._values[j] = mx * self.B._values[self.n_out*j + obs._values[t]]
     858                psi._values[N*t + j] = index
     859
     860        # Termination:
     861        mx, index = delta.max(index=True)
     862
     863        # Path (state sequence) backtracking:
     864        state_sequence._values[T-1] = index
     865        t = T-2
     866        while t >= 0:
     867            state_sequence._values[t] = psi._values[N*(t+1) + state_sequence._values[t+1]]
     868            t -= 1
     869
     870        return state_sequence, log(mx)
     871
     872
     873    cpdef _viterbi_scale(self, IntList obs):
     874        """
     875        Used internally to compute the viterbi path with rescaling.
     876
     877        INPUT:
     878
     879            - obs -- IntList
     880
     881        OUTPUT:
     882
     883            - IntList (most likely state sequence)
     884
     885            - log of probability that the observed sequence was
     886              produced by the Viterbi sequence of states.
     887
     888        EXAMPLES::
     889
     890            sage: m = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.9,0.1]], [[.5,.5],[0,1]], [.2,.8])
     891            sage: m._viterbi_scale(stats.IntList([1]*10))
     892            ([1, 0, 1, 0, 1, 0, 1, 0, 1, 0], -4.6371240950343733)
     893
     894        Long sequences should not overflow::
     895
     896            sage: m._viterbi_scale(stats.IntList([1]*1000))
     897            ([1, 0, 1, 0, 1 ... 0, 1, 0, 1, 0], -452.05188897345...)
     898        """
     899        # The algorithm is the same as _viterbi above, except
     900        # we take the logarithms of everything first, and add
     901        # instead of multiply.
     902        cdef Py_ssize_t t, T = obs._length
     903        cdef IntList state_sequence = IntList(T)
     904        if T == 0:
     905            return state_sequence, 0.0
     906        cdef int i, j, N = self.N
     907
     908        # delta[i] is the maximum of the probabilities over all
     909        # paths ending in state i.
     910        cdef TimeSeries delta = TimeSeries(N), delta_prev = TimeSeries(N)
     911
     912        # We view psi as an N x T matrix of ints. The quantity
     913        #          psi[N*t + j]
     914        # is a most probable hidden state at time t-1, given
     915        # that j is a most probable state at time j.
     916        cdef IntList psi = IntList(N * T)  # initialized to 0 by default
     917
     918        # Log Preprocessing
     919        cdef TimeSeries A = self.A.log()
     920        cdef TimeSeries B = self.B.log()
     921        cdef TimeSeries pi = self.pi.log()
     922
     923        # Initialization:
     924        for i in range(N):
     925            delta._values[i] = pi._values[i] + B._values[self.n_out*i + obs._values[0]]
     926
     927        # Recursion:
     928        cdef double mx, tmp, minus_inf = float('-inf')
     929        cdef int index
     930
     931        for t in range(1, T):
     932            delta_prev, delta = delta, delta_prev
     933            for j in range(N):
     934                # delta_t[j] = max_i(delta_{t-1}(i) a_{i,j}) * b_j(obs[t])
     935                mx = minus_inf
     936                index = -1
     937                for i in range(N):
     938                    tmp = delta_prev._values[i] + A._values[i*N+j]
     939                    if tmp > mx:
     940                        mx = tmp
     941                        index = i
     942                delta._values[j] = mx + B._values[self.n_out*j + obs._values[t]]
     943                psi._values[N*t + j] = index
     944
     945        # Termination:
     946        mx, index = delta.max(index=True)
     947
     948        # Path (state sequence) backtracking:
     949        state_sequence._values[T-1] = index
     950        t = T-2
     951        while t >= 0:
     952            state_sequence._values[t] = psi._values[N*(t+1) + state_sequence._values[t+1]]
     953            t -= 1
     954
     955        return state_sequence, mx
     956
     957    cdef TimeSeries _backward_scale_all(self, IntList obs, TimeSeries scale):
     958        r"""
     959        Return the scaled matrix of values `\beta_t(i)` that appear in
     960        the backtracking algorithm.  This function is used internally
     961        by the Baum-Welch algorithm.
     962
     963        The matrix is stored as a TimeSeries T, such that
     964        `\beta_t(i) = T[t*N + i]` where N is the number of states of
     965        the Hidden Markov Model.
     966
     967        The quantity beta_t(i) is the probability of observing the
     968        sequence obs[t+1:] assuming that the model is in state i at
     969        time t.
     970
     971        INPUT:
     972
     973            - obs -- IntList
     974            - scale -- series that is *changed* in place, so that
     975              after calling this function, scale[t] is value that is
     976              used to scale each of the `\beta_t(i)`.
     977
     978        OUTPUT:
     979
     980            - a TimeSeries of values beta_t(i).
     981            - the input object scale is modified
     982        """
     983        cdef Py_ssize_t t, T = obs._length
     984        cdef int N = self.N, i, j
     985        cdef double s
     986        cdef TimeSeries beta = TimeSeries(N*T, initialize=False)
     987
     988        # 1. Initialization
     989        for i in range(N):
     990            beta._values[(T-1)*N + i] = 1/scale._values[T-1]
     991
     992        # 2. Induction
     993        t = T-2
     994        while t >= 0:
     995            for i in range(N):
     996                s = 0
     997                for j in range(N):
     998                    s += self.A._values[i*N+j] * \
     999                         self.B._values[j*self.n_out+obs._values[t+1]] * beta._values[(t+1)*N+j]
     1000                beta._values[t*N + i] = s/scale._values[t]
     1001            t -= 1
     1002        return beta
     1003
     1004    cdef _forward_scale_all(self, IntList obs):
     1005        """
     1006        Return scaled values alpha_t(i), the sequence of scalings, and
     1007        the log probability.
     1008
     1009        INPUT:
     1010
     1011            - obs -- IntList
     1012
     1013        OUTPUT:
     1014
     1015            - TimeSeries alpha with alpha_t(i) = alpha[t*N + i]
     1016            - TimeSeries scale with scale[t] the scaling at step t
     1017            - float -- log_probability of the observation sequence
     1018              being produced by the model.
     1019        """
     1020        cdef Py_ssize_t i, j, t, T = len(obs)
     1021        cdef int N = self.N
     1022
     1023        # The running some of the log probabilities
     1024        cdef double log_probability = 0, sum, a
     1025
     1026        cdef TimeSeries alpha = TimeSeries(N*T, initialize=False)
     1027        cdef TimeSeries scale = TimeSeries(T, initialize=False)
     1028
     1029        # Initialization
     1030        sum = 0
     1031        for i in range(self.N):
     1032            a = self.pi._values[i] * self.B._values[i*self.n_out + obs._values[0]]
     1033            alpha._values[0*N + i] = a
     1034            sum += a
     1035
     1036        scale._values[0] = sum
     1037        log_probability = log(sum)
     1038        for i in range(self.N):
     1039            alpha._values[0*N + i] /= sum
     1040
     1041        # Induction
     1042        # The code below is just an optimized version of:
     1043        #     alpha = (alpha * A).pairwise_product(B[O[t+1]])
     1044        #     alpha = alpha.scale(1/alpha.sum())
     1045        # along with keeping track of the log of the scaling factor.
     1046        cdef double s
     1047        for t in range(1,T):
     1048            sum = 0
     1049            for j in range(N):
     1050                s = 0
     1051                for i in range(N):
     1052                    s += alpha._values[(t-1)*N + i] * self.A._values[i*N + j]
     1053                a = s * self.B._values[j*self.n_out + obs._values[t]]
     1054                alpha._values[t*N + j] = a
     1055                sum += a
     1056
     1057            log_probability += log(sum)
     1058            scale._values[t] = sum
     1059            for j in range(self.N):
     1060                alpha._values[t*N + j] /= sum
     1061
     1062        # Termination
     1063        return alpha, scale, log_probability
     1064
     1065    cdef TimeSeries _baum_welch_gamma(self, TimeSeries alpha, TimeSeries beta):
     1066        """
     1067        Used internally to compute the scaled quantity gamma_t(j)
     1068        appearing in the Baum-Welch reestimation algorithm.
     1069
     1070        The quantity gamma_t(j) is the (scaled!) probability of being
     1071        in state j at time t, given the observation sequence.
     1072
     1073        INPUT:
     1074
     1075            - alpha -- TimeSeries as output by the scaled forward algorithm
     1076            - beta -- TimeSeries as output by the scaled backward algorithm
     1077
     1078        OUTPUT:
     1079
     1080            - TimeSeries gamma such that gamma[t*N+j] is gamma_t(j).
     1081        """
     1082        cdef int j, N = self.N
     1083        cdef Py_ssize_t t, T = alpha._length//N
     1084        cdef TimeSeries gamma = TimeSeries(alpha._length, initialize=False)
     1085        cdef double denominator
     1086        _sig_on
     1087        for t in range(T):
     1088            denominator = 0
     1089            for j in range(N):
     1090                gamma._values[t*N + j] = alpha._values[t*N + j] * beta._values[t*N + j]
     1091                denominator += gamma._values[t*N + j]
     1092            for j in range(N):
     1093                gamma._values[t*N + j] /= denominator
     1094        _sig_off
     1095        return gamma
     1096
     1097    cdef TimeSeries _baum_welch_xi(self, TimeSeries alpha, TimeSeries beta, IntList obs):
     1098        """
     1099        Used internally to compute the scaled quantity xi_t(i,j)
     1100        appearing in the Baum-Welch reestimation algorithm.
     1101
     1102        INPUT:
     1103
     1104            - alpha -- TimeSeries as output by the scaled forward algorithm
     1105            - beta -- TimeSeries as output by the scaled backward algorithm
     1106            - obs -- IntList of observations
     1107
     1108        OUTPUT:
     1109
     1110            - TimeSeries xi such that xi[t*N*N + i*N + j] = xi_t(i,j).
     1111        """
     1112        cdef int i, j, N = self.N
     1113        cdef double sum
     1114        cdef Py_ssize_t t, T = alpha._length//N
     1115        cdef TimeSeries xi = TimeSeries(T*N*N, initialize=False)
     1116        _sig_on
     1117        for t in range(T-1):
     1118            sum = 0.0
     1119            for i in range(N):
     1120                for j in range(N):
     1121                    xi._values[t*N*N + i*N + j] = alpha._values[t*N + i]*beta._values[(t+1)*N + j]*\
     1122                       self.A._values[i*N + j] * self.B._values[j*self.n_out + obs._values[t+1]]
     1123                    sum += xi._values[t*N*N + i*N + j]
     1124            for i in range(N):
     1125                for j in range(N):
     1126                    xi._values[t*N*N + i*N + j] /= sum
     1127        _sig_off
     1128        return xi
     1129
     1130    def baum_welch(self, obs, int max_iter=100, double log_likelihood_cutoff=1e-4, bint fix_emissions=False):
     1131        """
     1132        Given an observation sequence obs, improve this HMM using the
     1133        Baum-Welch algorithm to increase the probability of observing obs.
     1134
     1135        INPUT:
     1136
     1137            - obs -- list of emissions
     1138
     1139            - max_iter -- integer (default: 100) maximum number
     1140              of Baum-Welch steps to take
     1141
     1142            - log_likehood_cutoff -- positive float (default: 1e-4);
     1143              the minimal improvement in likelihood with respect to
     1144              the last iteration required to continue. Relative value
     1145              to log likelihood.
     1146
     1147            - fix_emissions -- bool (default: False); if True, do not
     1148              change emissions when updating
     1149
     1150        OUTPUT:
     1151
     1152            - changes the model in places, and returns the log
     1153              likelihood and number of iterations.
     1154
     1155        EXAMPLES::
     1156
     1157            sage: m = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.9,0.1]], [[.5,.5],[0,1]], [.2,.8])
     1158            sage: m.baum_welch([1,0]*20, log_likelihood_cutoff=0)
     1159            (0.0, 4)
     1160            sage: m
     1161            Discrete Hidden Markov Model with 2 States and 2 Emissions
     1162            Transition matrix:
     1163            [1.35152697077e-51               1.0]
     1164            [              1.0               0.0]
     1165            Emission matrix:
     1166            [              1.0 6.46253713885e-52]
     1167            [              0.0               1.0]
     1168            Initial probabilities: [0.0000, 1.0000]
     1169
     1170        The following illustrates how Baum-Welch is only a local
     1171        optimizer, i.e., the above model is far more likely to produce
     1172        the sequence [1,0]*20 than the one we get below::
     1173
     1174            sage: m = hmm.DiscreteHiddenMarkovModel([[0.5,0.5],[0.5,0.5]], [[.5,.5],[.5,.5]], [.5,.5])
     1175            sage: m.baum_welch([1,0]*20, log_likelihood_cutoff=0)
     1176            (-27.725887222397784, 1)
     1177            sage: m
     1178            Discrete Hidden Markov Model with 2 States and 2 Emissions
     1179            Transition matrix:
     1180            [0.5 0.5]
     1181            [0.5 0.5]
     1182            Emission matrix:
     1183            [0.5 0.5]
     1184            [0.5 0.5]
     1185            Initial probabilities: [0.5000, 0.5000]
     1186
     1187        We illustrate fixing emissions::
     1188
     1189            sage: m = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.9,0.1]], [[.5,.5],[.2,.8]], [.2,.8])
     1190            sage: set_random_seed(0); v = m.sample(100)
     1191            sage: m.baum_welch(v,fix_emissions=True)
     1192            (-66.986308569187742, 100)
     1193            sage: m.emission_matrix()
     1194            [0.5 0.5]
     1195            [0.2 0.8]
     1196            sage: m = hmm.DiscreteHiddenMarkovModel([[0.1,0.9],[0.9,0.1]], [[.5,.5],[.2,.8]], [.2,.8])
     1197            sage: m.baum_welch(v)
     1198            (-66.7823606592935..., 100)
     1199            sage: m.emission_matrix()
     1200            [0.530308574863 0.469691425137]
     1201            [0.290977555017 0.709022444983]
     1202        """
     1203        if self._emission_symbols is not None:
     1204            obs = self._emission_symbols_to_IntList(obs)
     1205        elif not isinstance(obs, IntList):
     1206            obs = IntList(obs)
     1207        cdef IntList _obs = obs
     1208        cdef TimeSeries alpha, beta, scale, gamma, xi
     1209        cdef double log_probability, log_probability_prev, delta
     1210        cdef int i, j, k, N, n_iterations
     1211        cdef Py_ssize_t t, T
     1212        cdef double denominator_A, numerator_A, denominator_B, numerator_B
     1213
     1214        # Initialization
     1215        alpha, scale, log_probability = self._forward_scale_all(_obs)
     1216        beta = self._backward_scale_all(_obs, scale)
     1217        gamma = self._baum_welch_gamma(alpha, beta)
     1218        xi = self._baum_welch_xi(alpha, beta, _obs)
     1219        log_probability_prev = log_probability
     1220        N = self.N
     1221        n_iterations = 0
     1222        T = len(_obs)
     1223
     1224        # Re-estimation
     1225        while True:
     1226
     1227            # Reestimate frequency of state i in time t=0
     1228            for i in range(N):
     1229                self.pi._values[i] = gamma._values[0*N+i]
     1230
     1231            # Reestimate transition matrix and emissions probability in
     1232            # each state.
     1233            for i in range(N):
     1234                denominator_A = 0.0
     1235                for t in range(T-1):
     1236                    denominator_A += gamma._values[t*N+i]
     1237                for j in range(N):
     1238                    numerator_A = 0.0
     1239                    for t in range(T-1):
     1240                        numerator_A += xi._values[t*N*N+i*N+j]
     1241                    self.A._values[i*N+j] = numerator_A / denominator_A
     1242
     1243                if not fix_emissions:
     1244                    denominator_B = denominator_A + gamma._values[(T-1)*N + i]
     1245                    for k in range(self.n_out):
     1246                        numerator_B = 0.0
     1247                        for t in range(T):
     1248                            if _obs._values[t] == k:
     1249                                numerator_B += gamma._values[t*N + i]
     1250                        self.B._values[i*self.n_out + k] = numerator_B / denominator_B
     1251
     1252            # Initialization
     1253            alpha, scale, log_probability = self._forward_scale_all(_obs)
     1254            beta = self._backward_scale_all(_obs, scale)
     1255            gamma = self._baum_welch_gamma(alpha, beta)
     1256            xi = self._baum_welch_xi(alpha, beta, _obs)
     1257
     1258            # Compute the difference between the log probability of
     1259            # two iterations.
     1260            delta = log_probability - log_probability_prev
     1261            log_probability_prev = log_probability
     1262            n_iterations += 1
     1263
     1264            # If the log probability does not change by more than
     1265            # delta, then terminate
     1266            if delta <= log_likelihood_cutoff or n_iterations >= max_iter:
     1267                break
     1268
     1269        return log_probability, n_iterations
     1270
     1271
     1272# Keep this -- it's for backwards compatibility with the GHMM based implementation
     1273def unpickle_discrete_hmm_v0(A, B, pi, emission_symbols, name):
     1274    """
     1275    TESTS::
     1276
     1277        sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.0,1.0],[0.5,0.5]], [1,0])
     1278        sage: sage.stats.hmm.hmm.unpickle_discrete_hmm_v0(m.transition_matrix(), m.emission_matrix(), m.initial_probabilities(), ['a','b'], 'test model')
     1279        Discrete Hidden Markov Model with 2 States and 2 Emissions...
     1280    """
     1281    return DiscreteHiddenMarkovModel(A,B,pi,emission_symbols,normalize=False)
     1282
     1283def unpickle_discrete_hmm_v1(A, B, pi, n_out, emission_symbols, emission_symbols_dict):
     1284    """
     1285    TESTS::
     1286
     1287        sage: m = hmm.DiscreteHiddenMarkovModel([[0.4,0.6],[0.1,0.9]], [[0.0,1.0],[0.5,0.5]], [1,0],['a','b'])
     1288        sage: loads(dumps(m)) == m   # indirect test
     1289        True
     1290    """
     1291    cdef DiscreteHiddenMarkovModel m = PY_NEW(DiscreteHiddenMarkovModel)
     1292    m.A = A
     1293    m.B = B
     1294    m.pi = pi
     1295    m.n_out = n_out
     1296    m._emission_symbols = emission_symbols
     1297    m._emission_symbols_dict = emission_symbols_dict
     1298    return m
     1299
     1300

  • deleted file sage/stats/hmm/misc.pxi 

    diff --git a/sage/stats/hmm/misc.pxi b/sage/stats/hmm/misc.pxi
deleted file mode 100644
    + -  
    1 #############################################################################
    2 #       Copyright (C) 2008 William Stein <wstein@gmail.com>
    3 #  Distributed under the terms of the GNU General Public License (GPL),
    4 #  version 2 or any later version at your option.
    5 #  The full text of the GPL is available at:
    6 #                  http://www.gnu.org/licenses/
    7 #############################################################################
    8 
    9 cdef extern from "string.h":
    10     char *strcpy(char *s1, char *s2)
    11     void* memcpy(void* dst, void* src, size_t len)
    12 
    13 cdef double* to_double_array(v) except NULL:
    14     """
    15     Transform a Python list of floats to a C array of doubles.  The caller is
    16     responsible for deallocating the resulting memory.
    17    
    18     INPUT:
    19         v -- a list of objects coercible to floats
    20     OUTPUT:
    21         a newly allocated C array of doubles
    22     """
    23     cdef double x
    24     cdef double* w = <double*> safe_malloc(sizeof(double)*len(v))
    25     cdef Py_ssize_t i = 0
    26     for x in v:
    27         w[i] = x
    28         i += 1
    29     return w
    30 
    31 cdef int* to_int_array(v) except NULL:
    32     """
    33     Transform a Python list of ints to a C array of ints.  The caller is
    34     responsible for deallocating the resulting memory.
    35    
    36     INPUT:
    37         v -- a list of objects coercible to ints
    38     OUTPUT:
    39         a newly allocated C array of ints
    40     """
    41     cdef int x
    42     cdef int* w = <int*> safe_malloc(sizeof(int)*len(v))
    43     cdef Py_ssize_t i = 0
    44     for x in v:
    45         w[i] = x
    46         i += 1
    47     return w
    48 
    49 cdef void* safe_malloc(int bytes) except NULL:
    50     """
    51     malloc the given bytes of memory and check that the malloc
    52     succeeds -- if not raise a MemoryError.
    53 
    54     INPUT:
    55         bytes -- a nonnegatie integer
    56        
    57     OUTPUT:
    58         void pointer or raise a MemoryError.
    59     """
    60     cdef void* t = sage_malloc(bytes)
    61     if not t:
    62         raise MemoryError, "error allocating memory for Hidden Markov Model"
    63     return t
    64        

  • new file sage/stats/hmm/util.pxd 

    diff --git a/sage/stats/hmm/util.pxd b/sage/stats/hmm/util.pxd
new file mode 100644
    - +  
     1from sage.finance.time_series cimport TimeSeries
     2
     3cdef class HMM_Util:
     4    cpdef normalize_probability_TimeSeries(self, TimeSeries T, Py_ssize_t i, Py_ssize_t j)
     5    cpdef TimeSeries initial_probs_to_TimeSeries(self, pi, bint normalize)
     6    cpdef TimeSeries state_matrix_to_TimeSeries(self, A, int N, bint normalize)
     7   

  • new file sage/stats/hmm/util.pyx 

    diff --git a/sage/stats/hmm/util.pyx b/sage/stats/hmm/util.pyx
new file mode 100644
    - +  
     1"""
     2Hidden Markov Models -- Utility functions
     3
     4AUTHOR:
     5
     6   - William Stein, 2010-03
     7"""
     8
     9#############################################################################
     10#       Copyright (C) 2010 William Stein <wstein@gmail.com>
     11#  Distributed under the terms of the GNU General Public License (GPL) v2+.
     12#  The full text of the GPL is available at:
     13#                  http://www.gnu.org/licenses/
     14#############################################################################
     15   
     16include "../../ext/stdsage.pxi"
     17
     18from sage.matrix.all import is_Matrix
     19from sage.misc.flatten  import flatten
     20
     21cdef class HMM_Util:
     22    """
     23    A class used in order to share cdef's methods between different files.
     24    """
     25    cpdef normalize_probability_TimeSeries(self, TimeSeries T, Py_ssize_t i, Py_ssize_t j):
     26        """
     27        This function is used internally by the Hidden Markov Models code.
     28
     29        Replace entries of T[i:j] in place so that they are all
     30        nonnegative and sum to 1.  Negative entries are replaced by 0 and
     31        T[i:j] is then rescaled to ensure that the sum of the entries in
     32        each row is equal to 1.  If all entries are 0, replace them
     33        by 1/(j-i).
     34
     35        INPUT:
     36
     37            - T -- a TimeSeries
     38