[Depmix-commits] r472 - in pkg/depmixS4: R src

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

Author: ingmarvisser
Date: 2011-06-22 10:59:11 +0200 (Wed, 22 Jun 2011)
New Revision: 472

Modified:
   pkg/depmixS4/R/fb.R
   pkg/depmixS4/src/fb.cc
Log:
Working version of C-code for alpha, beta and xi variables (not optimezed yet).

Modified: pkg/depmixS4/R/fb.R
===================================================================
--- pkg/depmixS4/R/fb.R	2011-06-21 20:50:33 UTC (rev 471)
+++ pkg/depmixS4/R/fb.R	2011-06-22 08:59:11 UTC (rev 472)
@@ -55,9 +55,10 @@
 		
 		res$alpha <- matrix(res$alpha,nc=ns,byrow=TRUE)
 		res$beta <- matrix(res$beta,nc=ns,byrow=TRUE)
+		res$xi <- array(res$xi,dim=c(nt,ns,ns))
 		
-		res$gamma <- alpha*beta/sca
-		res$logLike <- -sum(log(sca))
+		res$gamma <- res$alpha*res$beta/res$sca
+		res$logLike <- -sum(log(sca))		
 		
 	} else {
 		
@@ -81,19 +82,11 @@
 			}
 			
 			beta[et[case],] <- 1*sca[et[case]] # initialize
-			
-# 			print(et[case])
-# 			print(round(beta[et[case],],6))
-			
+						
 			if(ntimes[case]>1) {
 				for(i in (et[case]-1):bt[case]) {
-# 					cat("t: ",i,"\n")
 					if(stationary) beta[i,] <-(B[i+1,]*beta[i+1,])%*%A[1,,]*sca[i]
 					else beta[i,] <-(B[i+1,]*beta[i+1,])%*%A[i,,]*sca[i]
-# 					print(round(sca[i],6))
-# 					print(round(beta[i,],6))
-#  					print(round(B[i+1,],6))
-#  					print(round(A[i,,],6))
 				}
 				
 				for(i in bt[case]:(et[case]-1)) {
Modified: pkg/depmixS4/src/fb.cc
===================================================================
--- pkg/depmixS4/src/fb.cc	2011-06-21 20:50:33 UTC (rev 471)
+++ pkg/depmixS4/src/fb.cc	2011-06-22 08:59:11 UTC (rev 472)
@@ -58,12 +58,27 @@
 			// 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
+			sca1 += alpha[(bt[cas]-1)*ns[0]+i];
+		}
+		sca[(bt[cas]-1)] = 1/sca1;
+		// scale alpha1 and copy to alpha
+		for(int i=0; i<ns[0]; i++) {
+			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
@@ -89,63 +104,52 @@
 			}
 		}
 		
+		// alphat without intervening matrices
+		
+		
 		// compute backward variables and xi
 		matrix betatp1(ns[0]);
 		matrix betat(ns[0]);
-// 		Rprintf("et[cas]: %d \n",et[cas]);		
-// 		Rprintf("sca et[cas]: %f \n",sca[et[cas]-1]);
 		// compute initial beta, ie for t=T (for each case)
 		for(int i=0; i<ns[0]; i++) {
 			betatp1(i+1) = sca[et[cas]-1];
 		}
-// 		betatp1.print();
 		for(int i=0; i<ns[0]; i++) {
 			beta[(et[cas]-1)*ns[0]+i] = betatp1(i+1);
 		}
  		if(ntimes[cas]>1) {
 			// loop from T-1 to 1 for each case
  			for(int t=(et[cas]-1); t>=bt[cas]; t--) {
-// 				Rprintf("t: %d \n", 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
-//  				if(stationary) beta[i,] <-(B[i+1,]*beta[i+1,])%*%A[1,,]*sca[i]
 				}
-//  				transpose(trans).print();
-//  				denst.print();
-// 				Rprintf("sca t: %f \n",sca[t-1]);
 				betat = trans*had(denst,betatp1)*sca[t-1];
-// 				betat.print();
 				// store betat somewhere
 				for(int i=0; i<ns[0]; i++) {
 					beta[(t-1)*ns[0]+i] = betat(i+1);
 				}
 				betatp1 = betat;
  			}
- 			
- 			
-//  			for(i in bt[case]:(et[case]-1)) {
-//  				if(stationary) xi[i,,] <- rep(alpha[i,],each=ns)*(B[i+1,]*beta[i+1,]*A[1,,])
-//  				else xi[i,,] <- rep(alpha[i,],each=ns)*(B[i+1,]*beta[i+1,]*A[i,,])
-// 			}
+			
+			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]);				
+				}
+			}			
  		}
-		
-		
-// 		R-code for the backwards/xi loop
-// 		if(ntimes[case]>1) {
-// 			for(i in (et[case]-1):bt[case]) {
-// 				if(stationary) beta[i,] <-(B[i+1,]*beta[i+1,])%*%A[1,,]*sca[i]
-// 				else beta[i,] <-(B[i+1,]*beta[i+1,])%*%A[i,,]*sca[i]
-// 			}
-// 			
-// 			for(i in bt[case]:(et[case]-1)) {
-// 				if(stationary) xi[i,,] <- rep(alpha[i,],each=ns)*(B[i+1,]*beta[i+1,]*A[1,,])
-// 				else xi[i,,] <- rep(alpha[i,],each=ns)*(B[i+1,]*beta[i+1,]*A[i,,])
-// 			}
-// 		}
-				
+						
 	} // end cases
 	
 }