[Dplr-commits] r676 - in branches/redfit: R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Sep 3 13:45:50 CEST 2013
Author: mvkorpel
Date: 2013-09-03 13:45:50 +0200 (Tue, 03 Sep 2013)
New Revision: 676
Modified:
branches/redfit/R/redfit.R
branches/redfit/src/redfit.c
Log:
* Removed parameter 'cbind' from C function spectr():
lookup now done in C code
* In C function rmtrend(), use R version of length():
theoretical support for long vectors while not requiring R >= 3.0.0
Modified: branches/redfit/R/redfit.R
===================================================================
--- branches/redfit/R/redfit.R 2013-08-30 20:48:31 UTC (rev 675)
+++ branches/redfit/R/redfit.R 2013-09-03 11:45:50 UTC (rev 676)
@@ -216,11 +216,10 @@
ia <- redfitInitArrays(t, freq, params)
## determine autospectrum of input data
dn50 <- as.numeric(n50)
- cbindfun <- match.fun("cbind")
lmfitfun <- tryCatch(match.fun(".lm.fit"),
error = function(...) match.fun("lm.fit"))
gxx <- .Call(dplR.spectr, t, x, np, ia[[1]], ia[[2]], ia[[3]], ia[[4]],
- nseg, nfreq, avgdt, freq, dn50, segskip, cbindfun, lmfitfun)
+ nseg, nfreq, avgdt, freq, dn50, segskip, lmfitfun)
## estimate data variance from autospectrum
varx <- df * sum(gxx)
## dplR: estimate lag-1 autocorrelation coefficient unless prescribed
@@ -244,7 +243,7 @@
grr[, i] <-
.Call(dplR.spectr, t, .Call(dplR.makear1, difft, np, tau), np,
ia[[1]], ia[[2]], ia[[3]], ia[[4]], nseg, nfreq, avgdt,
- freq, dn50, segskip, cbindfun, lmfitfun)
+ freq, dn50, segskip, lmfitfun)
## scale and sum red-noise spectra
varr1 <- df * sum(grr[, i])
grr[, i] <- varx / varr1 * grr[, i]
@@ -259,7 +258,7 @@
## setup AR(1) time series and estimate its spectrum
grr <- .Call(dplR.spectr, t, .Call(dplR.makear1, difft, np, tau),
np, ia[[1]], ia[[2]], ia[[3]], ia[[4]], nseg, nfreq,
- avgdt, freq, dn50, segskip, cbindfun, lmfitfun)
+ avgdt, freq, dn50, segskip, lmfitfun)
## scale and sum red-noise spectra
varr1 <- df * sum(grr)
grr <- varx / varr1 * grr
Modified: branches/redfit/src/redfit.c
===================================================================
--- branches/redfit/src/redfit.c 2013-08-30 20:48:31 UTC (rev 675)
+++ branches/redfit/src/redfit.c 2013-09-03 11:45:50 UTC (rev 676)
@@ -24,10 +24,10 @@
#include <complex.h>
SEXP seg50(SEXP k, SEXP nseg, SEXP segskip, SEXP np);
-void rmtrend(SEXP x, SEXP y, SEXP lmfit);
+void rmtrend(SEXP x, SEXP y, SEXP lengthfun, SEXP lmfit);
SEXP spectr(SEXP t, SEXP x, SEXP np, SEXP ww, SEXP tsin, SEXP tcos, SEXP wtau,
SEXP nseg, SEXP nfreq, SEXP avgdt, SEXP freq, SEXP n50,
- SEXP segskip, SEXP cbind, SEXP lmfit);
+ SEXP segskip, SEXP lmfit);
void ftfix(const double *xx, const double *tsamp, const size_t nxx,
const double *freq, const size_t nfreq, const double si,
const size_t lfreq, const double tzero, const double *tcos,
@@ -87,11 +87,13 @@
/* dplR: y <- lmfit(x, y)[["residuals"]]
*/
-void rmtrend(SEXP x, SEXP y, SEXP lmfit) {
+void rmtrend(SEXP x, SEXP y, SEXP lengthfun, SEXP lmfit) {
SEXP tmp, lmcall, lmres, lmnames, rduals;
+ SEXP sn, ncall;
+ PROTECT_INDEX ipx;
double *rdualsptr, *y_data;
- int i, nameslength;
- int n = 0;
+ size_t i, nameslength;
+ size_t n = 0;
Rboolean found = FALSE;
Rboolean mismatch = TRUE;
@@ -105,7 +107,14 @@
/* dplR: get residuals from the list given by lm.fit(x, y) */
lmnames = getAttrib(lmres, R_NamesSymbol);
- nameslength = length(lmnames);
+ PROTECT(tmp = ncall = allocList(2));
+ SET_TYPEOF(ncall, LANGSXP);
+ SETCAR(tmp, lengthfun); tmp = CDR(tmp);
+ SETCAR(tmp, lmnames);
+ PROTECT_WITH_INDEX(sn = eval(ncall, R_GlobalEnv), &ipx);
+ REPROTECT(sn = coerceVector(sn, REALSXP), ipx);
+ nameslength = (size_t) *REAL(sn);
+ UNPROTECT(2);
for (i = 0; i < nameslength; i++) {
if (strcmp(CHAR(STRING_ELT(lmnames, i)), "residuals") == 0) {
rduals = VECTOR_ELT(lmres, i);
@@ -115,10 +124,24 @@
}
}
+ /* dplR: compare length of y with length of residuals */
+ PROTECT(tmp = ncall = allocList(2));
+ SET_TYPEOF(ncall, LANGSXP);
+ SETCAR(tmp, lengthfun); tmp = CDR(tmp);
+ SETCAR(tmp, y);
+ PROTECT_WITH_INDEX(sn = eval(ncall, R_GlobalEnv), &ipx);
+ REPROTECT(sn = coerceVector(sn, REALSXP), ipx);
+ n = (size_t) *REAL(sn);
+ UNPROTECT(1);
if (found) {
- n = length(rduals);
- mismatch = n != length(y);
+ SETCAR(tmp, rduals);
+ PROTECT_WITH_INDEX(sn = eval(ncall, R_GlobalEnv), &ipx);
+ REPROTECT(sn = coerceVector(sn, REALSXP), ipx);
+ mismatch = n != (size_t) *REAL(sn);
+ UNPROTECT(1);
}
+ UNPROTECT(1);
+
y_data = REAL(y);
if (!mismatch) {
rdualsptr = REAL(rduals);
@@ -127,7 +150,6 @@
y_data[i] = rdualsptr[i];
}
} else {
- n = length(y);
for (i = 0; i < n; i++) {
y_data[i] = NA_REAL;
}
@@ -143,8 +165,8 @@
*/
SEXP spectr(SEXP t, SEXP x, SEXP np, SEXP ww, SEXP tsin, SEXP tcos, SEXP wtau,
SEXP nseg, SEXP nfreq, SEXP avgdt, SEXP freq, SEXP n50,
- SEXP segskip, SEXP cbind, SEXP lmfit) {
- SEXP gxx, twk, xwk, ftrx, ftix, tmp, cbindcall;
+ SEXP segskip, SEXP lmfit) {
+ SEXP gxx, twk, xwk, ftrx, ftix, tmp, cbindcall, lengthfun;
double dnseg, segskip_val, scal, np_val;
long double sumx, sqrt_nseg;
size_t i, j, nseg_val, nfreq_val, n50_val, segstart;
@@ -182,7 +204,7 @@
* though.*/
PROTECT(tmp = cbindcall = allocList(3));
SET_TYPEOF(cbindcall, LANGSXP);
- SETCAR(tmp, cbind); tmp = CDR(tmp);
+ SETCAR(tmp, install("cbind")); tmp = CDR(tmp);
SETCAR(tmp, ScalarReal(1.0)); tmp = CDR(tmp);
SETCAR(tmp, twk);
REPROTECT(twk = eval(cbindcall, R_GlobalEnv), pidx);
@@ -204,6 +226,7 @@
for (i = 0; i < nfreq_val; i++) {
gxx_data[i] = 0.0;
}
+ lengthfun = install("length");
for (i = 0; i < n50_val; i++) {
/* copy data of i'th segment into workspace */
segstart = (size_t) segfirst((double) i, segskip_val, np_val, dnseg);
@@ -212,7 +235,7 @@
xwk_data[j] = x_data[segstart + j];
}
/* detrend data */
- rmtrend(twk, xwk, lmfit);
+ rmtrend(twk, xwk, lengthfun, lmfit);
/* apply window to data */
sumx = 0.0L;
for (j = 0; j < nseg_val; j++) {
More information about the Dplr-commits
mailing list