noise under oscillations runs

R/gamma_beetles.R

@@ -68,20 +68,20 @@ gamma_example <- function(){
 
 
 	volume <- 100
-#	beetle_pars <- c(	b=5, ue= 0, ul = 0.001, up = 0, ua = 0.01, 
-#						ae = .13*k, al = .01*k, ap = .15*k,
-#						cle = .2, cap = .1, cae = 5, V=volume)
+	beetle_pars <- c(	b=5, ue= 0, ul = 0.001, up = 0, ua = 0.01, 
+						ae = .13*k, al = .01*k, ap = .15*k,
+						cle = .2, cap = .1, cae = 5, V=volume)
 ## 
-	beetle_pars <- c(	b=5, ue= 0, ul = 0.001, up = 0, ua = .001, 
-						ae = .1*k, al = .01*k, ap = .1*k,
-						cle = 1, cap = .4, cae = 1, V=volume)
+#	beetle_pars <- c(	b=5, ue= 0, ul = 0.001, up = 0, ua = .001, 
+#						ae = .1*k, al = .01*k, ap = .1*k,
+#						cle = 1, cap = .4, cae = 1, V=volume)
 
-	times <- seq(0,400,length=100)
+	times <- seq(0,4000,length=1000)
 	Xo <- numeric(adults)
-	Xo[1] <- 250
-	Xo[11] <- 50
-	Xo[21] <- 10
-	Xo[31] <- 50
+	Xo[1] <- 200
+#	Xo[11] <- 50
+#	Xo[21] <- 10
+#	Xo[31] <- 50
 #	Xo = rep(10, adults)
 
 	beetle_data <- linear_noise_approx	(Xo, times, beetle_pars, b_gamma, d_gamma, J_gamma, T_gamma, Omega=volume)
@@ -103,11 +103,11 @@ gamma_example <- function(){
 	data[,8] = beetle_data[,1+2*(3*k+1)]
 
 
-	png("poisson_noise.png", 800, 800)
+#	pdf("poisson_noise.pdf", 800, 800)
 	# Plot results
 	cols = c("yellow", "yellowgreen", "lightgreen", "darkgreen");
 	par(mfrow=c(2,1))
-	m<- max(data[,1:4])
+	m<- 100+max(data[,1:4])
 	plot(times, data[,1], col=cols[1], lwd=3, type='l', ylim=c(0,m), ylab="means", main="Gamma Aging" )
 	for (i in 2:4) {
 		lines(times, data[,i], col=cols[i], lwd=3)
@@ -122,17 +122,17 @@ gamma_example <- function(){
 
 
 
-	v<- max(data[,5:8])
-	plot(times, (data[,5]), col=cols[1], lwd=3, type='l', ylim=c(0,v) )
+	v<- max(2*data[,5:8])
+	plot(times, (data[,5]), col=cols[1], lwd=3, type='l', ylim=c(0,v), ylab="variance" )
 	for (i in 2:4) {
 		lines(times, (data[,4+i]), col=cols[i], lwd=3)
 	}
 
 
-	points(beetle_data[,1], sqrt(ibm$mv[[2,1]]), col="yellow")	
-	points(beetle_data[,1], sqrt(ibm$mv[[2,2]]), col="yellowgreen")	
-	points(beetle_data[,1], sqrt(ibm$mv[[2,3]]), col="lightgreen")	
-	points(beetle_data[,1], sqrt(ibm$mv[[2,4]]), col="darkgreen")
+	points(beetle_data[,1], (ibm$mv[[2,1]]), col="yellow")	
+	points(beetle_data[,1], (ibm$mv[[2,2]]), col="yellowgreen")	
+	points(beetle_data[,1], (ibm$mv[[2,3]]), col="lightgreen")	
+	points(beetle_data[,1], (ibm$mv[[2,4]]), col="darkgreen")
 
 	dev.off()
 

R/ind_based_models.R

@@ -60,6 +60,8 @@ gamma_beetles_ibm <- function(Xo = c(100,0,0,0),
 	n_states = 3*K+1;
 	max_time = 200;
 
+print(pars)
+
 	# start at beginning of each age class
 	inits = integer(n_states)
 	inits[1] = Xo[1]
@@ -88,7 +90,7 @@ gamma_beetles_ibm <- function(Xo = c(100,0,0,0),
 	}
 	moments <- sapply(1:4, calc_moments)
 
-	list(E = o[[8]], L = o[[9]], P = o[[10]], A = o[[11]], mv=moments, parameters = parameters, Xo = Xo, times = times)
+	list(E = o[[8]], L = o[[9]], P = o[[10]], A = o[[11]], mv=moments, parameters = parameters, Xo = o[[1]], times = times)
 }
 
 

src/Makevars

@@ -1,2 +1,5 @@
-PKG_CPPFLAGS=-Wall -fopenmp
-PKG_LIBS=-lm -lgsl -lgslcblas -lgomp 
+#PKG_CPPFLAGS=-Wall -fopenmp
+#PKG_LIBS=-lm -lgsl -lgslcblas -lgomp 
+
+PKG_CPPFLAGS=-Wall
+PKG_LIBS=-lm -lgsl -lgslcblas 

src/gamma_beetles.c

@@ -254,16 +254,17 @@ int ga(void)
 	const int K = 10;
 	int n_rates = 6*K+2;
 	int n_states = 3*K+1;
-	double max_time = 200;
+	double max_time = 400;
 	int samples = 50;
 	int ensembles = 2;
 
 	int * inits = (int *) calloc(n_states, sizeof(int));
 	inits[0] = 100;
-//                      {ae,  al,  ap, ue ,   ul,   up , ua , b, K,  cle,  cae,  cap, Vol
-	double parameters[13] = {1.3, 0.1, 1.5, .00, .001, .00, .003, 5, K,  .2,   0.5, .100, 100};
-
-	// outputs -- Doubles because the ensemble averaging will require it!  -- Wait, we're not ensemble averaging! 
+//                          {ae,  al,  ap,  ue,   ul,   up , ua , b, K,  cle,  cae,  cap, Vol
+//	double parameters[13] = {1.3, 0.1, 1.5, .00, .001, .00, .003, 5, K,  .2,   0.5, .100, 100};
+	double parameters[13] = {1.0, 0.1, 1.0, .00, .001, .00, .003, 5, K,  1.,   1.0, .400, 100};
+
+	// outputs -- each run concatinated  
 	int s1[50*2];
 	int s2[50*2];
 	int s3[50*2];
@@ -275,7 +276,7 @@ int ga(void)
 	int i;
 	for(i = 0; i < 100; i++)
 	{
-//		printf("%d %d %d %d\n", s1[i], s2[i], s3[i], s4[i]);
+		printf("%d %d %d %d\n", s1[i], s2[i], s3[i], s4[i]);
 	}
 
 	free(inits);