[Returnanalytics-commits] r3464 - in pkg/PerformanceAnalytics: . R sandbox src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jul 4 00:13:53 CEST 2014
Author: rossbennett34
Date: 2014-07-04 00:13:53 +0200 (Fri, 04 Jul 2014)
New Revision: 3464
Added:
pkg/PerformanceAnalytics/sandbox/testMM.R
pkg/PerformanceAnalytics/src/
pkg/PerformanceAnalytics/src/momentF.f90
Modified:
pkg/PerformanceAnalytics/NAMESPACE
pkg/PerformanceAnalytics/R/MultivariateMoments.R
Log:
Adding compiled M3 and M4 functions
Modified: pkg/PerformanceAnalytics/NAMESPACE
===================================================================
--- pkg/PerformanceAnalytics/NAMESPACE 2014-07-03 21:18:40 UTC (rev 3463)
+++ pkg/PerformanceAnalytics/NAMESPACE 2014-07-03 22:13:53 UTC (rev 3464)
@@ -234,3 +234,4 @@
importFrom(stats,sd)
importFrom(utils,packageDescription)
importFrom(zoo,rollapply)
+useDynLib(PerformanceAnalytics)
Modified: pkg/PerformanceAnalytics/R/MultivariateMoments.R
===================================================================
--- pkg/PerformanceAnalytics/R/MultivariateMoments.R 2014-07-03 21:18:40 UTC (rev 3463)
+++ pkg/PerformanceAnalytics/R/MultivariateMoments.R 2014-07-03 22:13:53 UTC (rev 3464)
@@ -16,7 +16,7 @@
###############################################################################
-M3.MM = function(R,...){
+M3.MM.old = function(R,...){
cAssets = ncol(R); T = nrow(R);
if(!hasArg(mu)) mu = apply(R,2,'mean') else mu=mu=list(...)$mu
M3 = matrix(rep(0,cAssets^3),nrow=cAssets,ncol=cAssets^2)
@@ -28,7 +28,35 @@
return( 1/T*M3 );
}
-M4.MM = function(R,...){
+#'@useDynLib PerformanceAnalytics
+M3.MM = function(R, ...){
+ if(!hasArg(mu)) mu = colMeans(R) else mu=list(...)$mu
+
+ # Pass variables directly from R to the Fortran subroutine
+ # Note that we also need to allocate the object that the Fortran subroutine
+ # returns
+ # It is ok to allocate objects in R and then pass to the Fortran subroutine
+ # https://stat.ethz.ch/pipermail/r-devel/2005-September/034570.html
+ x <- coredata(R)
+
+ nr <- NROW(x)
+ nc <- NCOL(x)
+
+ CC <- matrix(0, nc*nc, nc)
+ om <- matrix(0, nc, nc*nc)
+
+ out <- .Fortran("M3",
+ x=x,
+ mu=as.double(mu),
+ nr=as.integer(nr),
+ nc=as.integer(nc),
+ C=CC,
+ om=om,
+ PACKAGE="PerformanceAnalytics")$om
+ out
+}
+
+M4.MM.old = function(R,...){
cAssets = ncol(R); T = nrow(R);
if(!hasArg(mu)) mu = apply(R,2,'mean') else mu=list(...)$mu
M4 = matrix(rep(0,cAssets^4),nrow=cAssets,ncol=cAssets^3);
@@ -40,6 +68,33 @@
return( 1/T*M4 );
}
+#'@useDynLib PerformanceAnalytics
+M4.MM = function(R, ...){
+ if(!hasArg(mu)) mu = colMeans(R) else mu=list(...)$mu
+
+ # Pass variables directly from R to the Fortran subroutine
+ # Note that we also need to allocate the object that the Fortran subroutine
+ # returns
+ # It is ok to allocate objects in R and then pass to the Fortran subroutine
+ # https://stat.ethz.ch/pipermail/r-devel/2005-September/034570.html
+ x <- coredata(R)
+ nr <- NROW(x)
+ nc <- NCOL(x)
+
+ DD <- matrix(0, nc*nc*nc, nc)
+ om <- matrix(0, nc, nc*nc*nc)
+
+ out <- .Fortran("M4",
+ x=x,
+ mu=as.double(mu),
+ nr=as.integer(nr),
+ nc=as.integer(nc),
+ D=DD,
+ om=om,
+ PACKAGE="PerformanceAnalytics")$om
+ out
+}
+
multivariate_mean = function(w,mu){
return( t(w)%*%mu )
}
Added: pkg/PerformanceAnalytics/sandbox/testMM.R
===================================================================
--- pkg/PerformanceAnalytics/sandbox/testMM.R (rev 0)
+++ pkg/PerformanceAnalytics/sandbox/testMM.R 2014-07-03 22:13:53 UTC (rev 3464)
@@ -0,0 +1,36 @@
+library(microbenchmark)
+
+# Test for equality
+set.seed(123)
+x <- matrix(rnorm(1000 * 5), 1000, 5)
+
+all.equal(PerformanceAnalytics:::M3.MM(x), PerformanceAnalytics:::M3.MM.old(x))
+all.equal(PerformanceAnalytics:::M4.MM(x), PerformanceAnalytics:::M4.MM.old(x))
+
+data(edhec)
+all.equal(PerformanceAnalytics:::M3.MM(edhec), PerformanceAnalytics:::M3.MM.old(edhec))
+all.equal(PerformanceAnalytics:::M4.MM(edhec), PerformanceAnalytics:::M4.MM.old(edhec))
+
+# Do benchmarks with bigger data set
+m <- 100000
+n <- 10
+set.seed(21)
+r <- matrix(rnorm(m * n), m, n)
+
+benchM3 <- microbenchmark(
+ PerformanceAnalytics:::M3.MM(r),
+ PerformanceAnalytics:::M3.MM.old(r),
+ times=10
+)
+
+benchM3
+plot(benchM3, main="M3 Benchmarks")
+
+benchM4 <- microbenchmark(
+ PerformanceAnalytics:::M4.MM(r),
+ PerformanceAnalytics:::M4.MM.old(r),
+ times=10
+)
+
+benchM4
+plot(benchM4, main="M4 Benchmarks")
Added: pkg/PerformanceAnalytics/src/momentF.f90
===================================================================
--- pkg/PerformanceAnalytics/src/momentF.f90 (rev 0)
+++ pkg/PerformanceAnalytics/src/momentF.f90 2014-07-03 22:13:53 UTC (rev 3464)
@@ -0,0 +1,189 @@
+
+ subroutine asVecCov1(ia, n, oa)
+ ! compute oa = ia * ia**T and return as an array
+
+ ! ia : input array
+ ! n : length of input array
+ ! oa : output array of length n * n
+ implicit none
+
+ ! input/output variables
+ integer :: n
+ double precision :: ia(n), oa(n*n)
+
+ ! local variables
+ integer :: i, j, ii
+
+ ! element access of a matrix in column major ordering with 1-based index
+ ! oa(i + (j-1) * n) = ia(i) * ia(j)
+
+ ii = 1
+ do j=1,n
+ do i=1,n
+ oa(ii) = ia(i) * ia(j)
+ ii = ii + 1
+ end do
+ end do
+
+ end subroutine asVecCov1
+
+ subroutine asVec(im, nr, nc, oa)
+ ! convert matrix to array
+ ! reshape might be more efficient?
+
+ ! ia : input matrix
+ ! nr : number of rows of im
+ ! nc : number of columns of im
+ ! oa : output array of length nr * nc
+ implicit none
+
+ ! input/output variables
+ integer :: nr, nc
+ double precision :: im(nr,nc), oa(nr * nc)
+
+ ! local variables
+ integer :: i, j
+
+ ! element access of a matrix in column major ordering with 1-based index
+ ! oa(i + (j-1) * n) = im(i,j)
+
+ do j=1,nc
+ do i=1,nr
+ oa(i + (j-1) * nr) = im(i,j)
+ end do
+ end do
+
+ end subroutine asVec
+
+ subroutine M3(x, mu, nr, nc, C, om)
+ ! compute the third moment of x
+
+ ! x : input matrix of data
+ ! mu : vector of means to center x
+ ! nr : number of rows of x
+ ! nc : number of columns of x
+ ! C : temporary matrix of dimension (nc * nc, n)
+ ! om : output matrix of dimension om(nc, nc * nc)
+
+ implicit none
+
+ ! input/output variables
+ integer :: nr, nc
+ double precision :: x(nr, nc), mu(nc), om(nc, nc*nc), C(nc*nc, nc)
+
+ ! local variables
+ integer :: i
+ double precision :: alpha, beta
+ double precision :: centret(nc), tccr(nc * nc)
+
+ alpha = 1.d0
+ beta = 1.d0
+
+ do i=1,nr
+ centret = x(i,:) - mu
+ ! the output we care about here is tccr
+ ! tccr = as.vector(tcrossprod(centret))
+ call asVecCov1(centret, nc, tccr)
+
+ ! C = tccr * centret**T + C
+ ! (nc*nc x 1) * (1 x nc) = (nc*nc x nc)
+ ! C := alpha*op( A )*op( B ) + beta*C
+ ! DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
+ ! Note that TRANSB="N". B is just an array, so I specify the dimension
+ ! of B as 1 x nc with the arguments N, K, and LDB.
+ call DGEMM('N', 'N', nc*nc, nc, 1, alpha, tccr, nc*nc, centret, 1, beta, C, nc*nc)
+ end do
+
+ om = transpose(C) / DBLE(nr)
+
+ ! C := alpha*op( A )*op( B ) + beta*C
+ ! DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
+ ! TRANSA : 'N' or 'n', op( A ) = A; 'T' or 't', op( A ) = A**T
+ ! TRANSB : 'N' or 'n', op( B ) = A; 'T' or 't', op( B ) = A**T
+ ! M : number of rows of the matrix op( A ) and of the matrix C.
+ ! N : number of columns of the matrix op( B ) and
+ ! the number of columns of the matrix C.
+ ! K : number of columns of the matrix op( A ) and the number of
+ ! rows of the matrix op( B ).
+ ! ALPHA : specifies the scalar alpha
+ ! A : array of DIMENSION ( LDA, ka ), where ka is k when
+ ! TRANSA = 'N' or 'n', and is m otherwise.
+ ! LDA : the first dimension of A (i.e. number of rows)
+ ! B : array of DIMENSION ( LDB, kb ), where kb is n when
+ ! TRANSB = 'N' or 'n', and is k otherwise.
+ ! LDB : first dimension of B (i.e. number of rows of B)
+ ! BETA : specifies the scalar alpha
+ ! C : array of DIMENSION ( LDC, n )
+ ! LDC : first dimension of C
+
+
+ end subroutine M3
+
+ subroutine M4(x, mu, nr, nc, D, om)
+ ! compute the fourth moment of x
+
+ ! x : input matrix of data
+ ! mu : vector of means to center x
+ ! nr : number of rows of x
+ ! nc : number of columns of x
+ ! D : temporary matrix of dimension (nc * nc * nc, n)
+ ! om : output matrix of dimension om(nc, nc * nc * nc)
+
+ implicit none
+
+ ! input/output variables
+ integer :: nr, nc
+ double precision :: x(nr, nc), mu(nc), om(nc, nc*nc*nc), D(nc*nc*nc, nc)
+
+ ! local variables
+ integer :: i
+ double precision :: alpha, beta, beta1
+ double precision :: centret(nc), tccr(nc * nc), tccr2(nc * nc * nc), C(nc*nc, nc)
+
+ alpha = 1.d0
+ beta = 0.d0
+ beta1 = 1.d0
+
+ do i=1,nr
+ centret = x(i,:) - mu
+ ! the output we care about here is tccr
+ ! tccr is an nc^2 x 1 matrix (i.e. 1-d vector or array)
+ ! tccr = as.vector(tcrossprod(centret))
+ call asVecCov1(centret, nc, tccr)
+
+ ! C = tcrossprod(tccr, centret)
+ call DGEMM('N', 'N', nc*nc, nc, 1, alpha, tccr, nc*nc, centret, 1, beta, C, nc*nc)
+
+ ! tccr2 is an N^3 x 1 matrix (i.e. 1-d vector or array)
+ ! convert C to a N^3 array and assign to variable tccr2
+ call asVec(C, nc*nc, nc, tccr2)
+
+ ! D = tccr2 * centret**T + D
+ call DGEMM('N', 'N', nc*nc*nc, nc, 1, alpha, tccr2, nc*nc*nc, centret, 1, beta1, D, nc*nc*nc)
+ end do
+
+ om = transpose(D) / DBLE(nr)
+
+ ! C := alpha*op( A )*op( B ) + beta*C
+ ! DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
+ ! TRANSA : 'N' or 'n', op( A ) = A; 'T' or 't', op( A ) = A**T
+ ! TRANSB : 'N' or 'n', op( B ) = A; 'T' or 't', op( B ) = A**T
+ ! M : number of rows of the matrix op( A ) and of the matrix C.
+ ! N : number of columns of the matrix op( B ) and
+ ! the number of columns of the matrix C.
+ ! K : number of columns of the matrix op( A ) and the number of
+ ! rows of the matrix op( B ).
+ ! ALPHA : specifies the scalar alpha
+ ! A : array of DIMENSION ( LDA, ka ), where ka is k when
+ ! TRANSA = 'N' or 'n', and is m otherwise.
+ ! LDA : the first dimension of A (i.e. number of rows)
+ ! B : array of DIMENSION ( LDB, kb ), where kb is n when
+ ! TRANSB = 'N' or 'n', and is k otherwise.
+ ! LDB : first dimension of B (i.e. number of rows of B)
+ ! BETA : specifies the scalar alpha
+ ! C : array of DIMENSION ( LDC, n )
+ ! LDC : first dimension of C
+
+
+ end subroutine M4
+
\ No newline at end of file
More information about the Returnanalytics-commits
mailing list