[Depmix-commits] r473 - pkg/depmixS4/src

[noreply at r-forge.r-project.org]

Author: ingmarvisser
Date: 2011-06-22 11:24:11 +0200 (Wed, 22 Jun 2011)
New Revision: 473

Modified:
   pkg/depmixS4/src/fb.cc
Log:
Optimized computation of alpha (removed intermediate computations); working correctly.

Modified: pkg/depmixS4/src/fb.cc
===================================================================
--- pkg/depmixS4/src/fb.cc	2011-06-22 08:59:11 UTC (rev 472)
+++ pkg/depmixS4/src/fb.cc	2011-06-22 09:24:11 UTC (rev 473)
@@ -50,22 +50,22 @@
 	for(int cas=0; cas<nc[0]; cas++) {
 		// compute alpha1 for this case
 		double sca1=0.0;
-		matrix alpha1(ns[0],1);
-		matrix denst(ns[0]);
-		for(int i=0; i<ns[0]; i++) {
-			alpha1(i+1) = init[cas*ns[0]+i]*dens[(bt[cas]-1)*ns[0]+i];
-			denst(i+1) = dens[(bt[cas]-1)*ns[0]+i];
-			// compute scale for alpha1
-			sca1 += alpha1(i+1);
-		}
+// 		matrix alpha1(ns[0],1);
+ 		matrix denst(ns[0]);
+// 		for(int i=0; i<ns[0]; i++) {
+// 			alpha1(i+1) = init[cas*ns[0]+i]*dens[(bt[cas]-1)*ns[0]+i];
+// 			denst(i+1) = dens[(bt[cas]-1)*ns[0]+i];
+// 			// compute scale for alpha1
+// 			sca1 += alpha1(i+1);
+// 		}
+// 		
+// 		sca[(bt[cas]-1)] = 1/sca1;
+// 		// scale alpha1 and copy to alpha
+// 		for(int i=0; i<ns[0]; i++) {
+// 			alpha1(i+1) *= sca[(bt[cas]-1)];
+// 			alpha[(bt[cas]-1)*ns[0]+i] = alpha1(i+1);
+// 		}
 		
-		sca[(bt[cas]-1)] = 1/sca1;
-		// scale alpha1 and copy to alpha
-		for(int i=0; i<ns[0]; i++) {
-			alpha1(i+1) *= sca[(bt[cas]-1)];
-			alpha[(bt[cas]-1)*ns[0]+i] = alpha1(i+1);
-		}
-		
 		// compute alpha1 without intervening matrices
 		sca1 = 0.0;
 		for(int i=0; i<ns[0]; i++) {
@@ -79,34 +79,65 @@
 			alpha[(bt[cas]-1)*ns[0]+i] *= sca[(bt[cas]-1)];
 		}
 		
-		matrix alphat(ns[0],1);
-		matrix trans(ns[0],ns[0]);
-		// loop over ntimes[cas]>1
+ 		matrix alphat(ns[0],1);
+ 		matrix trans(ns[0],ns[0]);
+// 		// loop over ntimes[cas]>1
+// 		if(ntimes[cas]>0) {
+// 			for(int t=bt[cas]; t<et[cas]; t++) {
+// 				// get trans and dens values
+// 				for(int i=0; i<ns[0]; i++) {
+// 					for(int j=0; j<ns[0]; j++) {
+// 						trans(i+1,j+1) = trdens[(i*ns[0]+j)*nt[0]+t-1]; // A_t
+// 					}
+// 					denst(i+1) = dens[t*ns[0]+i]; //B_t+1
+// 				}					
+// 				// compute alphat
+// 				alphat = had(transpose(trans)*alpha1,denst);// A_t*alpha_t*B_t+1
+// 				// compute scale for t
+// 				sca[t] = 1/(alphat.msum());
+// 				// scale alphat values
+// 				for(int i=0; i<ns[0]; i++) {
+// 					alphat(i+1) *= sca[t];
+// 					alpha[t*ns[0]+i] = alphat(i+1);
+// 				}
+// 				alpha1 = alphat;
+// 			}
+// 		}
+		
+		// compute alphat without intervening matrices
+		
 		if(ntimes[cas]>0) {
 			for(int t=bt[cas]; t<et[cas]; t++) {
 				// get trans and dens values
-				for(int i=0; i<ns[0]; i++) {
-					for(int j=0; j<ns[0]; j++) {
-						trans(i+1,j+1) = trdens[(i*ns[0]+j)*nt[0]+t-1]; // A_t
+// 				for(int i=0; i<ns[0]; i++) {
+// 					for(int j=0; j<ns[0]; j++) {
+// 						trans(i+1,j+1) = trdens[(i*ns[0]+j)*nt[0]+t-1]; // A_t
+// 					}
+// 					denst(i+1) = dens[t*ns[0]+i]; //B_t+1
+// 				}					
+				// compute alphat
+				sca[t] = 0.0;
+				for(int j=0; j<ns[0]; j++) {
+					alpha[t*ns[0]+j] = 0.0;
+					for(int i=0; i<ns[0]; i++) {
+						alpha[t*ns[0]+j] += trdens[(i*ns[0]+j)*nt[0]+t-1]*alpha[(t-1)*ns[0]+i];
 					}
-					denst(i+1) = dens[t*ns[0]+i]; //B_t+1
-				}					
-				// compute alphat
-				alphat = had(transpose(trans)*alpha1,denst);// A_t*alpha_t*B_t+1
-				// compute scale for t
-				sca[t] = 1/(alphat.msum());
+					alpha[t*ns[0]+j] *= dens[t*ns[0]+j];
+					sca[t] += alpha[t*ns[0]+j];
+				}
+				sca[t] = 1/(sca[t]);				
 				// scale alphat values
 				for(int i=0; i<ns[0]; i++) {
-					alphat(i+1) *= sca[t];
-					alpha[t*ns[0]+i] = alphat(i+1);
+					alpha[t*ns[0]+i] *= sca[t];
 				}
-				alpha1 = alphat;
 			}
 		}
 		
-		// alphat without intervening matrices
 		
 		
+		
+		
+		
 		// compute backward variables and xi
 		matrix betatp1(ns[0]);
 		matrix betat(ns[0]);