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

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

Author: ingmarvisser
Date: 2011-06-22 12:06:22 +0200 (Wed, 22 Jun 2011)
New Revision: 474

Modified:
   pkg/depmixS4/src/fb.cc
Log:
Removed all intermediate variables in computation of alpha, beta and xi.

Modified: pkg/depmixS4/src/fb.cc
===================================================================
--- pkg/depmixS4/src/fb.cc	2011-06-22 09:24:11 UTC (rev 473)
+++ pkg/depmixS4/src/fb.cc	2011-06-22 10:06:22 UTC (rev 474)
@@ -48,26 +48,9 @@
 	// compute forward variables
 	// loop over cases
 	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);
-// 		}
-// 		
-// 		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++) {
 			alpha[(bt[cas]-1)*ns[0]+i] = init[cas*ns[0]+i]*dens[(bt[cas]-1)*ns[0]+i];
 			// compute scale for alpha1
@@ -79,46 +62,13 @@
 			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
-// 		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
-// 					}
-// 					denst(i+1) = dens[t*ns[0]+i]; //B_t+1
-// 				}					
+			for(int t=bt[cas]; t<et[cas]; t++) {				
 				// compute alphat
-				sca[t] = 0.0;
+// 				sca[t] = 0.0;
 				for(int j=0; j<ns[0]; j++) {
-					alpha[t*ns[0]+j] = 0.0;
+// 					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];
 					}
@@ -131,52 +81,33 @@
 					alpha[t*ns[0]+i] *= sca[t];
 				}
 			}
-		}
+		}	
 		
-		
-		
-		
-		
-		
-		// compute backward variables and xi
-		matrix betatp1(ns[0]);
-		matrix betat(ns[0]);
-		// compute initial beta, ie for t=T (for each case)
+		// compute backward variables
+		// initial backward vars for t=T
 		for(int i=0; i<ns[0]; i++) {
-			betatp1(i+1) = sca[et[cas]-1];
+			beta[(et[cas]-1)*ns[0]+i] = sca[et[cas]-1];
 		}
-		for(int i=0; i<ns[0]; i++) {
-			beta[(et[cas]-1)*ns[0]+i] = betatp1(i+1);
-		}
- 		if(ntimes[cas]>1) {
+		
+ 		// compute beta t-1
+  		if(ntimes[cas]>1) {			
 			// loop from T-1 to 1 for each case
- 			for(int t=(et[cas]-1); t>=bt[cas]; t--) {
+ 			for(int t=(et[cas]-1); t>=bt[cas]; t--) {				
 				for(int i=0; i<ns[0]; i++) {
+// 					beta[(t-1)*ns[0]+i] = 0.0;
 					for(int j=0; j<ns[0]; j++) {
-						trans(i+1,j+1) = trdens[(i*ns[0]+j)*nt[0]+t-1]; // A_t
+						beta[(t-1)*ns[0]+i] += trdens[(i*ns[0]+j)*nt[0]+t-1]*dens[t*ns[0]+j]*beta[t*ns[0]+j];
 					}
-  					denst(i+1) = dens[t*ns[0]+i]; //B_t+1
+					beta[(t-1)*ns[0]+i] *= sca[t-1];
 				}
-				betat = trans*had(denst,betatp1)*sca[t-1];
-				// store betat somewhere
-				for(int i=0; i<ns[0]; i++) {
-					beta[(t-1)*ns[0]+i] = betat(i+1);
-				}
-				betatp1 = betat;
  			}
 			
+			// compute xi 
 			for(int t=bt[cas]; t<et[cas]; t++) {
-// 				Rprintf("t: %d\n", t);
 				for(int i=0; i<ns[0]; i++) {
-// 					Rprintf("i: %d\n", i);
 					for(int j=0; j<ns[0]; j++) {
-// 						Rprintf("j: %d\n", j);
 						xi[(i*ns[0]+j)*nt[0]+t-1] = alpha[(t-1)*ns[0]+i]*trdens[(i*ns[0]+j)*nt[0]+t-1]*dens[t*ns[0]+j]*beta[t*ns[0]+j];
-// 						Rprintf("beta_t_j: %f \n", beta[t*ns[0]+j]);	
-// 						Rprintf("t, i, j: %d %d %d\n", t, i, j);
-// 						Rprintf("xi_t_i_j: %f \n", xi[t*ns[0]*ns[0]+i*ns[0]+j]);				
 					}
-// 					Rprintf("alpha_t_i: %f \n", alpha[(t-1)*ns[0]+i]);				
 				}
 			}			
  		}