--- /dev/null
+*.pdf filter=fat
+*.zip filter=fat
+*.tar.xz filter=fat
+*.png filter=fat
--- /dev/null
+[rsync]
+remote = gitfat@auder.net:~/files/valse
-#files generated by initialize.sh
-/.gitfat
-/.gitattributes
-
#ignore temporary files
*~
*.swp
*.Rproj*
.Rprofile
.Rbuildignore
+Rprof.out
#ignore R CMD build/check genrated files
/*.Rcheck/
*.o
*.so
*.exe
-
-#misc
-Rprof.out
-*.zip
-.Rproj.user
# VAriable seLection with mixtureS of modEls
-This code is the applied part of the PhD thesis of [Benjamin Gohehry](http://www.math.u-psud.fr/~goehry/).
+With [Emilie Devijver](http://ama.liglab.fr/~devijver/) and [Benjamin Gohehry](https://fr.linkedin.com/in/benjamin-goehry-275276124).
+
+Re-writing from [a similar project](https://github.com/yagu0/select) in MATLAB, which
+corresponds to applied parts of the PhD thesis of both co-authors.
+
+This git repo uses [git-fat](https://github.com/jedbrown/git-fat) to store binary objects. "git fat init && git fat pull" will get them.
## Description
-Apply git-fat patch
+n = 100; m = 70; p = 5
+X = matrix(runif(n*p, -10, 10), nrow=n)
+Y = matrix(runif(n*m, -5, 15), nrow=n)
-Continue debug: (message "cannot allocate vector 40gb" ?!
- X = matrix(runif(n*p, -10, 10), nrow=n)
- Y = matrix(runif(n*m, -5, 15), nrow=n)
- V1 = valse::valse(X, Y, fast=FALSE) #OK
- V2 = valse::valse(X, Y, fast=FALSE) #not OK... almost
+V1 = valse::valse(X, Y, fast=FALSE)
+Error in while (!pi2AllPositive) { :
+ missing value where TRUE/FALSE needed
+
+V2 = valse::valse(X, Y, fast=TRUE)
+list()
+Error in out[[ind_uniq[l]]] :
+ attempt to select less than one element in get1index
+
+==> Error identified: line 61 in initSmallEM.R, division by 0
+It occurs also for smallers values of n and m, e.g.: n = 20; m = 20; p = 3
#!/bin/sh
-./.git-fat/git-fat pull
-./.git-fat/git-fat push
-git submodule update --merge
+git-fat push
Package: valse
Title: Variable Selection With Mixture Of Models
-Date: 2016-12-01
+Date: 2020-01-11
Version: 0.1-0
Description: Two methods are implemented to cluster data with finite mixture
regression models. Those procedures deal with high-dimensional covariates and
(slope heuristic, BIC or AIC). Details of the procedure are provided in 'Model-
based clustering for high-dimensional data. Application to functional data' by
Emilie Devijver, published in Advances in Data Analysis and Clustering (2016).
-Author: Benjamin Auder <Benjamin.Auder@math.u-psud.fr> [aut,cre],
+Author: Benjamin Auder <benjamin.auder@universite-paris-saclay.fr> [aut,cre],
Emilie Devijver <Emilie.Devijver@kuleuven.be> [aut],
Benjamin Goehry <Benjamin.Goehry@math.u-psud.fr> [aut]
-Maintainer: Benjamin Auder <Benjamin.Auder@math.u-psud.fr>
+Maintainer: Benjamin Auder <benjamin.auder@universite-paris-saclay.fr>
Depends:
- R (>= 3.0.0)
+ R (>= 3.5.0)
Imports:
MASS,
parallel
testthat
URL: http://git.auder.net/?p=valse.git
License: MIT + file LICENSE
-RoxygenNote: 5.0.1
+RoxygenNote: 6.0.1
Collate:
'plot_valse.R'
'main.R'
-Copyright (c)
- 2014-2017, Benjamin Auder
- 2014-2017, Emilie Devijver
- 2016-2017, Benjamin Goehry
-
-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.
+YEAR: 2014-2020
+COPYRIGHT HOLDER: Benjamin Auder, Emilie Devijver, Benjamin Goehry
-#' EMGLLF
+#' EMGLLF
#'
#' Description de EMGLLF
#'
#' affec : ...
#'
#' @export
-EMGLLF <- function(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda,
+EMGLLF <- function(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda,
X, Y, eps, fast)
{
if (!fast)
{
# Function in R
- return(.EMGLLF_R(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda,
+ return(.EMGLLF_R(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda,
X, Y, eps))
}
p <- ncol(X) #nombre de covariables
m <- ncol(Y) #taille de Y (multivarié)
k <- length(piInit) #nombre de composantes dans le mélange
- .Call("EMGLLF", phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda,
- X, Y, eps, phi = double(p * m * k), rho = double(m * m * k), pi = double(k),
- llh = double(1), S = double(p * m * k), affec = integer(n), n, p, m, k,
+ .Call("EMGLLF", phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda,
+ X, Y, eps, phi = double(p * m * k), rho = double(m * m * k), pi = double(k),
+ llh = double(1), S = double(p * m * k), affec = integer(n), n, p, m, k,
PACKAGE = "valse")
}
# R version - slow but easy to read
-.EMGLLF_R <- function(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda,
+.EMGLLF_R <- function(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda,
X, Y, eps)
{
# Matrix dimensions
p <- ncol(X) #nombre de covariables
m <- ncol(Y) #taille de Y (multivarié)
k <- length(Pi) #nombre de composantes dans le mélange
- .Call("EMGrank", Pi, Rho, mini, maxi, X, Y, eps, rank, phi = double(p * m * k),
+ .Call("EMGrank", Pi, Rho, mini, maxi, X, Y, eps, as.integer(rank), phi = double(p * m * k),
LLF = double(1), n, p, m, k, PACKAGE = "valse")
}
-#' computeGridLambda
+#' computeGridLambda
#'
#' Construct the data-driven grid for the regularization parameters used for the Lasso estimator
#'
#' @return the grid of regularization parameters
#'
#' @export
-computeGridLambda <- function(phiInit, rhoInit, piInit, gamInit, X, Y, gamma, mini,
+computeGridLambda <- function(phiInit, rhoInit, piInit, gamInit, X, Y, gamma, mini,
maxi, eps, fast)
{
n <- nrow(X)
m <- ncol(Y)
k <- length(piInit)
- list_EMG <- EMGLLF(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda = 0,
+ list_EMG <- EMGLLF(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda = 0,
X, Y, eps, fast)
grid <- array(0, dim = c(p, m, k))
-#' constructionModelesLassoMLE
+#' constructionModelesLassoMLE
#'
#' Construct a collection of models with the Lasso-MLE procedure.
-#'
+#'
#' @param phiInit an initialization for phi, get by initSmallEM.R
#' @param rhoInit an initialization for rho, get by initSmallEM.R
#' @param piInit an initialization for pi, get by initSmallEM.R
#' @param ncores Number of cores, by default = 3
#' @param fast TRUE to use compiled C code, FALSE for R code only
#' @param verbose TRUE to show some execution traces
-#'
+#'
#' @return a list with several models, defined by phi, rho, pi, llh
#'
#' @export
-constructionModelesLassoMLE <- function(phiInit, rhoInit, piInit, gamInit, mini,
+constructionModelesLassoMLE <- function(phiInit, rhoInit, piInit, gamInit, mini,
maxi, gamma, X, Y, eps, S, ncores = 3, fast, verbose)
{
if (ncores > 1)
{
cl <- parallel::makeCluster(ncores, outfile = "")
- parallel::clusterExport(cl, envir = environment(), varlist = c("phiInit",
- "rhoInit", "gamInit", "mini", "maxi", "gamma", "X", "Y", "eps", "S",
+ parallel::clusterExport(cl, envir = environment(), varlist = c("phiInit",
+ "rhoInit", "gamInit", "mini", "maxi", "gamma", "X", "Y", "eps", "S",
"ncores", "fast", "verbose"))
}
# Individual model computation
computeAtLambda <- function(lambda)
{
- if (ncores > 1)
+ if (ncores > 1)
require("valse") #nodes start with an empty environment
- if (verbose)
+ if (verbose)
print(paste("Computations for lambda=", lambda))
n <- nrow(X)
sel.lambda <- S[[lambda]]$selected
# col.sel = which(colSums(sel.lambda)!=0) #if boolean matrix
col.sel <- which(sapply(sel.lambda, length) > 0) #if list of selected vars
- if (length(col.sel) == 0)
+ if (length(col.sel) == 0)
return(NULL)
# lambda == 0 because we compute the EMV: no penalization here
# log(piLambda[r]) + log(detRho[r]) - 0.5 *
# sum((Y[i, ] %*% rhoLambda[, , r] - X[i, ] %*% phiLambda[, , r])^2)
# })
- #
+ #
# #logGam <- logGam - max(logGam) #adjust without changing proportions -> change the LLH
# gam <- exp(logGam)
# norm_fact <- sum(gam)
lapply(1:length(S), computeAtLambda)
}
- if (ncores > 1)
+ if (ncores > 1)
parallel::stopCluster(cl)
out
#' @return a list with several models, defined by phi, rho, pi, llh
#'
#' @export
-constructionModelesLassoRank <- function(S, k, mini, maxi, X, Y, eps, rank.min, rank.max,
+constructionModelesLassoRank <- function(S, k, mini, maxi, X, Y, eps, rank.min, rank.max,
ncores, fast, verbose)
{
n <- nrow(X)
# (rank.max-rank.min)^(k-2) chaque chiffre, et on fait ça (rank.max-rank.min)^2
# fois ... Dans la dernière, on répète chaque chiffre une fois, et on fait ça
# (rank.min-rank.max)^(k-1) fois.
- RankLambda[, r] <- rep(rank.min + rep(0:(deltaRank - 1), deltaRank^(r - 1),
+ RankLambda[, r] <- rep(rank.min + rep(0:(deltaRank - 1), deltaRank^(r - 1),
each = deltaRank^(k - r)), each = L)
}
RankLambda[, k + 1] <- rep(1:L, times = Size)
if (ncores > 1)
{
cl <- parallel::makeCluster(ncores, outfile = "")
- parallel::clusterExport(cl, envir = environment(), varlist = c("A1", "Size",
- "Pi", "Rho", "mini", "maxi", "X", "Y", "eps", "Rank", "m", "phi", "ncores",
+ parallel::clusterExport(cl, envir = environment(), varlist = c("A1", "Size",
+ "Pi", "Rho", "mini", "maxi", "X", "Y", "eps", "Rank", "m", "phi", "ncores",
"verbose"))
}
{
lambdaIndex <- RankLambda[index, k + 1]
rankIndex <- RankLambda[index, 1:k]
- if (ncores > 1)
+ if (ncores > 1)
require("valse") #workers start with an empty environment
# 'relevant' will be the set of relevant columns
phi <- array(0, dim = c(p, m, k))
if (length(relevant) > 0)
{
- res <- EMGrank(S[[lambdaIndex]]$Pi, S[[lambdaIndex]]$Rho, mini, maxi,
+ res <- EMGrank(S[[lambdaIndex]]$Pi, S[[lambdaIndex]]$Rho, mini, maxi,
X[, relevant], Y, eps, rankIndex, fast)
llh <- c(res$LLF, sum(rankIndex * (length(relevant) - rankIndex + m)))
phi[relevant, , ] <- res$phi
lapply(seq_len(length(S) * Size), computeAtLambda)
}
- if (ncores > 1)
+ if (ncores > 1)
parallel::stopCluster(cl)
out
-#' generateXY
+#' generateXY
#'
#' Generate a sample of (X,Y) of size n
#'
class <- c(class, rep(i, sizePop[i]))
newBlockX <- MASS::mvrnorm(sizePop[i], meanX, covX)
X <- rbind(X, newBlockX)
- Y <- rbind(Y, t(apply(newBlockX, 1, function(row) MASS::mvrnorm(1, row %*%
+ Y <- rbind(Y, t(apply(newBlockX, 1, function(row) MASS::mvrnorm(1, row %*%
β[, , i], covY[, , i]))))
}
-#' initialization of the EM algorithm
+#' initialization of the EM algorithm
#'
#' @param k number of components
#' @param X matrix of covariates (of size n*p)
Z <- Zinit1[, repet]
Z_indice <- seq_len(n)[Z == r] #renvoit les indices où Z==r
if (length(Z_indice) == 1) {
- betaInit1[, , r, repet] <- MASS::ginv(crossprod(t(X[Z_indice, ]))) %*%
+ betaInit1[, , r, repet] <- MASS::ginv(crossprod(t(X[Z_indice, ]))) %*%
crossprod(t(X[Z_indice, ]), Y[Z_indice, ])
} else {
- betaInit1[, , r, repet] <- MASS::ginv(crossprod(X[Z_indice, ])) %*%
+ betaInit1[, , r, repet] <- MASS::ginv(crossprod(X[Z_indice, ])) %*%
crossprod(X[Z_indice, ], Y[Z_indice, ])
}
sigmaInit1[, , r, repet] <- diag(m)
{
dotProduct <- tcrossprod(Y[i, ] %*% rhoInit1[, , r, repet]
- X[i, ] %*% phiInit1[, , r, repet])
- Gam[i, r] <- piInit1[repet, r] *
+ Gam[i, r] <- piInit1[repet, r] *
det(rhoInit1[, , r, repet]) * exp(-0.5 * dotProduct)
}
sumGamI <- sum(Gam[i, ])
+ # TODO: next line is a division by zero if dotProduct is big
gamInit1[i, , repet] <- Gam[i, ]/sumGamI
}
-#' valse
+#' valse
#'
#' Main function
#'
#' @examples
#' #TODO: a few examples
#' @export
-valse <- function(X, Y, procedure = "LassoMLE", selecMod = "DDSE", gamma = 1, mini = 10,
- maxi = 50, eps = 1e-04, kmin = 2, kmax = 3, rank.min = 1, rank.max = 5, ncores_outer = 1,
+valse <- function(X, Y, procedure = "LassoMLE", selecMod = "DDSE", gamma = 1, mini = 10,
+ maxi = 50, eps = 1e-04, kmin = 2, kmax = 3, rank.min = 1, rank.max = 5, ncores_outer = 1,
ncores_inner = 1, thresh = 1e-08, grid_lambda = numeric(0), size_coll_mod = 10,
fast = TRUE, verbose = FALSE, plot = TRUE)
{
p <- ncol(X)
m <- ncol(Y)
- if (verbose)
+ if (verbose)
print("main loop: over all k and all lambda")
if (ncores_outer > 1) {
cl <- parallel::makeCluster(ncores_outer, outfile = "")
- parallel::clusterExport(cl = cl, envir = environment(), varlist = c("X",
- "Y", "procedure", "selecMod", "gamma", "mini", "maxi", "eps", "kmin",
- "kmax", "rank.min", "rank.max", "ncores_outer", "ncores_inner", "thresh",
+ parallel::clusterExport(cl = cl, envir = environment(), varlist = c("X",
+ "Y", "procedure", "selecMod", "gamma", "mini", "maxi", "eps", "kmin",
+ "kmax", "rank.min", "rank.max", "ncores_outer", "ncores_inner", "thresh",
"size_coll_mod", "verbose", "p", "m"))
}
# Compute models with k components
computeModels <- function(k)
{
- if (ncores_outer > 1)
+ if (ncores_outer > 1)
require("valse") #nodes start with an empty environment
- if (verbose)
+ if (verbose)
print(paste("Parameters initialization for k =", k))
# smallEM initializes parameters by k-means and regression model in each
# component, doing this 20 times, and keeping the values maximizing the
P <- initSmallEM(k, X, Y, fast)
if (length(grid_lambda) == 0)
{
- grid_lambda <- computeGridLambda(P$phiInit, P$rhoInit, P$piInit, P$gamInit,
+ grid_lambda <- computeGridLambda(P$phiInit, P$rhoInit, P$piInit, P$gamInit,
X, Y, gamma, mini, maxi, eps, fast)
}
- if (length(grid_lambda) > size_coll_mod)
+ if (length(grid_lambda) > size_coll_mod)
grid_lambda <- grid_lambda[seq(1, length(grid_lambda), length.out = size_coll_mod)]
- if (verbose)
+ if (verbose)
print("Compute relevant parameters")
# select variables according to each regularization parameter from the grid:
# S$selected corresponding to selected variables
- S <- selectVariables(P$phiInit, P$rhoInit, P$piInit, P$gamInit, mini, maxi,
+ S <- selectVariables(P$phiInit, P$rhoInit, P$piInit, P$gamInit, mini, maxi,
gamma, grid_lambda, X, Y, thresh, eps, ncores_inner, fast)
if (procedure == "LassoMLE") {
- if (verbose)
+ if (verbose)
print("run the procedure Lasso-MLE")
# compute parameter estimations, with the Maximum Likelihood Estimator,
# restricted on selected variables.
- models <- constructionModelesLassoMLE(P$phiInit, P$rhoInit, P$piInit,
+ models <- constructionModelesLassoMLE(P$phiInit, P$rhoInit, P$piInit,
P$gamInit, mini, maxi, gamma, X, Y, eps, S, ncores_inner, fast, verbose)
} else {
- if (verbose)
+ if (verbose)
print("run the procedure Lasso-Rank")
# compute parameter estimations, with the Low Rank Estimator, restricted on
# selected variables.
- models <- constructionModelesLassoRank(S, k, mini, maxi, X, Y, eps, rank.min,
+ models <- constructionModelesLassoRank(S, k, mini, maxi, X, Y, eps, rank.min,
rank.max, ncores_inner, fast, verbose)
}
# warning! Some models are NULL after running selectVariables
} else {
lapply(kmin:kmax, computeModels)
}
- if (ncores_outer > 1)
+ if (ncores_outer > 1)
parallel::stopCluster(cl)
if (!requireNamespace("capushe", quietly = TRUE))
# For a collection of models (same k, several lambda):
LLH <- sapply(models, function(model) model$llh[1])
k <- length(models[[1]]$pi)
- sumPen <- sapply(models, function(model) k * (dim(model$rho)[1] + sum(model$phi[,
+ sumPen <- sapply(models, function(model) k * (dim(model$rho)[1] + sum(model$phi[,
, 1] != 0) + 1) - 1)
- data.frame(model = paste(i, ".", seq_along(models), sep = ""), pen = sumPen/n,
+ data.frame(model = paste(i, ".", seq_along(models), sep = ""), pen = sumPen/n,
complexity = sumPen, contrast = -LLH)
}))
tableauRecap <- tableauRecap[which(tableauRecap[, 4] != Inf), ]
if (verbose == TRUE)
print(tableauRecap)
modSel <- capushe::capushe(tableauRecap, n)
- indModSel <- if (selecMod == "DDSE")
+ indModSel <- if (selecMod == "DDSE")
{
as.numeric(modSel@DDSE@model)
- } else if (selecMod == "Djump")
+ } else if (selecMod == "Djump")
{
as.numeric(modSel@Djump@model)
- } else if (selecMod == "BIC")
+ } else if (selecMod == "BIC")
{
modSel@BIC_capushe$model
- } else if (selecMod == "AIC")
+ } else if (selecMod == "AIC")
{
modSel@AIC_capushe$model
}
-#' Plot
+#' Plot
#'
#' It is a function which plots relevant parameters
#'
for (r in 1:K)
{
Melt <- melt(t((model$phi[, , r])))
- gReg[[r]] <- ggplot(data = Melt, aes(x = Var1, y = Var2, fill = value)) +
- geom_tile() + scale_fill_gradient2(low = "blue", high = "red", mid = "white",
+ gReg[[r]] <- ggplot(data = Melt, aes(x = Var1, y = Var2, fill = value)) +
+ geom_tile() + scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, space = "Lab") + ggtitle(paste("Regression matrices in cluster", r))
}
print(gReg)
Melt <- melt(t(model$phi[, , k1] - model$phi[, , k2]))
gDiff <- ggplot(data = Melt, aes(x = Var1, y = Var2, fill = value))
+ geom_tile()
- + scale_fill_gradient2(low = "blue", high = "red", mid = "white", midpoint = 0,
+ + scale_fill_gradient2(low = "blue", high = "red", mid = "white", midpoint = 0,
space = "Lab")
- + ggtitle(paste("Difference between regression matrices in cluster",
+ + ggtitle(paste("Difference between regression matrices in cluster",
k1, "and", k2))
print(gDiff)
}
matCov[, r] <- diag(model$rho[, , r])
MeltCov <- melt(matCov)
gCov <- ggplot(data = MeltCov, aes(x = Var1, y = Var2, fill = value)) + geom_tile()
- + scale_fill_gradient2(low = "blue", high = "red", mid = "white", midpoint = 0,
+ + scale_fill_gradient2(low = "blue", high = "red", mid = "white", midpoint = 0,
space = "Lab")
+ ggtitle("Covariance matrices")
print(gCov)
-#' selectVariables
+#' selectVariables
#'
#' It is a function which construct, for a given lambda, the sets of relevant variables.
#'
#'
#' @export
#'
-selectVariables <- function(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma,
+selectVariables <- function(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma,
glambda, X, Y, thresh = 1e-08, eps, ncores = 3, fast)
{
if (ncores > 1) {
cl <- parallel::makeCluster(ncores, outfile = "")
- parallel::clusterExport(cl = cl, varlist = c("phiInit", "rhoInit", "gamInit",
+ parallel::clusterExport(cl = cl, varlist = c("phiInit", "rhoInit", "gamInit",
"mini", "maxi", "glambda", "X", "Y", "thresh", "eps"), envir = environment())
}
# Computation for a fixed lambda
computeCoefs <- function(lambda)
{
- params <- EMGLLF(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda,
+ params <- EMGLLF(phiInit, rhoInit, piInit, gamInit, mini, maxi, gamma, lambda,
X, Y, eps, fast)
p <- ncol(X)
} else {
lapply(glambda, computeCoefs)
}
- if (ncores > 1)
+ if (ncores > 1)
parallel::stopCluster(cl)
-
+
print(out)
# Suppress models which are computed twice En fait, ca ca fait la comparaison de
# tous les parametres On veut juste supprimer ceux qui ont les memes variables
-# ...
+# Compute the determinant of a matrix, which can be 1x1 (scalar)
gdet <- function(M)
{
- if (is.matrix(M))
- return (det(M))
- return (M[1]) #numeric, double
+ ifelse(is.matrix(M), det(M), M[1])
}
SEXP maxi_,
SEXP X_,
SEXP Y_,
- SEXP tau_,
+ SEXP eps_,
SEXP rank_
) {
// Get matrices dimensions
// get scalar parameters
int mini = INTEGER_VALUE(mini_);
int maxi = INTEGER_VALUE(maxi_);
- double tau = NUMERIC_VALUE(tau_);
+ double eps = NUMERIC_VALUE(eps_);
// Get pointers from SEXP arrays ; WARNING: by columns !
double* Pi = REAL(Pi_);
// Call to EMGrank //
/////////////////////
- EMGrank_core(Pi, Rho, mini, maxi, X, Y, tau, rank,
+ EMGrank_core(Pi, Rho, mini, maxi, X, Y, eps, rank,
pPhi,pLLF,
n,p,m,k);
free(X2);
free(Y2);
free(sqNorm2);
-}\f
+}
--- /dev/null
+#$# git-fat 6ea97df8b04ef5722cc85b9f6c467cdee85155ff 43421
rm -f *.o ../pkg/src/sources/*.o
cclean: clean
- rm -f *.so $(TESTS)
+ rm -f *.so ../pkg/src/*.so $(TESTS)
.PHONY: all clean cclean
source("helper.R")
library(valse)
-generateRunSaveTest_EMGLLF = function(n=200, p=15, m=10, k=3, mini=5, maxi=10,
- gamma=1., lambda=0.5, eps=1e-6)
+generateRunSaveTest_EMGLLF = function(n=200, p=15, m=10, k=3, mini=5, maxi=10, gamma=1., lambda=0.5, eps=1e-6)
{
- testFolder = "data/"
+ testFolder = "./data/"
dir.create(testFolder, showWarnings=FALSE, mode="0755")
params = basicInitParameters(n, p, m, k)
source("helper.R")
library(valse)
-generateRunSaveTest_EMGrank = function(n=200, p=15, m=10, k=3, mini=5, maxi=10, gamma=1.0,
- rank = c(1,2,4))
+generateRunSaveTest_EMGrank = function(n=200, p=15, m=10, k=3, mini=5, maxi=10, gamma=1.0, rank = c(1,2,4))
{
eps = 1e-6
- pi = rep(1.0/k, k)
- rho = array(dim=c(m,m,k))
+ Pi = rep(1.0/k, k)
+ Rho = array(dim=c(m,m,k))
for(i in 1:k)
- rho[,,i] = diag(1,m)
+ Rho[,,i] = diag(1,m)
xy = generateXYdefault(n, p, m, k)
- testFolder = "../data/"
+ testFolder = "./data/"
dir.create(testFolder, showWarnings=FALSE, mode="0755")
#save inputs
- write.table(as.double(rho), paste(testFolder,"rho",sep=""),
+ write.table(as.double(Pi), paste(testFolder,"Pi",sep=""),
row.names=F, col.names=F)
- write.table(as.double(pi), paste(testFolder,"pi",sep=""),
+ write.table(as.double(Rho), paste(testFolder,"Rho",sep=""),
row.names=F, col.names=F)
write.table(as.integer(mini), paste(testFolder,"mini",sep=""),
row.names=F, col.names=F)
write.table(as.integer(c(n,p,m,k)), paste(testFolder,"dimensions",sep=""),
row.names=F, col.names=F)
- res = valse::EMGrank(pi,rho,mini,maxi,xy$X,xy$Y,eps,rank,fast=FALSE)
+ res = valse::EMGrank(Pi,Rho,mini,maxi,xy$X,xy$Y,eps,rank,fast=FALSE)
#save output
write.table(as.double(res$phi),paste(testFolder,"phi",sep=""),row.names=F,col.names=F)
////////////
// INPUTS //
- Real* rho = readArray_real("rho");
- Real* pi = readArray_real("pi");
+ Real* Rho = readArray_real("Rho");
+ Real* Pi = readArray_real("Pi");
int mini = read_int("mini");
int maxi = read_int("maxi");
Real* X = readArray_real("X");
Real* Y = readArray_real("Y");
- Real tau = read_real("tau");
+ Real eps = read_real("eps");
int* rank = readArray_int("rank");
////////////
//////////////////////////
// Main call to EMGrank //
- EMGrank_core(pi,rho,mini,maxi,X,Y,tau,rank,
+ EMGrank_core(Pi,Rho,mini,maxi,X,Y,eps,rank,
phi,LLF,
n,p,m,k);
//////////////////////////
- free(rho);
- free(pi);
+ free(Rho);
+ free(Pi);
free(X);
free(Y);
free(rank);
# NOTE: R typing is really terrible. as.double as.matrix ...and so on; don't remove.
#get inputs
-#TODO
+npmk = as.matrix(read.table(paste(testFolder,"dimensions",sep="")))
+n = npmk[1]; p=npmk[2]; m=npmk[3]; k=npmk[4]
+Pi = as.double(as.matrix(read.table(paste(testFolder,"Pi",sep="")))[,])
+Rho = array(as.double(as.matrix(read.table(paste(testFolder,"Rho",sep="")))), dim=c(m,m,k))
+mini = as.integer(as.matrix(read.table(paste(testFolder,"mini",sep="")))[1])
+maxi = as.integer(as.matrix(read.table(paste(testFolder,"maxi",sep="")))[1])
+X = matrix(as.double(as.matrix(read.table(paste(testFolder,"X",sep="")))), n,p)
+Y = matrix(as.double(as.matrix(read.table(paste(testFolder,"Y",sep="")))), n,m)
+eps = as.double(as.matrix(read.table(paste(testFolder,"eps",sep="")))[1])
+rank = as.double(as.matrix(read.table(paste(testFolder,"rank",sep="")))[,])
#get outputs
-#TODO
+phi = array(as.double(as.matrix(read.table(paste(testFolder,"phi",sep="")))), dim=c(p,m,k))
+LLF = as.double(as.matrix(read.table(paste(testFolder,"LLF",sep="")))[1])
-res = valse::EMGrank(...)
+res = valse::EMGrank(Pi,Rho,mini,maxi,X,Y,eps,rank,fast=TRUE)
#compare outputs
nd=7 #number of digits
-#TODO
+print( all(round(phi,nd) == round(res$phi,nd)) )
+print( all(round(LLF,nd) == round(res$LLF,nd)) )
projects / valse.git / commitdiff