ghmmunittests.py

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

        w = self.i_seq.getWeight(4)
        self.assertEqual(w,1.0)
        self.i_seq.setWeight(4,4.0)
        w = self.i_seq.getWeight(4)
        self.assertEqual(w,4.0)
        
        w2 = self.d_seq.getWeight(2)
        self.assertEqual(w2,1.0)
        self.d_seq.setWeight(2,7.0)
        w2 = self.d_seq.getWeight(2)
        self.assertEqual(w2,7.0)
        

    def testmerge(self):
        """Merging two SequenceSets   """
        #print"\n************************************ testmerge "
        wrong = 4  # wrong argument type to merge
        self.assertRaises(TypeError,self.i_seq.merge,wrong)

        mseq = ghmm.SequenceSet(self.i_alph,[ [1,4,0,4,5,3],[1,2,3,0] ])        
        self.i_seq.merge(mseq)
        self.assertEqual(len(self.i_seq),7)
        s = "\nNumber of sequences: 7\nSeq 0, length 5, weight 1.0:\n12345\nSeq 1, length 3, weight 1.0:\n030\nSeq 2, length 8, weight 1.0:\n43221111\nSeq 3, length 5, weight 1.0:\n00021\nSeq 4, length 6, weight 1.0:\n111111\nSeq 5, length 6, weight 1.0:\n140453\nSeq 6, length 4, weight 1.0:\n1230"
        self.assertEqual(self.i_seq.verboseStr(),s)

        #print self.i_seq

        d_mseq = ghmm.SequenceSet(self.d_alph,[ [7.5,4.0,1.2],[0.4,0.93,3.3,2.54] ])    
        self.d_seq.merge(d_mseq)
        self.assertEqual(len(self.d_seq),7)
        s2 = "\nNumber of sequences: 7\nSeq 0, length 5, weight 1.0:\n1.5 2.3 3.7 4.1 5.1 \nSeq 1, length 3, weight 1.0:\n0.0 3.1 0.7 \nSeq 2, length 8, weight 1.0:\n4.4 3.05 2.0 2.4 1.2 1.8 1.0 1.0 \nSeq 3, length 5, weight 1.0:\n0.4 0.1 0.33 2.7 1.345 \nSeq 4, length 6, weight 1.0:\n1.0 1.0 1.0 1.0 1.0 1.0 \nSeq 5, length 3, weight 1.0:\n7.5 4.0 1.2 \nSeq 6, length 4, weight 1.0:\n0.4 0.93 3.3 2.54 "
        self.assertEqual(self.d_seq.verboseStr(),s2)        
        
    
    def testgetsubset(self):
        #print"\n----------------- testgetsubset "
        i_subseq = self.i_seq.getSubset([2,1,3])
        s = "\nNumber of sequences: 3\nSeq 0, length 8, weight 1.0:\n43221111\nSeq 1, length 3, weight 1.0:\n030\nSeq 2, length 5, weight 1.0:\n00021"
        self.assertEqual(i_subseq.verboseStr(),s)
        self.assertEqual(len(i_subseq),3)
        self.assertEqual(i_subseq.sequenceLength(0),8)
        
        d_subseq = self.d_seq.getSubset([0,4])
        s2 = "\nNumber of sequences: 2\nSeq 0, length 5, weight 1.0:\n1.5 2.3 3.7 4.1 5.1 \nSeq 1, length 6, weight 1.0:\n1.0 1.0 1.0 1.0 1.0 1.0 "
        self.assertEqual(d_subseq.verboseStr(),s2)
        self.assertEqual(len(d_subseq),2)
        self.assertEqual(d_subseq.sequenceLength(0),5)
        
    def testwrite(self):
        #print"\n-----------------testwrite "
        self.i_seq.write("testdata/ghmmunittests_testwrite.seq") 
        self.d_seq.write("testdata/ghmmunittests_testwrite.seq") 

       
    def testlabelaccess(self):
        #print "\n ----------------- testlabelaccess"
        self.i_seq.getSeqLabel(2)   
        l = self.d_seq.getSeqLabel(3)
        self.assertEqual(l,-1)
        self.d_seq.setSeqLabel(3,8)
        l = self.d_seq.getSeqLabel(3)
        self.assertEqual(l,8)

    def testfilereading(self):
        #print "\n ----------------- testfilereading"
        dom = ghmm.IntegerRange(0,12)
        seqs = ghmm.SequenceSetOpen(dom, 'testdata/d_seq.sqd')
        seqs = ghmm.SequenceSetOpen(self.d_alph, 'testdata/test10.sqd')
        seqs = ghmm.SequenceSetOpen(ghmm.Float(), 'testdata/tiny.txt.sqd')

class HMMBaseClassTests(unittest.TestCase):
    def setUp(self):
        A   = [[0.3,0.3,0.4],[0.6,0.1,0.3],[1.0,0.0,0.0]]
        B   = [[0.0,0.5,0.5,0.0],[0.1,0.0,0.8,0.1], [0.0,0.0,0.0,0.0]]
        pi  = [1.0,0.0,0.0]
        tmp = ghmm.HMMFromMatrices(ghmm.DNA, ghmm.DiscreteDistribution(ghmm.DNA), A, B, pi)
        self.model = ghmm.HMM(ghmm.DNA, ghmm.DiscreteDistribution(ghmm.DNA), tmp.cmodel)

    def testpathPosteriorExeption(self):
        self.assertRaises(NotImplementedError, self.model.pathPosterior, [1,2], 34)

    def teststatePosteriorExeption(self):
        self.assertRaises(NotImplementedError, self.model.statePosterior, "sequence", "state", "time")

    def testposteriorExeption(self):
        self.assertRaises(NotImplementedError, self.model.posterior, "sequence")

    def testjoinedExeption(self):
        self.assertRaises(NotImplementedError, self.model.joined, "emissionSequence", "stateSequence")

    def testbaumWelchExeption(self):
        self.assertRaises(NotImplementedError, self.model.baumWelch, "trainingSequences", "nrSteps", "loglikelihoodCutoff")

    def testbaumWelchSetupExeption(self):
        self.assertRaises(NotImplementedError, self.model.baumWelchSetup, "trainingSequences", "nrSteps")

    def testbaumWelchStepExeption(self):
        self.assertRaises(NotImplementedError, self.model.baumWelchStep, "nrSteps", "loglikelihoodCutoff")
    
    def testbaumWelchDeleteExeption(self):
        self.assertRaises(NotImplementedError, self.model.baumWelchDelete)

    def teststateExeption(self):
        self.assertRaises(NotImplementedError, self.model.state, "stateLabel")

    def testsetEmissionExeption(self):
        self.assertRaises(NotImplementedError, self.model.setEmission, 1, "blah")

    def testasMatricesExeption(self):
        self.assertRaises(NotImplementedError, self.model.asMatrices)

    def testnormalizeExeption(self):
        self.assertRaises(NotImplementedError, self.model.normalize)

    def testrandomizeException(self):
        self.assertRaises(NotImplementedError, self.model.randomize, "noiseLevel")


class DiscreteEmissionHMMTests(unittest.TestCase):
    def setUp(self):
        log.debug("DiscreteEmissionHMMTests.setUp() -- begin")
        self.A = [[0.3,0.3,0.4],[0.6,0.1,0.3],[1.0,0.0,0.0]]
        self.B = [[0.0,0.5,0.5,0.0],[0.1,0.0,0.8,0.1], [0.0,0.0,0.0,0.0]]
        self.pi = [1.0,0.0,0.0]
        self.model = ghmm.HMMFromMatrices(ghmm.DNA,ghmm.DiscreteDistribution(ghmm.DNA), self.A, self.B, self.pi)
        log.debug("DiscreteEmissionHMMTests.setUp() -- end")

    def test__str__(self):
        log.debug("test__str__ -- begin")
        # we aren't interested in the output but the function should run fine
        str(self.model)

    def testAccessFunctions(self):
        log.debug("testAccessFunctions -- begin")

        self.assertEqual(self.model.N,3)
        self.assertEqual(self.model.M,4)
        
        pi = self.model.getInitial(2)
        self.assertEqual(pi,0)
        self.model.setInitial(2,0.5,fixProb=1)
        pi = self.model.getInitial(2)
        self.assertEqual(pi,0.5)
        
        trans = self.model.getTransition(0,1)
        self.assertEqual(trans, 0.3)
        self.model.setTransition(0,1,0.6)
        trans = self.model.getTransition(0,1)
        self.assertEqual(trans, 0.6)
        
        emission = self.model.getEmission(1)
        self.assertEqual(emission, [0.1, 0.0, 0.8, 0.1] )
        
        # introducing silent state
        self.model.setEmission(1,[0.0,0.0,0.0,0.0])
        emission = self.model.getEmission(1)
        self.assertEqual(emission,[0.0,0.0,0.0,0.0] ) 
        self.assertEqual(self.model.cmodel.model_type,260)
        self.assertEqual(self.model.getSilentFlag(1),1)
        
        # removing silent state
        self.model.setEmission(1,[0.2,0.2,0.2,0.4])
        emission = self.model.getEmission(1)
        self.assertEqual(emission,[0.2,0.2,0.2,0.4] )  
        self.assertEqual(self.model.cmodel.model_type,260)
        self.assertEqual(self.model.getSilentFlag(1),0)
        
        # removing last silent state
        self.model.setEmission(2,[0.25,0.25,0.25,0.25])
        emission = self.model.getEmission(2)
        self.assertEqual(emission,[0.25,0.25,0.25,0.25])  
        self.assertEqual(self.model.cmodel.model_type,256)
        self.assertEqual(self.model.getSilentFlag(2),0)
    
        # inserting silent state
        self.model.setEmission(2,[0.0,0.0,0.0,0.0])
        emission = self.model.getEmission(2)
        self.assertEqual(emission,[0.0,0.0,0.0,0.0])  
        self.assertEqual(self.model.cmodel.model_type,260)
        self.assertEqual(self.model.getSilentFlag(2),1)

    def testNewXML(self):
        log.debug("testNewXML -- begin")
        model = ghmm.HMMOpen('../doc/xml_example.xml')
    
    def getModel(self):
        A  = [[0.3, 0.6,0.1],[0.0, 0.5, 0.5],[0.0,0.0,1.0]]
        B  = [[0.5, 0.5],[0.5,0.5],[1.0,0.0]]
        pi = [1.0, 0.0, 0.0]
        return ghmm.HMMFromMatrices(ghmm.IntegerRange(0,2),
                                    ghmm.DiscreteDistribution(ghmm.IntegerRange(0,2)),
                                    A, B, pi)
        
    def testDel(self):
        """  test for explicit construction and destruction """
        log.debug("testDel -- begin")
        del self.model
        for i in range(100):
            mo = self.getModel()
    
    def testAsMatrices(self):
        log.debug("testAsMatrices -- begin")
        tA,tB,tpi = self.model.asMatrices()
        
        self.assertEqual(self.A,tA)
        self.assertEqual(self.B,tB)
        self.assertEqual(self.pi,tpi)        
            
    def testSample(self):
        log.debug("testSample -- begin")
        seq = self.model.sampleSingle(100, seed=3586662)
        #print seq
        seq2 = self.model.sample(10, 100, seed=3586662)


    def testBaumWelch(self):
        log.debug("testBaumWelch -- begin")
        seq = self.model.sample(100,100,seed=3586662)
        self.assertRaises(NotImplementedError, self.model.baumWelch, seq, 5, 0.01)
        
        self.model.setEmission(2,[0.25,0.25,0.25,0.25])
        self.assertEqual(self.model.cmodel.model_type & 4, 0)
        self.model.baumWelch(seq,5,0.01)
        self.model.baumWelch(seq)
        
    def testViterbi(self):
        log.debug("testViterbi -- begin")
        # Caution with an even number of consecutives 'g'
        # the model can produce two equal probable paths
        
        f = lambda x: round(x, 13)
        g = lambda x: map(f, x)

        seq = ghmm.EmissionSequence(ghmm.DNA, [ 'c','c','g','c','c','g','g','g','g','g','c','g','g','g','c' ])
                                               #[0,  2,  0,  1,  0,  2,  0,  1,  0,  1,  0,  1,  0,  1,  0]
                                               #[0,  2,  0,  1,  0,  2,  0,  1,  0,  1,  0,  1,  0,  1,  0, 1, 0]
                
        result = self.model.viterbi(seq)
        trueResult =  ([0, 2, 0, 1, 0, 2, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0], -19.698557965224637)
        
        self.assertEqual (result[0], trueResult[0])
        self.assertEqual (f(result[1]), f(trueResult[1]))
    
        
        seq2 = ghmm.SequenceSet(ghmm.DNA, [['c','c','g','c','c','g','g','g','g','g','c','g','g','g','c'],
                                           ['c','g','g','g','g','g','c','t','g','c','g','g','t','c','c'],
                                           ['g','g','c','g','c','c','g','c','c','c','c','g','g','g','t'],
                                           ['g','c','c','g','c','c','g','g','g','c','c','c','g','g','g'],
                                           ['g','g','g','g','c','g','c','a','g','g','c','c','g','g','g'],
                                           ['g','g','c','g','g','t','g','c','c','c','c','g','g','a','a'],
                                           ['c','g','g','g','c','c','t','g','c','g','c','g','a','g','g'],
                                           ['c','c','g','g','g','g','g','g','g','c','g','c','g','c','g'],
                                           ['g','g','g','g','c','c','g','g','g','c','g','c','g','g','g'],
                                           ['g','c','t','c','g','g','a','g','g','c','a','g','g','g','g']])
        
        
        path2 = self.model.viterbi(seq2)
        truePath2 = ([[0, 2, 0, 1, 0, 2, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0],
                      [0, 1, 0, 1, 0, 1, 0, 1, 0, 2, 0, 1, 0, 1, 0, 2, 0],
                      [0, 1, 0, 1, 0, 2, 0, 1, 0, 2, 0, 2, 0, 2, 0, 2, 0, 1, 0, 1],
                      [0, 2, 0, 2, 0, 1, 0, 2, 0, 1, 0, 1, 0, 2, 0, 2, 0, 1, 0, 1],
                      [0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 2, 0, 1, 0, 1],
                      [0, 1, 0, 1, 0, 1, 0, 2, 0, 2, 0, 2, 0, 2, 0, 1, 0, 1, 1],
                      [0, 1, 0, 1, 0, 2, 0, 1, 0, 2, 0, 1, 0, 2, 0, 1, 0, 1],
                      [0, 2, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1],
                      [0, 1, 0, 1, 0, 2, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1],
                      [0, 2, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1]],
                      [-19.698557965224637, -23.857441048584306, -22.77125127991475,
                       -20.691809738234902, -21.595677950110499, -27.846425095148575,
                       -24.262906156692473, -19.516236408430682, -19.5162364084307,
                       -25.754561033470189])
                
        self.assertEqual (path2[0], truePath2[0])    
        self.assertEqual (g(path2[1]), g(truePath2[1]))    
        

    def testLoglikelihood(self):
        log.debug("testLoglikelihood -- begin")
        seq = self.model.sampleSingle(100,seed=3586662)
        logp = self.model.loglikelihood(seq)
        self.assert_(logp-93.1053904716 < 10^-8, "Different results in loglikelihood ")
        log.debug("testLoglikelihood -- end")

    def testLogProb(self):
        log.debug("testLogProb -- begin")
        seq = self.model.sampleSingle(15,seed=3586662)
        path,logp = self.model.viterbi(seq)
        logp = self.model.joined(seq,path)
        self.assert_(logp - 22.4303246929 < 10^-8, "Different results in logprob ")
        
    def testFoBa(self):
        log.debug("testFoBa -- begin")
        #print self.model
        seq = ghmm.EmissionSequence(self.model.emissionDomain,
                                    ['g','g','g','c','t','g','g','c','g','g',
                                     'g','c','g','g','g','c','c','c','g','c',
                                     'g','g','c','c','g','c','g','c','c','c',
                                     'g','c','g','c','g','g','c','t','c','c'])

        (alpha,scale) = self.model.forward(seq)       

        f = lambda x: round(x, 13)
        g = lambda x: map(f, x)
        
        talpha = [[0.7142857142857143, 0.0, 0.28571428571428575], [0.43640897755610975, 0.29925187032418954, 0.26433915211970077], [0.50453092851201675, 0.22588976929475116, 0.26957930219323206], [0.7142857142857143, 0.0, 0.2857142857142857], [0.0, 0.76923076923076916, 0.23076923076923075], [0.61307901907356943, 0.10899182561307905, 0.27792915531335155], [0.46113149992850677, 0.27262761546160807, 0.26624088460988515], [0.7142857142857143, 0.0, 0.28571428571428575], [0.43640897755610975, 0.29925187032418954, 0.26433915211970077], [0.50453092851201675, 0.22588976929475116, 0.26957930219323206], [0.48775986167826385, 0.24395091819263898, 0.26828922012909723], [0.71428571428571419, 0.0, 0.2857142857142857], [0.43640897755610969, 0.29925187032418948, 0.26433915211970072], [0.50453092851201675, 0.22588976929475116, 0.26957930219323206], [0.48775986167826385, 0.24395091819263898, 0.26828922012909723], [0.71428571428571419, 0.0, 0.2857142857142857], [0.7142857142857143, 0.0, 0.2857142857142857], [0.7142857142857143, 0.0, 0.28571428571428575], [0.43640897755610975, 0.29925187032418954, 0.26433915211970077], [0.7142857142857143, 0.0, 0.28571428571428575], [0.43640897755610975, 0.29925187032418954, 0.26433915211970077], [0.50453092851201675, 0.22588976929475116, 0.26957930219323206], [0.7142857142857143, 0.0, 0.2857142857142857], [0.7142857142857143, 0.0, 0.28571428571428575], [0.43640897755610975, 0.29925187032418954, 0.26433915211970077], [0.7142857142857143, 0.0, 0.28571428571428575], [0.43640897755610975, 0.29925187032418954, 0.26433915211970077], [0.7142857142857143, 0.0, 0.28571428571428575], [0.7142857142857143, 0.0, 0.28571428571428575], [0.7142857142857143, 0.0, 0.28571428571428575], [0.43640897755610975, 0.29925187032418954, 0.26433915211970077], [0.7142857142857143, 0.0, 0.28571428571428575], [0.43640897755610975, 0.29925187032418954, 0.26433915211970077], [0.7142857142857143, 0.0, 0.28571428571428575], [0.43640897755610975, 0.29925187032418954, 0.26433915211970077], [0.50453092851201675, 0.22588976929475116, 0.26957930219323206], [0.7142857142857143, 0.0, 0.2857142857142857], [0.0, 0.76923076923076916, 0.23076923076923075], [0.71428571428571419, 0.0, 0.2857142857142857], [0.7142857142857143, 0.0, 0.2857142857142857]]