[Pomp-commits] r150 - pkg/src

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

Author: kingaa
Date: 2009-06-19 23:16:04 +0200 (Fri, 19 Jun 2009)
New Revision: 150

Modified:
   pkg/src/initstate.c
Log:
check for duplication of names between states and parameters
manage the pointer protection stack more sensibly

Modified: pkg/src/initstate.c
===================================================================
--- pkg/src/initstate.c	2009-06-17 16:49:11 UTC (rev 149)
+++ pkg/src/initstate.c	2009-06-19 21:16:04 UTC (rev 150)
@@ -10,8 +10,12 @@
 SEXP do_init_state (SEXP object, SEXP params, SEXP t0)
 {
   int nprotect = 0;
-  SEXP x, x1, par, fcall, fn, rho, covar, tcovar, covars;
-  int *dim, npar, nrep, nvar, index = 0;
+  SEXP x, x1, x2, par;
+  SEXP fcall, fn, rho;
+  SEXP covar, tcovar, covars;
+  SEXP paramnames, statenames, mindex;
+  int *dim, *midx;
+  int npar, nrep, nvar, index = 0;
   int xdim[2], j, k;
   double *p, *pp, *xp, *xpp;
 
@@ -20,7 +24,8 @@
   p = REAL(params);
 
   PROTECT(par = NEW_NUMERIC(npar)); nprotect++;
-  SET_NAMES(par,GET_ROWNAMES(GET_DIMNAMES(params)));
+  PROTECT(paramnames = GET_ROWNAMES(GET_DIMNAMES(params))); nprotect++;
+  SET_NAMES(par,paramnames);
   pp = REAL(par); 
 
   // extract the initializer function and its environment
@@ -54,27 +59,38 @@
 
   for (k = 0; k < npar; k++) pp[k] = p[k];
   PROTECT(x1 = eval(fcall,rho)); nprotect++; // do the call
-
-  if (!IS_NUMERIC(x1) || isNull(GET_NAMES(x1))) {
+  PROTECT(statenames = GET_NAMES(x1)); nprotect++;
+  
+  if (!IS_NUMERIC(x1) || isNull(statenames)) {
     UNPROTECT(nprotect);
     error("init.state error: user 'initializer' must return a named numeric vector");
   }
 
   nvar = LENGTH(x1);
   xp = REAL(x1);
+  PROTECT(mindex = match(paramnames,statenames,0)); nprotect++;
+  midx = INTEGER(mindex);
 
+  for (k = 0; k < nvar; k++) {
+    if (midx[k]!=0) {
+      UNPROTECT(nprotect);
+      error("a state variable and a parameter share a single name: '%s'",CHARACTER_DATA(STRING_ELT(statenames,k)));
+    }
+  }
+
   xdim[0] = nvar; xdim[1] = nrep;
   PROTECT(x = makearray(2,xdim)); nprotect++;
-  setrownames(x,GET_NAMES(x1),2);
+  setrownames(x,statenames,2);
   xpp = REAL(x);
 
   for (k = 0; k < nvar; k++) xpp[k] = xp[k];
 
   for (j = 1; j < nrep; j++) {
     for (k = 0; k < npar; k++) pp[k] = p[j*npar+k];
-    PROTECT(x1 = eval(fcall,rho)); nprotect++;
-    xp = REAL(x1);
+    PROTECT(x2 = eval(fcall,rho));
+    xp = REAL(x2);
     for (k = 0; k < nvar; k++) xpp[j*nvar+k] = xp[k];
+    UNPROTECT(1);
   }
 
   UNPROTECT(nprotect);