From: Benjamin Auder Date: Tue, 23 Jan 2018 01:08:08 +0000 (+0100) Subject: First commit X-Git-Url: https://git.auder.net/%7B%7B%20asset%28%27mixstore/css/static/img/doc/%7B%7B?a=commitdiff_plain;h=cbd88fe5729bf206a784238a2637aa60e697fcdc;p=morpheus.git First commit --- cbd88fe5729bf206a784238a2637aa60e697fcdc diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 0000000..2bc4d4d --- /dev/null +++ b/.gitattributes @@ -0,0 +1,6 @@ +*.pdf filter=fat +*.html filter=fat +*.zip filter=fat +*.tar.xz filter=fat +*.data filter=fat +*.png filter=fat diff --git a/.gitfat b/.gitfat new file mode 100644 index 0000000..ec2f583 --- /dev/null +++ b/.gitfat @@ -0,0 +1,2 @@ +[rsync] +remote = gitfat@auder.net:~/files/morpheus diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..83249ef --- /dev/null +++ b/.gitignore @@ -0,0 +1,23 @@ +#ignore roxygen2 generated files +NAMESPACE +*.Rd +!*-package.Rd + +#ignore temporary files +*~ +*.swp + +#ignore package for CRAN +/pkg-cran/ + +#ignore R session files +.Rhistory +.RData + +#ignore R CMD build/check genrated files +/*.Rcheck/ +/*.tar.gz + +#ignore generated object files +*.[oa] +*.so diff --git a/README.md b/README.md new file mode 100644 index 0000000..dbfd60a --- /dev/null +++ b/README.md @@ -0,0 +1,17 @@ +# Estimate parameters of mixtures of logistic regressions + +This code is one applied part in the PhD thesis of [Mor-Absa Loum](https://fr.linkedin.com/in/mor-absa-loum-3372aa73) + +## Description + +Mixture of lOgistic Regressions Parameters (H)Estimation with (U)Spectral methods. +The main methods take d-dimensional inputs + a vector of binary outputs, and return +parameters according to the GLMs mixture model (please see the package vignette). + +NOTE: greek unicode letters are used in the code, because it's much nicer to write λ, β +and Σ than lambda, beta and Sigma - and also more importantly because it works well :) +...However CRAN demands ASCII files - therefore the presence of to-cran.sh script. + +## Example + +Install the package, then ?multiRun is a possible starting point. diff --git a/TODO b/TODO new file mode 100644 index 0000000..63d0cb4 --- /dev/null +++ b/TODO @@ -0,0 +1,12 @@ +simus point initial chaque M-C --> 3 ou 5 ? multiples de mu ? perturbés ? +n petit, mu = mauvais point initial ? plusieurs ? combien ? recommandations ? +n grand : mu OK ? pour quel n en fonction de d ? + +Pour la vérification de M_3 indiquée par Francis, juste un point: +Si on note B la matrice qui contient les beta en colonne, et si on note A l'inverse de B, c'est +A M_3 (\rho_l) A^{T} qui est diagonale +(car M_3 (\rho_l) = B D_l B^{T} pour une matrice D_l diagonale). + +Partir du vrai beta + vrais param :: voir si on converge loin de beta (n diminue...) +--> partir du vrai mu, aussi. +betas sparses, bcp de 0 ?! (comment ça se comporte ?) diff --git a/doc/Tensor-Kuleshov2015.pdf b/doc/Tensor-Kuleshov2015.pdf new file mode 100644 index 0000000..c50c073 --- /dev/null +++ b/doc/Tensor-Kuleshov2015.pdf @@ -0,0 +1 @@ +#$# git-fat 16f72d1107b66f90f27dd19ae241fc1a528971a5 478827 diff --git a/patch_Bettina/FLXMRglm.R b/patch_Bettina/FLXMRglm.R new file mode 100644 index 0000000..6fab5f5 --- /dev/null +++ b/patch_Bettina/FLXMRglm.R @@ -0,0 +1,139 @@ +FLXMRglm <- function(formula=.~., family=gaussian, offset=NULL) +{ + if (is.character(family)) + family <- get(family, mode = "function", envir = parent.frame()) + if (is.function(family)) + family <- family() + if (is.null(family$family)) { + print(family) + stop("'family' not recognized") + } + glmrefit <- function(x, y, w) { + fit <- c(glm.fit(x, y, weights=w, offset=offset, + family=family), + list(call = sys.call(), offset = offset, + control = eval(formals(glm.fit)$control), + method = "weighted.glm.fit")) + fit$df.null <- sum(w) + fit$df.null - fit$df.residual - fit$rank + fit$df.residual <- sum(w) - fit$rank + fit$x <- x + fit + } + + z <- new("FLXMRglm", weighted=TRUE, formula=formula, + name=paste("FLXMRglm", family$family, sep=":"), offset = offset, + family=family$family, refit=glmrefit) + z@preproc.y <- function(x){ + if (ncol(x) > 1) + stop(paste("for the", family$family, "family y must be univariate")) + x + } + + if(family$family=="gaussian"){ + z@defineComponent <- function(para) { + predict <- function(x, ...) { + dotarg = list(...) + if("offset" %in% names(dotarg)) offset <- dotarg$offset + p <- x %*% para$coef + if (!is.null(offset)) p <- p + offset + family$linkinv(p) + } + + logLik <- function(x, y, ...) + dnorm(y, mean=predict(x, ...), sd=para$sigma, log=TRUE) + + new("FLXcomponent", + parameters=list(coef=para$coef, sigma=para$sigma), + logLik=logLik, predict=predict, + df=para$df) + } + + z@fit <- function(x, y, w, component){ + fit <- glm.fit(x, y, w=w, offset=offset, family = family) + z@defineComponent(para = list(coef = coef(fit), df = ncol(x)+1, + sigma = sqrt(sum(fit$weights * fit$residuals^2 / + mean(fit$weights))/ (nrow(x)-fit$rank)))) + } + } + else if(family$family=="binomial"){ + z@preproc.y <- function(x){ + if (ncol(x) != 2) + stop("for the binomial family, y must be a 2 column matrix\n", + "where col 1 is no. successes and col 2 is no. failures") + if (any(x < 0)) + stop("negative values are not allowed for the binomial family") + x + } + z@defineComponent <- function(para) { + predict <- function(x, ...) { + dotarg = list(...) + if("offset" %in% names(dotarg)) offset <- dotarg$offset + p <- x %*% para$coef + if (!is.null(offset)) p <- p + offset + family$linkinv(p) + } + logLik <- function(x, y, ...) + dbinom(y[,1], size=rowSums(y), prob=predict(x, ...), log=TRUE) + + new("FLXcomponent", + parameters=list(coef=para$coef), + logLik=logLik, predict=predict, + df=para$df) + } + + z@fit <- function(x, y, w, component){ + fit <- glm.fit(x, y, weights=w, family=family, offset=offset, start=component$coef) + z@defineComponent(para = list(coef = coef(fit), df = ncol(x))) + } + } + else if(family$family=="poisson"){ + z@defineComponent <- function(para) { + predict <- function(x, ...) { + dotarg = list(...) + if("offset" %in% names(dotarg)) offset <- dotarg$offset + p <- x %*% para$coef + if (!is.null(offset)) p <- p + offset + family$linkinv(p) + } + logLik <- function(x, y, ...) + dpois(y, lambda=predict(x, ...), log=TRUE) + + new("FLXcomponent", + parameters=list(coef=para$coef), + logLik=logLik, predict=predict, + df=para$df) + } + + z@fit <- function(x, y, w, component){ + fit <- glm.fit(x, y, weights=w, family=family, offset=offset, start=component$coef) + z@defineComponent(para = list(coef = coef(fit), df = ncol(x))) + } + } + else if(family$family=="Gamma"){ + z@defineComponent <- function(para) { + predict <- function(x, ...) { + dotarg = list(...) + if("offset" %in% names(dotarg)) offset <- dotarg$offset + p <- x %*% para$coef + if (!is.null(offset)) p <- p + offset + family$linkinv(p) + } + logLik <- function(x, y, ...) + dgamma(y, shape = para$shape, scale=predict(x, ...)/para$shape, log=TRUE) + + new("FLXcomponent", + parameters = list(coef = para$coef, shape = para$shape), + predict = predict, logLik = logLik, + df = para$df) + } + + z@fit <- function(x, y, w, component){ + fit <- glm.fit(x, y, weights=w, family=family, offset=offset, start=component$coef) + z@defineComponent(para = list(coef = coef(fit), df = ncol(x)+1, + shape = sum(fit$prior.weights)/fit$deviance)) + + } + } + else stop(paste("Unknown family", family)) + z +} diff --git a/patch_Bettina/code.R b/patch_Bettina/code.R new file mode 100644 index 0000000..829d82d --- /dev/null +++ b/patch_Bettina/code.R @@ -0,0 +1,25 @@ +library("flexmix") +data("tribolium", package = "flexmix") +set.seed(1234) +## uses FLXMRglm() from package. +tribMix <- initFlexmix(cbind(Remaining, Total - Remaining) ~ Species, + k = 2, nrep = 5, data = tribolium, + model = FLXMRglm(family = "binomial")) +parameters(tribMix); logLik(tribMix) +source("FLXMRglm.R") +set.seed(1234) +## uses FLXMRglm() from source file which allows to specify link. +tribMixNew <- initFlexmix(cbind(Remaining, Total - Remaining) ~ Species, + k = 2, nrep = 5, data = tribolium, + model = FLXMRglm(family = "binomial")) +parameters(tribMixNew); logLik(tribMixNew) +set.seed(1234) +tribMixlogit <- initFlexmix(cbind(Remaining, Total - Remaining) ~ Species, + k = 2, nrep = 5, data = tribolium, + model = FLXMRglm(family = binomial(link = logit))) +parameters(tribMixlogit); logLik(tribMixlogit) +set.seed(1234) +tribMixprobit <- initFlexmix(cbind(Remaining, Total - Remaining) ~ Species, + k = 2, nrep = 5, data = tribolium, + model = FLXMRglm(family = binomial(link = probit))) +parameters(tribMixprobit); logLik(tribMixprobit) diff --git a/pkg/DESCRIPTION b/pkg/DESCRIPTION new file mode 100644 index 0000000..eb5e6f6 --- /dev/null +++ b/pkg/DESCRIPTION @@ -0,0 +1,34 @@ +Package: morpheus +Title: Estimate Parameters of Mixtures of Logistic Regressions +Description: Mixture of lOgistic Regressions Parameters (H)Estimation with + (U)Spectral methods. The main methods take d-dimensional inputs and a vector + of binary outputs, and return parameters according to the GLMs mixture model + (please refer to the package vignette). +Version: 0.2-0 +Author: Benjamin Auder [aut,cre], + Mor-Absa Loum [aut] +Maintainer: Benjamin Auder +Depends: + R (>= 3.0.0), +Imports: + MASS, + jointDiag, + methods, + pracma +Suggests: + flexmix, + parallel, + testthat, + roxygen2, + tensor, + nloptr +License: MIT + file LICENSE +RoxygenNote: 5.0.1 +Collate: + 'utils.R' + 'A_NAMESPACE.R' + 'computeMu.R' + 'multiRun.R' + 'optimParams.R' + 'plot.R' + 'sampleIO.R' diff --git a/pkg/LICENSE b/pkg/LICENSE new file mode 100644 index 0000000..5d50df5 --- /dev/null +++ b/pkg/LICENSE @@ -0,0 +1,20 @@ +Copyright (c) 2016-2017, Mor-Absa Loum ; 2017, Benjamin Auder + +Permission is hereby granted, free of charge, to any person obtaining +a copy of this software and associated documentation files (the +"Software"), to deal in the Software without restriction, including +without limitation the rights to use, copy, modify, merge, publish, +distribute, sublicense, and/or sell copies of the Software, and to +permit persons to whom the Software is furnished to do so, subject to +the following conditions: + +The above copyright notice and this permission notice shall be +included in all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE +LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION +OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION +WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. diff --git a/pkg/R/A_NAMESPACE.R b/pkg/R/A_NAMESPACE.R new file mode 100644 index 0000000..ea2711e --- /dev/null +++ b/pkg/R/A_NAMESPACE.R @@ -0,0 +1,11 @@ +#' @include utils.R +#' +#' @useDynLib morpheus +#' +#' @importFrom jointDiag ajd +#' @importFrom stats rbinom rmultinom rnorm pnorm runif integrate +#' @importFrom graphics boxplot barplot hist par +#' @importFrom methods new +#' @importFrom pracma integral +#' +NULL diff --git a/pkg/R/computeMu.R b/pkg/R/computeMu.R new file mode 100644 index 0000000..87cc680 --- /dev/null +++ b/pkg/R/computeMu.R @@ -0,0 +1,92 @@ +#' Compute μ +#' +#' Estimate the normalized columns μ of the β matrix parameter in a mixture of +#' logistic regressions models, with a spectral method described in the package vignette. +#' +#' @param X Matrix of input data (size nxd) +#' @param Y Vector of binary outputs (size n) +#' @param optargs List of optional argument: +#' \itemize{ +#' \item 'jd_method', joint diagonalization method from the package jointDiag: +#' 'uwedge' (default) or 'jedi'. +#' \item 'jd_nvects', number of random vectors for joint-diagonalization +#' (or 0 for p=d, canonical basis by default) +#' \item 'M', moments of order 1,2,3: will be computed if not provided. +#' \item 'K', number of populations (estimated with ranks of M2 if not given) +#' } +#' +#' @return The estimated normalized parameters as columns of a matrix μ of size dxK +#' +#' @seealso \code{multiRun} to estimate statistics based on μ, +#' and \code{generateSampleIO} for I/O random generation. +#' +#' @examples +#' io = generateSampleIO(10000, 1/2, matrix(c(1,0,0,1),ncol=2), c(0,0), "probit") +#' μ = computeMu(io$X, io$Y, list(K=2)) #or just X and Y for estimated K +#' @export +computeMu = function(X, Y, optargs=list()) +{ + if (!is.matrix(X) || !is.numeric(X) || any(is.na(X))) + stop("X: real matrix, no NA") + n = nrow(X) + d = ncol(X) + if (!is.numeric(Y) || length(Y)!=n || any(Y!=0 & Y!=1)) + stop("Y: vector of 0 and 1, size nrow(X), no NA") + if (!is.list(optargs)) + stop("optargs: list") + + # Step 0: Obtain the empirically estimated moments tensor, estimate also K + M = if (is.null(optargs$M)) computeMoments(X,Y) else optargs$M + K = optargs$K + if (is.null(K)) + { + # TODO: improve this basic heuristic + Σ = svd(M[[2]])$d + large_ratio <- ( abs(Σ[-d] / Σ[-1]) > 3 ) + K <- if (any(large_ratio)) max(2, which.min(large_ratio)) else d + } + + # Step 1: generate a family of d matrices to joint-diagonalize to increase robustness + d = ncol(X) + fixed_design = FALSE + jd_nvects = ifelse(!is.null(optargs$jd_nvects), optargs$jd_nvects, 0) + if (jd_nvects == 0) + { + jd_nvects = d + fixed_design = TRUE + } + M2_t = array(dim=c(d,d,jd_nvects)) + for (i in seq_len(jd_nvects)) + { + rho = if (fixed_design) c(rep(0,i-1),1,rep(0,d-i)) else normalize( rnorm(d) ) + M2_t[,,i] = .T_I_I_w(M[[3]],rho) + } + + # Step 2: obtain factors u_i (and their inverse) from the joint diagonalisation of M2_t + jd_method = ifelse(!is.null(optargs$jd_method), optargs$jd_method, "uwedge") + V = + if (jd_nvects > 1) { + #NOTE: increasing itermax does not help to converge, thus we suppress warnings + suppressWarnings({jd = jointDiag::ajd(M2_t, method=jd_method)}) +# if (jd_method=="uwedge") jd$B else solve(jd$A) + if (jd_method=="uwedge") jd$B else MASS::ginv(jd$A) + } + else + eigen(M2_t[,,1])$vectors + + # Step 3: obtain final factors from joint diagonalisation of T(I,I,u_i) + M2_t = array(dim=c(d,d,K)) + for (i in seq_len(K)) + M2_t[,,i] = .T_I_I_w(M[[3]],V[,i]) + suppressWarnings({jd = jointDiag::ajd(M2_t, method=jd_method)}) +# U = if (jd_method=="uwedge") solve(jd$B) else jd$A + U = if (jd_method=="uwedge") MASS::ginv(jd$B) else jd$A + μ = normalize(U[,1:K]) + + # M1 also writes M1 = sum_k coeff_k * μ_k, where coeff_k >= 0 + # ==> search decomposition of vector M1 onto the (truncated) basis μ (of size dxK) + # This is a linear system μ %*% C = M1 with C of size K ==> C = psinv(μ) %*% M1 + C = MASS::ginv(μ) %*% M[[1]] + μ[,C < 0] = - μ[,C < 0] + μ +} diff --git a/pkg/R/multiRun.R b/pkg/R/multiRun.R new file mode 100644 index 0000000..0a5c843 --- /dev/null +++ b/pkg/R/multiRun.R @@ -0,0 +1,125 @@ +#' multiRun +#' +#' Estimate N times some parameters, outputs of some list of functions. This method is +#' thus very generic, allowing typically bootstrap or Monte-Carlo estimations of matrices +#' μ or β. Passing a list of functions opens the possibility to compare them on a fair +#' basis (exact same inputs). It's even possible to compare methods on some deterministic +#' design of experiments. +#' +#' @param fargs List of arguments for the estimation functions +#' @param estimParams List of nf function(s) to apply on fargs - shared signature +#' @param prepareArgs Prepare arguments for the functions inside estimParams +#' @param N Number of runs +#' @param ncores Number of cores for parallel runs (<=1: sequential) +#' @param verbose TRUE to indicate runs + methods numbers +#' +#' @return A list of nf aggregates of N results (matrices). +#' +#' @examples +#' \dontrun{ +#' β <- matrix(c(1,-2,3,1),ncol=2) +#' +#' # Bootstrap + computeMu, morpheus VS flexmix ; assumes fargs first 3 elts X,Y,K +#' io <- generateSampleIO(n=1000, p=1/2, β=β, b=c(0,0), "logit") +#' μ <- normalize(β) +#' res <- multiRun(list(X=io$X,Y=io$Y,optargs=list(K=2,jd_nvects=0)), list( +#' # morpheus +#' function(fargs) { +#' library(morpheus) +#' ind <- fargs$ind +#' computeMu(fargs$X[ind,],fargs$Y[ind],fargs$optargs) +#' }, +#' # flexmix +#' function(fargs) { +#' library(flexmix) +#' ind <- fargs$ind +#' K <- fargs$optargs$K +#' dat = as.data.frame( cbind(fargs$Y[ind],fargs$X[ind,]) ) +#' out = refit( flexmix( cbind(V1, 1 - V1) ~ 0+., data=dat, k=K, +#' model=FLXMRglm(family="binomial") ) ) +#' normalize( matrix(out@@coef[1:(ncol(fargs$X)*K)], ncol=K) ) +#' } ), +#' prepareArgs = function(fargs) { +#' fargs$ind <- sample(1:nrow(fargs$X),replace=TRUE) +#' fargs +#' }, N=10, ncores=3) +#' for (i in 1:2) +#' res[[i]] <- alignMatrices(res[[i]], ref=μ, ls_mode="exact") +#' +#' # Monte-Carlo + optimParams from X,Y, morpheus VS flexmix ; first args n,p,β,b +#' res <- multiRun(list(n=1000,p=1/2,β=β,b=c(0,0),optargs=list(link="logit")),list( +#' # morpheus +#' function(fargs) { +#' library(morpheus) +#' K <- fargs$optargs$K +#' μ <- computeMu(fargs$X, fargs$Y, fargs$optargs) +#' V <- list( p=rep(1/K,K-1), β=μ, b=c(0,0) ) +#' optimParams(V,fargs$optargs)$β +#' }, +#' # flexmix +#' function(fargs) { +#' library(flexmix) +#' K <- fargs$optargs$K +#' dat <- as.data.frame( cbind(fargs$Y,fargs$X) ) +#' out <- refit( flexmix( cbind(V1, 1 - V1) ~ 0+., data=dat, k=K, +#' model=FLXMRglm(family="binomial") ) ) +#' sapply( seq_len(K), function(i) as.double( out@@components[[1]][[i]][,1] ) ) +#' } ), +#' prepareArgs = function(fargs) { +#' library(morpheus) +#' io = generateSampleIO(fargs$n, fargs$p, fargs$β, fargs$b, fargs$optargs$link) +#' fargs$X = io$X +#' fargs$Y = io$Y +#' fargs$optargs$K = ncol(fargs$β) +#' fargs$optargs$M = computeMoments(io$X,io$Y) +#' fargs +#' }, N=10, ncores=3) +#' for (i in 1:2) +#' res[[i]] <- alignMatrices(res[[i]], ref=β, ls_mode="exact")} +#' @export +multiRun <- function(fargs, estimParams, + prepareArgs = function(x) x, N=10, ncores=3, agg=lapply, verbose=FALSE) +{ + if (!is.list(fargs)) + stop("fargs: list") + # No checks on fargs: supposedly done in estimParams[[i]]() + if (!is.list(estimParams)) + estimParams = list(estimParams) + # Verify that the provided parameters estimations are indeed functions + lapply(seq_along(estimParams), function(i) { + if (!is.function(estimParams[[i]])) + stop("estimParams: list of function(fargs)") + }) + if (!is.numeric(N) || N < 1) + stop("N: positive integer") + + estimParamAtIndex <- function(index) + { + fargs <- prepareArgs(fargs) + if (verbose) + cat("Run ",index,"\n") + lapply(seq_along(estimParams), function(i) { + if (verbose) + cat(" Method ",i,"\n") + out <- estimParams[[i]](fargs) + if (is.list(out)) + do.call(rbind, out) + else + out + }) + } + + if (ncores > 1) + { + cl = parallel::makeCluster(ncores, outfile="") + parallel::clusterExport(cl, c("fargs","verbose"), environment()) + list_res = parallel::clusterApplyLB(cl, 1:N, estimParamAtIndex) + parallel::stopCluster(cl) + } + else + list_res = lapply(1:N, estimParamAtIndex) + + # De-interlace results: output one list per function + nf <- length(estimParams) + lapply( seq_len(nf), function(i) lapply(seq_len(N), function(j) list_res[[j]][[i]]) ) +} diff --git a/pkg/R/optimParams.R b/pkg/R/optimParams.R new file mode 100644 index 0000000..9185efb --- /dev/null +++ b/pkg/R/optimParams.R @@ -0,0 +1,332 @@ +#' Optimize parameters +#' +#' Optimize the parameters of a mixture of logistic regressions model, possibly using +#' \code{mu <- computeMu(...)} as a partial starting point. +#' +#' @param K Number of populations. +#' @param link The link type, 'logit' or 'probit'. +#' @param optargs a list with optional arguments: +#' \itemize{ +#' \item 'M' : list of moments of order 1,2,3: will be computed if not provided. +#' \item 'X,Y' : input/output, mandatory if moments not given +#' \item 'exact': use exact formulas when available? +#' } +#' +#' @return An object 'op' of class OptimParams, initialized so that \code{op$run(x0)} +#' outputs the list of optimized parameters +#' \itemize{ +#' \item p: proportions, size K +#' \item β: regression matrix, size dxK +#' \item b: intercepts, size K +#' } +#' x0 is a vector containing respectively the K-1 first elements of p, then β by +#' columns, and finally b: \code{x0 = c(p[1:(K-1)],as.double(β),b)}. +#' +#' @seealso \code{multiRun} to estimate statistics based on β, and +#' \code{generateSampleIO} for I/O random generation. +#' +#' @examples +#' # Optimize parameters from estimated μ +#' io = generateSampleIO(10000, 1/2, matrix(c(1,-2,3,1),ncol=2), c(0,0), "logit") +#' μ = computeMu(io$X, io$Y, list(K=2)) +#' M <- computeMoments(io$X, io$Y) +#' o <- optimParams(2, "logit", list(M=M)) +#' x0 <- c(1/2, as.double(μ), c(0,0)) +#' par0 <- o$run(x0) +#' # Compare with another starting point +#' x1 <- c(1/2, 2*as.double(μ), c(0,0)) +#' par1 <- o$run(x1) +#' o$f( o$linArgs(par0) ) +#' o$f( o$linArgs(par1) ) +#' @export +optimParams = function(K, link=c("logit","probit"), optargs=list()) +{ + # Check arguments + link <- match.arg(link) + if (!is.list(optargs)) + stop("optargs: list") + if (!is.numeric(K) || K < 2) + stop("K: integer >= 2") + + M <- optargs$M + if (is.null(M)) + { + if (is.null(optargs$X) || is.null(optargs$Y)) + stop("If moments are not provided, X and Y are required") + M <- computeMoments(optargs$X,optargs$Y) + } + + # TODO: field?! + exactComp <<- optargs$exact + if (is.null(exactComp)) + exactComp <<- FALSE + + # Build and return optimization algorithm object + methods::new("OptimParams", "li"=link, "M1"=as.double(M[[1]]), + "M2"=as.double(M[[2]]), "M3"=as.double(M[[3]]), "K"=as.integer(K)) +} + +# Encapsulated optimization for p (proportions), β and b (regression parameters) +# +# @field li Link, 'logit' or 'probit' +# @field M1 Estimated first-order moment +# @field M2 Estimated second-order moment (flattened) +# @field M3 Estimated third-order moment (flattened) +# @field K Number of populations +# @field d Number of dimensions +# +setRefClass( + Class = "OptimParams", + + fields = list( + # Inputs + li = "character", #link 'logit' or 'probit' + M1 = "numeric", #order-1 moment (vector size d) + M2 = "numeric", #M2 easier to process as a vector + M3 = "numeric", #M3 easier to process as a vector + # Dimensions + K = "integer", + d = "integer" + ), + + methods = list( + initialize = function(...) + { + "Check args and initialize K, d" + + callSuper(...) + if (!hasArg("li") || !hasArg("M1") || !hasArg("M2") || !hasArg("M3") + || !hasArg("K")) + { + stop("Missing arguments") + } + + d <<- length(M1) + }, + + expArgs = function(x) + { + "Expand individual arguments from vector x" + + list( + # p: dimension K-1, need to be completed + "p" = c(x[1:(K-1)], 1-sum(x[1:(K-1)])), + "β" = matrix(x[K:(K+d*K-1)], ncol=K), + "b" = x[(K+d*K):(K+(d+1)*K-1)]) + }, + + linArgs = function(o) + { + " Linearize vectors+matrices into a vector x" + + c(o$p[1:(K-1)], as.double(o$β), o$b) + }, + + f = function(x) + { + "Sum of squares (Mi - hat_Mi)^2 where Mi is obtained from formula" + + P <- expArgs(x) + p <- P$p + β <- P$β + λ <- sqrt(colSums(β^2)) + b <- P$b + + # Tensorial products β^2 = β2 and β^3 = β3 must be computed from current β1 + β2 <- apply(β, 2, function(col) col %o% col) + β3 <- apply(β, 2, function(col) col %o% col %o% col) + + return( + sum( ( β %*% (p * .G(li,1,λ,b)) - M1 )^2 ) + + sum( ( β2 %*% (p * .G(li,2,λ,b)) - M2 )^2 ) + + sum( ( β3 %*% (p * .G(li,3,λ,b)) - M3 )^2 ) ) + }, + + grad_f = function(x) + { + "Gradient of f, dimension (K-1) + d*K + K = (d+2)*K - 1" + + P <- expArgs(x) + p <- P$p + β <- P$β + λ <- sqrt(colSums(β^2)) + μ <- sweep(β, 2, λ, '/') + b <- P$b + + # Tensorial products β^2 = β2 and β^3 = β3 must be computed from current β1 + β2 <- apply(β, 2, function(col) col %o% col) + β3 <- apply(β, 2, function(col) col %o% col %o% col) + + # Some precomputations + G1 = .G(li,1,λ,b) + G2 = .G(li,2,λ,b) + G3 = .G(li,3,λ,b) + G4 = .G(li,4,λ,b) + G5 = .G(li,5,λ,b) + + # (Mi - hat_Mi)^2 ' == (Mi - hat_Mi)' 2(Mi - hat_Mi) = Mi' Fi + F1 = as.double( 2 * ( β %*% (p * G1) - M1 ) ) + F2 = as.double( 2 * ( β2 %*% (p * G2) - M2 ) ) + F3 = as.double( 2 * ( β3 %*% (p * G3) - M3 ) ) + + km1 = 1:(K-1) + grad <- #gradient on p + t( sweep(as.matrix(β [,km1]), 2, G1[km1], '*') - G1[K] * β [,K] ) %*% F1 + + t( sweep(as.matrix(β2[,km1]), 2, G2[km1], '*') - G2[K] * β2[,K] ) %*% F2 + + t( sweep(as.matrix(β3[,km1]), 2, G3[km1], '*') - G3[K] * β3[,K] ) %*% F3 + + grad_β <- matrix(nrow=d, ncol=K) + for (i in 1:d) + { + # i determines the derivated matrix dβ[2,3] + + dβ_left <- sweep(β, 2, p * G3 * β[i,], '*') + dβ_right <- matrix(0, nrow=d, ncol=K) + block <- i + dβ_right[block,] <- dβ_right[block,] + 1 + dβ <- dβ_left + sweep(dβ_right, 2, p * G1, '*') + + dβ2_left <- sweep(β2, 2, p * G4 * β[i,], '*') + dβ2_right <- do.call( rbind, lapply(1:d, function(j) { + sweep(dβ_right, 2, β[j,], '*') + }) ) + block <- ((i-1)*d+1):(i*d) + dβ2_right[block,] <- dβ2_right[block,] + β + dβ2 <- dβ2_left + sweep(dβ2_right, 2, p * G2, '*') + + dβ3_left <- sweep(β3, 2, p * G5 * β[i,], '*') + dβ3_right <- do.call( rbind, lapply(1:d, function(j) { + sweep(dβ2_right, 2, β[j,], '*') + }) ) + block <- ((i-1)*d*d+1):(i*d*d) + dβ3_right[block,] <- dβ3_right[block,] + β2 + dβ3 <- dβ3_left + sweep(dβ3_right, 2, p * G3, '*') + + grad_β[i,] <- t(dβ) %*% F1 + t(dβ2) %*% F2 + t(dβ3) %*% F3 + } + grad <- c(grad, as.double(grad_β)) + + grad = c(grad, #gradient on b + t( sweep(β, 2, p * G2, '*') ) %*% F1 + + t( sweep(β2, 2, p * G3, '*') ) %*% F2 + + t( sweep(β3, 2, p * G4, '*') ) %*% F3 ) + + grad + }, + + run = function(x0) + { + "Run optimization from x0 with solver..." + + if (!is.numeric(x0) || any(is.na(x0)) || length(x0) != (d+2)*K-1 + || any(x0[1:(K-1)] < 0) || sum(x0[1:(K-1)]) > 1) + { + stop("x0: numeric vector, no NA, length (d+2)*K-1, sum(x0[1:(K-1) >= 0]) <= 1") + } + + op_res = constrOptim( x0, .self$f, .self$grad_f, + ui=cbind( + rbind( rep(-1,K-1), diag(K-1) ), + matrix(0, nrow=K, ncol=(d+1)*K) ), + ci=c(-1,rep(0,K-1)) ) + + expArgs(op_res$par) + } + ) +) + +# Compute vectorial E[g^{(order)}(<β,x> + b)] with x~N(0,Id) (integral in R^d) +# = E[g^{(order)}(z)] with z~N(b,diag(λ)) +# +# @param link Link, 'logit' or 'probit' +# @param order Order of derivative +# @param λ Norm of columns of β +# @param b Intercept +# +.G <- function(link, order, λ, b) +{ + # NOTE: weird "integral divergent" error on inputs: + # link="probit"; order=2; λ=c(531.8099,586.8893,523.5816); b=c(-118.512674,-3.488020,2.109969) + # Switch to pracma package for that (but it seems slow...) + + if (exactComp && link == "probit") + { + # Use exact computations + sapply( seq_along(λ), function(k) { + .exactProbitIntegral(order, λ[k], b[k]) + }) + } + + else + { + # Numerical integration + sapply( seq_along(λ), function(k) { + res <- NULL + tryCatch({ + # Fast code, may fail: + res <- stats::integrate( + function(z) .deriv[[link]][[order]](λ[k]*z+b[k]) * exp(-z^2/2) / sqrt(2*pi), + lower=-Inf, upper=Inf )$value + }, error = function(e) { + # Robust slow code, no fails observed: + sink("/dev/null") #pracma package has some useless printed outputs... + res <- pracma::integral( + function(z) .deriv[[link]][[order]](λ[k]*z+b[k]) * exp(-z^2/2) / sqrt(2*pi), + xmin=-Inf, xmax=Inf, method="Kronrod") + sink() + }) + res + }) + } +} + +# TODO: check these computations (wrong atm) +.exactProbitIntegral <- function(order, λ, b) +{ + c1 = (1/sqrt(2*pi)) * exp( -.5 * b/((λ^2+1)^2) ) + if (order == 1) + return (c1) + c2 = b - λ^2 / (λ^2+1) + if (order == 2) + return (c1 * c2) + if (order == 3) + return (c1 * (λ^2 - 1 + c2^2)) + if (order == 4) + return ( (c1*c2/((λ^2+1)^2)) * (-λ^4*((b+1)^2+1) - + 2*λ^3 + λ^2*(2-2*b*(b-1)) + 6*λ + 3 - b^2) ) + if (order == 5) #only remaining case... + return ( c1 * (3*λ^4+c2^4+6*c1^2*(λ^2-1) - 6*λ^2 + 6) ) +} + +# Derivatives list: g^(k)(x) for links 'logit' and 'probit' +# +.deriv <- list( + "probit"=list( + # 'probit' derivatives list; + # TODO: exact values for the integral E[g^(k)(λz+b)] + function(x) exp(-x^2/2)/(sqrt(2*pi)), #g' + function(x) exp(-x^2/2)/(sqrt(2*pi)) * -x, #g'' + function(x) exp(-x^2/2)/(sqrt(2*pi)) * ( x^2 - 1), #g^(3) + function(x) exp(-x^2/2)/(sqrt(2*pi)) * (-x^3 + 3*x), #g^(4) + function(x) exp(-x^2/2)/(sqrt(2*pi)) * ( x^4 - 6*x^2 + 3) #g^(5) + ), + "logit"=list( + # Sigmoid derivatives list, obtained with http://www.derivative-calculator.net/ + # @seealso http://www.ece.uc.edu/~aminai/papers/minai_sigmoids_NN93.pdf + function(x) {e=exp(x); .zin(e /(e+1)^2)}, #g' + function(x) {e=exp(x); .zin(e*(-e + 1) /(e+1)^3)}, #g'' + function(x) {e=exp(x); .zin(e*( e^2 - 4*e + 1) /(e+1)^4)}, #g^(3) + function(x) {e=exp(x); .zin(e*(-e^3 + 11*e^2 - 11*e + 1) /(e+1)^5)}, #g^(4) + function(x) {e=exp(x); .zin(e*( e^4 - 26*e^3 + 66*e^2 - 26*e + 1)/(e+1)^6)} #g^(5) + ) +) + +# Utility for integration: "[return] zero if [argument is] NaN" (Inf / Inf divs) +# +# @param x Ratio of polynoms of exponentials, as in .S[[i]] +# +.zin <- function(x) +{ + x[is.nan(x)] <- 0. + x +} diff --git a/pkg/R/plot.R b/pkg/R/plot.R new file mode 100644 index 0000000..f94f19a --- /dev/null +++ b/pkg/R/plot.R @@ -0,0 +1,139 @@ +# extractParam +# +# Extract successive values of a projection of the parameter(s) +# +# @inheritParams plotHist +# +extractParam <- function(mr, x=1, y=1) +{ + # Obtain L vectors where L = number of res lists in mr + lapply( mr, function(mr_list) { + sapply(mr_list, function(m) m[x,y]) + } ) +} + +#' plotHist +#' +#' Plot histogram +#' +#' @param mr Output of multiRun(), list of lists of functions results +#' @param x Row index of the element inside the aggregated parameter +#' @param y Colomn index of the element inside the aggregated parameter +#' +#' @examples +#' \dontrun{ +#' β <- matrix(c(1,-2,3,1),ncol=2) +#' mr <- multiRun(...) #see bootstrap example in ?multiRun : return lists of mu_hat +#' μ <- normalize(β) +#' for (i in 1:2) +#' mr[[i]] <- alignMatrices(res[[i]], ref=μ, ls_mode="exact") +#' plotHist(mr, 2, 1) #second row, first column} +#' @export +plotHist <- function(mr, x, y) +{ + params <- extractParam(mr, x, y) + L = length(params) + # Plot histograms side by side + par(mfrow=c(1,L), cex.axis=1.5, cex.lab=1.5, mar=c(4.7,5,1,1)) + for (i in 1:L) + hist(params[[i]], breaks=40, freq=FALSE, xlab="Parameter value", ylab="Density") +} + +#' plotBox +#' +#' Draw boxplot +#' +#' @inheritParams plotHist +#' +#' @examples +#' #See example in ?plotHist +#' @export +plotBox <- function(mr, x, y) +{ + params <- extractParam(mr, x, y) + L = length(params) + # Plot boxplots side by side + par(mfrow=c(1,L), cex.axis=1.5, cex.lab=1.5, mar=c(4.7,5,1,1)) + for (i in 1:L) + boxplot(params[[i]], ylab="Parameter value") +} + +#' plotCoefs +#' +#' Draw coefs estimations + standard deviations +#' +#' @inheritParams plotHist +#' @param params True value of parameters matrix +#' +#' @examples +#' #See example in ?plotHist +#' @export +plotCoefs <- function(mr, params) +{ + nf <- length(mr) + L <- nrow(mr[[1]][[1]]) + K <- ncol(mr[[1]][[1]]) + + params_hat <- vector("list", nf) + stdev <- vector("list", nf) + for (i in 1:nf) + { + params_hat[[i]] <- matrix(nrow=L, ncol=K) + stdev[[i]] <- matrix(nrow=L, ncol=K) + } + for (x in 1:L) + { + for (y in 1:K) + { + estims <- extractParam(mr, x, y) + for (i in 1:nf) + { + params_hat[[i]][x,y] <- mean(estims[[i]]) +# stdev[[i]][x,y] <- sqrt( mean( (estims[[i]] - params[x,y])^2 ) ) + # HACK remove extreme quantile in estims[[i]] before computing sd() + stdev[[i]][x,y] <- sd( estims[[i]] [ estims[[i]] < max(estims[[i]]) & estims[[i]] > min(estims[[i]]) ] ) + } + } + } + + par(mfrow=c(1,nf), cex.axis=1.5, cex.lab=1.5, mar=c(4.7,5,1,1)) + params <- as.double(params) + o <- order(params) + for (i in 1:nf) + { + avg_param <- as.double(params_hat[[i]]) + std_param <- as.double(stdev[[i]]) + matplot(cbind(params[o],avg_param[o],avg_param[o]+std_param[o],avg_param[o]-std_param[o]), + col=c(2,1,1,1), lty=c(1,1,2,2), type="l", lwd=2, xlab="param", ylab="value") + } + + #print(o) #not returning o to avoid weird Jupyter issue... (TODO:) +} + +#' plotQn +#' +#' Draw 3D map of objective function values +#' +#' @param N Number of starting points +#' @param β Regression matrix (target) +#' @param link Link function (logit or probit) +#' +#' @export +plotQn <- function(N, n, p, β, b, link) +{ + d <- nrow(β) + K <- ncol(β) + io <- generateSampleIO(n, p, β, b, link) + op <- optimParams(K, link, list(X=io$X, Y=io$Y)) + # N random starting points gaussian (TODO: around true β?) + res <- matrix(nrow=d*K+1, ncol=N) + for (i in seq_len(N)) + { + β_init <- rnorm(d*K) + par <- op$run( c(rep(1/K,K-1), β_init, rep(0,K)) ) + par <- op$linArgs(par) + Qn <- op$f(par) + res[,i] = c(Qn, par[K:(K+d*K-1)]) + } + res #TODO: plot this, not just return it... +} diff --git a/pkg/R/sampleIO.R b/pkg/R/sampleIO.R new file mode 100644 index 0000000..5e45837 --- /dev/null +++ b/pkg/R/sampleIO.R @@ -0,0 +1,71 @@ +#' Generate sample inputs-outputs +#' +#' Generate input matrix X of size nxd and binary output of size n, where Y is subdivided +#' into K groups of proportions p. Inside one group, the probability law P(Y=1) is +#' described by the corresponding column parameter in the matrix β + intercept b. +#' +#' @param n Number of individuals +#' @param p Vector of K-1 populations relative proportions (sum <= 1) +#' @param β Vectors of model parameters for each population, of size dxK +#' @param b Vector of intercept values (use rep(0,K) for no intercept) +#' @param link Link type; "logit" or "probit" +#' +#' @return A list with +#' \itemize{ +#' \item{X: the input matrix (size nxd)} +#' \item{Y: the output vector (size n)} +#' \item{index: the population index (in 1:K) for each row in X} +#' } +#' +#' @export +generateSampleIO = function(n, p, β, b, link) +{ + # Check arguments + tryCatch({n = as.integer(n)}, error=function(e) stop("Cannot convert n to integer")) + if (length(n) > 1) + warning("n is a vector but should be scalar: only first element used") + if (n <= 0) + stop("n: positive integer") + if (!is.matrix(β) || !is.numeric(β) || any(is.na(β))) + stop("β: real matrix, no NAs") + K = ncol(β) + if (!is.numeric(p) || length(p)!=K-1 || any(is.na(p)) || any(p<0) || sum(p) > 1) + stop("p: positive vector of size K-1, no NA, sum<=1") + p <- c(p, 1-sum(p)) + if (!is.numeric(b) || length(b)!=K || any(is.na(b))) + stop("b: real vector of size K, no NA") + + #random generation of the size of each population in X~Y (unordered) + classes = rmultinom(1, n, p) + + d = nrow(β) + zero_mean = rep(0,d) + id_sigma = diag(rep(1,d)) + # Always consider an intercept (use b=0 for none) + d = d + 1 + β = rbind(β, b) + X = matrix(nrow=0, ncol=d) + Y = c() + index = c() + for (i in 1:ncol(β)) + { + index = c(index, rep(i, classes[i])) + newXblock = cbind( MASS::mvrnorm(classes[i], zero_mean, id_sigma), 1 ) + arg_link = newXblock%*%β[,i] + probas = + if (link == "logit") + { + e_arg_link = exp(arg_link) + e_arg_link / (1 + e_arg_link) + } + else #"probit" + pnorm(arg_link) + probas[is.nan(probas)] = 1 #overflow of exp(x) + X = rbind(X, newXblock) + Y = c( Y, vapply(probas, function(p) (rbinom(1,1,p)), 1) ) + } + shuffle = sample(n) + # Returned X should not contain an intercept column (it's an argument of estimation + # methods) + list("X"=X[shuffle,-d], "Y"=Y[shuffle], "index"=index[shuffle]) +} diff --git a/pkg/R/utils.R b/pkg/R/utils.R new file mode 100644 index 0000000..6d1c361 --- /dev/null +++ b/pkg/R/utils.R @@ -0,0 +1,169 @@ +#' normalize +#' +#' Normalize a vector or a matrix (by columns), using euclidian norm +#' +#' @param X Vector or matrix to be normalized +#' +#' @return The normalized matrix (1 column if X is a vector) +#' +#' @export +normalize = function(X) +{ + X = as.matrix(X) + norm2 = sqrt( colSums(X^2) ) + sweep(X, 2, norm2, '/') +} + +# Computes a tensor-vector product +# +# @param Te third-order tensor (size dxdxd) +# @param w vector of size d +# +# @return Matrix of size dxd +# +.T_I_I_w = function(Te, w) +{ + d = length(w) + Ma = matrix(0,nrow=d,ncol=d) + for (j in 1:d) + Ma = Ma + w[j] * Te[,,j] + Ma +} + +# Computes the second-order empirical moment between input X and output Y +# +# @param X matrix of covariates (of size n*d) +# @param Y vector of responses (of size n) +# +# @return Matrix of size dxd +# +.Moments_M2 = function(X, Y) +{ + n = nrow(X) + d = ncol(X) + M2 = matrix(0,nrow=d,ncol=d) + matrix( .C("Moments_M2", X=as.double(X), Y=as.double(Y), pn=as.integer(n), + pd=as.integer(d), M2=as.double(M2), PACKAGE="morpheus")$M2, nrow=d, ncol=d) +} + +# Computes the third-order empirical moment between input X and output Y +# +# @param X matrix of covariates (of size n*d) +# @param Y vector of responses (of size n) +# +# @return Array of size dxdxd +# +.Moments_M3 = function(X, Y) +{ + n = nrow(X) + d = ncol(X) + M3 = array(0,dim=c(d,d,d)) + array( .C("Moments_M3", X=as.double(X), Y=as.double(Y), pn=as.integer(n), + pd=as.integer(d), M3=as.double(M3), PACKAGE="morpheus")$M3, dim=c(d,d,d) ) +} + +#' computeMoments +#' +#' Compute cross-moments of order 1,2,3 from X,Y +#' +#' @inheritParams computeMu +#' +#' @return A list L where L[[i]] is the i-th cross-moment +#' +#' @export +computeMoments = function(X, Y) + list( colMeans(Y * X), .Moments_M2(X,Y), .Moments_M3(X,Y) ) + +# Find the optimal assignment (permutation) between two sets (minimize cost) +# +# @param distances The distances matrix, in columns (distances[i,j] is distance between i +# and j) +# +# @return A permutation minimizing cost +# +.hungarianAlgorithm = function(distances) +{ + n = nrow(distances) + .C("hungarianAlgorithm", distances=as.double(distances), pn=as.integer(n), + assignment=integer(n), PACKAGE="morpheus")$assignment +} + +#' alignMatrices +#' +#' Align a set of parameters matrices, with potential permutations. +#' +#' @param Ms A list of matrices, all of same size DxK +#' @param ref Either a reference matrix or "mean" to align on empirical mean +#' @param ls_mode How to compute the labels assignment: "exact" for exact algorithm +#' (default, but might be time-consuming, complexity is O(K^3) ), or "approx1", or +#' "approx2" to apply a greedy matching algorithm (heuristic) which for each column in +#' reference (resp. in current row) compare to all unassigned columns in current row +#' (resp. in reference) +#' +#' @return The aligned list (of matrices), of same size as Ms +#' +#' @export +alignMatrices = function(Ms, ref, ls_mode) +{ + if (!is.matrix(ref) && ref != "mean") + stop("ref: matrix or 'mean'") + if (!ls_mode %in% c("exact","approx1","approx2")) + stop("ls_mode in {'exact','approx1','approx2'}") + + K <- ncol(Ms[[1]]) + if (is.character(ref)) #ref=="mean" + m_sum = Ms[[1]] + L <- length(Ms) + for (i in ifelse(is.character(ref),2,1):L) + { + m_ref = if (is.character(ref)) m_sum / (i-1) else ref + m = Ms[[i]] #shorthand + + if (ls_mode == "exact") + { + #distances[i,j] = distance between m column i and ref column j + distances = apply( m_ref, 2, function(col) ( sqrt(colSums((m-col)^2)) ) ) + assignment = .hungarianAlgorithm(distances) + col <- m[,assignment] + if (is.list(Ms)) Ms[[i]] <- col else Ms[,,i] <- col + } + else + { + # Greedy matching: + # approx1: li[[i]][,j] is assigned to m[,k] minimizing dist(li[[i]][,j],m[,k']) + # approx2: m[,j] is assigned to li[[i]][,k] minimizing dist(m[,j],li[[i]][,k']) + available_indices = 1:K + for (j in 1:K) + { + distances = + if (ls_mode == "approx1") + { + apply(as.matrix(m[,available_indices]), 2, + function(col) ( sqrt(sum((col - m_ref[,j])^2)) ) ) + } + else #approx2 + { + apply(as.matrix(m_ref[,available_indices]), 2, + function(col) ( sqrt(sum((col - m[,j])^2)) ) ) + } + indMin = which.min(distances) + if (ls_mode == "approx1") + { + col <- m[ , available_indices[indMin] ] + if (is.list(Ms)) Ms[[i]][,j] <- col else Ms[,j,i] <- col + } + else #approx2 + { + col <- available_indices[indMin] + if (is.list(Ms)) Ms[[i]][,col] <- m[,j] else Ms[,col,i] <- m[,j] + } + available_indices = available_indices[-indMin] + } + } + + # Update current sum with "label-switched" li[[i]] + if (is.character(ref)) #ref=="mean" + m_sum = m_sum + Ms[[i]] + } + Ms +} diff --git a/pkg/man/morpheus-package.Rd b/pkg/man/morpheus-package.Rd new file mode 100644 index 0000000..e8181aa --- /dev/null +++ b/pkg/man/morpheus-package.Rd @@ -0,0 +1,47 @@ +\name{morpheus-package} +\alias{morpheus-package} +\alias{morpheus} +\docType{package} + +\title{ + \packageTitle{morpheus} +} + +\description{ + \packageDescription{morpheus} +} + +\details{ + The package devtools should be useful in development stage, since we rely on testthat for + unit tests, and roxygen2 for documentation. knitr is used to generate the package vignette. + Concerning the other suggested packages: + \itemize{ + \item{tensor is used for comparing to some reference functions initially coded in R; + it should not be required in further package versions;} + \item{jointDiag allows to solve a joint diagonalization problem, providing a more + robust solution compared to a single diagonalization;} + \item{parallel (generally) permits to run the bootstrap method faster.} + } + + The three main functions are located in R/main.R: + \itemize{ + \item{getParamsDirs_ref: reference method to estimate parameters directions;} + \item{getParamsDirs: method of choice to estimate parameters directions, using a + spectral decomposition of inputs/outputs;} + \item{getBootstrapParams: run getParamsDirs on B bootstrap replicates.} + } +} + +\author{ + \packageAuthor{morpheus} + + Maintainer: \packageMaintainer{morpheus} +} + +%\references{ +% TODO: Literature or other references for background information +%} + +%\examples{ +% TODO: simple examples of the most important functions +%} diff --git a/pkg/src/functions.c b/pkg/src/functions.c new file mode 100644 index 0000000..0d632d3 --- /dev/null +++ b/pkg/src/functions.c @@ -0,0 +1,55 @@ +#include + +//index matrix (by columns) +int mi(int i, int j, int d1, int d2) +{ + return j*d1+i; +} + +//index 3-tensor (by columns, matrices ordered by last dim) +int ti(int i, int j, int k, int d1, int d2, int d3) +{ + return k*d1*d2 + j*d1 + i; +} + +// Emprical cross-moment of order 2 between X size nxd and Y size n +void Moments_M2(double* X, double* Y, int* pn, int* pd, double* M2) +{ + int n=*pn, d=*pd; + //double* M2 = (double*)calloc(d*d,sizeof(double)); + + // M2 = E[Y*X^*2] - E[Y*e^*2] = E[Y (X^*2 - I)] + for (int j=0; j +#include + +#define HUNGARIAN_NOT_ASSIGNED 0 +#define HUNGARIAN_ASSIGNED 1 + +#define HUNGARIAN_MODE_MINIMIZE_COST 0 +#define HUNGARIAN_MODE_MAXIMIZE_UTIL 1 + +#define INF (1e32) +#define verbose (0) + +typedef struct { + int num_rows; + int num_cols; + double** cost; + int** assignment; +} hungarian_problem_t; + +int hungarian_imax(int a, int b) { + return (anum_rows = rows; + p->num_cols = cols; + + p->cost = (double**)calloc(rows,sizeof(double*)); + p->assignment = (int**)calloc(rows,sizeof(int*)); + + for(i=0; inum_rows; i++) { + p->cost[i] = (double*)calloc(cols,sizeof(double)); + p->assignment[i] = (int*)calloc(cols,sizeof(int)); + for(j=0; jnum_cols; j++) { + p->cost[i][j] = (i < org_rows && j < org_cols) ? cost_matrix[i][j] : 0.; + p->assignment[i][j] = 0; + + if (max_cost < p->cost[i][j]) + max_cost = p->cost[i][j]; + } + } + + if (mode == HUNGARIAN_MODE_MAXIMIZE_UTIL) { + for(i=0; inum_rows; i++) { + for(j=0; jnum_cols; j++) + p->cost[i][j] = max_cost - p->cost[i][j]; + } + } + else if (mode == HUNGARIAN_MODE_MINIMIZE_COST) { + // nothing to do + } +// else +// fprintf(stderr,"%s: unknown mode. Mode was set to \ +// HUNGARIAN_MODE_MINIMIZE_COST !\n", __FUNCTION__); + + return rows; +} + +/** Free the memory allocated by init. **/ +void hungarian_free(hungarian_problem_t* p) { + int i; + for(i=0; inum_rows; i++) { + free(p->cost[i]); + free(p->assignment[i]); + } + free(p->cost); + free(p->assignment); + p->cost = NULL; + p->assignment = NULL; +} + +/** This method computes the optimal assignment. **/ +void hungarian_solve(hungarian_problem_t* p) +{ + int i, j, m, n, k, l, t, q, unmatched; + double cost, s; + int* col_mate; + int* row_mate; + int* parent_row; + int* unchosen_row; + double* row_dec; + double* col_inc; + double* slack; + int* slack_row; + + cost = 0.; + m =p->num_rows; + n =p->num_cols; + + col_mate = (int*)calloc(p->num_rows,sizeof(int)); + unchosen_row = (int*)calloc(p->num_rows,sizeof(int)); + row_dec = (double*)calloc(p->num_rows,sizeof(double)); + slack_row = (int*)calloc(p->num_rows,sizeof(int)); + + row_mate = (int*)calloc(p->num_cols,sizeof(int)); + parent_row = (int*)calloc(p->num_cols,sizeof(int)); + col_inc = (double*)calloc(p->num_cols,sizeof(double)); + slack = (double*)calloc(p->num_cols,sizeof(double)); + + for (i=0; inum_rows; i++) { + col_mate[i]=0; + unchosen_row[i]=0; + row_dec[i]=0.; + slack_row[i]=0; + } + for (j=0; jnum_cols; j++) { + row_mate[j]=0; + parent_row[j] = 0; + col_inc[j]=0.; + slack[j]=0.; + } + + for (i=0; inum_rows; ++i) + for (j=0; jnum_cols; ++j) + p->assignment[i][j]=HUNGARIAN_NOT_ASSIGNED; + + // Begin subtract column minima in order to start with lots of zeroes 12 + for (l=0; lcost[0][l]; + for (k=1; kcost[k][l]cost[k][l]; + cost+=s; + if (s!=0.) + for (k=0; kcost[k][l]-=s; + } + // End subtract column minima in order to start with lots of zeroes 12 + + // Begin initial state 16 + t=0; + for (l=0; lcost[k][0]; + for (l=1; lcost[k][l]cost[k][l]; + row_dec[k]=s; + for (l=0; lcost[k][l] && row_mate[l]<0) + { + col_mate[k]=l; + row_mate[l]=k; + goto row_done; + } + col_mate[k]= -1; + unchosen_row[t++]=k; +row_done: + ; + } + // End initial state 16 + + // Begin Hungarian algorithm 18 + if (t==0) + goto done; + unmatched=t; + while (1) + { + q=0; + while (1) + { + while (q 0.) + { + double del=p->cost[k][l]-s+col_inc[l]; + if (del0. && slack[l]0.) + { + slack[l]-=s; + if (slack[l]==0.) + { + // Begin look at a new zero 22 + k=slack_row[l]; + if (row_mate[l]<0) + { + for (j=l+1; jcost[k][l]cost[k][l]!=row_dec[k]-col_inc[l]) + exit(0); + } + k=0; + for (l=0; lm) + exit(0); + // End doublecheck the solution 23 +*/ // End Hungarian algorithm 18 + + for (i=0; iassignment[i][col_mate[i]]=HUNGARIAN_ASSIGNED; + /*TRACE("%d - %d\n", i, col_mate[i]);*/ + } + for (k=0; kcost[k][l]-row_dec[k]+col_inc[l]);*/ + p->cost[k][l]=p->cost[k][l]-row_dec[k]+col_inc[l]; + } + /*TRACE("\n");*/ + } + for (i=0; i