test_ghmm.py

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

from ghmm import *
################################################################################
#
#       This file is part of the General Hidden Markov Model Library,
#       GHMM version 0.8_beta1, see http://ghmm.org
#
#       file:    test_ghmm.py
#       authors: Benjamin Georgi
#
#       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
#
#
#       This file is version $Revision: 1791 $
#                       from $Date: 2006-11-15 15:09:43 +0100 (Wed, 15 Nov 2006) $
#             last change by $Author: grunau $.
#
################################################################################
print "GHMM Tests:"

print "\n*** EmissionSequence / SequenceSet ***"

alph = IntegerRange(0,7)

print "- Discrete case"
print "\nlist input EmssionSequence:"
s = EmissionSequence(alph,[1,2,0,0,0,3,4])
print s

#print "list(EmissionSequence)"
#l = list(s)
#print l

print "\n\nlist input SequenceSet:"
s2 = SequenceSet(alph,[ [1,2,3,4,5],[0,3,0],[4,3,2,2,1,1,1,1], [0,0,0,2,1],[1,1,1,1,1,1] ])
print s2

print "\nwriting SequenceSet to file."
s2.write("test_ghmm_d.seq")

print "\nSequenceSet.getSubset"
s4 = s2.getSubset([0,2,4])
print s4

print "\npointer input EmissionSequence"
s3 = s2[1] # call to SequenceSet.__getitem__
print s3


print "- Continous case"

alph = Float()

print "\nlist input EmssionSequence:"
s = EmissionSequence(alph,[1.3,2.1,0.8,0.1,0.03,3.6,43.3])
print s

#print "list(EmissionSequence)"
#l = list(s)
#print l

print "\n\nlist input SequenceSet:"
s2 = SequenceSet(alph,[ [1.5,2.3,3.7,4.1,5.1],[0.0,3.1,0.7],[4.4,3.05,2.0,2.4,1.2,1.8,1.0,1.0], [0.4,0.1,0.33,2.7,1.345],[1.0,1.0,1.0,1.0,1.0,1.0] ])
print s2

print "\nwriting SequenceSet to file."
s2.write("test_ghmm_c.seq")

print "\nSequenceSet.getSubset"
s4 = s2.getSubset([1,3,0])
print s4

print "\npointer input EmissionSequence"
s3 = s2[1] # call to SequenceSet.__getitem__
print s3



print "\n*** Discrete Emission Model ***"

m = HMMFromMatrices(DNA,DiscreteDistribution(DNA),
                       [[0.3,0.3,0.4],[0.6,0.1,0.3],[1.0,0.0,0.0]],
                       [[0.0,0.5,0.5,0.0],[0.1,0.0,0.8,0.1], [0.0,0.0,0.0,0.0]],
                       [1.0,0,0])

print m

trans = m.getTransition(0,1)
print "a[0,1] = " + str(trans)

emission = m.getEmission(1)
print emission

print "\nSample:"					   
s4 = m.sample(4,15)
print str(s4) + "\n"

print "\nSampleSingle :"					   
s5 = m.sampleSingle(10)
print str(s5) + "\n"

print "\nMerging two sequences:"
s4.merge(s5)
print s4

print "\nTraining model"
m.baumWelch(s5)
print m

print "\nViterbi:"
path = m.viterbi(s5)
print str(path) + "\n"

print "\nForward"
logp1 = m.loglikelihood(s5)
print "logp = " + str(logp1) + "\n"

print "\nlogprob:"
logp2 = m.logprob(s5,[0, 2, 0, 2, 0, 1, 0, 1, 0, 1, 0, 1])
print logp2, " -> " + str(2.71**logp2)

print "\nForward matrices"
(alpha,scale) = m.forward(s5)
print "alpha:\n" + str(alpha) + "\n"
print "scale = " + str(scale) + "\n"	

beta = m.backward(s5,scale)
print "beta = \n " + str(beta) + "\n"


print "\n\n\n *** Gaussian Model ***"
F = Float()            		   
m2 = HMMFromMatrices(F,GaussianDistribution(F),
                         [[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,0])

m4 = HMMFromMatrices(F,GaussianDistribution(F),
                         [[0.0,1.0,0],[0.5,0.0,0.5],[0.3,0.3,0.4]],
                         [[0.0,1.3],[-1.0,0.1], [1.0,0.6]],
                         [1.0,0,0])

print m2
trans = m2.getTransition(2,0)
print "a[2,2] = " + str(trans)
                         
print "\nSample:"
cs1 = m2.sample(4,15)                         
print str(cs1) + "\n"

print "\nSampleSingle:"
cs2 = m2.sampleSingle(10)                         
print str(cs2) + "\n"

print "\nget subset"
cs3 = cs1.getSubset([0,2])
print cs3

print "\nViterbi"
spath = m2.viterbi(cs1)
print str(spath) + "\n"

print "\nForward"
logp = m2.loglikelihood(cs1)    
print "logp = " + str(logp) + "\n"

print "\nForward matrices"
(salpha,sscale) = m2.forward(cs2)
print "alpha:\n" + str(salpha) + "\n"
print "scale = " + str(sscale) + "\n"	

print "\nBackward matrix"
beta = m2.backward(cs2,sscale)
print "beta = \n " + str(beta) + "\n"

print "Reading SequenceSet from .sqd file"
l = SequenceSetOpen(F,"seq_test.sqd")
print l

print "Model distances (continous):"
d = m2.distance(m4,1000)
print "distance= " + str(d)


print "*** Normalizing  discrete model***"
m3 = HMMFromMatrices(DNA,DiscreteDistribution(DNA),
                     [[10.0,10.0,10.0],[0.0,0.0,100.0],[25.0,25.0,50.0]],
                     [[10.0,0.0,010.0,0.0],[0.0,3.5,3.5,0.0], [5.0,5.0,5.0,5.0]],
                     [1.0,0,0])

#print m3

m3.normalize()

print "\nModel distance (discrete):"
d = m3.distance(m,1000) 
print "distance = " + str(d)

#print m3

print "\nWriting to file:"
m.write("er.log")
m3.write("er.log")
m2.write("er2.log")

mList = [m,m2,m3]
HMMwriteList("er.log",mList)


print "\n\n\n *** GaussianMixture Model ***"
F = Float()            		   
m_mix = HMMFromMatrices(F,GaussianMixtureDistribution(F),
                         [[0.0,1.0,0],[0.5,0.0,0.5],[0.3,0.3,0.0]],
                         
                         [  [ [10.0,40.0], [1,1], [0.5,0.5]  ],
                            [ [-10.0,-0.5], [2,2], [0.5,0.5] ],
                            [ [1.0,0.2], [3,3],  [0.5,0.5]  ]
                         ],

                         [1.0,0,0])

print m_mix


trans = m2.getTransition(2,0)
print "a[2,0] = " + str(trans)
                         
print "\nSample:"
mseq1 = m_mix.sample(4,15)                         
print str(mseq1) + "\n"

print "\nSampleSingle:"
mseq2 = m_mix.sampleSingle(10)                         
print str(mseq2) + "\n"

print "\nViterbi"
mixpath = m_mix.viterbi(mseq2)
print str(mixpath) + "\n"


print "\nForward"
logp = m_mix.loglikelihood(mseq1)    
print "logp = " + str(logp) + "\n"

print "\nForward matrices"
(salpha,sscale) = m_mix.forward(mseq2)
print "alpha:\n" + str(salpha) + "\n"
print "scale = " + str(sscale) + "\n"	

print "\nBackward matrices"
beta = m_mix.backward(mseq2,sscale)
print "beta = \n " + str(beta) + "\n"



#print "*** Reading XML-HMM from file."
#HMMOpen = HMMOpenFactory(GHMM_FILETYPE_XML)
#m4 = HMMOpen('/project/algorithmics/Sopra/python/simple2.xml')
#print m4