projects
/
morpheus.git
/ blobdiff
commit
grep
author
committer
pickaxe
?
search:
re
summary
|
shortlog
|
log
|
commit
|
commitdiff
|
tree
raw
|
inline
| side by side
Allow to optimize from a given matrix W
[morpheus.git]
/
pkg
/
R
/
optimParams.R
diff --git
a/pkg/R/optimParams.R
b/pkg/R/optimParams.R
index
f42571d
..
8dfe3dc
100644
(file)
--- 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 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
#'
#' @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
#' 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)))
{
# 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,
# 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)
}
# 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 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(
# @field W Weights matrix (initialized at identity)
#
setRefClass(
@@
-91,6
+93,7
@@
setRefClass(
K = "integer",
n = "integer",
d = "integer",
K = "integer",
n = "integer",
d = "integer",
+ nc = "integer",
# Weights matrix (generalized least square)
W = "matrix"
),
# Weights matrix (generalized least square)
W = "matrix"
),
@@
-102,7
+105,7
@@
setRefClass(
callSuper(...)
if (!hasArg("X") || !hasArg("Y") || !hasArg("K")
callSuper(...)
if (!hasArg("X") || !hasArg("Y") || !hasArg("K")
- || !hasArg("li") || !hasArg("Mhat"))
+ || !hasArg("li") || !hasArg("Mhat")
|| !hasArg("nc")
)
{
stop("Missing arguments")
}
{
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),
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),
+ pn
c=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)
},
W=as.double(W), PACKAGE="morpheus")$W, nrow=dd, ncol=dd )
MASS::ginv(Omega)
},
@@
-250,7
+253,8
@@
setRefClass(
res
},
res
},
- run = function(θ0)
+ # userW allows to bypass the W optimization by giving a W matrix
+ run = function(θ0, userW=NULL)
{
"Run optimization from θ0 with solver..."
{
"Run optimization from θ0 with solver..."
@@
-276,13
+280,14
@@
setRefClass(
stop("θ0$b: length K, no NA")
# (Re)Set W to identity, to allow several run from the same object
stop("θ0$b: length K, no NA")
# (Re)Set W to identity, to allow several run from the same object
- W <<-
diag(d+d^2+d^3)
+ W <<-
if (is.null(userW)) diag(d+d^2+d^3) else userW
- loopMax <- 2 #TODO: loopMax = 3 ? Seems not improving...
+ #NOTE: loopMax = 3 seems to not improve the final results.
+ loopMax <- ifelse(is.null(userW), 2, 1)
x_init <- linArgs(θ0)
for (loop in 1:loopMax)
{
x_init <- linArgs(θ0)
for (loop in 1:loopMax)
{
- op_res
=
constrOptim( x_init, .self$f, .self$grad_f,
+ op_res
<-
constrOptim( x_init, .self$f, .self$grad_f,
ui=cbind(
rbind( rep(-1,K-1), diag(K-1) ),
matrix(0, nrow=K, ncol=(d+1)*K) ),
ui=cbind(
rbind( rep(-1,K-1), diag(K-1) ),
matrix(0, nrow=K, ncol=(d+1)*K) ),