[R] Aggregate matrix in a 2 by 2 manor
Jeff Newmiller
jdnewmil at dcn.davis.ca.us
Sun Jul 31 08:13:04 CEST 2016
If you don't need all that FUN flexibility, you can get this done way
faster with the aperm and colMeans functions:
tst <- matrix( seq.int( 1440 * 360 )
, ncol = 1440
, nrow = 360
)
tst.small <- matrix( seq.int( 8 * 4 )
, ncol = 8
, nrow = 4
)
aggregate.nx.ny.expand.grid <- function( dta, nx = 2, ny = 2, FUN = mean, ... )
{
ilon <- seq( 1, ncol( dta ), nx )
ilat <- seq( 1, nrow( dta ), ny )
cells <- as.matrix( expand.grid( ilat, ilon ) )
blocks <- apply( cells
, 1
, function( x ) dta[ x[ 1 ]:( x[ 1 ] + 1 ), x[ 2 ]:( x[ 2 ] + 1 ) ] )
block.means <- colMeans( blocks )
matrix( block.means
, nrow( dta ) / ny
, ncol( dta ) / nx
)
}
aggregate.nx.ny.array.apply <- function( dta, nx = 2, ny = 2, FUN = mean, ... ) {
a <- array( dta
, dim = c( ny
, nrow( dta ) %/% ny
, nx
, ncol( dta ) %/% nx
)
)
apply( a, c( 2, 4 ), FUN, ... )
}
aggregate.nx.ny.array.aperm.mean <- function( dta, nx = 2, ny = 2, ... ) {
# number of rows in result
nnr <- nrow( dta ) %/% ny
# number of columns in result
nnc <- ncol( dta ) %/% nx
# number of values to take mean of
nxny <- nx * ny
# describe existing layout of values in dta as 4-d array
a1 <- array( dta, dim = c( ny, nnr, nx, nnc ) )
# swap data in dimensions 2 and 3
a2 <- aperm( a1, c( 1, 3, 2, 4 ) )
# treat first two dimensions as column vectors, remaining as columns
a3 <- matrix( a2, nrow = nxny )
# fast calculation of column means
v <- colMeans( a3, ... )
# reframe result vector as a matrix
matrix( v, ncol = nnc )
}
aggregate.nx.ny.array.aperm.apply <- function( dta, nx = 2, ny = 2, FUN = mean, ... ) {
# number of rows in result
nnr <- nrow( dta ) %/% ny
# number of columns in result
nnc <- ncol( dta ) %/% nx
# number of values to apply FUN to
nxny <- nx * ny
# describe existing layout of values in dta as 4-d array
a1 <- array( dta, dim = c( ny, nnr, nx, nnc ) )
# swap data in dimensions 2 and 3
a2 <- aperm( a1, c( 1, 3, 2, 4 ) )
# treat first two dimensions as column vectors, remaining as columns
a3 <- matrix( a2, nrow = nxny )
# apply FUN to column vectors
v <- apply( a3, 2, FUN = FUN, ... )
matrix( v, ncol = nnc )
}
test1 <- aggregate.nx.ny.expand.grid( tst )
test2 <- aggregate.nx.ny.array.apply( tst )
test3 <- aggregate.nx.ny.array.aperm.mean( tst )
test4 <- aggregate.nx.ny.array.aperm.apply( tst )
library(microbenchmark)
microbenchmark(
aggregate.nx.ny.expand.grid( tst, 2, 2, mean, na.rm = TRUE )
, aggregate.nx.ny.array.apply( tst, 2, 2, mean, na.rm = TRUE )
, aggregate.nx.ny.array.aperm.mean( tst, 2, 2, na.rm = TRUE )
, aggregate.nx.ny.array.aperm.apply( tst, 2, 2, mean, na.rm = TRUE )
)
#Unit: milliseconds
# expr min
# aggregate.nx.ny.expand.grid(tst, 2, 2, mean, na.rm = TRUE) 628.528322
# aggregate.nx.ny.array.apply(tst, 2, 2, mean, na.rm = TRUE) 846.883314
# aggregate.nx.ny.array.aperm.mean(tst, 2, 2, na.rm = TRUE) 8.904369
# aggregate.nx.ny.array.aperm.apply(tst, 2, 2, mean, na.rm = TRUE) 619.691851
# lq mean median uq max neval cld
# 675.470967 916.39630 778.54090 873.9754 2452.695 100 b
# 920.831966 1126.94691 1000.33830 1094.9233 3412.639 100 c
# 9.191747 21.98528 10.30099 15.9169 158.687 100 a
# 733.246331 936.73359 757.58383 844.2016 2824.557 100 b
On Sat, 30 Jul 2016, Jeff Newmiller wrote:
> For the record, the array.apply code can be fixed as below, but then it is slower than the expand.grid version.
>
> aggregate.nx.ny.array.apply <- function(dta,nx=2,ny=2, FUN=mean,...)
> {
> a <- array(dta, dim = c(ny, nrow( dta ) %/% ny, nx, ncol( dta ) %/% nx))
> apply( a, c(2, 4), FUN, ... )
> }
>
> --
> Sent from my phone. Please excuse my brevity.
>
> On July 30, 2016 11:06:16 AM PDT, "Anthoni, Peter (IMK)" <peter.anthoni at kit.edu> wrote:
>> Hi all,
>>
>> thanks for the suggestions, I did some timing tests, see below.
>> Unfortunately the aggregate.nx.ny.array.apply, does not produce the
>> expected result.
>> So the fastest seems to be the aggregate.nx.ny.expand.grid, though the
>> double for loop is not that much slower.
>>
>> many thanks
>> Peter
>>
>>> tst=matrix(1:(1440*360),ncol=1440,nrow=360)
>>> system.time( {for(i in 1:10)
>> tst_2x2=aggregate.nx.ny.forloop(tst,2,2,mean,na.rm=T)})
>> user system elapsed
>> 11.227 0.073 11.371
>>> system.time( {for(i in 1:10)
>> tst_2x2=aggregate.nx.ny.interaction(tst,2,2,mean,na.rm=T)})
>> user system elapsed
>> 26.354 0.475 26.880
>>> system.time( {for(i in 1:10)
>> tst_2x2=aggregate.nx.ny.expand.grid(tst,2,2,mean,na.rm=T)})
>> user system elapsed
>> 9.683 0.055 9.763
>>> system.time( {for(i in 1:10)
>> tst_2x2=aggregate.nx.ny.array.apply(tst,2,2,mean,na.rm=T)})
>> user system elapsed
>> 7.693 0.055 7.800
>>
>>> tst.small=matrix(1:(8*4),ncol=8,nrow=4)
>>> aggregate.nx.ny.forloop = function(data,nx=2,ny=2, FUN=mean,...)
>> + {
>> + nlon=nrow(data)
>> + nlat=ncol(data)
>> + newdata=matrix(NA,nrow=nlon/nx,ncol=nlat/ny)
>> + dim(newdata)
>> + for(ilon in seq(1,nlon,nx)) {
>> + for(ilat in seq(1,nlat,ny)) {
>> + ilon_new=1+(ilon-1)/nx
>> + ilat_new=1+(ilat-1)/ny
>> + newdata[ilon_new,ilat_new] = FUN(data[ilon+0:1,ilat+0:1],...)
>> + }
>> + }
>> + newdata
>> + }
>>> aggregate.nx.ny.forloop(tst.small)
>> [,1] [,2] [,3] [,4]
>> [1,] 3.5 11.5 19.5 27.5
>> [2,] 5.5 13.5 21.5 29.5
>>>
>>> aggregate.nx.ny.interaction = function(data,nx=2,ny=2, FUN=mean,...)
>> + {
>> +
>> + nlon=nrow(data)
>> + nlat=ncol(data)
>> + newdata=matrix(NA,nrow=nlon/nx,ncol=nlat/ny)
>> + newdata[] <- tapply( data, interaction( (row(data)+1) %/% 2,
>> (col(data)+1) %/% 2 ), FUN, ...)
>> + newdata
>> + }
>>> aggregate.nx.ny.interaction(tst.small)
>> [,1] [,2] [,3] [,4]
>> [1,] 3.5 11.5 19.5 27.5
>> [2,] 5.5 13.5 21.5 29.5
>>>
>>> aggregate.nx.ny.expand.grid = function(data,nx=2,ny=2, FUN=mean,...)
>> + {
>> + ilon <- seq(1,ncol(data),nx)
>> + ilat <- seq(1,nrow(data),ny)
>> + cells <- as.matrix(expand.grid(ilat, ilon))
>> + blocks <- apply(cells, 1, function(x)
>> data[x[1]:(x[1]+1),x[2]:(x[2]+1)])
>> + block.means <- colMeans(blocks)
>> + matrix(block.means, nrow(data)/ny, ncol(data)/nx)
>> + }
>>> aggregate.nx.ny.expand.grid(tst.small)
>> [,1] [,2] [,3] [,4]
>> [1,] 3.5 11.5 19.5 27.5
>> [2,] 5.5 13.5 21.5 29.5
>>>
>>> aggregate.nx.ny.array.apply = function(data,nx=2,ny=2, FUN=mean,...)
>> {
>> + a <- array(data, dim = c(ny, nrow( data ) %/% ny, ncol( data ) %/%
>> nx))
>> + apply( a, c(2, 3), FUN, ... )
>> + }
>>> aggregate.nx.ny.array.apply(tst.small)
>> [,1] [,2] [,3] [,4]
>> [1,] 1.5 5.5 9.5 13.5
>> [2,] 3.5 7.5 11.5 15.5
>>
>>
>>
>>> On 28 Jul 2016, at 00:26, David Winsemius <dwinsemius at comcast.net>
>> wrote:
>>>
>>>
>>>> On Jul 27, 2016, at 12:02 PM, Jeff Newmiller
>> <jdnewmil at dcn.davis.ca.us> wrote:
>>>>
>>>> An alternative (more compact, not necessarily faster, because apply
>> is still a for loop inside):
>>>>
>>>> f <- function( m, nx, ny ) {
>>>> # redefine the dimensions of my
>>>> a <- array( m
>>>> , dim = c( ny
>>>> , nrow( m ) %/% ny
>>>> , ncol( m ) %/% nx )
>>>> )
>>>> # apply mean over dim 1
>>>> apply( a, c( 2, 3 ), FUN=mean )
>>>> }
>>>> f( tst, nx, ny )
>>>
>>> Here's an apparently loopless strategy, although I suspect the code
>> for interaction (and maybe tapply as well?) uses a loop.
>>>
>>>
>>> tst_2X2 <- matrix(NA, ,ncol=4,nrow=2)
>>>
>>> tst_2x2[] <- tapply( tst, interaction( (row(tst)+1) %/% 2,
>> (col(tst)+1) %/% 2 ), mean)
>>>
>>> tst_2x2
>>>
>>> [,1] [,2] [,3] [,4]
>>> [1,] 3.5 11.5 19.5 27.5
>>> [2,] 5.5 13.5 21.5 29.5
>>>
>>> --
>>> David.
>>>
>>>
>>>>
>>>> --
>>>> Sent from my phone. Please excuse my brevity.
>>>>
>>>> On July 27, 2016 9:08:32 AM PDT, David L Carlson <dcarlson at tamu.edu>
>> wrote:
>>>>> This should be faster. It uses apply() across the blocks.
>>>>>
>>>>>> ilon <- seq(1,8,nx)
>>>>>> ilat <- seq(1,4,ny)
>>>>>> cells <- as.matrix(expand.grid(ilat, ilon))
>>>>>> blocks <- apply(cells, 1, function(x) tst[x[1]:(x[1]+1),
>>>>> x[2]:(x[2]+1)])
>>>>>> block.means <- colMeans(blocks)
>>>>>> tst_2x2 <- matrix(block.means, 2, 4)
>>>>>> tst_2x2
>>>>> [,1] [,2] [,3] [,4]
>>>>> [1,] 3.5 11.5 19.5 27.5
>>>>> [2,] 5.5 13.5 21.5 29.5
>>>>>
>>>>> -------------------------------------
>>>>> David L Carlson
>>>>> Department of Anthropology
>>>>> Texas A&M University
>>>>> College Station, TX 77840-4352
>>>>>
>>>>>
>>>>>
>>>>> -----Original Message-----
>>>>> From: R-help [mailto:r-help-bounces at r-poject.org] On Behalf Of
>> Anthoni,
>>>>> Peter (IMK)
>>>>> Sent: Wednesday, July 27, 2016 6:14 AM
>>>>> To: r-help at r-project.org
>>>>> Subject: [R] Aggregate matrix in a 2 by 2 manor
>>>>>
>>>>> Hi all,
>>>>>
>>>>> I need to aggregate some matrix data (1440x720) to a lower
>> dimension
>>>>> (720x360) for lots of years and variables
>>>>>
>>>>> I can do double for loop, but that will be slow. Anybody know a
>> quicker
>>>>> way?
>>>>>
>>>>> here an example with a smaller matrix size:
>>>>>
>>>>> tst=matrix(1:(8*4),ncol=8,nrow=4)
>>>>> tst_2x2=matrix(NA,ncol=4,nrow=2)
>>>>> nx=2
>>>>> ny=2
>>>>> for(ilon in seq(1,8,nx)) {
>>>>> for (ilat in seq(1,4,ny)) {
>>>>> ilon_2x2=1+(ilon-1)/nx
>>>>> ilat_2x2=1+(ilat-1)/ny
>>>>> tst_2x2[ilat_2x2,ilon_2x2] = mean(tst[ilat+0:1,ilon+0:1])
>>>>> }
>>>>> }
>>>>>
>>>>> tst
>>>>> tst_2x2
>>>>>
>>>>>> tst
>>>>> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
>>>>> [1,] 1 5 9 13 17 21 25 29
>>>>> [2,] 2 6 10 14 18 22 26 30
>>>>> [3,] 3 7 11 15 19 23 27 31
>>>>> [4,] 4 8 12 16 20 24 28 32
>>>>>
>>>>>> tst_2x2
>>>>> [,1] [,2] [,3] [,4]
>>>>> [1,] 3.5 11.5 19.5 27.5
>>>>> [2,] 5.5 13.5 21.5 29.5
>>>>>
>>>>>
>>>>> I though a cast to 3d-array might do the trick and apply over the
>> new
>>>>> dimension, but that does not work, since it casts the data along
>> the
>>>>> row.
>>>>>> matrix(apply(array(tst,dim=c(nx,ny,8)),3,mean),nrow=nrow(tst)/ny)
>>>>> [,1] [,2] [,3] [,4]
>>>>> [1,] 2.5 10.5 18.5 26.5
>>>>> [2,] 6.5 14.5 22.5 30.5
>>>>>
>>>>>
>>>>> cheers
>>>>> Peter
>>>>>
>>>>> ______________________________________________
>>>>> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
>>>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>>>> PLEASE do read the posting guide
>>>>> http://www.R-project.org/posting-guide.html
>>>>> and provide commented, minimal, self-contained, reproducible code.
>>>>>
>>>>> ______________________________________________
>>>>> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
>>>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>>>> PLEASE do read the posting guide
>>>>> http://www.R-project.org/posting-guide.html
>>>>> and provide commented, minimal, self-contained, reproducible code.
>>>>
>>>> ______________________________________________
>>>> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
>>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>>>> and provide commented, minimal, self-contained, reproducible code.
>>>
>>> David Winsemius
>>> Alameda, CA, USA
>>>
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
---------------------------------------------------------------------------
Jeff Newmiller The ..... ..... Go Live...
DCN:<jdnewmil at dcn.davis.ca.us> Basics: ##.#. ##.#. Live Go...
Live: OO#.. Dead: OO#.. Playing
Research Engineer (Solar/Batteries O.O#. #.O#. with
/Software/Embedded Controllers) .OO#. .OO#. rocks...1k
More information about the R-help
mailing list