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)13cViterbi::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()32cViterbi::~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)42void 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