1 /**************************************************************
2 *** RHmm package
3 ***
4 *** File: RHmm.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
13 #ifdef _RDLL_
14
15
16
17 BEG_EXTERN_C
RBaumWelch(SEXP theParamHMM,SEXP theYt,SEXP theParamBW)18 DECL_DLL_EXPORT SEXP RBaumWelch (SEXP theParamHMM, SEXP theYt, SEXP theParamBW)
19 {
20 initEnum myTypeInit ;
21 distrDefinitionEnum myDistrType ;
22 uint myNbClasses, myDimObs, myNbMixt, myNbProba ;
23 cRUtil myRUtil ;
24 // Retrieve parameters for given HMM
25 myRUtil.GetValSexp(theParamHMM, eNClasses, myNbClasses) ;
26 myRUtil.GetValSexp(theParamHMM, eObsDim, myDimObs) ;
27 myRUtil.GetValSexp(theParamHMM, eNMixt, myNbMixt) ;
28 myRUtil.GetValSexp(theParamHMM, eNProba, myNbProba) ;
29
30 char myStr[255] ;
31 char *myString = (char *)myStr ;
32 myRUtil.GetValSexp(theParamHMM, eDistrType, myString) ;
33 myDistrType = eUnknownDistr ;
34 if (strcmp(myString, "NORMAL") == 0)
35 { if (myDimObs == 1)
36 myDistrType = eNormalDistr ;
37 else
38 myDistrType = eMultiNormalDistr ;
39 }
40 else
41 { if (strcmp(myString, "DISCRETE") == 0)
42 myDistrType = eDiscreteDistr ;
43 else
44 { if (strcmp(myString, "MIXTURE") == 0)
45 if (myDimObs==1)
46 myDistrType = eMixtUniNormalDistr ;
47 else
48 myDistrType = eMixtMultiNormalDistr ;
49 }
50 }
51
52 myRUtil.GetValSexp(theParamBW, eInitType, (char *)myString) ;
53 if (strcmp(myString, "RANDOM") == 0)
54 myTypeInit = eRandom ;
55 else
56 { if (strcmp(myString, "KMEANS") == 0)
57 myTypeInit = eKMeans ;
58 else
59 if (strcmp(myString, "USER") == 0)
60 myTypeInit = eUser ;
61 }
62
63 uint myNbIterMax ;
64 myRUtil.GetValSexp(theParamBW, eNMaxIter, myNbIterMax) ;
65 double myTol ;
66 myRUtil.GetValSexp(theParamBW, eTol, myTol) ;
67 uint myVerbose ;
68 myRUtil.GetValSexp(theParamBW, eVerbose, myVerbose) ;
69 uint myNbIterInit ;
70 myRUtil.GetValSexp(theParamBW, eNInitIter, myNbIterInit) ;
71 uint myNbIterMaxInit ;
72 myRUtil.GetValSexp(theParamBW, eNMaxIterinit, myNbIterMaxInit) ;
73
74 uint myNbSample = length(theYt) ;
75
76 cDVector* myY = new cDVector[myNbSample] ;
77
78 for (register uint n = 0 ; n < myNbSample ; n++)
79 {
80 SEXP myAux ;
81 myRUtil.GetValSexp(theYt, n, myAux) ;
82 myY[n].ReAlloc(length(myAux), REAL(myAux)) ;
83 }
84
85 cBaumWelchInParam myParamEntree = cBaumWelchInParam(myNbSample, myDimObs, myY, myDistrType, myNbClasses, myNbMixt, myNbProba) ;
86
87 myParamEntree.mNMaxIter = myNbIterMax ;
88 myParamEntree.mTol = myTol ;
89 myParamEntree.mVerbose = myVerbose ;
90 myParamEntree.mNInitIter = myNbIterInit ;
91 myParamEntree.mNMaxIterInit = myNbIterMaxInit ;
92 myParamEntree.mInitType = myTypeInit ;
93
94 cHmmFit myParamSortie = cHmmFit(myParamEntree) ;
95
96 if (myTypeInit == eUser)
97 { cHmm myHMM = cHmm(myParamEntree) ;
98 SEXP myAux1 ;
99 myRUtil.GetValSexp(theParamBW, eInitPoint, myAux1) ;
100 myRUtil.GetVectSexp(myAux1, 0, myHMM.mInitProba) ;
101 myRUtil.GetMatSexp(myAux1, 1, myHMM.mTransMatVector[0]) ; /* FIXME */
102 if (myHMM.mTransMatVector.size() > 1)
103 warning("Time-inhomogeneous Markov chain not supported yet for BaumWelch algorithm.");
104
105 SEXP myAux ;
106 myRUtil.GetValSexp(myAux1, 2, myAux) ; // $distribution
107 switch (myDistrType)
108 { case eNormalDistr :
109 { cUnivariateNormal* myParam = dynamic_cast<cUnivariateNormal *>(myHMM.mDistrParam) ;
110 myRUtil.GetVectSexp(myAux, 3, myParam->mMean) ;
111 myRUtil.GetVectSexp(myAux, 4, myParam->mVar) ;
112 }
113 break ;
114 case eMultiNormalDistr :
115 { cMultivariateNormal* myParam = dynamic_cast<cMultivariateNormal *>(myHMM.mDistrParam) ;
116 myRUtil.GetListVectSexp(myAux, 3, myNbClasses, myParam->mMean) ;
117 myRUtil.GetListMatSexp(myAux, 4, myNbClasses, myParam->mCov) ;
118 }
119 break ;
120 case eMixtUniNormalDistr :
121 { cMixtUnivariateNormal* myParam = dynamic_cast<cMixtUnivariateNormal *>(myHMM.mDistrParam) ;
122 myRUtil.GetListVectSexp(myAux, 4, myNbClasses, myParam->mMean) ;
123 myRUtil.GetListVectSexp(myAux, 5, myNbClasses, myParam->mVar) ;
124 myRUtil.GetListVectSexp(myAux, 6, myNbClasses, myParam->mp) ;
125 }
126 break ;
127 case eMixtMultiNormalDistr :
128 { cMixtMultivariateNormal* myParam = dynamic_cast<cMixtMultivariateNormal *>(myHMM.mDistrParam) ;
129 myRUtil.GetListListVectSexp(myAux, 4, myNbClasses, myNbMixt, myParam->mMean) ;
130 myRUtil.GetListListMatSexp(myAux, 5, myNbClasses, myNbMixt, myParam->mCov) ;
131 myRUtil.GetListVectSexp(myAux, 6, myNbClasses, myParam->mp) ;
132 }
133 break ;
134 case eDiscreteDistr :
135 { cDiscrete* myParam = dynamic_cast<cDiscrete *>(myHMM.mDistrParam) ;
136 myRUtil.GetEmissionSexp(myAux, 3, myParam->mProbaMatVector);
137
138 if (myParam->mProbaMatVector.size() > 1)
139 warning("Variable discrete emission probabilities not supported yet for BaumWelch algorithm.");
140 }
141 break ;
142 }
143 myParamSortie.CopyHmm(myHMM) ;
144 }
145 else
146 myParamSortie.BaumWelchAlgoInit(myParamEntree) ;
147 myParamSortie.BaumWelchAlgo(myParamEntree) ;
148
149 for (register uint n = 0 ; n < myNbSample ; n++)
150 myY[n].Delete() ;
151 delete[] myY ;
152
153
154 SEXP myRes,
155 myAux[6] ;
156
157 myRUtil.SetVectSexp(myParamSortie.mInitProba, myAux[0]) ;
158 myRUtil.SetMatSexp(myParamSortie.mTransMatVector[0], myAux[1]) ; // TODO: Do it for the real list
159
160 switch (myDistrType)
161 { case eNormalDistr :
162 { cUnivariateNormal* myParam = dynamic_cast<cUnivariateNormal *>(myParamSortie.mDistrParam) ;
163 myRUtil.SetVectSexp(myParam->mMean, myAux[2]) ;
164 myRUtil.SetVectSexp(myParam->mVar, myAux[3]) ;
165 }
166 break ;
167 case eMultiNormalDistr :
168 { cMultivariateNormal* myParam = dynamic_cast<cMultivariateNormal *>(myParamSortie.mDistrParam) ;
169 myRUtil.SetListVectSexp(myParam->mMean, myNbClasses, myAux[2]) ;
170 myRUtil.SetListMatSexp(myParam->mCov, myNbClasses, myAux[3]) ;
171 }
172 break ;
173 case eDiscreteDistr :
174 { cDiscrete* myParam = dynamic_cast<cDiscrete *>(myParamSortie.mDistrParam) ;
175 myRUtil.SetListVectSexp(myParam->mProbaMatVector[0], myAux[2]); // TODO: Do it for the real list
176 }
177 break ;
178 case eMixtUniNormalDistr :
179 { cMixtUnivariateNormal* myParam = dynamic_cast<cMixtUnivariateNormal *>(myParamSortie.mDistrParam) ;
180 myRUtil.SetListVectSexp(myParam->mMean, myNbClasses, myAux[2]) ;
181 myRUtil.SetListVectSexp(myParam->mVar, myNbClasses, myAux[3]) ;
182 myRUtil.SetListVectSexp(myParam->mp, myNbClasses,myAux[4]) ;
183 }
184 break ;
185 case eMixtMultiNormalDistr :
186 { cMixtMultivariateNormal* myParam = dynamic_cast<cMixtMultivariateNormal *>(myParamSortie.mDistrParam) ;
187 myRUtil.SetListListVectSexp(myParam->mMean, myNbClasses, myNbMixt, myAux[2]) ;
188 myRUtil.SetListListMatSexp(myParam->mCov, myNbClasses, myNbMixt, myAux[3]) ;
189 myRUtil.SetListVectSexp(myParam->mp, myNbClasses, myAux[4]) ;
190 }
191 break ;
192 default :
193 break ;
194 }
195
196 PROTECT(myAux[5] = allocVector(REALSXP, 5)) ;
197 REAL(myAux[5])[0] = myParamSortie.mLLH ;
198 REAL(myAux[5])[1] = myParamSortie.mBic ;
199 REAL(myAux[5])[2] = (double)( myParamSortie.mNIter) ;
200 REAL(myAux[5])[3] = myParamSortie.mTol ;
201 REAL(myAux[5])[4] = myParamSortie.mAic ;
202
203 uint myNAlloc ;
204 switch (myDistrType)
205 { case eNormalDistr :
206 case eMultiNormalDistr :
207 myNAlloc = 5 ;
208 break ;
209 case eMixtUniNormalDistr :
210 case eMixtMultiNormalDistr :
211 myNAlloc = 6 ;
212 break ;
213 case eDiscreteDistr :
214 myNAlloc = 4 ;
215 break ;
216 case eUnknownDistr :
217 default :
218 myNAlloc = 0 ;
219 break ;
220 }
221 PROTECT(myRes = allocVector(VECSXP, myNAlloc)) ;
222
223 for (register uint i = 0 ; i < 3 ; i++)
224 SET_VECTOR_ELT(myRes, i, myAux[i]) ;
225 for ( register uint i = 3 ; i < myNAlloc-1 ; i++)
226 SET_VECTOR_ELT(myRes, i, myAux[i]) ;
227 SET_VECTOR_ELT(myRes, myNAlloc-1, myAux[5]) ;
228
229 UNPROTECT(2) ;
230 myRUtil.EndProtect() ;
231 return(myRes) ;
232 }
233 END_EXTERN_C
234
235 BEG_EXTERN_C
RViterbi(SEXP theHMM,SEXP theYt)236 DECL_DLL_EXPORT SEXP RViterbi(SEXP theHMM, SEXP theYt)
237 {
238 distrDefinitionEnum myDistrType ;
239 uint myDimObs=1, myNbClasses, myNbProba=0, myNbMixt=0 ;
240 cRUtil myRUtil ;
241
242 SEXP myDistSEXP ;
243 myRUtil.GetValSexp(theHMM, fDistr, myDistSEXP) ; // Loi de proba
244 char myString[255] ;
245 char *myStr = (char *)myString ;
246 myRUtil.GetValSexp(myDistSEXP, gType, myStr) ;
247 myRUtil.GetValSexp(myDistSEXP, gNClasses, myNbClasses) ;
248 if (strcmp(myStr, "NORMAL") == 0)
249 { myRUtil.GetValSexp(myDistSEXP, 2, myDimObs) ;
250 if (myDimObs == 1)
251 myDistrType = eNormalDistr ;
252 else
253 myDistrType = eMultiNormalDistr ;
254 }
255 else
256 { if (strcmp(myStr, "DISCRETE") == 0)
257 { myDistrType = eDiscreteDistr ;
258 myRUtil.GetValSexp(myDistSEXP, 2, myNbProba) ;
259 }
260 else
261 { if (strcmp(myStr, "MIXTURE") == 0)
262 { myRUtil.GetValSexp(myDistSEXP, 3, myDimObs) ;
263 if (myDimObs == 1)
264 myDistrType = eMixtUniNormalDistr ;
265 else
266 myDistrType = eMixtMultiNormalDistr ;
267 myRUtil.GetValSexp(myDistSEXP, 2, myNbMixt) ;
268 }
269 }
270 }
271 uint myNbSample = length(theYt) ;
272 uint* myT = new uint[myNbSample] ;
273 //double **myY ;
274 cDVector* myY = new cDVector[myNbSample] ;
275 for (register uint n = 0 ; n < myNbSample ; n++)
276 {
277 SEXP myAux ;
278 myRUtil.GetValSexp(theYt, n, myAux) ;
279 myT[n] = length(myAux) / myDimObs ;
280 myY[n].ReAlloc(myT[n]*myDimObs) ;
281 myY[n]= REAL(myAux) ;
282 }
283
284 cHmm myHMM = cHmm(myDistrType, myNbClasses, myDimObs, myNbMixt, myNbProba) ;
285
286 myRUtil.GetVectSexp(theHMM, fInitProba, myHMM.mInitProba) ;
287 myRUtil.GetMatListSexp(theHMM, fTransMat, myHMM.mTransMatVector) ;
288
289 switch (myDistrType)
290 { case eNormalDistr :
291 {
292 cUnivariateNormal* myLoi = (cUnivariateNormal *)(myHMM.mDistrParam) ;
293 myRUtil.GetVectSexp(myDistSEXP, 3, myLoi->mMean) ;
294 myRUtil.GetVectSexp(myDistSEXP, 4, myLoi->mVar) ;
295 }
296 break ;
297 case eMultiNormalDistr :
298 {
299 cMultivariateNormal* myLoi = (cMultivariateNormal *)(myHMM.mDistrParam) ;
300 myRUtil.GetListVectSexp(myDistSEXP, 3, myNbClasses, myLoi->mMean) ;
301 myRUtil.GetListMatSexp(myDistSEXP, 4, myNbClasses, myLoi->mCov) ;
302 }
303 break ;
304 case eMixtUniNormalDistr :
305 {
306 cMixtUnivariateNormal* myParam = (cMixtUnivariateNormal *)(myHMM.mDistrParam) ;
307 myRUtil.GetListVectSexp(myDistSEXP, 4, myNbClasses, myParam->mMean) ;
308 myRUtil.GetListVectSexp(myDistSEXP, 5, myNbClasses, myParam->mVar) ;
309 myRUtil.GetListVectSexp(myDistSEXP, 6, myNbClasses, myParam->mp) ;
310 }
311 break ;
312 case eMixtMultiNormalDistr :
313 {
314 cMixtMultivariateNormal* myParam = (cMixtMultivariateNormal *)(myHMM.mDistrParam) ;
315 myRUtil.GetListListVectSexp(myDistSEXP, 4, myNbClasses, myNbMixt, myParam->mMean) ;
316 myRUtil.GetListListMatSexp(myDistSEXP, 5, myNbClasses, myNbMixt, myParam->mCov) ;
317 myRUtil.GetListVectSexp(myDistSEXP, 6, myNbClasses, myParam->mp) ;
318 }
319 break ;
320
321 case eDiscreteDistr :
322 {
323 cDiscrete* myParam = (cDiscrete *)(myHMM.mDistrParam) ;
324 myRUtil.GetEmissionSexp(myDistSEXP, 3, myParam->mProbaMatVector);
325 }
326 break ;
327 case eUnknownDistr :
328 default :
329 break ;
330 }
331
332 cInParam myParamEntree(myNbSample, myDimObs, myY) ;
333 myParamEntree.mDimObs = myDimObs ;
334 myParamEntree.mNMixt = myNbMixt ;
335 myParamEntree.mNProba = myNbProba ;
336 myParamEntree.mNClass = myNbClasses ;
337 myParamEntree.mDistrType = myDistrType ;
338 cViterbi myViterbi = cViterbi(myParamEntree) ;
339 myViterbi.ViterbiPath(myParamEntree, myHMM) ;
340
341 SEXP myAux[2] ;
342 myRUtil.SetListVectSexp(myViterbi.mSeq, myNbSample, myT, myAux[0]) ;
343 myRUtil.SetListValSexp(myViterbi.mLogProb, myAux[1]) ;
344
345 SEXP myRes ;
346 PROTECT(myRes = allocVector(VECSXP, 2)) ;
347 for (register uint i = 0 ; i < 2 ; i++)
348 SET_VECTOR_ELT(myRes, i, myAux[i]) ;
349 myRUtil.EndProtect() ;
350 UNPROTECT(1) ;
351 return(myRes) ;
352 }
353 END_EXTERN_C
354
355 BEG_EXTERN_C
Rforwardbackward(SEXP theHMM,SEXP theYt,SEXP theLogData)356 DECL_DLL_EXPORT SEXP Rforwardbackward(SEXP theHMM, SEXP theYt, SEXP theLogData)
357 {
358 distrDefinitionEnum myDistrType ;
359 uint myDimObs=1,
360 myNbClasses,
361 myNbProba=0,
362 myNbMixt=0 ;
363 cRUtil myRUtil ;
364
365
366 SEXP myDistSEXP ;
367
368 int myValBool = *INTEGER(theLogData) ;
369
370 bool myLogData = (myValBool != 0) ;
371
372 myRUtil.GetValSexp(theHMM, fDistr, myDistSEXP) ; // Loi de proba
373 char myString[255] ;
374 char *myStr = (char *)myString ;
375 myRUtil.GetValSexp(myDistSEXP, gType, myStr) ;
376 myRUtil.GetValSexp(myDistSEXP, gNClasses, myNbClasses) ;
377 if (strcmp(myStr, "NORMAL") == 0)
378 {myRUtil.GetValSexp(myDistSEXP, 2, myDimObs) ;
379 if (myDimObs == 1)
380 myDistrType = eNormalDistr ;
381 else
382 myDistrType = eMultiNormalDistr ;
383 }
384 else
385 {if (strcmp(myStr, "DISCRETE") == 0)
386 { myDistrType = eDiscreteDistr ;
387 myRUtil.GetValSexp(myDistSEXP, 2, myNbProba) ;
388 }
389 else
390 { if (strcmp(myStr, "MIXTURE") == 0)
391 { myRUtil.GetValSexp(myDistSEXP, 2, myNbMixt) ;
392 myRUtil.GetValSexp(myDistSEXP, 3, myDimObs) ;
393 if (myDimObs == 1)
394 myDistrType = eMixtUniNormalDistr ;
395 else
396 myDistrType = eMixtMultiNormalDistr ;
397 }
398 }
399 }
400 uint myNbSample = length(theYt) ;
401 uint* myT = new uint[myNbSample] ;
402
403 cDVector* myY = new cDVector[myNbSample] ;
404
405 for (register uint n = 0 ; n < myNbSample ; n++)
406 {
407 SEXP myAux ;
408 myRUtil.GetValSexp(theYt, n, myAux) ;
409 myT[n] = length(myAux) / myDimObs ;
410 myY[n].ReAlloc(myT[n]*myDimObs) ;
411 myY[n]= REAL(myAux) ;
412 }
413
414 cHmm myHMM = cHmm(myDistrType, myNbClasses, myDimObs, myNbMixt, myNbProba) ;
415 myRUtil.GetVectSexp(theHMM, fInitProba, myHMM.mInitProba) ;
416 myRUtil.GetMatListSexp(theHMM, fTransMat, myHMM.mTransMatVector) ;
417
418 switch (myDistrType)
419 { case eNormalDistr :
420 {
421 cUnivariateNormal* myLoi = (cUnivariateNormal *)(myHMM.mDistrParam) ;
422 myRUtil.GetVectSexp(myDistSEXP, 3, myLoi->mMean) ;
423 myRUtil.GetVectSexp(myDistSEXP, 4, myLoi->mVar) ;
424 }
425 break ;
426 case eMultiNormalDistr :
427 {
428 cMultivariateNormal* myLoi = (cMultivariateNormal *)(myHMM.mDistrParam) ;
429 myRUtil.GetListVectSexp(myDistSEXP, 3, myNbClasses, myLoi->mMean) ;
430 myRUtil.GetListMatSexp(myDistSEXP, 4, myNbClasses, myLoi->mCov) ;
431 }
432 break ;
433 case eMixtUniNormalDistr :
434 {
435 cMixtUnivariateNormal* myParam = (cMixtUnivariateNormal *)(myHMM.mDistrParam) ;
436 myRUtil.GetListVectSexp(myDistSEXP, 4, myNbClasses, myParam->mMean) ;
437 myRUtil.GetListVectSexp(myDistSEXP, 5, myNbClasses, myParam->mVar) ;
438 myRUtil.GetListVectSexp(myDistSEXP, 6, myNbClasses, myParam->mp) ;
439 }
440 break ;
441
442 case eMixtMultiNormalDistr :
443 {
444 cMixtMultivariateNormal* myParam = (cMixtMultivariateNormal *)(myHMM.mDistrParam) ;
445 myRUtil.GetListListVectSexp(myDistSEXP, 4, myNbClasses, myNbMixt, myParam->mMean) ;
446 myRUtil.GetListListMatSexp(myDistSEXP, 5, myNbClasses, myNbMixt, myParam->mCov) ;
447 myRUtil.GetListVectSexp(myDistSEXP, 6, myNbClasses, myParam->mp) ;
448 }
449 break ;
450
451 case eDiscreteDistr :
452 {
453 cDiscrete* myParam = (cDiscrete *)(myHMM.mDistrParam) ;
454 myRUtil.GetEmissionSexp(myDistSEXP, 3, myParam->mProbaMatVector);
455 }
456 break ;
457 case eUnknownDistr :
458 default :
459 break ;
460 }
461
462 cDMatrix* myProbaCond = new cDMatrix[myNbSample] ;
463
464 for (register uint n = 0 ; n < myNbSample ; n++)
465 myProbaCond[n].ReAlloc(myNbClasses, myT[n]) ;
466
467 myHMM.mDistrParam->ComputeCondProba(myY, myNbSample, myProbaCond) ;
468
469 cBaumWelch myBaumWelch=cBaumWelch(myNbSample, myT, myNbClasses) ;
470
471 myBaumWelch.OutForwardBackward(myProbaCond, myHMM, myLogData) ;
472
473 for (register uint n = 0 ; n < myNbSample ; n++)
474 { myProbaCond[n].Delete() ;
475 myY[n].Delete() ;
476 }
477
478 delete [] myY ;
479
480 delete [] myProbaCond ;
481
482 SEXP myAux[7] ;
483
484
485
486 myRUtil.SetListMatSexp(myBaumWelch.mAlpha, myNbSample,myAux[0]) ;
487 myRUtil.SetListMatSexp(myBaumWelch.mBeta, myNbSample, myAux[1]) ;
488 myRUtil.SetListMatSexp(myBaumWelch.mDelta, myNbSample, myAux[2]) ;
489 myRUtil.SetListMatSexp(myBaumWelch.mGamma, myNbSample, myAux[3]) ;
490 myRUtil.SetListListMatSexp(myBaumWelch.mXsi, myNbSample, myT, myAux[4]) ;
491 myRUtil.SetListVectSexp(myBaumWelch.mRho, myNbSample, myAux[5]) ;
492 myRUtil.SetListValSexp(myBaumWelch.mLogVrais, myAux[6]) ;
493
494 delete [] myT ;
495 SEXP myRes ;
496 PROTECT(myRes = allocVector(VECSXP, 7)) ;
497 for (register int i = 0 ; i < 7 ; i++)
498 SET_VECTOR_ELT(myRes, i, myAux[i]) ;
499 myRUtil.EndProtect() ;
500
501 UNPROTECT(1) ;
502 return(myRes) ;
503 }
504 END_EXTERN_C
505
506 BEG_EXTERN_C
RLogforwardbackward(SEXP theHMM,SEXP theYt)507 DECL_DLL_EXPORT SEXP RLogforwardbackward(SEXP theHMM, SEXP theYt)
508 {
509 distrDefinitionEnum myDistrType ;
510 uint myDimObs=1,
511 myNbClasses,
512 myNbProba=0,
513 myNbMixt=0 ;
514 cRUtil myRUtil ;
515
516
517 SEXP myDistSEXP ;
518
519 myRUtil.GetValSexp(theHMM, fDistr, myDistSEXP) ; // Loi de proba
520 char myString[255] ;
521 char *myStr = (char *)myString ;
522 myRUtil.GetValSexp(myDistSEXP, gType, myStr) ;
523 myRUtil.GetValSexp(myDistSEXP, gNClasses, myNbClasses) ;
524 if (strcmp(myStr, "NORMAL") == 0)
525 { myRUtil.GetValSexp(myDistSEXP, 2, myDimObs) ;
526 if (myDimObs == 1)
527 myDistrType = eNormalDistr ;
528 else
529 myDistrType = eMultiNormalDistr ;
530 }
531 else
532 { if (strcmp(myStr, "DISCRETE") == 0)
533 { myDistrType = eDiscreteDistr ;
534 myRUtil.GetValSexp(myDistSEXP, 2, myNbProba) ;
535 }
536 else
537 { if (strcmp(myStr, "MIXTURE") == 0)
538 { myRUtil.GetValSexp(myDistSEXP, 2, myNbMixt) ;
539 myRUtil.GetValSexp(myDistSEXP, 3, myDimObs) ;
540 if (myDimObs == 1)
541 myDistrType = eMixtUniNormalDistr ;
542 else
543 myDistrType = eMixtMultiNormalDistr ;
544 }
545 }
546 }
547 uint myNbSample = length(theYt) ;
548 uint* myT = new uint[myNbSample] ;
549
550 cDVector* myY = new cDVector[myNbSample] ;
551
552
553 for (register uint n = 0 ; n < myNbSample ; n++)
554 { SEXP myAux ;
555 myRUtil.GetValSexp(theYt, n, myAux) ;
556 myT[n] = length(myAux) / myDimObs ;
557 myY[n].ReAlloc(myT[n]*myDimObs) ;
558 myY[n]= REAL(myAux) ;
559 }
560
561 cHmm myHMM = cHmm(myDistrType, myNbClasses, myDimObs, myNbMixt, myNbProba) ;
562 myRUtil.GetVectSexp(theHMM, fInitProba, myHMM.mInitProba) ;
563 myRUtil.GetMatListSexp(theHMM, fTransMat, myHMM.mTransMatVector) ;
564
565 switch (myDistrType)
566 { case eNormalDistr :
567 { cUnivariateNormal* myLoi = (cUnivariateNormal *)(myHMM.mDistrParam) ;
568 myRUtil.GetVectSexp(myDistSEXP, 3, myLoi->mMean) ;
569 myRUtil.GetVectSexp(myDistSEXP, 4, myLoi->mVar) ;
570 }
571 break ;
572 case eMultiNormalDistr :
573 { cMultivariateNormal* myLoi = (cMultivariateNormal *)(myHMM.mDistrParam) ;
574 myRUtil.GetListVectSexp(myDistSEXP, 3, myNbClasses, myLoi->mMean) ;
575 myRUtil.GetListMatSexp(myDistSEXP, 4, myNbClasses, myLoi->mCov) ;
576 }
577 break ;
578 case eMixtUniNormalDistr :
579 { cMixtUnivariateNormal* myParam = (cMixtUnivariateNormal *)(myHMM.mDistrParam) ;
580 myRUtil.GetListVectSexp(myDistSEXP, 4, myNbClasses, myParam->mMean) ;
581 myRUtil.GetListVectSexp(myDistSEXP, 5, myNbClasses, myParam->mVar) ;
582 myRUtil.GetListVectSexp(myDistSEXP, 6, myNbClasses, myParam->mp) ;
583 }
584 break ;
585
586 case eMixtMultiNormalDistr :
587 { cMixtMultivariateNormal* myParam = (cMixtMultivariateNormal *)(myHMM.mDistrParam) ;
588 myRUtil.GetListListVectSexp(myDistSEXP, 4, myNbClasses, myNbMixt, myParam->mMean) ;
589 myRUtil.GetListListMatSexp(myDistSEXP, 5, myNbClasses, myNbMixt, myParam->mCov) ;
590 myRUtil.GetListVectSexp(myDistSEXP, 6, myNbClasses, myParam->mp) ;
591 }
592 break ;
593
594 case eDiscreteDistr :
595 { cDiscrete* myParam = (cDiscrete *)(myHMM.mDistrParam) ;
596 myRUtil.GetEmissionSexp(myDistSEXP, 3, myParam->mProbaMatVector);
597 }
598 break ;
599 case eUnknownDistr :
600 default :
601 break ;
602 }
603
604 cDMatrix* myProbaCond = new cDMatrix[myNbSample] ;
605
606 for (register uint n = 0 ; n < myNbSample ; n++)
607 myProbaCond[n].ReAlloc(myNbClasses, myT[n]) ;
608
609 myHMM.mDistrParam->ComputeCondProba(myY, myNbSample, myProbaCond) ;
610
611 cLogBaumWelch myLogBaumWelch=cLogBaumWelch(myNbSample, myT, myNbClasses) ;
612 myLogBaumWelch.LogForwardBackward(myProbaCond, myHMM) ;
613
614 for (register uint n = 0 ; n < myNbSample ; n++)
615 { myProbaCond[n].Delete() ;
616 myY[n].Delete() ;
617 }
618
619 delete [] myY ;
620
621 delete [] myProbaCond ;
622
623 SEXP myAux[6] ;
624 uint* myLigne = new uint[myNbSample] ;
625
626
627 for (register uint n = 0 ; n < myNbSample ; n++)
628 myLigne[n] = myNbClasses ;
629
630 myRUtil.SetListMatSexp(myLogBaumWelch.mLogAlpha, myNbSample,myAux[0]) ;
631 myRUtil.SetListMatSexp(myLogBaumWelch.mLogBeta, myNbSample, myAux[1]) ;
632 myRUtil.SetListMatSexp(myLogBaumWelch.mLogGamma, myNbSample, myAux[2]) ;
633 myRUtil.SetListListMatSexp(myLogBaumWelch.mLogXsi, myNbSample, myT, myAux[3]) ;
634 myRUtil.SetListVectSexp(myLogBaumWelch.mLogRho, myNbSample, myAux[4]) ;
635 myRUtil.SetListValSexp(myLogBaumWelch.mLogVrais, myAux[5]) ;
636
637 delete [] myLigne ;
638 delete [] myT ;
639 SEXP myRes ;
640 PROTECT(myRes = allocVector(VECSXP, 6)) ;
641 for (register int i = 0 ; i < 6 ; i++)
642 SET_VECTOR_ELT(myRes, i, myAux[i]) ;
643 myRUtil.EndProtect() ;
644
645 UNPROTECT(1) ;
646 return(myRes) ;
647 }
648 END_EXTERN_C
649
650 BEG_EXTERN_C
RComputeCov(SEXP theHMM,SEXP theYt)651 DECL_DLL_EXPORT SEXP RComputeCov(SEXP theHMM, SEXP theYt)
652 {
653
654 distrDefinitionEnum myDistrType ;
655 uint myDimObs=1,
656 myNClass,
657 myNProba=0,
658 myNMixt=0 ;
659 cRUtil myRUtil ;
660
661
662
663 SEXP myDistSEXP ;
664
665 myRUtil.GetValSexp(theHMM, fDistr, myDistSEXP) ; // Loi de proba
666 char myString[255] ;
667 char *myStr = (char *)myString ;
668 myRUtil.GetValSexp(myDistSEXP, gType, myStr) ;
669 myRUtil.GetValSexp(myDistSEXP, gNClasses, myNClass) ;
670 if (strcmp(myStr, "NORMAL") == 0)
671 { myRUtil.GetValSexp(myDistSEXP, 2, myDimObs) ;
672 if (myDimObs == 1)
673 myDistrType = eNormalDistr ;
674 else
675 myDistrType = eMultiNormalDistr ;
676 }
677 else
678 { if (strcmp(myStr, "DISCRETE") == 0)
679 { myDistrType = eDiscreteDistr ;
680 myRUtil.GetValSexp(myDistSEXP, 2, myNProba) ;
681 }
682 else
683 { if (strcmp(myStr, "MIXTURE") == 0)
684 { myRUtil.GetValSexp(myDistSEXP, 2, myNMixt) ;
685 myRUtil.GetValSexp(myDistSEXP, 3, myDimObs) ;
686 if (myDimObs == 1)
687 myDistrType = eMixtUniNormalDistr ;
688 else
689 myDistrType = eMixtMultiNormalDistr ;
690 }
691 }
692 }
693
694
695 uint myNSample = length(theYt) ;
696 uint* myT = new uint[myNSample] ;
697
698 cDVector* myY = new cDVector[myNSample] ;
699
700
701 for (register uint n = 0 ; n < myNSample ; n++)
702 { SEXP myAux ;
703 myRUtil.GetValSexp(theYt, n, myAux) ;
704 myT[n] = length(myAux) / myDimObs ;
705 myY[n].ReAlloc(myT[n]*myDimObs) ;
706 myY[n]= REAL(myAux) ;
707 }
708
709 cHmm myHmm = cHmm(myDistrType, myNClass, myDimObs, myNMixt, myNProba) ;
710 myRUtil.GetVectSexp(theHMM, fInitProba, myHmm.mInitProba) ;
711 myRUtil.GetMatListSexp(theHMM, fTransMat, myHmm.mTransMatVector) ;
712
713 switch (myDistrType)
714 { case eNormalDistr :
715 { cUnivariateNormal* myLoi = (cUnivariateNormal *)(myHmm.mDistrParam) ;
716 myRUtil.GetVectSexp(myDistSEXP, 3, myLoi->mMean) ;
717 myRUtil.GetVectSexp(myDistSEXP, 4, myLoi->mVar) ;
718 }
719 break ;
720 case eMultiNormalDistr :
721 { cMultivariateNormal* myLoi = (cMultivariateNormal *)(myHmm.mDistrParam) ;
722 myRUtil.GetListVectSexp(myDistSEXP, 3, myNClass, myLoi->mMean) ;
723 myRUtil.GetListMatSexp(myDistSEXP, 4, myNClass, myLoi->mCov) ;
724 }
725 break ;
726 case eMixtUniNormalDistr :
727 { cMixtUnivariateNormal* myParam = (cMixtUnivariateNormal *)(myHmm.mDistrParam) ;
728 myRUtil.GetListVectSexp(myDistSEXP, 4, myNClass, myParam->mMean) ;
729 myRUtil.GetListVectSexp(myDistSEXP, 5, myNClass, myParam->mVar) ;
730 myRUtil.GetListVectSexp(myDistSEXP, 6, myNClass, myParam->mp) ;
731 }
732 break ;
733
734 case eMixtMultiNormalDistr :
735 { cMixtMultivariateNormal* myParam = (cMixtMultivariateNormal *)(myHmm.mDistrParam) ;
736 myRUtil.GetListListVectSexp(myDistSEXP, 4, myNClass, myNMixt, myParam->mMean) ;
737 myRUtil.GetListListMatSexp(myDistSEXP, 5, myNClass, myNMixt, myParam->mCov) ;
738 myRUtil.GetListVectSexp(myDistSEXP, 6, myNClass, myParam->mp) ;
739 }
740 break ;
741
742 case eDiscreteDistr :
743 { cDiscrete* myParam = (cDiscrete *)(myHmm.mDistrParam) ;
744 myRUtil.GetEmissionSexp(myDistSEXP, 3, myParam->mProbaMatVector);
745 }
746 break ;
747 case eUnknownDistr :
748 default :
749 break ;
750 }
751
752 cInParam myInParam(myNSample, myDimObs, myY, myDistrType, myNClass, myNMixt, myNProba) ;
753 uint myNFreeParam = myHmm.GetNFreeParam() ;
754 cDerivative myDerivative(myInParam, myNFreeParam) ;
755
756 myDerivative.ComputeDerivative(myHmm, myInParam) ;
757
758 cDMatrix myCov ;
759 myDerivative.ComputeCov(myHmm, myCov) ;
760
761
762 for (register uint n = 0 ; n < myNSample ; n++)
763 { myY[n].Delete() ;
764
765 }
766
767 delete [] myY ;
768
769 delete [] myT;
770
771 SEXP myRes ;
772 myRUtil.SetMatSexp(myCov, myRes) ;
773 myRUtil.EndProtect() ;
774
775 return(myRes) ;
776 }
777 END_EXTERN_C
778
779 BEG_EXTERN_C
RScoreAndInformation(SEXP theHMM,SEXP theYt)780 DECL_DLL_EXPORT SEXP RScoreAndInformation(SEXP theHMM, SEXP theYt)
781 {
782
783 distrDefinitionEnum myDistrType ;
784 uint myDimObs=1,
785 myNClass,
786 myNProba=0,
787 myNMixt=0 ;
788 cRUtil myRUtil ;
789
790
791
792 SEXP myDistSEXP ;
793
794 myRUtil.GetValSexp(theHMM, fDistr, myDistSEXP) ; // Loi de proba
795 char myString[255] ;
796 char *myStr = (char *)myString ;
797 myRUtil.GetValSexp(myDistSEXP, gType, myStr) ;
798 myRUtil.GetValSexp(myDistSEXP, gNClasses, myNClass) ;
799 if (strcmp(myStr, "NORMAL") == 0)
800 { myRUtil.GetValSexp(myDistSEXP, 2, myDimObs) ;
801 if (myDimObs == 1)
802 myDistrType = eNormalDistr ;
803 else
804 myDistrType = eMultiNormalDistr ;
805 }
806 else
807 { if (strcmp(myStr, "DISCRETE") == 0)
808 { myDistrType = eDiscreteDistr ;
809 myRUtil.GetValSexp(myDistSEXP, 2, myNProba) ;
810 }
811 else
812 { if (strcmp(myStr, "MIXTURE") == 0)
813 { myRUtil.GetValSexp(myDistSEXP, 2, myNMixt) ;
814 myRUtil.GetValSexp(myDistSEXP, 3, myDimObs) ;
815 if (myDimObs == 1)
816 myDistrType = eMixtUniNormalDistr ;
817 else
818 myDistrType = eMixtMultiNormalDistr ;
819 }
820 }
821 }
822
823
824 uint myNSample = length(theYt) ;
825 uint* myT = new uint[myNSample] ;
826
827 cDVector* myY = new cDVector[myNSample] ;
828
829
830 for (register uint n = 0 ; n < myNSample ; n++)
831 { SEXP myAux ;
832 myRUtil.GetValSexp(theYt, n, myAux) ;
833 myT[n] = length(myAux) / myDimObs ;
834 myY[n].ReAlloc(myT[n]*myDimObs) ;
835 myY[n]= REAL(myAux) ;
836 }
837
838 cHmm myHmm = cHmm(myDistrType, myNClass, myDimObs, myNMixt, myNProba) ;
839 myRUtil.GetVectSexp(theHMM, fInitProba, myHmm.mInitProba) ;
840 myRUtil.GetMatListSexp(theHMM, fTransMat, myHmm.mTransMatVector) ;
841
842 switch (myDistrType)
843 { case eNormalDistr :
844 { cUnivariateNormal* myLoi = (cUnivariateNormal *)(myHmm.mDistrParam) ;
845 myRUtil.GetVectSexp(myDistSEXP, 3, myLoi->mMean) ;
846 myRUtil.GetVectSexp(myDistSEXP, 4, myLoi->mVar) ;
847 }
848 break ;
849 case eMultiNormalDistr :
850 { cMultivariateNormal* myLoi = (cMultivariateNormal *)(myHmm.mDistrParam) ;
851 myRUtil.GetListVectSexp(myDistSEXP, 3, myNClass, myLoi->mMean) ;
852 myRUtil.GetListMatSexp(myDistSEXP, 4, myNClass, myLoi->mCov) ;
853 }
854 break ;
855 case eMixtUniNormalDistr :
856 { cMixtUnivariateNormal* myParam = (cMixtUnivariateNormal *)(myHmm.mDistrParam) ;
857 myRUtil.GetListVectSexp(myDistSEXP, 4, myNClass, myParam->mMean) ;
858 myRUtil.GetListVectSexp(myDistSEXP, 5, myNClass, myParam->mVar) ;
859 myRUtil.GetListVectSexp(myDistSEXP, 6, myNClass, myParam->mp) ;
860 }
861 break ;
862
863 case eMixtMultiNormalDistr :
864 { cMixtMultivariateNormal* myParam = (cMixtMultivariateNormal *)(myHmm.mDistrParam) ;
865 myRUtil.GetListListVectSexp(myDistSEXP, 4, myNClass, myNMixt, myParam->mMean) ;
866 myRUtil.GetListListMatSexp(myDistSEXP, 5, myNClass, myNMixt, myParam->mCov) ;
867 myRUtil.GetListVectSexp(myDistSEXP, 6, myNClass, myParam->mp) ;
868 }
869 break ;
870
871 case eDiscreteDistr :
872 { cDiscrete* myParam = (cDiscrete *)(myHmm.mDistrParam) ;
873 myRUtil.GetEmissionSexp(myDistSEXP, 3, myParam->mProbaMatVector);
874 }
875 break ;
876 case eUnknownDistr :
877 default :
878 break ;
879 }
880
881 cInParam myInParam(myNSample, myDimObs, myY, myDistrType, myNClass, myNMixt, myNProba) ;
882 uint myNFreeParam = myHmm.GetNFreeParam() ;
883 cDerivative myDerivative(myInParam, myNFreeParam) ;
884
885 myDerivative.ComputeDerivative(myHmm, myInParam) ;
886 cDVector myScore(myNFreeParam) ;
887 cDMatrix myInformation(myNFreeParam, myNFreeParam) ;
888 myDerivative.ComputeScoreAndInformation(myScore, myInformation) ;
889
890 for (register uint n = 0 ; n < myNSample ; n++)
891 { myY[n].Delete() ;
892
893 }
894
895 delete [] myY ;
896
897 delete [] myT;
898
899
900 SEXP myAux[2] ;
901 myRUtil.SetVectSexp(myScore, myAux[0]) ;
902 myRUtil.SetMatSexp(myInformation, myAux[1]) ;
903 SEXP myRes ;
904 PROTECT(myRes = allocVector(VECSXP, 2)) ;
905 for (register int i = 0 ; i < 2 ; i++)
906 SET_VECTOR_ELT(myRes, i, myAux[i]) ;
907 myRUtil.EndProtect() ;
908
909 UNPROTECT(1) ;
910 return(myRes) ;}
911 END_EXTERN_C
912
913 #endif //_RDLL_
914