[Pomp-commits] r351 - in pkg: R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Sep 29 23:06:10 CEST 2010
Author: kingaa
Date: 2010-09-29 23:06:10 +0200 (Wed, 29 Sep 2010)
New Revision: 351
Added:
pkg/src/pomp_mat.h
Removed:
pkg/src/mat.c
Modified:
pkg/R/basic-probes.R
pkg/src/pomp_internal.h
pkg/src/probe_marginal.c
pkg/src/probe_nlar.c
Log:
- remove 'mat.c' in favor of inline functions for regression, defined in the new header 'pomp_mat.h'
Modified: pkg/R/basic-probes.R
===================================================================
--- pkg/R/basic-probes.R 2010-09-29 17:49:45 UTC (rev 350)
+++ pkg/R/basic-probes.R 2010-09-29 21:06:10 UTC (rev 351)
@@ -108,18 +108,6 @@
}
}
-probe.marginal <- function (var, ref, order = 3, diff = 1, transform = identity) {
- if (length(var)>1) stop(sQuote("probe.marginal")," is a univariate probe")
- transform <- match.fun(transform)
- setup <- .Call(probe_marginal_setup,transform(ref),order,diff)
- function (y) .Call(
- probe_marginal_solve,
- x=transform(y[var,]),
- setup=setup,
- diff=diff
- )
-}
-
probe.acf <- function (var, lag.max, type = c("covariance", "correlation"), transform = identity) {
type <- match.arg(type)
transform <- match.fun(transform)
@@ -132,7 +120,20 @@
)
}
+probe.marginal <- function (var, ref, order = 3, diff = 1, transform = identity) {
+ if (length(var)>1) stop(sQuote("probe.marginal")," is a univariate probe")
+ transform <- match.fun(transform)
+ setup <- .Call(probe_marginal_setup,transform(ref),order,diff)
+ function (y) .Call(
+ probe_marginal_solve,
+ x=transform(y[var,,drop=TRUE]),
+ setup=setup,
+ diff=diff
+ )
+}
+
probe.nlar <- function (var, lags, powers, transform = identity) {
+ if (length(var)>1) stop(sQuote("probe.nlar")," is a univariate probe")
transform <- match.fun(transform)
if (any(lags<1)||any(powers<1))
stop(sQuote("lags")," and ",sQuote("powers")," must be positive integers")
@@ -147,7 +148,7 @@
powers <- as.integer(powers)
function (y) .Call(
probe_nlar,
- x=transform(y[var,,drop=FALSE]),
+ x=transform(y[var,,drop=TRUE]),
lags=lags,
powers=powers
)
Deleted: pkg/src/mat.c
===================================================================
--- pkg/src/mat.c 2010-09-29 17:49:45 UTC (rev 350)
+++ pkg/src/mat.c 2010-09-29 21:06:10 UTC (rev 351)
@@ -1,55 +0,0 @@
-// -*- mode: C++; -*-
-// some functions for linear algebra written by Simon N. Wood.
-// codes have been modified by AAK to suit his own peccadilloes
-
-#include "pomp_internal.h"
-#include <R_ext/Linpack.h>
-#include <R_ext/Lapack.h>
-
-void pomp_backsolve (double *R, int *r, int *c, double *B, double *C, int *bc) {
- double x;
- int i, j, k;
- for (j = 0; j < *bc; j++) { // work across columns of B & C
- for (i = *c-1; i >= 0; i--) { // work up each column of B & C
- for (k = i+1, x = 0.0; k < *c; k++)
- x += R[i+*r*k]*C[k+*c*j];
- C[i+*c*j] = (B[i+*c*j]-x)/R[i+*r*i];
- }
- }
-}
-
-void pomp_qr (double *x, int *r, int *c, int *pivot, double *tau) {
- 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++;
- {
- 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;
- int lda, info, lwork = -1;
- double work1;
-
- 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++;
- {
- double work[lwork];
- // actual call
- F77_NAME(dormqr)(&side,&trans,r,c,k,a,&lda,tau,b,r,work,&lwork,&info);
- }
-}
Modified: pkg/src/pomp_internal.h
===================================================================
--- pkg/src/pomp_internal.h 2010-09-29 17:49:45 UTC (rev 350)
+++ pkg/src/pomp_internal.h 2010-09-29 21:06:10 UTC (rev 351)
@@ -149,8 +149,4 @@
#endif
-void pomp_backsolve(double *, int *, int *, double *, double *, int *);
-void pomp_qr(double *, int *, int *, int *, double *);
-void pomp_qrqy(double *, double *, double *, int *, int *, int *, int *, int *);
-
#endif
Added: pkg/src/pomp_mat.h
===================================================================
--- pkg/src/pomp_mat.h (rev 0)
+++ pkg/src/pomp_mat.h 2010-09-29 21:06:10 UTC (rev 351)
@@ -0,0 +1,47 @@
+#ifndef _POMP_MAT_H_
+#define _POMP_MAT_H_
+
+#include "pomp_internal.h"
+#include <R_ext/Lapack.h>
+
+static R_INLINE void pomp_backsolve (double *a, int m, int n, double *x, int incx) {
+ // Level 2 BLAS triangular-matrix solver
+ // DTRSV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX)
+ F77_NAME(dtrsv)("U","N","N",&n,a,&m,x,&incx);
+}
+
+static R_INLINE void pomp_qr (double *a, int m, int n, int *pivot, double *tau) {
+ int info, j, lwork = -1;
+ double work1;
+ for (j = 0; j < n; j++) pivot[j] = 0; // zero out the pivots (assumed by DGEQP3)
+ // LAPACK QR decomposition routine
+ // DGEQP3(M,N,A,LDA,JPVT,TAU,WORK,LWORK,INFO)
+ F77_NAME(dgeqp3)(&m,&n,a,&m,pivot,tau,&work1,&lwork,&info); // workspace query
+ lwork = (int) ceil(work1);
+ {
+ double work[lwork];
+ F77_NAME(dgeqp3)(&m,&n,a,&m,pivot,tau,work,&lwork,&info); // actual call
+ }
+ for (j = 0; j < n; j++) pivot[j]--;
+}
+
+static R_INLINE void pomp_qrqy (double *c, double *a, double *tau, int m, int n, int k, int left, int tp) {
+ int lda, info, lwork = -1;
+ char side, trans;
+ double work1;
+
+ side = (left) ? 'L' : 'R';
+ lda = (left) ? m : n;
+ trans = (tp) ? 'T' : 'N';
+
+ // workspace query
+ // DORMQR(SIDE,TRANS,M,N,K,A,LDA,TAU,C,LDC,WORK,LWORK,INFO)
+ F77_NAME(dormqr)(&side,&trans,&m,&n,&k,a,&lda,tau,c,&m,&work1,&lwork,&info);
+ lwork = (int) ceil(work1);
+ { // actual call
+ double work[lwork];
+ F77_NAME(dormqr)(&side,&trans,&m,&n,&k,a,&lda,tau,c,&m,work,&lwork,&info);
+ }
+}
+
+#endif
Property changes on: pkg/src/pomp_mat.h
___________________________________________________________________
Added: svn:eol-style
+ native
Modified: pkg/src/probe_marginal.c
===================================================================
--- pkg/src/probe_marginal.c 2010-09-29 17:49:45 UTC (rev 350)
+++ pkg/src/probe_marginal.c 2010-09-29 21:06:10 UTC (rev 351)
@@ -1,9 +1,8 @@
// -*- mode: C++ -*-
#include "pomp_internal.h"
+#include "pomp_mat.h"
#include <stdio.h>
-#include <R_ext/Linpack.h>
-#include <R_ext/Lapack.h>
static void order_reg_solve (double *beta, double *x, double *mm, double *tau, int *pivot, int n, int np, int diff);
static void order_reg_model_matrix (double *z, double *X, double *tau, int *pivot, int n, int np, int diff);
@@ -111,15 +110,14 @@
for (j = 0; j < nx; j++) X2[j] = X1[j]*z[j];
// QR decompose the model matrix
- for (i = 0; i < np; i++) pivot[i] = 0;
- pomp_qr(X,&nx,&np,pivot,tau);
+ pomp_qr(X,nx,np,pivot,tau);
}
// thanks to Simon N. Wood for the original version of the following code
static void order_reg_solve (double *beta, double *x, double *mm, double *tau, int *pivot, int n, int np, int diff) {
int nx, one = 1;
- double xx, coef[np];
+ double xx;
int i, j;
for (i = 0, nx = n; i < diff; i++) { // differencing loop
@@ -137,10 +135,10 @@
R_qsort(x,1,nx);
// solve R b = Q'x for b
- pomp_qrqy(x,mm,tau,&nx,&one,&np,&one,&one); // y <- Q'y
- pomp_backsolve(mm,&nx,&np,x,coef,&one); // b <- R^{-1} Q'y
+ pomp_qrqy(x,mm,tau,nx,1,np,1,1); // y <- Q'y
+ pomp_backsolve(mm,nx,np,x,1); // y <- R^{-1} Q'y
- // unpivot beta
- for (i = 0; i < np; i++) beta[i] = coef[pivot[i]];
+ // unpivot and store the coefficients in beta
+ for (i = 0; i < np; i++) beta[pivot[i]] = x[i];
}
Modified: pkg/src/probe_nlar.c
===================================================================
--- pkg/src/probe_nlar.c 2010-09-29 17:49:45 UTC (rev 350)
+++ pkg/src/probe_nlar.c 2010-09-29 21:06:10 UTC (rev 351)
@@ -1,9 +1,10 @@
// -*- mode: C++; -*-
#include "pomp_internal.h"
+#include "pomp_mat.h"
#include <stdio.h>
-static void pomp_nlar(double *beta, double *y, int n, int nterms, int *lag, int *power, double *X);
+static void poly_nlar_fit(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;
@@ -20,7 +21,7 @@
PROTECT(beta = NEW_NUMERIC(nterms)); nprotect++;
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);
+ poly_nlar_fit(REAL(beta),REAL(y),n,nterms,INTEGER(lags),INTEGER(powers),mm);
PROTECT(beta_names = NEW_STRING(nterms)); nprotect++;
for (k = 0; k < nterms; k++) {
@@ -36,8 +37,8 @@
// Code for polynomial auto-regression
// 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, double *X) {
+static void poly_nlar_fit (double *beta, double *y, int n,
+ 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.
@@ -114,18 +115,17 @@
// response vector is now length nx
{
- double tau[nterms], b[nterms];
+ double tau[nterms];
int pivot[nterms];
- // QR decompose the model matrix
- for (i = 0; i < nterms; i++) pivot[i] = 0;
- pomp_qr(X,&nx,&nterms,pivot,tau);
- // solve R b = Q'y for b
- pomp_qrqy(yp,X,tau,&nx,&one,&nterms,&one,&one); // y <- Q'y
- pomp_backsolve(X,&nx,&nterms,yp,b,&one); // b <- R^{-1} Q'y
+ // first QR decompose the model matrix
+ pomp_qr(X,nx,nterms,pivot,tau);
+ // then solve R b = Q'y for b
+ pomp_qrqy(yp,X,tau,nx,1,nterms,1,1); // y <- Q'y
+ pomp_backsolve(X,nx,nterms,yp,1); // y <- R^{-1} Q'y
- // store b in kth column of beta
- for (i = 0; i < nterms; i++) beta[pivot[i]] = b[i];
+ // unpivot and store coefficients in beta
+ for (i = 0; i < nterms; i++) beta[pivot[i]] = yp[i];
}
More information about the pomp-commits
mailing list