[Pomp-commits] r817 - in pkg/pomp: . R src

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

Author: kingaa
Date: 2013-01-15 16:26:18 +0100 (Tue, 15 Jan 2013)
New Revision: 817

Added:
   pkg/pomp/src/mif.c
Modified:
   pkg/pomp/NAMESPACE
   pkg/pomp/R/mif.R
   pkg/pomp/src/pomp_internal.h
Log:
- put mif update into C


Modified: pkg/pomp/NAMESPACE
===================================================================
--- pkg/pomp/NAMESPACE	2013-01-15 10:28:51 UTC (rev 816)
+++ pkg/pomp/NAMESPACE	2013-01-15 15:26:18 UTC (rev 817)
@@ -10,6 +10,7 @@
           lookup_in_table,
           SSA_simulator,
           R_Euler_Multinom,D_Euler_Multinom,R_GammaWN,
+          mif_update,
           pfilter_computations,
           simulation_computations,
           iterate_map,traj_transp_and_copy,

Modified: pkg/pomp/R/mif.R
===================================================================
--- pkg/pomp/R/mif.R	2013-01-15 10:28:51 UTC (rev 816)
+++ pkg/pomp/R/mif.R	2013-01-15 15:26:18 UTC (rev 817)
@@ -298,10 +298,7 @@
     switch(
            method,
            mif={              # original Ionides et al. (2006) average
-             v <- pfp at pred.var[pars,,drop=FALSE]
-             v1 <- cool.sched$gamma * (1 + var.factor^2) * sigma[pars]^2
-             theta.hat <- cbind(theta[pars],pfp at filter.mean[pars,,drop=FALSE])
-             theta[pars] <- theta[pars] + colSums(apply(theta.hat,1,diff)/t(v))*v1
+             theta <- .Call(mif_update,pfp,theta,cool.sched$gamma,var.factor,sigma,pars)
            },
            unweighted={                 # unweighted average
              theta[pars] <- rowMeans(pfp at filter.mean[pars,,drop=FALSE])

Added: pkg/pomp/src/mif.c
===================================================================
--- pkg/pomp/src/mif.c	                        (rev 0)
+++ pkg/pomp/src/mif.c	2013-01-15 15:26:18 UTC (rev 817)
@@ -0,0 +1,49 @@
+// -*- C++ -*-
+
+#include "pomp_internal.h"
+#include <Rdefines.h>
+
+// implements the Ionides et al. (2006) MIF update rule
+
+SEXP mif_update (SEXP pfp, SEXP theta, SEXP gamma, SEXP varfactor, 
+		 SEXP sigma, SEXP pars)
+{
+  int nprotect = 0;
+  double *v, *m1, *m2;
+  double scal, sig, grad;
+  int npar, ntimes, nfm, npv;
+  SEXP FM, PV, newtheta;
+  int *sidx, *thidx, *midx, *vidx, *dim;
+  int i, j;
+
+  sig = *(REAL(varfactor));
+  scal = *(REAL(gamma))*(1+sig*sig);
+
+  PROTECT(FM = GET_SLOT(pfp,install("filter.mean"))); nprotect++;
+  PROTECT(PV = GET_SLOT(pfp,install("pred.var"))); nprotect++;
+
+  npar = LENGTH(pars);
+  dim = INTEGER(GET_DIM(FM)); nfm = dim[0]; ntimes = dim[1];
+  dim = INTEGER(GET_DIM(PV)); npv = dim[0];
+
+  sidx = INTEGER(PROTECT(MATCHNAMES(sigma,pars))); nprotect++;
+  thidx = INTEGER(PROTECT(MATCHNAMES(theta,pars))); nprotect++;
+  midx = INTEGER(PROTECT(MATCHROWNAMES(FM,pars))); nprotect++;
+  vidx = INTEGER(PROTECT(MATCHROWNAMES(PV,pars))); nprotect++;
+
+  PROTECT(newtheta = duplicate(theta)); nprotect++;
+
+  for (i = 0; i < npar; i++) {
+    sig = REAL(sigma)[sidx[i]];
+    m1 = REAL(theta)+thidx[i];
+    m2 = REAL(FM)+midx[i];
+    v = REAL(PV)+vidx[i];
+    grad = (*m2-*m1)/(*v);
+    for (j = 1, m1 = m2, m2 += nfm, v += npv; j < ntimes; j++, m1 = m2, m2 += nfm, v += npv) 
+      grad += (*m2-*m1)/(*v);
+    REAL(newtheta)[thidx[i]] += scal*sig*sig*grad;
+  }
+
+  UNPROTECT(nprotect);
+  return newtheta;
+}

Modified: pkg/pomp/src/pomp_internal.h
===================================================================
--- pkg/pomp/src/pomp_internal.h	2013-01-15 10:28:51 UTC (rev 816)
+++ pkg/pomp/src/pomp_internal.h	2013-01-15 15:26:18 UTC (rev 817)
@@ -10,6 +10,7 @@
 
 #include "pomp.h"
 
+# define MATCHNAMES(X,N) (matchnames(GET_NAMES(X),(N)))
 # define MATCHROWNAMES(X,N) (matchnames(GET_ROWNAMES(GET_DIMNAMES(X)),(N)))
 # define MATCHCOLNAMES(X,N) (matchnames(GET_COLNAMES(GET_DIMNAMES(X)),(N)))
 # define MATCH_CHAR_TO_ROWNAMES(X,N,A) (match_char_to_names(GET_ROWNAMES(GET_DIMNAMES(X)),(N),(A)))