1 /**************************************************************
2 *** RHmm package
3 ***
4 *** File: cDerivative.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
cDerivative(uint theNSample,uint * theT,uint theNClass,uint theNFreeParam)13 cDerivative::cDerivative(uint theNSample, uint* theT, uint theNClass, uint theNFreeParam)
14 {
15 mvNFreeParam = theNFreeParam ;
16 mvNClass = theNClass ;
17 mvNSample = theNSample ;
18 mvT = new uint[theNSample] ;
19
20 mPsi = new cDVector**[theNSample] ;
21 mOmega = new cDMatrix**[theNSample] ;
22 mScore = new cDVector[theNSample] ;
23 mInformation = new cDMatrix[theNSample] ;
24
25 for (register uint n = 0 ; n < theNSample ; n++)
26 { mPsi[n] = new cDVector*[mvNClass] ;
27 mOmega[n] = new cDMatrix*[mvNClass] ;
28 mScore[n].ReAlloc(mvNFreeParam) ;
29 mvT[n] = theT[n] ;
30 mInformation[n].ReAlloc(mvNFreeParam, mvNFreeParam) ;
31 for (register uint i = 0 ; i < mvNClass ; i++)
32 { mPsi[n][i] = new cDVector[theT[n]] ;
33 mOmega[n][i] = new cDMatrix[theT[n]] ;
34 for (register uint t = 0 ; t < mvT[n] ; t++)
35 { mPsi[n][i][t].ReAlloc(mvNFreeParam, 0.0) ;
36 mOmega[n][i][t].ReAlloc(mvNFreeParam, mvNFreeParam) ;
37 }
38 }
39 }
40 MESS_CREAT("cDerivative")
41 }
42
cDerivative(const cInParam & theInParam,uint theNFreeParam)43 cDerivative::cDerivative(const cInParam &theInParam, uint theNFreeParam)
44 {
45 mvNFreeParam = theNFreeParam ;
46 mvNClass = theInParam.mNClass ;
47 mvNSample = theInParam.mNSample ;
48 mvT = new uint[mvNSample] ;
49
50 mPsi = new cDVector**[mvNSample] ;
51 mOmega = new cDMatrix**[mvNSample] ;
52 mScore = new cDVector[mvNSample] ;
53 mInformation = new cDMatrix[mvNSample] ;
54
55 for (register uint n = 0 ; n < mvNSample ; n++)
56 {
57 uint myT = theInParam.mY[n].GetSize() / theInParam.mDimObs ;
58 mPsi[n] = new cDVector*[mvNClass] ;
59 mOmega[n] = new cDMatrix*[mvNClass] ;
60 mScore[n].ReAlloc(mvNFreeParam) ;
61 mvT[n] = myT ;
62 mInformation[n].ReAlloc(mvNFreeParam, mvNFreeParam) ;
63 for (register uint j = 0 ; j < mvNClass ; j++)
64 { mPsi[n][j] = new cDVector[myT] ;
65 mOmega[n][j] = new cDMatrix[myT] ;
66 for (register uint t = 0 ; t < myT ; t++)
67 { mPsi[n][j][t].ReAlloc(mvNFreeParam) ;
68 mOmega[n][j][t].ReAlloc(mvNFreeParam, mvNFreeParam) ;
69 }
70 }
71 }
72 MESS_CREAT("cDerivative")
73 }
74
~cDerivative()75 cDerivative::~cDerivative()
76 {
77 for (register uint n = 0 ; n < mvNSample ; n++)
78 { for (register uint j = 0 ; j < mvNClass ; j++)
79 { for (register uint t = 0 ; t < mvT[n] ; t++)
80 { mPsi[n][j][t].Delete() ;
81 mOmega[n][j][t].Delete() ;
82 }
83 delete [] mPsi[n][j] ;
84 delete [] mOmega[n][j] ;
85 }
86 delete [] mPsi[n] ;
87 delete [] mOmega[n] ;
88 mScore[n].Delete() ;
89 mInformation[n].Delete() ;
90 }
91 delete [] mPsi ;
92 delete [] mOmega ;
93 delete [] mScore ;
94 delete [] mInformation ;
95 MESS_DESTR("cDerivative")
96 }
97
ComputeDerivative(cHmm & theHmm,cInParam & theInParam)98 void cDerivative::ComputeDerivative(cHmm& theHmm, cInParam& theInParam)
99 {
100 cDMatrix* myCondProba = new cDMatrix[theInParam.mNSample] ;
101
102 for (register uint n = 0 ; n < mvNSample ; n++)
103 {
104 uint mySize = theInParam.mY[n].mSize/theInParam.mDimObs ;
105 myCondProba[n].ReAlloc(theInParam.mNClass, mySize) ;
106 }
107 theHmm.mDistrParam->ComputeCondProba(theInParam.mY, mvNSample, myCondProba) ;
108
109 cDMatrix* myLambda = new cDMatrix[theInParam.mNSample] ;
110 cDVector* mySumLambda = new cDVector[theInParam.mNSample] ;
111
112 // CALCUL DES LAMBDAs
113 for (register uint n = 0 ; n < mvNSample ; n++)
114 {
115 mySumLambda[n].ReAlloc(mvT[n]) ;
116 myLambda[n].ReAlloc(mvNClass, mvT[n]) ;
117
118 // t = 0
119 mySumLambda[n][0] = 0.0 ;
120 for (register uint j = 0 ; j < mvNClass ; j++)
121 {
122 myLambda[n][j][0] = theHmm.mInitProba[j] * myCondProba[n][j][0] ;
123 mySumLambda[n][0] += myLambda[n][j][0] ;
124 }
125
126 // t > 0
127 uint myT = theInParam.mY[n].GetSize() / theInParam.mDimObs ;
128 for (register uint t = 1 ; t < myT ; t++)
129 { mySumLambda[n][t] = 0.0 ;
130 for (register uint j = 0 ; j < mvNClass ; j++)
131 {
132 double myAux = 0.0 ;
133 for (register uint i = 0 ; i < mvNClass ; i++)
134 myAux += myLambda[n][i][t-1]*theHmm.mTransMatVector[t][i][j] ;
135 myLambda[n][j][t] = myAux * myCondProba[n][j][t]/mySumLambda[n][t-1] ;
136 mySumLambda[n][t] += myLambda[n][j][t] ;
137 }
138 }
139 }
140
141 cDVector* myGradInitProb = new cDVector[mvNClass] ;
142 cDVector** myGradTransMat = new cDVector*[mvNClass] ;
143 cDVector** myGradCondProba = new cDVector*[mvNClass] ;
144 cDMatrix** myHessCondProba = new cDMatrix*[mvNClass] ;
145
146 for (register uint n = 0 ; n < mvNClass ; n++)
147 { myGradInitProb[n].ReAlloc(mvNFreeParam, 0.0L) ;
148 myGradTransMat[n] = new cDVector[mvNClass] ;
149 for (register uint i = 0 ; i < mvNClass ; i++)
150 myGradTransMat[n][i].ReAlloc(mvNFreeParam, 0.0L) ;
151 }
152
153 /* D�riv�es probabilit�s initiales */
154 uint myNFreeClass = mvNClass - 1 ;
155 for (register uint s = 0 ; s < myNFreeClass ; s++)
156 { myGradInitProb[s][s] = 1.0L ;
157 myGradInitProb[myNFreeClass][s] = -1.0L ;
158 }
159
160 /* D�riv�es matrice de transition */
161 uint myBegIndex = myNFreeClass ;
162 for (register uint i = 0 ; i < mvNClass ; i++)
163 { for (register uint j = 0 ; j < myNFreeClass ; j++)
164 { myGradTransMat[i][j][j+myBegIndex] = 1.0L ;
165 myGradTransMat[i][myNFreeClass][j+myBegIndex] = -1.0L ;
166 }
167 myBegIndex += myNFreeClass ;
168 }
169
170 for (register uint n = 0 ; n < mvNSample ; n++) // Boucle sur le nombre d'�chantillon
171 {
172 for (register uint j = 0 ; j < mvNClass ; j++)
173 { myGradCondProba[j] = new cDVector[mvT[n]] ;
174 myHessCondProba[j] = new cDMatrix[mvT[n]] ;
175 for (register uint t = 0 ; t < mvT[n] ; t++)
176 { myGradCondProba[j][t].ReAlloc(mvNFreeParam) ;
177 myHessCondProba[j][t].ReAlloc(mvNFreeParam, mvNFreeParam) ;
178 }
179 }
180 theHmm.mDistrParam->ComputeDerivative(theInParam.mY[n], myGradCondProba, myHessCondProba) ;
181 /* t = 0 */
182
183 for (register uint j = 0 ; j < mvNClass ; j++)
184 {
185 mPsi[n][j][0] = theHmm.mInitProba[j] * myGradCondProba[j][0] + myCondProba[n][j][0] * myGradInitProb[j] ;
186 mOmega[n][j][0] = theHmm.mInitProba[j] * myHessCondProba[j][0] + myGradCondProba[j][0] * Transpose(myGradInitProb[j])
187 + myGradInitProb[j] * Transpose(myGradCondProba[j][0]) ;
188 }
189
190 /* t > 0 */
191 uint myT = mvT[n] ;
192
193 for (register uint t = 1 ; t < myT ; t++)
194 { for (register uint j = 0 ; j < mvNClass ; j++)
195 { mPsi[n][j][t] = 0.0 ;
196 mOmega[n][j][t] = 0.0 ;
197 for ( register uint i = 0 ; i < mvNClass ; i++)
198 {
199 cDVector myVect1 = myCondProba[n][j][t] * theHmm.mTransMatVector[t][i][j] * mPsi[n][i][t-1] ;
200 cDVector myVect2 = myLambda[n][i][t-1] * theHmm.mTransMatVector[t][i][j] * myGradCondProba[j][t] ;
201 cDVector myVect3 = myLambda[n][i][t-1] * myCondProba[n][j][t] * myGradTransMat[i][j] ;
202
203 mPsi[n][j][t] += myVect1 + myVect2 + myVect3 ;
204
205 cDMatrix myMat1 = myCondProba[n][j][t] * theHmm.mTransMatVector[0][i][j] * mOmega[n][i][t-1] ;
206 cDMatrix myMat2 = theHmm.mTransMatVector[t][i][j] * (mPsi[n][i][t-1] * Transpose(myGradCondProba[j][t]) + myGradCondProba[j][t] * Transpose(mPsi[n][i][t-1])) ;
207 cDMatrix myMat3 = myCondProba[n][j][t] * (mPsi[n][i][t-1] * Transpose(myGradTransMat[i][j]) + myGradTransMat[i][j] * Transpose(mPsi[n][i][t-1])) ;
208 cDMatrix myMat4 = myLambda[n][i][t-1] * (myGradCondProba[j][t] * Transpose(myGradTransMat[i][j]) + myGradTransMat[i][j] * Transpose(myGradCondProba[j][t])) ;
209 cDMatrix myMat5 = myLambda[n][i][t-1] * theHmm.mTransMatVector[t][i][j] * myHessCondProba[j][t] ;
210 mOmega[n][j][t] += myMat1 + myMat2 + myMat3 + myMat4 + myMat5 ;
211 }
212 mPsi[n][j][t] /= mySumLambda[n][t-1] ;
213 mOmega[n][j][t] /= mySumLambda[n][t-1] ;
214 }
215 } /* for t */
216 mScore[n] = 0.0 ;
217 mInformation[n] = 0.0 ;
218 for ( register uint j = 0 ; j < mvNClass ; j++)
219 { mScore[n] += mPsi[n][j][myT-1] ;
220 mInformation[n] -= mOmega[n][j][myT-1] ;
221 }
222 mScore[n] /= mySumLambda[n][myT-1] ;
223 mInformation[n] /= mySumLambda[n][myT-1] ;
224 mInformation[n] += mScore[n] * Transpose(mScore[n]) ;
225
226 for (register uint p = 0 ; p < mvNClass ; p++)
227 { for (register uint t = 0 ; t < mvT[n] ; t++)
228 { myGradCondProba[p][t].Delete() ;
229 myHessCondProba[p][t].Delete() ;
230 }
231 delete [] myGradCondProba[p] ;
232 delete [] myHessCondProba[p] ;
233 }
234
235 } /* for n */
236
237 for (register uint n = 0 ; n < mvNClass ; n++)
238 { myGradInitProb[n].Delete() ;
239 for (register uint j = 0 ; j < mvNClass ; j++)
240 myGradTransMat[n][j].Delete() ;
241 delete [] myGradTransMat[n] ;
242 }
243
244 for (register uint n = 0 ; n < mvNSample ; n++)
245 { myLambda[n].Delete() ;
246 mySumLambda[n].Delete() ;
247 }
248
249 delete [] myGradCondProba ;
250 delete [] myHessCondProba ;
251 delete [] myGradInitProb ;
252 delete [] myGradTransMat ;
253 delete [] myLambda ;
254 delete [] mySumLambda ;
255 }
256
ComputeScoreAndInformation(cDVector & theScore,cDMatrix & theInformation)257 void cDerivative::ComputeScoreAndInformation(cDVector& theScore, cDMatrix& theInformation)
258 {
259 theScore = 0.0 ;
260 theInformation = 0.0 ;
261 uint myT = 0 ;
262 for (register uint n = 0 ; n < mvNSample ; n++)
263 { myT += mvT[n] ;
264 theScore += mvT[n]*mScore[n] ;
265 theInformation += mvT[n]*mInformation[n] ;
266 }
267 theScore /= myT ;
268 theInformation /= myT ;
269 }
270
ComputeCov(cHmm & theHmm,cDMatrix & theCov)271 void cDerivative::ComputeCov(cHmm& theHmm, cDMatrix& theCov)
272 {
273 uint myNParam = theHmm.GetNParam() ;
274 uint* myPlace = new uint[myNParam] ;
275 cDVector myScore(mvNFreeParam) ;
276 cDMatrix myInformation(mvNFreeParam, mvNFreeParam) ;
277 ComputeScoreAndInformation(myScore, myInformation) ;
278
279 theCov = Inv(myInformation) ;
280
281 uint myNFreeClass = mvNClass - 1 ;
282 uint myFreeIndex, myIndexCour, mySizeCour ;
283 cDVector myU(mvNFreeParam) ;
284
285 // initProb
286 myU = 0.0 ;
287 for (register uint i = 0 ; i < myNFreeClass ; i++)
288 myU[i] = -1.0 ;
289 theCov = AddOneVariable(theCov, myU) ;
290 mySizeCour = mvNFreeParam + 1 ;
291 // transMat
292 uint myBeg = 0 ;
293 for (register uint n = 0 ; n < mvNClass ; n++)
294 {
295 myU.ReAlloc(mySizeCour, 0.0) ;
296 myBeg += myNFreeClass ;
297 for (register uint i = myBeg ; i < myBeg+myNFreeClass ; i++)
298 myU[i] = -1.0 ;
299 theCov = AddOneVariable(theCov, myU) ;
300 mySizeCour++ ;
301 }
302 // Distribution
303
304 theHmm.mDistrParam->ComputeCov(theCov) ;
305
306 // Sorting
307 uint myNumParam = theHmm.GetNParam() ;
308 cDVector myParamNum(mvNFreeParam) ;
309 cDVector myNFreeClassVect(myNFreeClass) ;
310 uint myIndCour ;
311 uint myNextInd = mvNFreeParam ;
312 cDVector myResNum ;
313
314 for (register uint j = 0 ; j < mvNFreeParam ; j++)
315 myParamNum[j] = j ;
316 // initProba
317 myIndCour = 0 ;
318 GetSubVector(myParamNum, myIndCour, myNFreeClass, myNFreeClassVect) ;
319 myResNum = cat(myNFreeClassVect, (double)myNextInd) ;
320 myNextInd++ ;
321 // transMat
322 for (register uint i = 0 ; i < mvNClass ; i++)
323 { myIndCour += myNFreeClass ;
324 GetSubVector(myParamNum, myIndCour, myNFreeClass, myNFreeClassVect) ;
325 myResNum = cat(myResNum, myNFreeClassVect) ;
326 myResNum = cat(myResNum, (double)myNextInd) ;
327 myNextInd++ ;
328 }
329 // distribution
330 myIndCour += myNFreeClass ;
331 cDVector myParamNumDistr ;
332 GetSubVector(myParamNum, myIndCour, mvNFreeParam-myIndCour, myParamNumDistr) ;
333 cDVector myParamNumDistrAll = theHmm.mDistrParam->GetDistrNumParam(myParamNumDistr, myNextInd) ;
334 myResNum = cat(myResNum, myParamNumDistrAll) ;
335
336 cDMatrix myCov = theCov ;
337 for (register uint i = 0 ; i < myNParam ; i++)
338 for (register uint j = 0 ; j < myNParam ; j++)
339 theCov[i][j] = myCov[(int)myResNum[i]][(int)myResNum[j]] ;
340 }
341