[Pomp-commits] r350 - pkg/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Sep 29 19:49:45 CEST 2010
Author: kingaa
Date: 2010-09-29 19:49:45 +0200 (Wed, 29 Sep 2010)
New Revision: 350
Modified:
pkg/src/euler.c
pkg/src/mat.c
pkg/src/pfilter.c
pkg/src/probe.c
pkg/src/probe_acf.c
pkg/src/probe_marginal.c
pkg/src/probe_nlar.c
Log:
- remove the last usages of dynamic allocation from 'mat.c'
- changes in the compiled codes to eliminate pedantic compiler warnings
Modified: pkg/src/euler.c
===================================================================
--- pkg/src/euler.c 2010-09-29 17:47:21 UTC (rev 349)
+++ pkg/src/euler.c 2010-09-29 17:49:45 UTC (rev 350)
@@ -19,7 +19,7 @@
int covdim = ndim[5];
int nzero = ndim[6];
double covar_fn[covdim];
- int j, k, p, step, neuler;
+ int j, k, p, step, neuler = 0;
double dt, tol;
struct lookup_table covariate_table = {covlen, covdim, 0, time_table, covar_table};
Modified: pkg/src/mat.c
===================================================================
--- pkg/src/mat.c 2010-09-29 17:47:21 UTC (rev 349)
+++ pkg/src/mat.c 2010-09-29 17:49:45 UTC (rev 350)
@@ -2,11 +2,9 @@
// some functions for linear algebra written by Simon N. Wood.
// codes have been modified by AAK to suit his own peccadilloes
-#include <math.h>
-#include <R.h>
+#include "pomp_internal.h"
#include <R_ext/Linpack.h>
#include <R_ext/Lapack.h>
-#include "pomp_internal.h"
void pomp_backsolve (double *R, int *r, int *c, double *B, double *C, int *bc) {
double x;
@@ -21,40 +19,37 @@
}
void pomp_qr (double *x, int *r, int *c, int *pivot, double *tau) {
- int info, lwork = -1, i;
- double work1, *work;
+ int info, j, lwork = -1;
+ double work1;
// workspace query
F77_NAME(dgeqp3)(r,c,x,r,pivot,tau,&work1,&lwork,&info);
lwork = (int) floor(work1);
- if ((work1-lwork) >0.5) lwork++;
- work = (double *) Calloc(lwork,double);
- // actual call
- F77_NAME(dgeqp3)(r,c,x,r,pivot,tau,work,&lwork,&info);
- Free(work);
- for (i = 0; i < *c; i++) pivot[i]--;
+ if ((work1-lwork) > 0.5) lwork++;
+ {
+ double work[lwork];
+ // actual call
+ F77_NAME(dgeqp3)(r,c,x,r,pivot,tau,work,&lwork,&info);
+ }
+ for (j = 0; j < *c; j++) pivot[j]--;
// ... for 'tis C in which we work and not the 'cursed Fortran...
}
void pomp_qrqy (double *b, double *a, double *tau, int *r, int *c, int *k, int *left, int *tp) {
- char side, trans='N';
- int lda, lwork = -1, info;
- double *work, work1;
+ char side, trans;
+ int lda, info, lwork = -1;
+ double work1;
- if (! *left) {
- side='R';
- lda = *c;
- } else {
- side = 'L';
- lda= *r;
- }
- if (*tp) trans='T';
+ side = (*left) ? 'L' : 'R';
+ lda = (*left) ? *r : *c;
+ trans = (*tp) ? 'T' : 'N';
// workspace query
F77_NAME(dormqr)(&side,&trans,r,c,k,a,&lda,tau,b,r,&work1,&lwork,&info);
lwork = (int) floor(work1);
if ((work1-lwork) > 0.5) lwork++;
- work = (double *) Calloc(lwork,double);
- // actual call
- F77_NAME(dormqr)(&side,&trans,r,c,k,a,&lda,tau,b,r,work,&lwork,&info);
- Free(work);
+ {
+ double work[lwork];
+ // actual call
+ F77_NAME(dormqr)(&side,&trans,r,c,k,a,&lda,tau,b,r,work,&lwork,&info);
+ }
}
Modified: pkg/src/pfilter.c
===================================================================
--- pkg/src/pfilter.c 2010-09-29 17:47:21 UTC (rev 349)
+++ pkg/src/pfilter.c 2010-09-29 17:49:45 UTC (rev 350)
@@ -16,13 +16,12 @@
SEXP filtmean, SEXP weights, SEXP tol)
{
int nprotect = 0;
- SEXP pm, pv, fm, ess, fail, loglik;
+ SEXP pm = R_NilValue, pv = R_NilValue, fm = R_NilValue, ess, fail, loglik;
SEXP retval, retvalnames;
- double *xpm, *xpv, *xfm, *xw;
+ double *xpm = 0, *xpv = 0, *xfm = 0, *xw = 0, *xx = 0, *xp = 0;
SEXP dimX, Pnames, pindex;
- double *xx, *xp;
int *dim, *pidx, lv;
- int nvars, npars, nrw, nreps, offset, nlost;
+ int nvars, npars = 0, nrw = 0, nreps, offset, nlost;
int do_rw, do_pm, do_pv, do_fm, all_fail = 0;
double sum, sumsq, vsq, ws, w, toler;
int j, k;
Modified: pkg/src/probe.c
===================================================================
--- pkg/src/probe.c 2010-09-29 17:47:21 UTC (rev 349)
+++ pkg/src/probe.c 2010-09-29 17:49:45 UTC (rev 350)
@@ -33,7 +33,7 @@
int nprobe, nsims, nvars, ntimes, nvals;
int xdim[2];
double *xp, *yp;
- int p, s, i, j, k, len0, len = 0;
+ int p, s, i, j, k, len0 = 0, len = 0;
PROTECT(nsim = AS_INTEGER(nsim)); nprotect++;
if ((LENGTH(nsim)>1) || (INTEGER(nsim)[0]<=0))
@@ -79,6 +79,8 @@
for (p = 0, k = 0; p < nprobe; p++, k += len) { // loop over probes
+ R_CheckUserInterrupt();
+
for (s = 0; s < nsims; s++) { // loop over simulations
// copy the data from y[,s,-1] to x[,]
Modified: pkg/src/probe_acf.c
===================================================================
--- pkg/src/probe_acf.c 2010-09-29 17:47:21 UTC (rev 349)
+++ pkg/src/probe_acf.c 2010-09-29 17:49:45 UTC (rev 350)
@@ -8,25 +8,24 @@
SEXP probe_acf (SEXP x, SEXP lag_max, SEXP corr) {
int nprotect = 0;
SEXP acf, acf_names, cacf;
- SEXP X_names, Dim;
+ SEXP X, X_names;
int maxlag, correlation, nvars, n;
int j, k, l;
double *p, *p1;
- char tmp[BUFSIZ];
+ char tmp[BUFSIZ], *nm;
maxlag = *(INTEGER(AS_INTEGER(lag_max))); // maximum lag
correlation = *(INTEGER(AS_INTEGER(corr))); // correlation, or covariance?
- PROTECT(Dim = GET_DIM(x)); nprotect++;
- nvars = INTEGER(Dim)[0]; // nvars = # of variables
- n = INTEGER(Dim)[1]; // n = # of observations
- // depending on how this function is called, it may be necessary to duplicate x in the next line
- PROTECT(x = AS_NUMERIC(x)); nprotect++;
+ nvars = INTEGER(GET_DIM(x))[0]; // nvars = # of variables
+ n = INTEGER(GET_DIM(x))[1]; // n = # of observations
+
+ PROTECT(X = duplicate(AS_NUMERIC(x))); nprotect++;
PROTECT(X_names = GET_ROWNAMES(GET_DIMNAMES(x))); nprotect++;
PROTECT(acf = NEW_NUMERIC((maxlag+1)*nvars)); nprotect++;
- pomp_acf(REAL(acf),REAL(x),n,nvars,maxlag,correlation);
+ pomp_acf(REAL(acf),REAL(X),n,nvars,maxlag,correlation);
if (correlation) {
PROTECT(cacf = NEW_NUMERIC(maxlag*nvars)); nprotect++;
@@ -34,7 +33,8 @@
for (j = 0, l = 0, p = REAL(cacf), p1 = REAL(acf)+1; j < nvars; j++, p1++) {
for (k = 1; k <= maxlag; k++, p++, p1++, l++) {
*p = *p1; // copy correlations from acf into cacf
- snprintf(tmp,BUFSIZ,"acf.%ld.%s",k,CHARACTER_DATA(STRING_ELT(X_names,j)));
+ nm = (char *) CHARACTER_DATA(STRING_ELT(X_names,j));
+ snprintf(tmp,BUFSIZ,"acf.%d.%s",k,nm);
SET_STRING_ELT(acf_names,l,mkChar(tmp));
}
}
@@ -43,7 +43,8 @@
PROTECT(acf_names = NEW_STRING((maxlag+1)*nvars)); nprotect++;
for (j = 0, l = 0; j < nvars; j++) {
for (k = 0; k <= maxlag; k++, l++) {
- snprintf(tmp,BUFSIZ,"acf.%ld.%s",k,CHARACTER_DATA(STRING_ELT(X_names,j)));
+ nm = (char *) CHARACTER_DATA(STRING_ELT(X_names,j));
+ snprintf(tmp,BUFSIZ,"acf.%d.%s",k,nm);
SET_STRING_ELT(acf_names,l,mkChar(tmp));
}
}
@@ -58,41 +59,42 @@
// vectorized routine for ACF calculation
// thanks to Simon N. Wood for the original version of this code
// modifications due to AAK
+// note that the behavior of this ACF is slightly different from that of R's 'acf' function
+// we first center the series and then compute means of products
static void pomp_acf (double *acf, double *x, int n, int nvars, int maxlag, int correlation) {
- double sum, *p, *p0, *p1, *p2;
+ double xx, *p, *p0, *p1, *p2;
int j, k, lag, ct;
// first center each row
for (j = 0, p = x; j < nvars; j++, p++) {
- for (k = 0, p0 = p, sum = 0, ct = 0; k < n; p0 += nvars, k++) {
+ for (k = 0, p0 = p, xx = 0, ct = 0; k < n; p0 += nvars, k++) {
if (R_FINITE(*p0)) {
- sum += *p0;
+ xx += *p0;
ct++;
}
}
if (ct < 1) error("series %ld has no data",j);
- sum /= ((double) ct); // mean of x[j,]
+ xx /= ct; // mean of x[j,]
for (k = 0, p0 = p; k < n; p0 += nvars, k++)
- if (R_FINITE(*p0)) *p0 -= sum;
+ if (R_FINITE(*p0)) *p0 -= xx;
}
// compute covariances
for (j = 0, p0 = x, p = acf; j < nvars; j++, p0++) { // loop over series
for (lag = 0; lag <= maxlag; lag++, p++) { // loop over lags
- for (k = 0, ct = 0, sum = 0, p1 = p0, p2 = p0+lag*nvars; k < n-lag; k++, p1 += nvars, p2 += nvars)
+ for (k = 0, ct = 0, xx = 0, p1 = p0, p2 = p0+lag*nvars; k < n-lag; k++, p1 += nvars, p2 += nvars)
if (R_FINITE(*p1) && R_FINITE(*p2)) {
- sum += (*p1)*(*p2);
+ xx += (*p1)*(*p2);
ct++;
}
- *p = (ct > 0) ? sum/((double) ct) : R_NaReal;
- //// *p = (ct > 0) ? sum/((double) n) : R_NaReal; // strangely, this is apparently what R's 'acf' function does
+ *p = (ct > 0) ? xx/ct : R_NaReal;
}
}
// convert to correlations if desired
if (correlation)
for (j = 0, p = acf; j < nvars; j++, p += maxlag+1)
- for (lag = 0, p0 = p, sum = *p; lag <= maxlag; lag++, p0++)
- *p0 /= sum;
+ for (lag = 0, p0 = p, xx = *p; lag <= maxlag; lag++, p0++)
+ *p0 /= xx;
}
Modified: pkg/src/probe_marginal.c
===================================================================
--- pkg/src/probe_marginal.c 2010-09-29 17:47:21 UTC (rev 349)
+++ pkg/src/probe_marginal.c 2010-09-29 17:49:45 UTC (rev 350)
@@ -44,7 +44,6 @@
SEXP probe_marginal_solve (SEXP x, SEXP setup, SEXP diff) {
int nprotect = 0;
- SEXP dimx, dimm;
SEXP X, mm, tau, pivot, beta, beta_names;
int n, nx, np, df;
int i;
@@ -58,9 +57,8 @@
PROTECT(tau = VECTOR_ELT(setup,1)); nprotect++; // diagonals
PROTECT(pivot = VECTOR_ELT(setup,2)); nprotect++; // pivots
- PROTECT(dimm = GET_DIM(mm)); nprotect++;
- nx = INTEGER(dimm)[0]; // nx = number of rows in model matrix
- np = INTEGER(dimm)[1]; // np = order of polynomial
+ nx = INTEGER(GET_DIM(mm))[0]; // nx = number of rows in model matrix
+ np = INTEGER(GET_DIM(mm))[1]; // np = order of polynomial
if (n-df != nx) error("length of 'ref' must equal length of data");
PROTECT(X = duplicate(AS_NUMERIC(x))); nprotect++;
@@ -68,7 +66,7 @@
PROTECT(beta = NEW_NUMERIC(np)); nprotect++;
PROTECT(beta_names = NEW_STRING(np)); nprotect++;
for (i = 0; i < np; i++) {
- snprintf(tmp,BUFSIZ,"marg.%ld",i+1);
+ snprintf(tmp,BUFSIZ,"marg.%d",i+1);
SET_STRING_ELT(beta_names,i,mkChar(tmp));
}
SET_NAMES(beta,beta_names);
Modified: pkg/src/probe_nlar.c
===================================================================
--- pkg/src/probe_nlar.c 2010-09-29 17:47:21 UTC (rev 349)
+++ pkg/src/probe_nlar.c 2010-09-29 17:49:45 UTC (rev 350)
@@ -3,13 +3,14 @@
#include "pomp_internal.h"
#include <stdio.h>
-static void pomp_nlar(double *beta, double *y, int n, int nterms, int *lag, int *power);
+static void pomp_nlar(double *beta, double *y, int n, int nterms, int *lag, int *power, double *X);
SEXP probe_nlar (SEXP x, SEXP lags, SEXP powers) {
int nprotect = 0;
SEXP y, beta, beta_names;
int n, nterms;
int k;
+ double *mm;
char tmp[BUFSIZ];
n = LENGTH(x); // n = # of observations
@@ -17,12 +18,13 @@
PROTECT(y = duplicate(AS_NUMERIC(x))); nprotect++;
PROTECT(beta = NEW_NUMERIC(nterms)); nprotect++;
-
- pomp_nlar(REAL(beta),REAL(y),n,nterms,INTEGER(lags),INTEGER(powers));
+ mm = (double *) R_alloc(n*nterms,sizeof(double)); // storage for the model matrix
+ pomp_nlar(REAL(beta),REAL(y),n,nterms,INTEGER(lags),INTEGER(powers),mm);
+
PROTECT(beta_names = NEW_STRING(nterms)); nprotect++;
for (k = 0; k < nterms; k++) {
- snprintf(tmp,BUFSIZ,"nlar.%ld^%ld",INTEGER(lags)[k],INTEGER(powers)[k]);
+ snprintf(tmp,BUFSIZ,"nlar.%d^%d",INTEGER(lags)[k],INTEGER(powers)[k]);
SET_STRING_ELT(beta_names,k,mkChar(tmp));
}
SET_NAMES(beta,beta_names);
@@ -35,7 +37,7 @@
// The original version of the following code is due to Simon N. Wood.
// Modifications by AAK.
static void pomp_nlar(double *beta, double *y, int n,
- int nterms, int *lag, int *power) {
+ int nterms, int *lag, int *power, double *X) {
// 'x' is an n vector of data.
// 'nterms' gives the number of terms on the rhs of the autoregression.
// 'lag[i]' gives the lag of the ith term on the rhs.
@@ -81,21 +83,19 @@
} else { // data not all the same
- double *X, *Xp;
+ double *Xp;
int finite[ny];
// test for NA rows in model matrix and response vector
for (nx = 0, yp = y+maxlag, j = 0; j < ny; j++) {
finite[j] = (R_FINITE(yp[j])) ? 1 : 0; // finite response?
- for (i = 0; i < nterms; i++) {
+ for (i = 0; i < nterms; i++)
finite[j] = (R_FINITE(yp[j-lag[i]])) ? finite[j] : 0; // finite in model matrix row j, column i
- }
if (finite[j]) nx++;
}
// nx is now the number of non-NA rows in the model matrix
// build the model matrix, omitting NA rows
- X = (double *) Calloc(nx*nterms,double);
for (Xp = X, i = 0; i < nterms; i++) { // work through the terms
for (j = 0; j < ny; j++) {
if (finite[j]) {
@@ -109,11 +109,8 @@
// X is now the nx by nterms model matrix
// drop the NA rows from the response data
- for (i = 0, j = 0; j < ny; j++) {
- if (finite[j]) { // keep this row
- yp[i++] = yp[j];
- }
- }
+ for (i = 0, j = 0; j < ny; j++)
+ if (finite[j]) yp[i++] = yp[j]; // keep this row
// response vector is now length nx
{
@@ -132,8 +129,6 @@
}
- Free(X);
-
}
}
More information about the pomp-commits
mailing list