5 Appendix: Code
5.1 R Code
# select
aplSelect <- function(a, x, drop = FALSE) {
sa <- aplShape(a)
ra <- aplRank(a)
sz <- sapply(x, length)
z <- array(0, sz)
nz <- prod(sz)
for (i in 1:nz) {
ivec <- aplEncode(i, sz)
jvec <- vector()
for (j in 1:ra)
jvec <- c(jvec, x[[j]][ivec[j]])
z[i] <- a[aplDecode(jvec, sa)]
}
if (drop)
return(drop(z))
else
return(z)
}
# drop
aplDrop <- function(a, x, drop = FALSE) {
sa <- aplShape(a)
ra <- aplRank(a)
y <- as.list(rep(0, ra))
for (i in 1:ra) {
ss <- sa[i]
xx <- x[i]
sx <- ss + xx
if (xx >= 0)
y[[i]] <- (xx + 1):ss
if (xx < 0)
y[[i]] <- 1:sx
}
return(aplSelect(a, y, drop))
}
# take
aplTake <- function(a, x, drop = FALSE) {
sa <- aplShape(a)
ra <- aplRank(a)
y <- as.list(rep(0, ra))
for (i in 1:ra) {
ss <- sa[i]
xx <- x[i]
sx <- ss + xx
if (xx > 0)
y[[i]] <- 1:xx
if (xx < 0)
y[[i]] <- (sx + 1):ss
}
return(aplSelect(a, y, drop))
}
# reduce vector
aplRDV <- function(x, f = "+") {
if (length(x) == 0)
return(x)
s <- x[1]
if (length(x) == 1)
return(s)
for (i in 2:length(x))
s <- match.fun(f)(s, x[i])
return(s)
}
# scan vector
aplSCV <- function(x, f = "+") {
if (length(x) <= 1)
return(x)
return(sapply(1:length(x), function(i)
aplRDV(x[1:i], f)))
}
# inner product vector
aplIPV <- function(x, y, f = "*", g = "+") {
if (length(x) != length(y))
stop("Incorrect vector length")
if (length(x) == 0)
return(x)
z <- match.fun(f)(x, y)
return(aplRDV(z, g))
}
# expand vector
aplEXV <- function(x, y) {
z <- rep(0, length(y))
m <- which(y == TRUE)
if (length(m) != length(x))
stop("Incorrect vector length")
z[m] <- x
return(z)
}
# expand
aplExpand <- function(x, y, axis = 1) {
if (is.vector(x))
return(aplEXV(x, y))
d <- dim(x)
m <- which(y == TRUE)
n <- length (y)
e <- d
e[axis] <- n
if (length(m) != d[axis])
stop("Incorrect dimension length")
z <- array(0, e)
for (i in 1:prod(d)) {
k <- aplEncode (i, d)
k[axis] <- m[k[axis]]
z[aplDecode (k, e)] <- x[i]
}
return (z)
}
# compress/replicate vector
aplCRV <- function(x, y) {
n <- aplShape(x)
m <- aplShape(y)
if (m == 1)
y <- rep(y, n)
if (length(y) != n)
stop("Length Error")
z <- vector()
for (i in 1:n)
z <- c(z, rep(x[i], y[i]))
return(z)
}
# compress/replicate
aplReplicate <- function(x, y, k = aplRank (y)) {
if (is.vector(x))
return(aplCRV(x, y))
sx <- aplShape(x)
sy <- aplShape(y)
sk <- sx[k]
if (max(sy) == 1)
y <- rep(y, sk)
if (length(y) != sk)
stop("Length Error")
sz <- sx
sz[k] <- sum(y)
nz <- prod(sz)
gg <- aplCRV(1:sk, y)
z <- array(0, sz)
for (i in 1:nz) {
jvec <- aplEncode(i, sz)
jvec[k] <- gg[jvec[k]]
z[i] <- x[aplDecode(jvec, sx)]
}
return(z)
}
# rotate vector
aplRTV <- function(a, k) {
n <- aplShape(a)
if (k > 0)
return(c(a[-(1:k)], a[1:k]))
if (k < 0)
return(c(a[(n + k + 1):n], a[1:(n + k)]))
return(a)
}
# rotate
aplRotate <- function(a, b, axis = aplRank (a)) {
if (is.vector(a))
return(aplRTV(a, b))
sa <- aplShape(a)
sx <- aplShape(b)
if (max(sx) == 1)
b <- array(b, sa[-axis])
if (!identical(sa[-axis], aplShape(b)))
stop("Index Error")
z <- array(0, sa)
sz <- sa
nz <- prod(sz)
sk <- sz[axis]
for (i in 1:nz) {
ivec <- aplEncode(i, sz)
xx <- b[aplDecode(ivec[-axis], sx)]
ak <- rep(0, sk)
for (j in 1:sk) {
jvec <- ivec
jvec[axis] <- j
ak[j] <- a[aplDecode(jvec, sz)]
}
bk <- aplRTV(ak, xx)
for (j in 1:sk) {
jvec <- ivec
jvec[axis] <- j
z[aplDecode(jvec, sz)] <- bk[j]
}
}
return(z)
}
# transpose -- will be overwritten by the C version
aplTranspose <- function(a, x = rev(1:aplRank(a))) {
sa <- aplShape(a)
ra <- aplRank(a)
if (length(x) != ra)
stop("Length Error")
rz <- max(x)
sz <- rep(0, rz)
for (i in 1:rz)
sz[i] <- min(sa[which(x == i)])
nz <- prod(sz)
z <- array(0, sz)
for (i in 1:nz)
z[i] <- a[aplDecode(aplEncode(i, sz)[x], sa)]
return(z)
}
# representation -- will be overwritten by the C version
aplEncode <- function(rrr, base) {
b <- c(1, butLast(cumprod(base)))
r <- rep(0, length(b))
s <- rrr - 1
for (j in length(base):1) {
r[j] <- s %/% b[j]
s <- s - r[j] * b[j]
}
return(1 + r)
}
# base value -- will be overwritten by the C version
aplDecode <- function(ind, base) {
b <- c(1, butLast(cumprod(base)))
return(1 + sum(b * (ind - 1)))
}
# get
aplGet <- function(a, cell) {
dims <- dim(a)
n <- length(dims)
b <- 0
if (any(cell > dims) || any(cell < 1))
stop("No such cell")
return(a[aplDecode(cell, dims)])
}
# set
aplSet <- function(a, b, cell) {
dims <- dim(a)
n <- length(dims)
if (any(cell > dims) || any(cell < 1))
stop("No such cell")
a[aplDecode(cell, dims)] <- b
return(a)
}
# join
aplJoin <- function(a, b, axis = 1) {
if (is.vector(a) && is.vector(b))
return(c(a, b))
sa <- aplShape(a)
sb <- aplShape(b)
ra <- aplRank(a)
rb <- aplRank(b)
if (ra != rb)
stop("Rank error in aplJoin")
if (!identical(sa[-axis], sb[-axis]))
stop("Shape error")
sz <- sa
sz[axis] <- sz[axis] + sb[axis]
nz <- prod(sz)
u <- unit(axis, ra)
z <- array(0, sz)
for (i in 1:nz) {
ivec <- aplEncode(i, sz)
if (ivec[axis] <= sa[axis])
z[i] <- a[aplDecode(ivec, sa)]
else
z[i] <- b[aplDecode(ivec - sa[axis] * u, sb)]
}
return(z)
}
# ravel
aplRavel <- function(a) {
as.vector(a)
}
# outer product
aplOuterProduct <- function(x, y, f = "*") {
return(outer(x, y, f))
}
# shape
aplShape <- function(a) {
if (is.vector(a))
return(length(a))
return(dim(a))
}
# rank
aplRank <- function(a) {
aplShape(aplShape(a))
}
# reshape
aplReshape <- function(a, d) {
return(array(a, d))
}
# reduce -- will be overwritten by C version
aplReduce <- function(a, k), f = "+") {
if (is.vector(a))
return(aplRDV(a, f))
ff <- if (is.function(f))
f
else
match.fun(f)
sa <- aplShape(a)
ra <- aplRank(a)
na <- prod(sa)
sz <- sa[(1:ra)[-k]]
z <- array(0, sz)
nz <- prod(sz)
ind <- rep(0, nz)
for (i in 1:na) {
ivec <- aplEncode(i, sa)
jind <- aplDecode(ivec[-k], sz)
if (ind[jind] == 0) {
z[jind] <- a[i]
ind[jind] <- 1
}
else
z[jind] <- ff(z[jind], a[i])
}
return(z)
}
# scan -- will be overwritten by the C version
aplScan <- function(a, k = aplRank(a), f = "+") {
if (is.vector(a))
return(aplSCV(a, f))
ff <- if (is.function(f))
f
else
match.fun(f)
sa <- aplShape(a)
ra <- aplRank(a)
sk <- sa[k]
u <- unit(k, ra)
na <- prod(sa)
z <- a
for (i in 1:na) {
ivec <- aplEncode(i, sa)
sk <- ivec[k]
if (sk == 1)
z[i] <- a[i]
else
z[i] <- ff(z[aplDecode(ivec - u, sa)], a[i])
}
return(z)
}
# inner product -- will be overwritten by the C version
aplInnerProduct <- function(a, b, f = "*", g = "+") {
sa <- aplShape(a)
sb <- aplShape(b)
ra <- aplRank(a)
rb <- aplRank(b)
ia <- 1:(ra - 1)
ib <- (ra - 1) + (1:(rb - 1))
ff <- match.fun(f)
gg <- match.fun(g)
ns <- last(sa)
nt <- first(sb)
if (ns != nt)
stop("Incompatible array dimensions")
sz <- c(butLast(sa), butFirst(sb))
nz <- prod(sz)
z <- array(0, sz)
for (i in 1:nz) {
ivec <- aplEncode(i, sz)
for (j in 1:ns) {
aa <- a[aplDecode(c(ivec[ia], j), sa)]
bb <- b[aplDecode(c(j, ivec[ib]), sb)]
tt <- ff(aa, bb)
if (j == 1)
z[i] <- tt
else
z[i] <- gg(z[i], tt)
}
}
return(z)
}
# member of
aplMemberOf <- function(a, b) {
sa <- aplShape(a)
sb <- aplShape(b)
na <- prod(sa)
nb <- prod(sb)
z <- array (0, sa)
for (i in 1:na) {
aa <- a[i]
for (j in 1:nb)
if (aa == b[j])
z[i] <- 1
}
return(z)
}
# utilities below
first <- function(x) {
return(x[1])
}
butFirst <- function(x) {
return(x[-1])
}
last <- function(x) {
return(x[length(x)])
}
butLast <- function(x) {
return(x[-length(x)])
}
unit <- function(i, n) {
ifelse(i == (1:n), 1, 0)
}5.2 R Glue for .C()
aplDecode <-
function(cell, dims) {
n <- length(dims)
if (any(cell > dims) || any(cell < 1))
stop("No such cell")
.C("aplDecodeC ",
as.integer(cell),
as.integer(dims),
as.integer(n),
as.integer(1))[[4]]
}
aplEncode <-
function(ind, dims) {
n <- length(dims)
cell <- integer(n)
if ((ind < 1) || (ind > prod(dims)))
stop("No such cell")
.C("aplEncodeC ",
as.integer(cell),
as.integer(dims),
as.integer(n),
as.integer(ind))[[1]]
}
aplInnerProduct <-
function(a, b, f = "*", g = "+") {
sa <- aplShape(a)
sb <- aplShape(b)
ra <- aplRank(a)
rb <- aplRank(b)
ia <- 1:(ra - 1)
ib <- (ra - 1) + (1:(rb - 1))
ff <- match.fun(f)
gg <- match.fun(g)
ns <- last(sa)
nt <- first(sb)
if (ns != nt)
stop("Incompatible array dimensions")
sz <- c(butLast(sa), butFirst(sb))
nz <- prod(sz)
z <- array(0, sz)
rz <- aplRank(z)
res <-
.C(
"aplInnerProductC",
list(ff, gg),
as.double(a),
as.double(b),
as.integer(sa),
as.integer(ra),
as.integer(sb),
as.integer(rb),
as.integer(sz),
as.integer(rz),
as.integer(nz),
as.integer(ns),
as.double(z)
)
return(array(res[[12]], sz))
}
aplReduce <-
function(a, k, f = "+") {
if (is.vector(a))
return(aplRDV(a, f))
ff <- if (is.function(f))
f
else
match.fun(f)
sa <- aplShape(a)
ra <- aplRank(a)
na <- prod(sa)
sz <- sa[(1:ra)[-k]]
z <- array(0, sz)
rz <- aplRank(z)
nz <- prod(sz)
z <-
.C(
"aplReduceC",
list(ff),
as.double(a),
as.integer(k),
as.integer(length(k)),
as.integer(na),
as.integer(sa),
as.integer(ra),
as.integer(nz),
as.integer(sz),
as.integer(rz),
as.double(z)
)
return(array(z[[11]], sz))
}
aplScan <-
function(a, k, f = "+") {
if (is.vector(a))
return(aplSCV(a, f))
ff <- if (is.function(f))
f
else
match.fun(f)
sa <- aplShape(a)
ra <- aplRank(a)
sk <- sa[k]
u <- unit(k, ra)
na <- prod(sa)
z <- a
res <- .C(
"aplScanC",
list(ff),
as.double(a),
as.integer(k),
as.integer(na),
as.integer(sa),
as.integer(ra),
as.double(z)
)
return(array(res[[7]], sa))
}
aplSelect <-
function(a, x, drop = FALSE) {
sa <- aplShape(a)
ra <- aplRank(a)
sz <- sapply(x, length)
z <- array(0, sz)
rz <- aplRank(z)
nz <- prod(sz)
z <-
array(
.C(
"aplSelectC",
as.double(a),
as.integer(sa),
as.integer(ra),
lapply(x, as.integer),
as.double(z),
as.integer(sz),
as.integer(rz),
as.integer(nz)
)[[5]],
sz
)
if (drop)
return(drop(z))
else
return(z)
}
aplTranspose <-
function(a, x = rev(1:aplRank(a))) {
sa <- aplShape(a)
ra <- aplRank(a)
na <- prod(sa)
if (length(x) != ra)
stop("Length Error")
rz <- max(x)
sz <- rep(0, rz)
for (i in 1:rz)
sz[i] <- min(sa[which(x == i)])
nz <- prod(sz)
z <- array(0, sz)
array(
.C(
"aplTransposeC",
as.double(a),
as.integer(x),
as.integer(sa),
as.integer(ra),
as.integer(na),
as.integer(sz),
as.integer(rz),
as.integer(nz),
as.double(z)
)[[9]],
sz
)
}5.3 R Glue for .Call()
aplDecode <- function(cell, dims) {
if (length(cell) != length(dims)) {
stop("Dimension error")
}
if (any(cell > dims) || any (cell < 1)) {
stop("No such cell")
}
.Call("APLDECODE", as.integer(cell), as.integer(dims))
}
aplEncode <- function(ind, dims) {
if (length(ind) > 1) {
stop ("Dimension error")
}
if ((ind < 1) || (ind > prod(dims))) {
stop ("No such cell")
}
.Call("APLENCODE", as.integer(ind), as.integer(dims))
}
aplInnerProduct <- function(a, b, f = "*", g = "+") {
sa <- aplShape(a)
ra <- aplRank(a)
ia <- 1:(ra - 1)
sb <-
aplShape(b)
rb <- aplRank(b)
ib <- (ra - 1) + (1:(rb - 1))
if (!is.function(f)) {
f <- match.fun(f)
}
if (!is.function(g)) {
g <- match.fun(g)
}
env <- new.env()
environment(f) <- env
environment(g) <- env
ns <- last(sa)
nt <- first(sb)
if (ns != nt) {
stop("Incompatible array dimensions")
}
sz <- c(butLast(sa), butFirst(sb))
z <- .Call(
"APLINNERPRODUCT",
f,
g,
as.double(a),
as.double(b),
as.integer(sa),
as.integer(sb),
as.integer(sz),
as.integer(ns),
env
)
return(array(z, sz))
}
aplReduce <- function(a, k, f = "+") {
if (is.vector(a)) {
return(aplRDV(a, f))
}
nk <- length(k)
sa <- aplShape(a)
sz <- sa[(1:length(sa))[-k]]
if (!is.function(f)) {
f <- match.fun(f)
}
env <- new.env()
environment(f) <- env
z <- .Call("APLREDUCE",
f,
as.double(a),
as.integer(k),
as.integer(sa),
as.integer(sz),
env)
return(array(z, sz))
}
aplScan <- function(a, k = aplRank (a), f = "+") {
if (is.vector(a)) {
return(aplSCV(a, f))
}
if (!is.function(f)) {
f <- match.fun(f)
}
env <- new.env()
environment(f) <- env
sa <- aplShape(a)
ra <- aplRank(a)
z <- .Call("APLSCAN",
f,
as.double(a),
as.integer(k),
as.integer(sa),
as.integer(ra),
env)
return(array(z, sa))
}
aplSelect <- function (a, x, drop = TRUE) {
dima = aplShape(a)
if (length(dima) != length(x)) {
stop("Dimension error")
}
z <- .Call("APLSELECT",
as.double(a),
as.integer(dima),
lapply(x, as.integer))
z <- array(z, sapply(x, length))
if (drop) {
return(drop(z))
}
return(z)
}
aplTranspose <- function(a, x = rev(1:aplRank(a))) {
sa <- aplShape(a)
if (length(x) != length(sa)) {
stop("Length Error")
}
rz <- max(x)
sz <- rep(0, rz)
for (i in 1:rz) {
sz[i] <- min(sa[which(x == i)])
}
z <- .Call(
"APLTRANSPOSE",
as.double(a),
as.integer(x),
as.integer(sa),
as.integer(sz),
as.integer(rz)
)
return(array(z, sz))
}5.4 C Code using .C()
#include <R.h>
#include <Rinternals.h>
#include <Rinterface.h>
static char* f2_char;
static char* g2_char;
void aplDecodeC(int* cell, int* dims, int* n, int* ind) {
int i, aux = 1; *ind = 1;
for (i = 0; i < *n; i++) {
*ind += aux * (cell[i] - 1);
aux *= dims[i];
}
}
void aplEncodeC(int* cell, int* dims, int* n, int* ind) {
int i, aux = *ind, nn = *n, pdim = 1;
for (i = 0; i < nn - 1; i++)
pdim *= dims[i];
for (i = nn - 1; i > 0; i--) {
cell[i] = (aux - 1)/pdim;
aux -= pdim*cell[i];
pdim /= dims[i-1];
cell[i] += 1;
}
cell[0] = aux;
}
void aplSelectC(double *a, int *sa, int *ra, SEXP *list, double *z, int *sz, int *rz, int *nz)
{
int i, j, k;
int ivec[*rz], jvec[*ra];
for (i = 0; i < *nz; i++) {
k = i + 1;
(void) aplEncodeC (ivec,sz,rz,&k);
for (j = 0; j < *ra; j++) {
SEXP vec = list[j];
jvec[j] = INTEGER(vec)[ivec[j]-1];
}
k = 1;
(void) aplDecodeC (jvec,sa,ra,&k);
z[i] = a[k - 1];
}
}
void aplTransposeC(double *a, int *x, int *sa, int *ra, int *na, int *sz, int *rz, int *nz, double *z)
{
int i, j, r, ivec[*rz], jvec[*ra];
for (i = 0; i < *nz; i++){
r = i + 1;
(void) aplEncodeC(ivec,sz,rz,&r);
for (j = 0; j < *ra; j++)
jvec[j] = ivec[x[j]-1];
(void) aplDecodeC(jvec,sa,ra,&r);
z[i] = a[r - 1];
}
}
void aplReduceC(char** funclist, double *a, int *k, int *nk, int *na, int *sa, int *ra, int *nz, int *sz, int *rz, double *z) {
double f2_glue(double, double);
double f2_comp(double (*)(),double,double);
int i, j, u, v, r, m, kk = (*ra) - (*nk), ivec[*ra], kvec[kk], ind[*nz];
for (i = 0; i < *nz; i++) ind[i] = 0;
f2_char = funclist[0];
for (i = 0; i < *na; i++){
r = i + 1;
(void) aplEncodeC(ivec,sa,ra,&r);
u = 0;
for (j = 0; j < *ra; j++) {
r = 0;
for (v = 0; v < *nk; v++) {
if (j == (k[v] - 1)) r = 1;
}
if (r == 0) {
kvec[u] = ivec[j];
u += 1;
}
}
(void) aplDecodeC(kvec,sz,rz,&m);
if (ind[m - 1] == 0) {
z[m - 1] = a[i];
ind[m - 1] = 1;
}
else
z[m - 1] = f2_comp((double(*)())f2_glue,z[m - 1],a[i]);
}
}
void aplScanC(char** funclist,double *a, int *k, int *na, int *sa, int *ra, double *z)
{
double f2_glue(double, double);
double f2_comp(double (*)(),double,double);
int i, r, sk, l, ivec[*ra];
l = *k - 1;
f2_char = funclist[0];
for (i = 0; i < *na; i++){
r = i + 1;
(void) aplEncodeC(ivec,sa,ra,&r);
sk = ivec[l];
if (sk == 1) z[i] = a[i];
else {
ivec[l] -= 1;
(void) aplDecodeC(ivec,sa,ra,&r);
z[i] = f2_comp((double(*)())f2_glue,z[r - 1],a[i]);
}
}
}
void aplInnerProductC(char** funclist, double *a, double *b, int *sa, int *ra, int *sb, int *rb, int *sz, int *rz, int *nz, int *ns, double *z) {
double f2_glue(double, double);
double g2_glue(double, double);
double f2_comp(double (*)(),double,double);
double g2_comp(double (*)(),double,double);
double t; int i, j, r, k, l, u, ivec[*rz], jvec[*ra], kvec[*rb];
f2_char = funclist[0]; g2_char= funclist[1]; k = l = 0;
for (i = 0; i < *nz; i++) {
r = i + 1;
(void) aplEncodeC(ivec,sz,rz,&r);
for (j = 0; j < *ns; j++) {
for (u = 0; u < *ra - 1; u++)
jvec[u] = ivec[u];
jvec[*ra - 1] = j + 1;
(void) aplDecodeC(jvec,sa,ra,&k);
for (u = 1; u < *rb; u++)
kvec[u] = ivec[*ra + u - 2];
kvec[0] = j + 1;
(void) aplDecodeC(kvec,sb,rb,&l);
t = f2_comp((double(*)())f2_glue,a[k-1],b[l-1]);
if (j == 0) z[i] = t;
else z[i] = g2_comp((double(*)())g2_glue,t,z[i]);
}
}
}
double f2_glue(double x, double y){
char *modes[2], *values[1];
void *args[2];
double xx[1], yy[1], *result;
long lengths[2], nargs = (long)2, nvals = (long)1;
lengths[0] = lengths[1] = (long)1;
nargs = (long)2;
args[0] = (void *)xx; xx[0] = x;
args[1] = (void *)yy; yy[0] = y;
modes[0] = modes[1] = "double";
call_R(f2_char,nargs,args,modes,lengths,0,nvals,values);
result = (double*)values[0];
return(result[0]);
}
double g2_glue(double x, double y){
char *modes[2], *values[1];
void *args[2];
double xx[1], yy[1], *result;
long lengths[2], nargs = (long)2, nvals = (long)1;
lengths[0] = lengths[1] = (long)1;
nargs = (long)2;
args[0] = (void *)xx; xx[0] = x;
args[1] = (void *)yy; yy[0] = y;
modes[0] = modes[1] = "double";
call_R(g2_char,nargs,args,modes,lengths,0,nvals,values);
result = (double*)values[0];
return(result[0]);
}
double f2_comp(double (*f)(), double x, double y) {
return f(x,y);
}
double g2_comp(double (*g)(), double x, double y) {
return g(x,y);
}5.5 C Code using .Call()
#include <R.h>
#include <Rinternals.h>
SEXP APLDECODE( SEXP, SEXP );
SEXP APLENCODE( SEXP, SEXP );
SEXP APLSELECT( SEXP, SEXP, SEXP );
SEXP APLTRANSPOSE( SEXP, SEXP, SEXP, SEXP, SEXP );
SEXP APLSCAN( SEXP, SEXP, SEXP, SEXP, SEXP );
SEXP APLREDUCE( SEXP, SEXP, SEXP, SEXP, SEXP, SEXP );
SEXP APLINNERPRODUCT( SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP );
SEXP
APLDECODE( SEXP cell, SEXP dims )
{
int aux = 1, n = length(dims);
SEXP ind;
PROTECT( ind = allocVector( INTSXP, 1 ) );
INTEGER( ind )[0] = 1;
for( int i = 0; i < n; i++ ) {
INTEGER( ind )[0] += aux * ( INTEGER( cell )[i] - 1 );
aux *= INTEGER(dims)[i];
}
UNPROTECT( 1 );
return (ind);
}
SEXP
APLENCODE( SEXP ind, SEXP dims )
{
int n = length(dims), aux = INTEGER(ind)[0], pdim = 1;
SEXP cell;
PROTECT( cell = allocVector( INTSXP, n ) );
for ( int i = 0; i < n - 1; i++ )
pdim *= INTEGER( dims )[i];
for ( int i = n - 1; i > 0; i-- ){
INTEGER( cell )[i] = ( aux - 1 ) / pdim;
aux -= pdim * INTEGER( cell )[i];
pdim /= INTEGER( dims )[i - 1];
INTEGER( cell )[i] += 1;
}
INTEGER( cell )[0] = aux;
UNPROTECT( 1 );
return cell;
}
SEXP
APLSELECT( SEXP a, SEXP dima, SEXP list ) {
int r = length (dima), lz = 1, dimzi, nProtect = 0;
SEXP dimz, itel, cell, czll, nind, z;
PROTECT(dimz = allocVector (INTSXP, r)); nProtect++;
PROTECT(cell = allocVector (INTSXP, r)); nProtect++;
PROTECT(czll = allocVector (INTSXP, r)); nProtect++;
PROTECT(itel = allocVector (INTSXP, 1)); nProtect++;
PROTECT(nind = allocVector (INTSXP, 1)); nProtect++;
for( int i = 0; i < r; i++ ){
dimzi = length( VECTOR_ELT( list, i ) );
INTEGER( dimz )[i] = dimzi;
lz *= dimzi;
}
PROTECT( z = allocVector( REALSXP, lz ) );
nProtect++;
for( int i = 0; i < lz; i++ ){
INTEGER(itel)[0] = i + 1;
cell = APLENCODE (itel, dimz);
for (int j = 0; j < r; j++) {
INTEGER (czll)[j] = INTEGER (VECTOR_ELT (list, j))[INTEGER(cell)[j] - 1];
}
nind = APLDECODE( czll, dima );
REAL(z)[i] = REAL(a)[INTEGER(nind)[0] - 1];
}
UNPROTECT(nProtect);
return(z);
}
SEXP
APLTRANSPOSE(SEXP a, SEXP x, SEXP sa, SEXP sz, SEXP rz)
{
int i, j, na=1, nz=1, ra = length (sa), lsz = length( sz ), nProtected=0;
SEXP ivec, jvec, z, itel, nind;
for( i=0;i<ra ;i++){ na *= INTEGER( sa )[i]; }
for( i=0;i<lsz;i++){ nz *= INTEGER( sz )[i]; }
PROTECT(itel = allocVector( INTSXP, 1 ) ); ++nProtected;
PROTECT(nind = allocVector( INTSXP, 1 ) ); ++nProtected;
PROTECT(ivec = allocVector( INTSXP, INTEGER(rz)[0] ) ); ++nProtected;
PROTECT(jvec = allocVector( INTSXP, ra ) ); ++nProtected;
PROTECT(z = allocVector( REALSXP, nz ) ); ++nProtected;
for( i = 0; i < nz; i++ ){
INTEGER( itel )[0] = i + 1;
ivec = APLENCODE( itel, sz );
for( j = 0; j < ra; j++ ){
INTEGER( jvec )[j] = INTEGER( ivec )[INTEGER( x )[j] - 1];
}
nind = APLDECODE( jvec, sa );
REAL( z )[i] = REAL( a )[INTEGER(nind)[0] - 1];
}
UNPROTECT( nProtected );
return z;
}
SEXP
APLREDUCE(SEXP f, SEXP a, SEXP k, SEXP sa, SEXP sz, SEXP env)
{
int i, j, u, v, r, kk, na=1, nz=1, nProtected = 0;
int nk = length( k ), ra = length( sa ), rz = length( sz );
SEXP ivec, kvec, ind, z,itel, nind, R_fcall= R_NilValue;
SEXP Z = R_NilValue, A = R_NilValue;
for( i=0;i<ra;i++){ na *= INTEGER( sa )[i]; }
for( i=0;i<rz;i++){ nz *= INTEGER( sz )[i]; }
kk = ra-nk;
PROTECT( R_fcall= lang3(f, R_NilValue, R_NilValue) ); ++nProtected;
PROTECT( Z = allocVector( REALSXP, 1 ) ); ++nProtected;
PROTECT( A = allocVector( REALSXP, 1 ) ); ++nProtected;
PROTECT( itel = allocVector( INTSXP, 1 ) ); ++nProtected;
PROTECT( ivec = allocVector( INTSXP, ra ) ); ++nProtected;
PROTECT( kvec = allocVector( INTSXP, kk ) ); ++nProtected;
PROTECT( ind = allocVector( INTSXP, nz ) ); ++nProtected;
PROTECT( z = allocVector( REALSXP, nz ) ); ++nProtected;
for( i = 0; i < nz; i++ ){
INTEGER( ind )[i] = 0;
}
for( i = 0; i < na; i++ ){
INTEGER( itel )[0] = i + 1;
ivec = APLENCODE( itel, sa );
u = 0;
for( j = 0; j < ra; j++ ){
r = 0;
for( v = 0; v < nk; v++ ){
if( j == ( INTEGER( k )[v] - 1 ) ) r = 1;
}
if( r == 0 ){
INTEGER( kvec )[u] = INTEGER( ivec )[j];
u += 1;
}
}
nind = APLDECODE( kvec, sz );
if ( INTEGER( ind )[INTEGER( nind )[0] - 1] == 0 ) {
REAL( z )[INTEGER( nind )[0] - 1] = REAL( a )[i];
INTEGER( ind )[INTEGER( nind )[0] - 1] = 1;
} else {
REAL( Z )[0] = REAL( z )[INTEGER( nind )[0] - 1];
REAL( A )[0] = REAL( a )[i];
SETCADR( R_fcall, Z );
SETCADDR( R_fcall, A );
REAL( z )[INTEGER( nind )[0] - 1] = REAL( eval( R_fcall, env ) )[0];
}
}
UNPROTECT( nProtected );
return z;
}
SEXP
APLSCAN( SEXP f, SEXP a, SEXP k, SEXP sa, SEXP ra, SEXP env )
{
int i, sk, l, na=1, ka = INTEGER( ra )[0],nProtected=0;
for( i=0;i<ka ;i++ ){ na *= INTEGER( sa )[i]; }
SEXP ivec, z, itel, nind, Z, A,R_fcall=R_NilValue;
PROTECT( R_fcall = lang3( f, R_NilValue, R_NilValue ) ); ++nProtected;
PROTECT( itel = allocVector( INTSXP, 1 ) ); ++nProtected;
PROTECT( nind = allocVector( INTSXP, 1 ) ); ++nProtected;
PROTECT( Z = allocVector( REALSXP, 1 ) ); ++nProtected;
PROTECT( A = allocVector( REALSXP, 1 ) ); ++nProtected;
PROTECT( ivec = allocVector( INTSXP, ka ) ); ++nProtected;
PROTECT( z = allocVector( REALSXP, na ) ); ++nProtected;
l = INTEGER( k )[0] - 1;
for( i = 0; i < na; i++ ){
INTEGER( itel )[0] = i + 1;
ivec = APLENCODE( itel, sa );
sk = INTEGER( ivec )[l];
if( sk == 1 ){
REAL( z )[i] = REAL( a )[i];
}else{
INTEGER( ivec )[l] -= 1;
nind = APLDECODE( ivec, sa );
REAL( Z )[0]=REAL( z )[INTEGER( nind )[0]-1];
REAL( A )[0]=REAL( a )[i];
SETCADR( R_fcall, Z );
SETCADDR( R_fcall, A );
REAL( z )[i] = REAL( eval( R_fcall, env ) )[0];
}
}
UNPROTECT( nProtected );
return z;
}
SEXP
APLINNERPRODUCT(SEXP f, SEXP g, SEXP a, SEXP b, SEXP sa, SEXP sb, SEXP sz, SEXP ns, SEXP env)
{
int i, j, u, nz = 1, nProtected = 0;
int ra = length( sa ),rb = length( sb ), rz = length( sz );
SEXP ivec, jvec, kvec, z, A, B, Z, itel, k, l, t;
SEXP R_fcall = R_NilValue, R_gcall = R_NilValue;
for( i = 0; i < rz; i++ ){ nz *= INTEGER( sz )[i]; }
PROTECT(R_fcall = lang3( f, R_NilValue, R_NilValue ) ); ++nProtected;
PROTECT(R_gcall = lang3( g, R_NilValue, R_NilValue ) ); ++nProtected;
PROTECT( itel = allocVector( INTSXP, 1 ) ); ++nProtected;
PROTECT( k = allocVector( INTSXP, 1 ) ); ++nProtected;
PROTECT( l = allocVector( INTSXP, 1 ) ); ++nProtected;
PROTECT( ivec = allocVector( INTSXP, rz ) ); ++nProtected;
PROTECT( jvec = allocVector( INTSXP, ra ) ); ++nProtected;
PROTECT( kvec = allocVector( INTSXP, rb ) ); ++nProtected;
PROTECT( z = allocVector( REALSXP, nz ) ); ++nProtected;
PROTECT( t = allocVector( REALSXP, 1 ) ); ++nProtected;
PROTECT( Z = allocVector( REALSXP, 1 ) ); ++nProtected;
PROTECT( B = allocVector( REALSXP, 1 ) ); ++nProtected;
PROTECT( A = allocVector( REALSXP, 1 ) ); ++nProtected;
for( i = 0; i < nz; i++ ){
INTEGER( itel )[0] = i + 1;
ivec = APLENCODE( itel, sz );
for( j = 0; j < INTEGER( ns )[0]; j++ ){
for( u = 0; u < ra - 1; u++ ){
INTEGER( jvec )[u] = INTEGER( ivec )[u];
}
INTEGER( jvec )[ra - 1] = j + 1;
k = APLDECODE( jvec, sa );
for( u = 1; u < rb; u++ ){
INTEGER( kvec )[u] = INTEGER( ivec )[ra + u - 2];
}
INTEGER( kvec )[0] = j + 1;
l = APLDECODE( kvec,sb );
REAL( A )[0] = REAL( a )[INTEGER( k )[0]-1];
REAL( B )[0] = REAL( b )[INTEGER( l )[0]-1];
SETCADR( R_fcall, A );
SETCADDR( R_fcall, B );
REAL( t )[0] = REAL( eval( R_fcall, env ) )[0];
if( j == 0 ){
REAL( z )[i] = REAL( t )[0];
} else {
REAL( Z )[0] = REAL( z )[i];
SETCADR( R_gcall, t );
SETCADDR( R_gcall, Z );
REAL( z )[i] = REAL( eval( R_fcall, env ) )[0];
}
}
}
UNPROTECT( nProtected );
return z;
}5.6 Rcpp code
5.6.1 hpp
#ifndef _apl_RCPP_H
#define _apl_RCPP_H
/* The Rcpp header takes care of everything you should need. */
#include <Rcpp.h>
/*
* note : RcppExport is an alias to `extern "C"` defined by Rcpp.
*
* It gives C calling convention to the rcpp_hello_world function so that
* it can be called from .Call in R. Otherwise, the C++ compiler mangles the
* name of the function and .Call can't find it.
*
* It is only useful to use RcppExport when the function is intended to be called
* by .Call. See the thread http://thread.gmane.org/gmane.comp.lang.r.rcpp/649/focus=672
* on Rcpp-devel for a misuse of RcppExport
*/
//RcppExport SEXP rcpp_hello_world() ;
RcppExport SEXP APLDECODErcpp( SEXP, SEXP );
RcppExport SEXP APLENCODErcpp( SEXP, SEXP );
RcppExport SEXP APLSELECTrcpp( SEXP, SEXP, SEXP );
RcppExport SEXP APLTRANSPOSErcpp( SEXP, SEXP, SEXP, SEXP, SEXP );
RcppExport SEXP APLSCANrcpp( SEXP, SEXP, SEXP, SEXP, SEXP );
RcppExport SEXP APLREDUCErcpp( SEXP, SEXP, SEXP, SEXP, SEXP, SEXP );
RcppExport SEXP APLINNERPRODUCTrcpp( SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP );
#endif5.6.2 cpp
//
#include "aplRcpp.hpp"
using namespace Rcpp ;
SEXP
APLDECODErcpp( SEXP cell_r, SEXP dims_r )
{
IntegerVector cell( cell_r );
IntegerVector dims( dims_r );
int n = dims.size(), aux = 1;
IntegerVector ind( 1 ); ind[0] = 1;
for( int i = 0; i < n; i++ ) {
ind[0] += aux * ( cell[i] - 1 );
aux *= dims[i];
}
return wrap( ind );
}
SEXP
APLENCODErcpp( SEXP ind_r, SEXP dims_r )
{
IntegerVector dims( dims_r );
IntegerVector ind( ind_r );
int n = dims.size(), aux = ind[0], pdim = 1;
IntegerVector cell( n );
for( int i = 0; i < n - 1; i++ ){
pdim *= dims[i];
}
for( int i = n - 1; i > 0; i-- ){
cell[i] = ( aux - 1 ) / pdim;
aux -= pdim * cell[i];
pdim /= dims[i - 1];
cell[i] += 1;
}
cell[0] = aux;
return wrap( cell );
}
SEXP
APLSELECTrcpp( SEXP a_r, SEXP dima_r, SEXP list_r )
{
IntegerVector dima( dima_r );
NumericVector a( a_r );
int r = dima.size(), lz = 1, dimzi;
List list( list_r );
IntegerVector dimz( r );
IntegerVector cell( r );
IntegerVector czll( r );
IntegerVector itel( 1 );
IntegerVector nind( 1 );
for( int i = 0; i < r; i++ ){
dimzi = as<IntegerVector>( list[i] ).size();
dimz[i] = dimzi;
lz *= dimzi;
}
NumericVector z(lz);
for( int i = 0; i < lz; i++ ){
itel[0] = i + 1;
cell = APLENCODErcpp( wrap( itel ), wrap( dimz ) );
for ( int j = 0; j < r; j++ ) {
IntegerVector listj= as<IntegerVector>( list[j] );
czll[j] = listj[cell[j] - 1];
}
nind = APLDECODErcpp( wrap( czll ), wrap( dima ) );
z[i] = a[nind[0] - 1];
}
return wrap( z );
}
SEXP
APLTRANSPOSErcpp( SEXP a_r, SEXP x_r, SEXP sa_r, SEXP sz_r, SEXP rz_r )
{
IntegerVector sa( sa_r ), sz( sz_r );
int na = 1, nz = 1, ra = sa.size(), lsz = sz.size();
IntegerVector rz( rz_r ), x( x_r ), a( a_r );
for( int i = 0; i < ra - 1; i++ ){ na *= sa[i]; }
for( int i = 0; i < lsz - 1; i++ ){ nz *= sz[i]; }
IntegerVector itel( 1 ), nind( 1 ), ivec( rz[0] ), jvec( ra );
NumericVector z( nz );
for( int i = 0; i < nz; i++ ){
itel[0] = i + 1;
ivec = APLENCODErcpp( wrap( itel ), wrap( sz ) );
for( int j = 0; j < ra; j++ ){
jvec[j] = ivec[x[j] - 1];
}
nind = APLDECODErcpp( wrap( jvec ), wrap( sa ) );
z[i] = a[nind[0] - 1];
}
return wrap( z );
}
SEXP
APLREDUCErcpp( SEXP f_r, SEXP a_r, SEXP k_r, SEXP sa_r, SEXP sz_r, SEXP env_r )
{
int u, r, kk, na = 1, nz = 1, nProtected=0;
IntegerVector k( k_r ), sa( sa_r ), sz( sz_r );
NumericVector a( a_r );
int nk = k.size(), ra = sa.size(), rz = sz.size();
SEXP R_fcall= R_NilValue;
for( int i = 0; i < ra; i++ ){ na *= sa[i]; }
for( int i = 0; i < rz; i++ ){ nz *= sz[i]; }
kk = ra - nk;
PROTECT( R_fcall= Rf_lang3(f_r, R_NilValue, R_NilValue) ); ++nProtected;
IntegerVector itel( 1 ), nind( 1 ), ind( nz ), ivec( ra ), kvec( kk );
NumericVector z( nz );
for( int i = 0; i < na; i++ ){
itel[0] = i + 1;
ivec = APLENCODErcpp( wrap( itel ), wrap( sa ) );
u = 0;
for( int j = 0; j < ra; j++ ){
r = 0;
for( int v = 0; v < nk; v++ ){
if( j == ( k[v] - 1 ) ) r = 1;
}
if( r == 0 ){
kvec[u] = ivec[j];
u += 1;
}
}
nind = APLDECODErcpp( wrap( kvec ), wrap( sz ) );
if ( ind[nind[0] - 1] == 0 ) {
z[nind[0] - 1] = a[i];
ind[nind[0] - 1] = 1;
} else {
SETCADR( R_fcall, wrap( z[nind[0] - 1] ) );
SETCADDR( R_fcall, wrap( a[i] ) );
z[nind[0] - 1] = REAL( Rf_eval( R_fcall, env_r ) )[0];
}
}
UNPROTECT( nProtected );
return wrap( z );
}
SEXP
APLSCANrcpp( SEXP f_r, SEXP a_r, SEXP k_r, SEXP sa_r, SEXP env_r )
{
IntegerVector k( k_r ) , sa( sa_r );
NumericVector a( a_r );
int sk, l, na=1, ra = sa.size(), nProtected=0;
for( int i = 0; i < ra; i++ ){ na *= sa[i]; }
SEXP R_fcall=R_NilValue;
PROTECT( R_fcall = Rf_lang3( f_r, R_NilValue, R_NilValue ) ); ++nProtected;
IntegerVector itel( 1 ), nind( 1 ), ivec( ra );
NumericVector z( na );
l = k[0] - 1;
for( int i = 0; i < na; i++ ){
itel[0] = i + 1;
ivec = APLENCODErcpp( itel, sa );
sk = ivec[l];
if( sk == 1 ){
z[i] = a[i];
} else {
ivec[l] -= 1;
nind = APLDECODErcpp( wrap( ivec ), wrap( sa ) );
SETCADR( R_fcall, wrap( z[nind[0]-1] ) );
SETCADDR( R_fcall, wrap( a[i] ) );
z[i] = REAL( Rf_eval( R_fcall, env_r ) )[0];
}
}
UNPROTECT( nProtected );
return wrap( z );
}
SEXP
APLINNERPRODUCTrcpp(SEXP f_r, SEXP g_r, SEXP a_r, SEXP b_r, SEXP sa_r, SEXP sb_r, SEXP sz_r, SEXP ns_r, SEXP env_r)
{
int nz = 1, nProtected = 0;
IntegerVector sa( sa_r ) , sb( sb_r ) , sz( sz_r ), ns( ns_r );
NumericVector a( a_r ), b( b_r );
int ra = sa.size(),rb = sb.size(), rz = sz.size();
SEXP R_fcall = R_NilValue, R_gcall = R_NilValue;
for( int i = 0; i < rz; i++ ){ nz *= sz[i]; }
PROTECT( R_fcall = Rf_lang3( f_r, R_NilValue, R_NilValue ) ); ++nProtected;
PROTECT( R_gcall = Rf_lang3( g_r, R_NilValue, R_NilValue ) ); ++nProtected;
IntegerVector itel( 1 ), k( 1 ), l( 1 ), ivec( rz ), jvec( ra ), kvec( rb );
NumericVector z( nz ), t( 1 ), Z( 1 ), B( 1 ), A( 1 );
for( int i = 0; i < nz; i++ ){
itel[0] = i + 1;
ivec = APLENCODErcpp( itel, sz );
for( int j = 0; j < ns[0]; j++ ){
for( int u = 0; u < ra - 1; u++ ){
jvec[u] = ivec[u];
}
jvec[ra - 1] = j + 1;
k = APLDECODErcpp( jvec, sa );
for( int u = 1; u < rb; u++ ){
kvec[u] = ivec[ra + u - 2];
}
kvec[0] = j + 1;
l = APLDECODErcpp( kvec,sb );
SETCADR( R_fcall, wrap( a[k[0]-1]) );
SETCADDR( R_fcall, wrap( b[l[0]-1] ));
t[0] = REAL( Rf_eval( R_fcall, env_r ) )[0];
if( j == 0 ){
z[i] = t[0];
} else {
SETCADR( R_gcall, wrap( t[0] ) );
SETCADDR( R_gcall, wrap( z[i] ) );
z[i] = REAL( Rf_eval( R_gcall, env_r ) )[0];
}
}
}
UNPROTECT( nProtected );
return wrap( z );
}