ghmm.py

General Hidden Markov Model Library 一个通用的隐马尔科夫模型的C代码库

#!/usr/bin/env python2.3
################################################################################
#
#       This file is part of the General Hidden Markov Model Library,
#       GHMM version 0.8_beta1, see http://ghmm.org
#
#       file:    ghmm.py
#       authors: Benjamin Georgi, Wasinee Rungsarityotin, Alexander Schliep,
#                Janne Grunau
#
#       Copyright (C) 1998-2004 Alexander Schliep
#       Copyright (C) 1998-2001 ZAIK/ZPR, Universitaet zu Koeln
#       Copyright (C) 2002-2004 Max-Planck-Institut fuer Molekulare Genetik,
#                               Berlin
#
#       Contact: [e-mail]
#
#       This library is free software; you can redistribute it and/or
#       modify it under the terms of the GNU Library General Public
#       License as published by the Free Software Foundation; either
#       version 2 of the License, or (at your option) any later version.
#
#       This library is distributed in the hope that it will be useful,
#       but WITHOUT ANY WARRANTY; without even the implied warranty of
#       MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
#       Library General Public License for more details.
#
#       You should have received a copy of the GNU Library General Public
#       License along with this library; if not, write to the Free
#       Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
#
#
#
################################################################################

"""
 
The Design of ghmm.py

HMMs are stochastic models which encode a probability density over
sequences of symbols. These symbols can be discrete letters (A,C,G and
T for DNA; 1,2,3,4,5,6 for dice), real numbers (weather measurement
over time: temperature) or vectors of either or the combination
thereof (weather again: temperature, pressure, percipitation).

Note: We will always talk about emissions, emission sequence and so
forth when we refer to the sequence of symbols. Another name
for the same object is observation resp. observation sequence.


The objects one has to deal with in HMM modelling are the following

1) The domain the emissions come from: the EmissionDomain. Domain
is to be understood mathematically and to encompass both discrete,
  finite alphabets and fields such as the real numbers or intervals
   of the reals.

   For technical reasons there can be two representations of an
   emission symbol: an external and an internal. The external
   representation is the view of the application using ghmm.py. The
   internal one is what is used in both ghmm.py and the ghmm
   C-library. Representations can coincide, but this is not
   guaranteed. Discrete alphabets of size k are represented as
   [0,1,2,...,k-1] internally.  It is the domain objects job to
   provide a mapping between representations in both directions.

   NOTE: Do not make assumptions about the internal
   representations. It might change.


2) Every domain has to afford a distribution, which is usually
   parameterized. A distribution associated with a domain
   should allow us to compute $\Prob[x| distribution parameters]$
   efficiently.

   The distribution defines the *type* of distribution which
   we will use to model emissions in *every state* of the HMM.
   The *type* of distribution will be identical for all states,
   their *parameterizations* will differ from state to state.


3) We will consider a Sequence of emissions from the same emission
   domain and very often sets of such sequences: SequenceSet


4) The HMM: The HMM consists of two major components: A Markov chain
   over states (implemented as a weighted directed graph with
   adjacency and inverse-adjacency lists) and the emission
   distributions per-state. For reasons of efficiency the HMM itself
   is *static*, as far as the topology of the underlying Markov chain
   (and obviously the EmissionDomain) are concerned. You cannot add or
   delete transitions in an HMM.

   Transition probabilities and the parameters of the per-state
   emission distributions can be easily modified. Particularly,
   Baum-Welch reestimation is supported.  While a transition cannot be
   deleted from the graph, you can set the transition probability to
   zero, which has the same effect from the theoretical point of
   view. However, the corresponding edge in the graph is still
   traversed in the computation.

   States in HMMs are referred to by their integer index. State sequences
   are simply list of integers.

   If you want to store application specific data for each state you
   have to do it yourself.

   Subclasses of HMM implement specific types of HMM. The type depends
   on the EmissionDomain, the Distribution used, the specific
   extensions to the 'standard' HMMs and so forth


5) HMMFactory: This provides a way of constucting HMMs. Classes derived
   from HMMFactory allow to read HMMs from files, construct them explicitly
   from, for a discrete alphabet, transition matrix, emission matrix and prior
   or serve as the basis for GUI-based model building.

   There are several ways of using the HMMFactory.

   Static construction:

   HMMOpen(fileName) # Calls an object of type HMMOpen instantiated in ghmm
   HMMOpen(fileName, type=HMM.FILE_XML)
   HMMFromMatrices(emission_domain, distribution, A, B, pi) # B is a list of distribution parameters

   Examples:

   hmm = HMMOpen('some-hmm.xml')


   hmm.bla ....



"""

import ghmmwrapper
import ghmmhelper
import modhmmer
import re
import StringIO
import copy
import math
import sys
import os
import logging
from string import join
from textwrap import fill

# Initialize logging to stderr
#logging.basicConfig(format="%(asctime)s %(filename)s:%(lineno)d %(levelname)-5s - %(message)s")
log = logging.getLogger("GHMM")

# creating StreamHandler to stderr
hdlr = logging.StreamHandler(sys.stderr)

# setting message format
#fmt = logging.Formatter("%(name)s %(asctime)s %(filename)s:%(lineno)d %(levelname)s %(thread)-5s - %(message)s")
fmt = logging.Formatter("%(name)s %(filename)s:%(lineno)d  - %(message)s")
hdlr.setFormatter(fmt)

# adding handler to logger object
log.addHandler(hdlr)

# Set the minimal severity of a message to be shown. The levels in
# increasing severity are: DEBUG, INFO, WARNING, ERROR, CRITICAL

log.setLevel(logging.ERROR)
log.info( " I'm the ghmm in "+ __file__)

c_log = [log.critical, log.error, log.warning, log.info, log.debug]
def logwrapper(level, message):
    c_log[level](message)

ghmmwrapper.set_pylogging(logwrapper)

# Initialize global random number generator by system time
ghmmwrapper.ghmm_rng_init()
ghmmwrapper.time_seed()


#-------------------------------------------------------------------------------
#- Exceptions ------------------------------------------------------------------

class GHMMError(Exception):
    """Base class for exceptions in this module."""
    def __init__(self, message):
        self.message = message
    def __str__(self):
        return repr(self.message)

class UnknownInputType(GHMMError):
    def __init__(self,message):
       self.message = message
    def __str__(self):
        return repr(self.message)


class NoValidCDataType(GHMMError):
    def __init__(self,message):
        self.message = message
    def __str__(self):
        return repr(self.message)


class badCPointer(GHMMError):
    def __init__(self,message):
        self.message = message
    def __str__(self):
        return repr(self.message)


class SequenceCannotBeBuild(GHMMError):
    def __init__(self,message):
        self.message = message
    def __str__(self):
        return repr(self.message)

class InvalidModelParameters(GHMMError):
    def __init__(self,message):
        self.message = message
    def __str__(self):
        return repr(self.message)

class GHMMOutOfDomain(GHMMError):
    def __init__(self,message):
        self.message = message
    def __str__(self):
        return repr(self.message)

class UnsupportedFeature(GHMMError):
    def __init__(self,message):
       self.message = message
    def __str__(self):
        return repr(self.message)
    
class WrongFileType(GHMMError):
    def __init__(self,message):
        self.message = message
    def __str__(self):
        return repr(self.message)

class ParseFileError(GHMMError):
    def __init__(self,message):
        self.message = message
    def __str__(self):
        return repr(self.message)

#-------------------------------------------------------------------------------
#- constants -------------------------------------------------------------------
kNotSpecified            = ghmmwrapper.kNotSpecified
kLeftRight               = ghmmwrapper.kLeftRight
kSilentStates            = ghmmwrapper.kSilentStates
kTiedEmissions           = ghmmwrapper.kTiedEmissions
kHigherOrderEmissions    = ghmmwrapper.kHigherOrderEmissions
kBackgroundDistributions = ghmmwrapper.kBackgroundDistributions
kLabeledStates           = ghmmwrapper.kLabeledStates
kTransitionClasses       = ghmmwrapper.kTransitionClasses
kDiscreteHMM             = ghmmwrapper.kDiscreteHMM
kContinuousHMM           = ghmmwrapper.kContinuousHMM
kPairHMM                 = ghmmwrapper.kPairHMM
types = {
    kLeftRight:'kLeftRight',
    kSilentStates:'kSilentStates',
    kTiedEmissions:'kTiedEmissions',
    kHigherOrderEmissions:'kHigherOrderEmissions',
    kBackgroundDistributions:'kBackgroundDistributions',
    kLabeledStates:'kLabeledStates',
    kTransitionClasses:'kTransitionClasses',
    kDiscreteHMM:'kDiscreteHMM',
    kContinuousHMM:'kContinuousHMM',
    kPairHMM:'kPairHMM',
    }
#-------------------------------------------------------------------------------
#- EmissionDomain and derived  -------------------------------------------------
class EmissionDomain(object):
    """ Abstract base class for emissions produced by an HMM.

        There can be two representations for emissions:
        1) An internal, used in ghmm.py and the ghmm C-library
        2) An external, used in your particular application

        Example:
        The underlying library represents symbols from a finite,
        discrete domain as integers (see Alphabet).

        EmissionDomain is the identity mapping
    """

    def internal(self, emission):
        """ Given a emission return the internal representation
        """
        return emission


    def internalSequence(self, emissionSequence):
        """ Given a emissionSequence return the internal representation
        """
        return emissionSequence


    def external(self, internal):
        """ Given an internal representation return the
            external representation
        """
        return internal

    def externalSequence(self, internalSequence):
        """ Given a sequence with the internal representation return the
            external representation
        """
        return internalSequence


    def isAdmissable(self, emission):
        """ Check whether emission is admissable (contained in) the domain
            raises GHMMOutOfDomain else
        """
        return None


class Alphabet(EmissionDomain):
    """ Discrete, finite alphabet

    """
    def __init__(self, listOfCharacters, calphabet = None):
        """
        Creates an alphabet out of a listOfCharacters
        @param listOfCharacters: a list of strings (single characters most of
        the time), ints, or other objects that can be used as dictionary keys
        for a mapping of the external sequences to the internal representation
        Can alternatively be an C alphabet_s struct
        
        Note: Alphabets should be considered as imutable. That means the
        listOfCharacters and the mapping should never be touched after
        construction.
        """
        self.index = {} # Which index belongs to which character

        if calphabet is None:
            self.listOfCharacters = listOfCharacters