now have parallel (openmp) support compiled into R via C code

R/crowley.R

@@ -31,17 +31,28 @@ T_crowley <- function(t,y,p){
 
 
 crowley_example <- function(){
+
 	pop <- 10000
 	crowley_parameters <- c(b1=.11/pop, b2=.6/pop, d1=0.1, d2=.1, c1=0.1/pop, c2=4/pop, K=pop)
 	#crowley_parameters <- c(b1=.2/pop, b2=.6/pop, d1=0.1, d2=.1, c1=0.1/pop, c2=.2/pop, K=pop)
-	times <- seq(0,2000,length=1000)
-	Xo <- c(595, 4550)
+	times <- seq(0,1000,length=500)
+	Xo <- c(505, 4550)
+	ibm <- crowley_ibm(Xo = Xo, par=crowley_parameters, time=times,reps=1)
+
+
+	crowley_sim <- linear_noise_approx(
+				Xo, times, crowley_parameters, 
+				b_crowley, d_crowley, J_crowley,
+				T_crowley, Omega=pop)
 
-	crowley_sim <- linear_noise_approx(Xo, times, crowley_parameters, 
-										b_crowley, d_crowley, J_crowley,
-										T_crowley, Omega=pop)
+	ens <- sfLapply(1:10, 
+			function(i){
+				ibm <- crowley_ibm(Xo = Xo, 
+					par=crowley_parameters, 
+					time=times)
+				list(ibm$x1, ibm$x2)
+		})
 
-	ibm <- crowley_ibm(Xo = Xo, par=crowley_parameters, time=times)
 
 	#png("crowley_noise.png")
 	par(mfrow=c(2,1))

R/ind_based_models.R

@@ -3,16 +3,19 @@
 # @section DESCRIPTION wrapper to the C code containing the gillespie simulation for individual-based models.  
 
 # Pars = {x, y, bx, by, dx, dy, cx, cy, K} 
-crowley_ibm <- function(Xo = c(500,4500), parameters=c(0.11, .6, .1, .1, .1, 4, 10000), times = seq(0,1000,length=500)){
-	N <- length(times)
+crowley_ibm <- function(Xo = c(500,4500), parameters=c(0.11/10000, .6/10000, .1, .1, .1/10000, 4/10000, 10000), times = seq(0,1000,length=500), reps=1){
+	samples <- length(times)
+	N <- reps*(1+samples)
 	maxtime <- max(times)
 	pars <- c(Xo, parameters)
-	o <- .C("crowley", double(N), double(N), as.double(pars), as.integer(N), as.double(maxtime) )
+	o <- .C("crowley", double(N), double(N), as.double(pars), as.integer(samples), as.integer(reps), as.double(maxtime) )
 	list(x1 = o[[1]], x2 = o[[2]], parameters = parameters, Xo = Xo)
 }
 
-# Pars = {E, L, P, A, b, ue, ul, up, ua, ae, al, ap, cle, cap, cae} 
 
+
+
+# Pars = {E, L, P, A, b, ue, ul, up, ua, ae, al, ap, cle, cap, cae} 
 beetles_ibm <- function(Xo = c(100,0,0,0), 
 						parameters= c(5., 0, 0.001, 0, 0.003, 1/3.8, 1/(20.2-3.8), 1/(25.5-20.2), 0.01, 0.004, 0.01),
 						times = seq(0,1000,length=500)

src/beetles.c

@@ -168,7 +168,7 @@ void * beetles_reset(const void * inits)
  * Currently this function isn't very api like, as it is 
  * closely tied to the logical structure of the record data structure 
  * which itself could be abstracted more */
-void beetles_fixed_interval(const double t, const void * mypars, void * myrecord)
+void beetles_fixed_interval(const double t, void * mypars, void * myrecord, int rep)
 {
 	record * my_record = (record *) myrecord;
 	if (t > my_record->i * my_record->t_step) 

src/crowley.c

@@ -20,20 +20,19 @@
 
 /** record class for Crowley model */
 typedef struct {
-	int i;
 	double t_step;
+	size_t N;
 	double * s1;
 	double * s2;
 } record;
 
-record* c_record_alloc(size_t N, double maxtime)
+record* c_record_alloc(size_t N, size_t replicates, double maxtime)
 {
 	record* myrecord = (record*) malloc(sizeof(record));
 	myrecord->t_step = maxtime/N;
-	printf("%g\n", maxtime/N);
-	myrecord->i = 0;
-	myrecord->s1 = (double*) calloc(N+1, sizeof(double));
-	myrecord->s2 = (double*) calloc(N+1, sizeof(double));
+	myrecord->N = N;
+	myrecord->s1 = (double*) calloc(replicates*(N+1), sizeof(double));
+	myrecord->s2 = (double*) calloc(replicates*(N+1), sizeof(double));
 	return myrecord;
 }
 
@@ -45,6 +44,20 @@ void c_record_free(record * myrecord)
 }
 
 
+/** Execute all tasks that occur on fixed interval schedule  */
+void crowley_fixed_interval(const double t, void * mypars, void * myrecord, int rep)
+{
+	record * my_record = (record *) myrecord;
+	double * s = (double *) mypars;
+	if (t > s[9]*my_record->t_step) 
+	{
+		my_record->s1[(int) s[9]+rep*my_record->N] = s[0]; 
+		my_record->s2[(int) s[9]+rep*my_record->N] = s[1]; 
+		//printf("%g %g %g\n", t, s[0], s[1]);
+		s[9] += 1; //increment sample counter
+	}
+}
+
 
 /*		   0  1  2   3   4   5   6   7   8
  * Pars = {x, y, bx, by, dx, dy, cx, cy, K} */
@@ -125,31 +138,16 @@ double d2_out(void * ss)
 void * crowley_reset(const void * inits)
 {
 	const double * ss = (const double *) inits;
-	double * s = (double *) calloc(9, sizeof(double));
+	double * s = (double *) calloc(10, sizeof(double));
 	int i;
 	for(i=0;i<9;i++) s[i] = ss[i];
+	s[9] = 0; // sampling counter
 	return s;
 }
 
-/** Execute all tasks that occur on fixed interval schedule 
- * Currently this function isn't very api like, as it is 
- * closely tied to the logical structure of the record data structure 
- * which itself could be abstracted more */
-void crowley_fixed_interval(const double t, const void * mypars, void * myrecord)
+void crowley(double* s1, double* s2, double* inits, int* n_samples, int* replicates, double * maxtime)
 {
-	record * my_record = (record *) myrecord;
-	if (t > my_record->i * my_record->t_step) 
-	{
-		double * s = (double *) mypars;
-		my_record->s1[my_record->i] = s[0]; 
-		my_record->s2[my_record->i] = s[1]; 
-//		printf("%g %g %g\n", t, s[0], s[1]);
-		++my_record->i;
-	}
-}
 
-void crowley(double* s1, double* s2, double* inits, int* n_samples, double * maxtime)
-{
 	/** Create a list of all event functions and their associated outcomes 
 	 *  Order doesn't matter, but make sure outcomes are paired with the 
 	 *  appropriate function! */
@@ -168,36 +166,43 @@ void crowley(double* s1, double* s2, double* inits, int* n_samples, double * max
 	FIXED fixed_interval_fn = &crowley_fixed_interval;
 
 	const size_t n_event_types = 4;
-	const size_t ensembles = 1;
-//	const size_t maxtime = 5000;
 
-	record * my_record = c_record_alloc(*n_samples, *maxtime);
-	gillespie(rate_fn, outcome, n_event_types, inits, my_record, *maxtime, ensembles, reset_fn, fixed_interval_fn);
+	record * my_record = c_record_alloc(*n_samples, *replicates, *maxtime);
+	gillespie(rate_fn, outcome, n_event_types, inits, my_record, *maxtime, (size_t) *replicates, reset_fn, fixed_interval_fn);
 
 	int i;
-	for(i = 0; i< *n_samples; i++)
+	for(i = 0; i< (*n_samples)*(*replicates); i++)
 	{ 
 		s1[i] = my_record->s1[i]; 
 		s2[i] = my_record->s2[i];
 	}
+
 	c_record_free(my_record);
-//	euler(my_pars, MAX_TIME, stderr);
-//	gslode(my_pars, MAX_TIME, theory);
 }
 
 
 
 int crow(void)
 {
-	int n_samples = 100;
-	double s1[100];
-	double s2[100];
+	const int n_samples = 10;
+	const int replicates = 2;
+	double s1[n_samples*replicates];
+	double s2[n_samples*replicates];
+
+	int n = n_samples;
+	int reps = replicates;
 	//			0  1  2   3   4   5   6   7   8
 	//	Pars = {x, y, bx, by, dx, dy, cx, cy, K} 
 	double c_inits[9] = {500, 4500, 0.11/10000, .6/10000, .1, .1, .1/10000, 4/10000, 10000};
-	double maxtime = 5000; 
+	double maxtime = 1; 
+
+	crowley(s1, s2, c_inits, &n, &reps, &maxtime);
+
+	int i;
+	for(i = 0; i< n_samples * replicates; i++)
+		printf("%g, %g\n", s1[i], s2[i] );
+
 
-	crowley(s1, s2, c_inits, &n_samples, &maxtime);
 	printf("m1 = %g sd1 = %g\nm2 = %g sd2 = %g\n",
 			gsl_stats_mean(s1, 1, n_samples), 
 			sqrt(gsl_stats_variance(s1, 1, n_samples)),

src/gillespie.c

@@ -92,12 +92,10 @@ gillespie(	const event_fn * rate_fn,
 	 * the parallel region.  Currently this uses the pars_alloc function 
 	 * to allocate space in user-defined way without assuming the function type */
 	#pragma omp parallel shared(rng, rate_fn, outcome, my_record, inits, reset_fn, fixed_interval_fn) \
-	private(lambda, t, tmp, i, check, rates_data)
+	private(lambda, t, tmp, i, check, rates_data,l)
 	{
 		/* The vector to store cumulative sum of rates */
 		rates_data = (double *) calloc (n_event_types,sizeof(double) );
-		/* allocates parameters based on user defined function */
-//		my_pars = pars_cpy(mypars);
 
 		/* Loop over ensembles, will be parallelized if compiled with -fopenmp */
 		#pragma omp for
@@ -115,7 +113,7 @@ gillespie(	const event_fn * rate_fn,
 				t += gsl_ran_exponential(rng, 1/lambda);
 
 				/* Events such as sampling occur at regular intervals */
-				fixed_interval_fn(t, my_pars, my_record);
+				fixed_interval_fn(t, my_pars, my_record, l);
 
 				/* Determine if event is a birth or death */
 				tmp = gsl_rng_uniform(rng);
@@ -128,9 +126,10 @@ gillespie(	const event_fn * rate_fn,
 			if(check) break;
 
 		   }/* end evolution */
+			free(my_pars);
 		}/* end ensembles */
 		free(rates_data);
-	}
+	} /* end parallel */
 	gsl_rng_free(rng);
 }
 

src/gillespie.h

@@ -11,11 +11,11 @@
 
 typedef double (* event_fn)(void * my_pars);
 typedef void * (* RESET)(const void * inits);
-typedef void (* FIXED)(const double t, const void * my_pars, void * my_record);
+typedef void (* FIXED)(const double t, void * my_pars, void * my_record, int rep);
 
 void * reset(void * inits);
 
-void fixed_interval_tasks(const double t, const void * my_pars, void * my_record);
+//void fixed_interval_tasks(const double t, const void * my_pars, void * my_record);
 
 
 

src/main.c

@@ -1,40 +1,19 @@
-#include "warning_signals.h"
-#define SAMPLE_TIME 1
-#define ENSEMBLES 1
-#define SAMPLE_FREQ 1.
-#define MAX_TIME 500
-#define START_POLLUTING 450
-#define S 500 
 void correlation(void);
-
 int tribol(void);
 int crow(void);
 int beetle(void);
+int ws(void);
 
 int main(void){
 
-/* Test Warning Signals */
-
-	double  time[S], a[S], means[S], vars[S], skews[S], ar1[S], arN[S];
-	double sample_freq = SAMPLE_FREQ, start_polluting = START_POLLUTING, pollute_rate = 1, pollute_increment = 1;
-	int sample_time = SAMPLE_TIME, max_time = MAX_TIME, n_ensembles = ENSEMBLES;
-	warning_signals(time, a, means, vars, skews, ar1, arN, &sample_time, 
-					&sample_freq, &max_time, &n_ensembles, &start_polluting,
-					&pollute_rate, &pollute_increment);
-	int i;
-	for(i=0; i<S; i++)
-		printf("%g %g %g\n", time[i], a[i], means[i] );
-
-
+	/* Test warning signals */
+//	ws();
 	/* Test correlation (for warning_signals) */
-	correlation();
-
+//	correlation();
 	/* Test beetles "main"*/
-	beetle();
-
+//	beetle();
 	/* Test crowley "main" */
 	crow();
-
 	/* Test tribolium "main" */
 //	tribol();
 	return 0;

src/makefile

@@ -4,34 +4,31 @@
 CCP = g++
 CC = gcc
 C_OMP = -fopenmp 
-CFLAGS = -Wall -lgsl -lgslcblas -lm -O3 -fpic
-BEETLES = tribolium.cpp kde.cpp 
+CFLAGS = -Wall -lgsl -lgslcblas -lm -O3 -fpic -g
+TRIBOL = tribolium.cpp kde.cpp 
 WARNING = warning_signals.c odeintegrators.c gillespie.c correlation.c saddle_node.c record.c pars.c main.c 
-CROWLEY = crowley.c gillespie.c 
-STAGES = beetles.c gillespie.c
+CROWLEY = crowley.c gillespie.c main.c
+BEETLES = beetles.c gillespie.c main.c
+
+build:
+
+warning: 
+	$(CC) $(CFLAGS) $(C_OMP) -o warning.exe $(WARNING)
+
+crowley:
+	$(CC) $(CFLAGS) $(C_OMP) -o crowley.exe $(CROWLEY)
+
+beetles:
+	$(CC) $(CFLAGS) $(C_OMP) -o beetles.exe $(BEETLES)
+
+tribol:
+	$(CCP) $(CFLAGS) $(C_OMP) -o tribol.exe $(TRIBOL)
+
+clean:
+	rm -f *.o *.so
+
 
-build: 
 
-#build: 
-#	$(CC) $(CFLAGS) -o warning.exe $(WARNING)
-#	$(CC) $(CFLAGS) -o crowley.exe $(CROWLEY)
-#	$(CCP) $(CFLAGS) -o beetles.exe $(BEETLES)
-#
-#stages:
-#	$(CC) $(CFLAGS) -o stages.exe $(STAGES)
-#
-#warning: $(WARNING)
-#	$(CC) $(CFLAGS) -o warning.exe $(WARNING)
-#
-#crowley: $(CROWLEY)
-#	$(CC) $(CFLAGS) -o crowley.exe $(CROWLEY)
-#
-#beetles: $(BEETLES)
-#	$(CCP) $(CFLAGS) -o beetles.exe $(BEETLES)
-#
-#clean:
-#	rm -f *.o *.so
-#
 ### run make profile, then execute program, then run gprof -a $(EXE) gmon.out
 #profile:  
 #	$(CC) $(CFLAGS) -g -pg -o $(EXE) $(FILES)

src/warning_signals.c

@@ -44,7 +44,7 @@ void reorder(double * data, const int i, const size_t windowsize){
  * Currently this function isn't very api like, as it is 
  * closely tied to the logical structure of the record data structure 
  * which itself could be abstracted more */
-void warning_fixed_interval(const double t, const void * mypars, void * myrecord)
+void warning_fixed_interval(const double t, void * mypars, void * myrecord, int rep)
 {
 	pars * my_pars = (pars *) mypars;
 	record * my_record = (record *) myrecord;
@@ -192,4 +192,25 @@ void warning_signals(
 	record_free(my_record);
 }
 
+/* ws "main" function */
+int ws(void)
+{
+	const size_t S = 500;
+	double  time[S], a[S], means[S], vars[S], skews[S], ar1[S], arN[S];
+	double sample_freq = 1;
+	double start_polluting = 450;
+	double pollute_rate = 1; 
+	double pollute_increment = 1;
+	int sample_time = 1; 
+	int max_time = 500;
+	int n_ensembles = 1;
+	warning_signals(time, a, means, vars, skews, ar1, arN, &sample_time, 
+					&sample_freq, &max_time, &n_ensembles, &start_polluting,
+					&pollute_rate, &pollute_increment);
+	int i;
+	for(i=0; i<S; i++)
+		printf("%g %g %g\n", time[i], a[i], means[i] );
+
+	return 0;
+}