context("Apply a function to a dimension of an array")
# list of acts: (name: act for arrApply; value: if NULL, the same  act is for apply(),
#   if list, its fields are: rf -> name for equivalent r function, args -> list
#   of ... in apply(), argc -> ... for arrApply())
set.seed(7)
n=3
v=rnorm(n)
m=matrix(rnorm(n*(n+1L)), n, n+1L)
d3=seq(n, n+2L)
ar3d=array(rnorm(prod(d3)), dim=d3)
d4=seq(n, n+3L)
ar4d=array(rnorm(prod(d4)), dim=d4)
vp=lapply(d4, rnorm)
p=sort(runif(3L))
lacts=list("sum"=NULL, "prod"=NULL, "all"=NULL, "any"=NULL, "min"=NULL, "max"=NULL, "mean"=NULL, "median"=NULL, "sd"=NULL, "var"=NULL, "cumsum"=NULL, "cumprod"=NULL, "diff"=NULL,
   # translated acts
   norm=list(rf="norm", argr=list(type='2'), argc=list(p=2)),
   trapz=list(rf=function(v) {n=length(v); return(sum(v)-0.5*(v[1]+v[n]))}),
   normalise=list(rf=function(v) v/norm(v, '2'), argc=list(p=2)),
   multv=list(rf=function(v, vv) v*vv, argr=quote(list(vv=vp[[idim]])), argc=quote(list(v=vp[[idim]]))), 
   divv=list(rf=function(v, vv) v/vv, argr=quote(list(vv=vp[[idim]])), argc=quote(list(v=vp[[idim]]))), 
   addv=list(rf=function(v, vv) v+vv, argr=quote(list(vv=vp[[idim]])), argc=quote(list(v=vp[[idim]]))), 
   subv=list(rf=function(v, vv) v-vv, argr=quote(list(vv=vp[[idim]])), argc=quote(list(v=vp[[idim]]))),
   conv=list(rf=function(v, vv) convolve(v, rev(vv), type="open"), argr=quote(list(vv=vp[[idim]])), argc=quote(list(v=vp[[idim]]))),
   quantile=list(rf=function(v, vv) quantile(v, vv, type=5), argr=list(vv=p), argc=list(p=p))
)
test_ar=function(ar, tol=1.e-14, acts=lacts, ndi=seq_along(dim(ar))) {
    # compare arrApply() to translated r call
    vec=FALSE
    if (length(ndi) == 0) {
        # we have a vector
        ndi=1
        vec=TRUE
    }
    for (act in names(acts)) {
        ract=acts[[act]]
        rfu=if (is.null(ract)) act else ract$rf
        for (idim in ndi) {
            argc=if (is.language(ract$argc)) eval(ract$argc) else ract$argc
            r1=do.call(arrApply, c(list(ar, idim, act), argc))
            if (!vec && length(dim(r1)) == length(dim(ar))) {
                # permute to the same order as in apply
                r1=aperm(r1, c(idim, ndi[-idim]))
            }
            argr=if (is.language(ract$argr)) eval(ract$argr) else ract$argr
            if (vec) {
                r2=suppressWarnings(do.call(apply, c(list(as.matrix(ar), 2, rfu), argr)))
            } else {
                r2=suppressWarnings(do.call(apply, c(list(ar, ndi[-idim], rfu), argr)))
            }
            expect_equal(as.numeric(r1), as.numeric(r2), tolerance=tol, scale=1, info=sprintf("'%s' on idim=%d in dims=(%s)", act, idim, paste(if (vec) length(ar) else dim(ar), collapse=", ")))
        }
    }
}
test_that("arrApply on a vector", {
    test_ar(v)
})
test_that("arrApply on a matrix", {
    test_ar(m)
})
test_that("arrApply on an array 3D", {
    test_ar(ar3d)
})
test_that("arrApply on an array 4D", {
    test_ar(ar4d)
})
test_that("dim preserving",
    expect_equal(dim(arrApply(ar3d, 3L, "sum")), dim(ar3d)[-3L])
)
test_that("dimnames preserving", {
    d3nm=lapply(seq_along(d3), function(i) paste0(letters[seq_len(d3[i])], "_", i))
    dimnames(ar3d)=d3nm
    expect_equal(dimnames(arrApply(ar3d, 3L, "sum")), d3nm[-3L])
    expect_equal(dimnames(arrApply(ar3d, 3L, "addv", v=double(d3[3L]))), d3nm)
    # one of dimnames is NULL
    d3nm[2L]=list(NULL)
    dimnames(ar3d)=d3nm
    expect_equal(dimnames(arrApply(ar3d, 3L, "addv", v=double(d3[3L]))), d3nm)
})