From 5af71d43f3f2dba21c6667939fcff88923af3b7b Mon Sep 17 00:00:00 2001 From: Benjamin Auder Date: Tue, 14 Jan 2020 18:24:36 +0100 Subject: [PATCH] Add a number_of_cores parameter for OpenMP // in Compute_Omega --- pkg/DESCRIPTION | 2 +- pkg/R/optimParams.R | 11 +++++++---- pkg/src/Makevars | 7 ++++--- pkg/src/functions.c | 23 ++++++++++++++++++----- pkg/src/morpheus_init.c | 4 ++-- reports/accuracy.R | 2 +- reports/multistart.R | 2 +- reports/run_time_cl.sh | 20 ++++++++++++-------- reports/timings.R | 4 ++-- 9 files changed, 48 insertions(+), 27 deletions(-) diff --git a/pkg/DESCRIPTION b/pkg/DESCRIPTION index bd92ed5..2a20240 100644 --- a/pkg/DESCRIPTION +++ b/pkg/DESCRIPTION @@ -6,7 +6,7 @@ Description: Mixture of logistic regressions parameters (H)estimation with (General Linear Model). For more details see chapter 3 in the PhD thesis of Mor-Absa Loum: , available here . -Version: 1.0-0 +Version: 1.0-2 Author: Benjamin Auder [aut,cre], Mor-Absa Loum [aut] Maintainer: Benjamin Auder diff --git a/pkg/R/optimParams.R b/pkg/R/optimParams.R index f42571d..411e010 100644 --- a/pkg/R/optimParams.R +++ b/pkg/R/optimParams.R @@ -5,6 +5,7 @@ #' @param K Number of populations. #' @param link The link type, 'logit' or 'probit'. #' @param M the empirical cross-moments between X and Y (optional) +#' @param nc Number of cores (default: 0 to use all) #' #' @return An object 'op' of class OptimParams, initialized so that #' \code{op$run(θ0)} outputs the list of optimized parameters @@ -36,7 +37,7 @@ #' o$f( o$linArgs(par1) )} #' #' @export -optimParams <- function(X, Y, K, link=c("logit","probit"), M=NULL) +optimParams <- function(X, Y, K, link=c("logit","probit"), M=NULL, nc=0) { # Check arguments if (!is.matrix(X) || any(is.na(X))) @@ -61,7 +62,7 @@ optimParams <- function(X, Y, K, link=c("logit","probit"), M=NULL) # Build and return optimization algorithm object methods::new("OptimParams", "li"=link, "X"=X, - "Y"=as.integer(Y), "K"=as.integer(K), "Mhat"=as.double(M)) + "Y"=as.integer(Y), "K"=as.integer(K), "Mhat"=as.double(M), "nc"=as.integer(nc)) } # Encapsulated optimization for p (proportions), β and b (regression parameters) @@ -76,6 +77,7 @@ optimParams <- function(X, Y, K, link=c("logit","probit"), M=NULL) # @field K Number of populations # @field n Number of sample points # @field d Number of dimensions +# @field nc Number of cores (OpenMP //) # @field W Weights matrix (initialized at identity) # setRefClass( @@ -91,6 +93,7 @@ setRefClass( K = "integer", n = "integer", d = "integer", + nc = "integer", # Weights matrix (generalized least square) W = "matrix" ), @@ -102,7 +105,7 @@ setRefClass( callSuper(...) if (!hasArg("X") || !hasArg("Y") || !hasArg("K") - || !hasArg("li") || !hasArg("Mhat")) + || !hasArg("li") || !hasArg("Mhat") || !hasArg("nc")) { stop("Missing arguments") } @@ -141,7 +144,7 @@ setRefClass( M <- Moments(θ) Omega <- matrix( .C("Compute_Omega", X=as.double(X), Y=as.integer(Y), M=as.double(M), - pn=as.integer(n), pd=as.integer(d), + pnc=as.integer(nc), pn=as.integer(n), pd=as.integer(d), W=as.double(W), PACKAGE="morpheus")$W, nrow=dd, ncol=dd ) MASS::ginv(Omega) }, diff --git a/pkg/src/Makevars b/pkg/src/Makevars index 1f1960a..ee13498 100644 --- a/pkg/src/Makevars +++ b/pkg/src/Makevars @@ -1,4 +1,5 @@ -# Override R weird defaults: -# https://stackoverflow.com/questions/23414448/r-makevars-file-to-overwrite-r-cmds-default-g-options +PKG_CFLAGS = $(SHLIB_OPENMP_CFLAGS) +PKG_LIBS = $(SHLIB_OPENMP_CFLAGS) + #%.o: %.c -# $(CC) -std=c11 -I"$(R_INCLUDE_DIR)" -DNDEBUG -fPIC -O2 -c $< -o $@ +# $(CC) -std=c11 -I"$(R_INCLUDE_DIR)" -DNDEBUG -fPIC -fopenmp -O2 -c $< -o $@ diff --git a/pkg/src/functions.c b/pkg/src/functions.c index b634c98..94ea12d 100644 --- a/pkg/src/functions.c +++ b/pkg/src/functions.c @@ -1,4 +1,5 @@ #include +#include // Index matrix (by columns) #define mi(i, j, d1, d2) (j*d1 + i) @@ -50,9 +51,9 @@ void Moments_M3(double* X, double* Y, int* pn, int* pd, double* M3) // W = 1/N sum( t(g(Zi,theta)) g(Zi,theta) ) // with g(Zi, theta) = i-th contribution to all moments (size dim) - real moments -void Compute_Omega(double* X, int* Y, double* M, int* pn, int* pd, double* W) +void Compute_Omega(double* X, int* Y, double* M, int* pnc, int* pn, int* pd, double* W) { - int n=*pn, d=*pd; + int nc=*pnc, n=*pn, d=*pd; int dim = d + d*d + d*d*d; //double* W = (double*)malloc(dim*dim*sizeof(double)); @@ -63,6 +64,8 @@ void Compute_Omega(double* X, int* Y, double* M, int* pn, int* pd, double* W) W[j*dim+k] = 0.0; } double* g = (double*)malloc(dim*sizeof(double)); + omp_set_num_threads(nc >= 1 ? nc : omp_get_num_procs()); + #pragma omp parallel for for (int i=0; i=0; j--) { - for (int k=0; k=0; k--) + W[baseIdx+k] += gj * g[k]; } } + // Normalize W: x 1/n + for (int j=0; jout 2>&1 +R --slave --args N=$N nc=$nc link=$link out 2>&1 diff --git a/reports/timings.R b/reports/timings.R index 704fd82..6d26bc0 100644 --- a/reports/timings.R +++ b/reports/timings.R @@ -16,8 +16,8 @@ ourOptim <- function(X, Y, K, link) { M <- computeMoments(X, Y) mu <- computeMu(X, Y, list(K=K,M=M)) - x_init = c(1/2, as.double(mu), c(0,0)) - optimParams(K, link, list(M=M))$run(x_init) + x_init = list(p=rep(1/K,K-1), beta=mu, b=rep(0,K)) + optimParams(X, Y, K, link, M, 1)$run(x_init) NULL } -- 2.44.0