1 /**************************************************************
2  *** RHmm package
3  ***
4  *** File: cViterbi.cpp
5  ***
6  *** Author: Ollivier TARAMASCO <Ollivier.Taramasco@imag.fr>
7  *** Author: Sebastian BAUER <sebastian.bauer@charite.de>
8  ***
9  **************************************************************/
10 
11 #include "StdAfxRHmm.h"
12 
cViterbi(cInParam & theInParam)13 cViterbi::cViterbi(cInParam& theInParam)
14 {       MESS_CREAT("cViterbi")
15         if (theInParam.mNSample > 0)
16         {
17         register uint n ;
18                 mSeq = new uint*[theInParam.mNSample] ;
19                 for (n = 0 ; n < theInParam.mNSample ; n++)
20                         mSeq[n] = new uint[theInParam.mY[n].mSize] ;
21 
22                 mLogProb.ReAlloc(theInParam.mNSample) ;
23                 for (n = 0 ; n < theInParam.mNSample ; n++)
24                         mLogProb[n] = -1e100 ;
25         }
26         else
27         {       mSeq = NULL ;
28                 mLogProb.Delete() ;
29         }
30 }
31 
~cViterbi()32 cViterbi::~cViterbi()
33 {       MESS_DESTR("cViterbi")
34         if (mLogProb.mSize > 0)
35         {       for (register uint n = 0 ; n < mLogProb.mSize ; n++)
36                         delete [] mSeq[n] ;
37                 delete [] mSeq ;
38                 mLogProb.Delete() ;
39         }
40 }
41 
ViterbiPath(cInParam & theInParam,cHmm & theHMM)42 void cViterbi::ViterbiPath(cInParam &theInParam, cHmm &theHMM)
43 {
44 int             myIndAux                                                ;
45 double  myAux,
46                 myAux1                                                  ;
47 uint    myNSample = theInParam.mNSample ;
48 
49 cDMatrix* myProbaCond = new cDMatrix[myNSample] ;
50         for (register uint n = 0 ; n < myNSample ; n++)
51         {
52         uint mySize = theInParam.mY[n].mSize/theInParam.mDimObs ;
53                 myProbaCond[n].ReAlloc(theInParam.mNClass, mySize) ;
54         }
55 
56 cDVector* myDelta = new cDVector[theInParam.mNClass] ;
57 int** myPsi = new int*[theInParam.mNClass] ;
58         theHMM.mDistrParam->ComputeCondProba(theInParam.mY, myNSample, myProbaCond) ;
59         for (register uint n = 0 ; n < myNSample ; n++)
60         {
61         uint mySize = theInParam.mY[n].mSize/theInParam.mDimObs ;
62                 for (register uint j = 0 ; j < theInParam.mNClass ; j++)
63                 {       myPsi[j] = new int[mySize] ;
64                         myDelta[j].ReAlloc(mySize) ;
65                 }
66         // Initialization
67                 for (register uint i = 0 ; i < theInParam.mNClass  ; i++)
68                 {       myDelta[i][0] = log(theHMM.mInitProba[i]) + log(myProbaCond[n][i][0]) ;
69                         myPsi[i][0] = 0 ;
70                 }
71 
72         // Recursion
73                 for (register int t = 0 ; t < (int)mySize -1 ; t++)
74                 {
75                         for (register uint j = 0 ; j < theInParam.mNClass  ; j++)
76                         {       myAux = myDelta[0][t] + log(theHMM.mTransMatVector[t][0][j]) ;
77                                 myIndAux = 0 ;
78                                 for (register uint i = 1 ; i < theInParam.mNClass ; i++)
79                                 {
80                                         if ((myAux1 = myDelta[i][t] + log(theHMM.mTransMatVector[t][i][j])) > myAux)
81                                         {       myAux = myAux1 ;
82                                                 myIndAux = i ;
83                                         }
84                                 }
85                                 myDelta[j][t+1] = myAux + log(myProbaCond[n][j][t+1]) ;
86                                 myPsi[j][t+1] = myIndAux ;
87                         }
88                 }
89         // Terminaison
90                 mLogProb[n] = myDelta[0][mySize-1] ;
91                 mSeq[n][mySize-1] = 0 ;
92                 for (register uint i = 1 ; i < theInParam.mNClass ; i++)
93                 {       if (myDelta[i][mySize-1] > mLogProb[n])
94                         {       mLogProb[n] = myDelta[i][mySize-1] ;
95                                 mSeq[n][mySize-1] = i ;
96                         }
97                 }
98 
99                 for (register int t =  (int)(mySize-2) ; t >= 0 ; t--)
100                         mSeq[n][t] = myPsi[mSeq[n][t+1]][t+1] ;
101 
102 /*              for (register uint j = 0 ; j < theInParam.mNClass ; j++)
103                 {       myPsi[j] = new int[theInParam.mY[n].mSize] ;
104                         myDelta[j].ReAlloc(theInParam.mY[n].mSize) ;
105                 }
106 */
107                 for (register uint j = 0 ; j < theInParam.mNClass ; j++)
108                 {
109                         delete [] myPsi[j] ;
110                         myDelta[j].Delete() ;
111                 }
112         }
113         for (register uint n = 0 ; n < myNSample ; n++)
114                 myProbaCond[n].Delete() ;
115 //      delete myPsi ;
116 //      delete myDelta ;
117 //      delete myProbaCond ;
118 }
119 
120