apl.R

language: R
license: GPL 2

Code for Snippet:

                
#
#   apl package
#
#   Copyright (C) 2008  Jan de Leeuw <[email protected]>
#   UCLA Department of Statistics, Box 951554, Los Angeles, CA 90095-1554
#   
#   This program is free software; you can redistribute it and/or modify
#   it under the terms of the GNU General Public License as published by
#   the Free Software Foundation; either version 2 of the License, or
#   (at your option) any later version.
#
#   This program is distributed in the hope that it will be useful,
#   but WITHOUT ANY WARRANTY; without even the implied warranty of
#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#   GNU General Public License for more details.
#   
#   You should have received a copy of the GNU General Public License
#   along with this program; if not, write to the Free Software
#   Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
#
###################################################################
#
# version 0.0.1, 2008-08-16,        first incomplete release
# version 0.1.0, 2008-08-18,        switched to index computation
# version 0.2.0, 2008-08-19,        added take, drop, compress, scan
# version 0.3.0, 2008-08-19,        added transpose and rotate
# version 0.4.0, 2008-11-22,        decode and encode using C
# version 0.4.1, 2008-11-23,        added aplGet and aplSet
# version 0.4.2, 2008-11-26,        aplSelect using C
# version 0.4.3, 2008-11-28,        aplReduce using C
# version 0.4.4, 2008-11-28,        aplTranspose using C
# version 0.4.5, 2008-12-07,        aplScan using C
# version 0.5.0, 2008-12-08,        aplInnerProduct using C
#
# To do:
# -- Remainder to C
# -- Generalize inner product
# -- Do not allow vectors
#
 
if (!exists("USE_C")) USE_C<-TRUE
 
if (USE_C) dyn.load("apl.so")
 
# just to be precise the shape of an array a is
# dim(a), the shape of a vector is length(a),
# the rank of a is the length of its shape.
 
# a is an m-dimensional, x is a list of length
# m. Each x[[i]] is a non-empty subvector 
# of 1:dim(a)[i]. aplSelect() returns
# the slice a[x[[1]],...,x[[m]]]
 
# select  -- R version
 
#
# Note the following
# USE_C<-FALSE
# source("apl.R")
# a<-array(1:100000,rep(10,5))
# system.time(for (i in 1:100) aplSelect(a,list(1:5,1:5,1:5,1:5,1:5))) 
#   user  system elapsed 
# 23.211   0.042  23.260 
# USE_C<-TRUE
# source("apl.R")
# system.time(for (i in 1:100) aplSelect(a,list(1:5,1:5,1:5,1:5,1:5))) 
#   user  system elapsed 
#  0.444   0.284   0.742 
#
 
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)
if (length(m) != length(x))
    stop("Incorrect vector length")
z[which(y)]<-x
return(z)
}
 
#  expand
 
aplExpand<-function(x,y,axis=1) {
    if (is.vector(x)) return(aplEXV(x,y))
    d<-dim(x); m<-which(y); e<-d; e[axis]<-m
    if (m != d[axis])
        stop("Incorrect dimension length")
    z<-array(0,e)
    apply(z,(1:n)[-axis],function(i) z[i]<-x[i])
}
 
#  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) {
    if (is.vector(x)) return(aplCRV(x,y))
    sx<-aplShape(x); sy<-aplShape(y); sk<-sx[k]
    if (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,x,k) {
    if (is.vector(a)) return(aplRTV(a,k))
    sa<-aplShape(a); sx<-aplShape(x)
    if (sx == 1) x<-array(x,sa[-k])
    if (!identical(sa[-k],aplShape(x)))
        stop("Index Error")
    z<-array(0,sa); sz<-sa; nz<-prod(sz); sk<-sz[k]
    for (i in 1:nz) {
        ivec<-aplEncode(i,sz)
        xx<-x[aplDecode(ivec[-k],sx)]
        ak<-rep(0,sk)
        for (j in 1:sk) {       
            jvec<-ivec; jvec[k]<-j
            ak[j]<-a[aplDecode(jvec,sz)]
            }       
        bk<-aplRTV(ak,xx)
        for (j in 1:sk) {       
            jvec<-ivec; jvec[k]<-j
            z[aplDecode(jvec,sz)]<-bk[j]
            }       
    }
return(z)
}
 
#  transpose
#
#  this has a somewhat tortuous relationship
#  with R's aperm. aplTranspose(a) and
#  aperm(a) are always the same. If x is a permutation
#  of 1,2,..,aplRank(a), then aperm(a,x) is
#  equal to aplTranspose(a,order(x)). If x is not
#  a permutation, then aperm is undefined.
#  If x has aplRank(a) elements equal to
#  one of 1,2,..,m with each of 1,2,..,m
#  occurring a least once, then aplTranspose(a,x)
#  has rank m. For obvious reasons dyadic
#  transpose is not used a great deal. For
#  permutations we could consequently define
#  aplTranspose(a,x) as aperm(a,order(x)).
#
#
#  Note the following
#  a<-array(1:100000,rep(10,5))
#  USE_C<-FALSE
#  source("apl.R")
#  system.time(aplTranspose(a))
#   user  system elapsed 
#  5.889   0.025   5.916 
#  USE_C<-TRUE
#  source("apl.R")
#  system.time(aplTranspose(a))
#   user  system elapsed 
#  0.035   0.007   0.041 
#  
 
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,k) {
    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[-k],sb[-k]))
        stop("Shape error")
    sz<-sa; sz[k]<-sz[k]+sb[k]; nz<-prod(sz); u<-unit(k,ra)
    z<-array(0,sz)
    for (i in 1:nz) {
        ivec<-aplEncode(i,sz)
        if (ivec[k] <- sa[k]) z[i]<-a[aplDecode(ivec,sa)]
            else z[i]<-b[aplDecode(ivec-sa[k]*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
 
# in most cases aplReduce could be computed by using
# apply. in aplReduce k are the dimensions along
# which a is reduced, in apply it's the dimensions
# along which a is NOT reduced. also in aplReduce
# f is a binary function, while in apply f is a 
# vector-reducing function such as sum or prod or max
# this means we can expect apply to be faster.
 
# Note the following
# a<-array(1:100000,rep(10,5))
# USE_C<-TRUE
# source("apl.R")
# system.time(aplReduce(a,1,max))
#  user  system elapsed 
# 0.134   0.009   0.142 
# USE_C<-FALSE
# source("apl.R")
# system.time(aplReduce(a,1,max))
#  user  system elapsed 
# 6.434   0.005   6.440 
# system.time(apply(a,2:5,max))
#   user  system elapsed 
#  0.086   0.001   0.088 
#
 
 
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 -- along a single dimension
 
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
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 -- along a single dimension
 
#
# Note the following (this is on the MacBook Air)
# USE_C<-TRUE
# source("apl.R")
# a<-array(1:10000,c(10,10,100))
# b<-array(1:10000,c(100,10,10))
# system.time(aplInnerProduct(a,b))
#   user  system elapsed 
#  3.060   0.010   3.072 
# USE_C<-FALSE
# source("apl.R")
# system.time(aplInnerProduct(a,b))
#   user  system elapsed 
# 59.286   0.274  59.927 
#
 
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) {
    if (!identical(typeof(a),typeof(b)))
        warning("Arguments of different types")
    arrTest(a); arrTest(b)
    sa<-aplShape(a); sb<-aplShape(b)
    na<-prod(sa); nb<-prod(sb)
    z<-array(0,sa)
    for (i in 1:na) {
        z[i]<-0; aa<-a[i]
        for (j in 1:nb)
            if (identical(aa,b[j])) z[i]<-1
        }
return(z)   
}
 
# .C versions
 
if (USE_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]]
}
 
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)
}
 
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))
}
 
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))
}
 
}
 
# 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)
 
arrTest<-function(x) if (!is.array(x)) stop("Functions in the apl package take array arguments")
 
# first lib
 
.First.lib <- function(lib, pkg) {
  library.dynam("apl", pkg, lib)
}
 
comments powered by Disqus

Info

Link to this snippet:


Download to Code Collector

To use the direct link to your snippet on CodeCollector.net either copy the html from the above section or drag the Download to Code Collector to where you would like to use it.

More Info:

Times Viewed: 668
Date Added: 2013-03-07 17:47:36
Last Modified: 2013-04-17 22:12:20

Web Analytics