From cbd88fe5729bf206a784238a2637aa60e697fcdc Mon Sep 17 00:00:00 2001
From: Benjamin Auder <benjamin.auder@somewhere>
Date: Tue, 23 Jan 2018 02:08:08 +0100
Subject: [PATCH] First commit

---
 .gitattributes                               |   6 +
 .gitfat                                      |   2 +
 .gitignore                                   |  23 ++
 README.md                                    |  17 +
 TODO                                         |  12 +
 doc/Tensor-Kuleshov2015.pdf                  |   1 +
 patch_Bettina/FLXMRglm.R                     | 139 +++++++
 patch_Bettina/code.R                         |  25 ++
 pkg/DESCRIPTION                              |  34 ++
 pkg/LICENSE                                  |  20 +
 pkg/R/A_NAMESPACE.R                          |  11 +
 pkg/R/computeMu.R                            |  92 +++++
 pkg/R/multiRun.R                             | 125 ++++++
 pkg/R/optimParams.R                          | 332 ++++++++++++++++
 pkg/R/plot.R                                 | 139 +++++++
 pkg/R/sampleIO.R                             |  71 ++++
 pkg/R/utils.R                                | 169 ++++++++
 pkg/man/morpheus-package.Rd                  |  47 +++
 pkg/src/functions.c                          |  55 +++
 pkg/src/hungarian.c                          | 398 +++++++++++++++++++
 pkg/tests/testthat.R                         |   5 +
 pkg/tests/testthat/test-Moments_M2.R         |  43 ++
 pkg/tests/testthat/test-Moments_M3.R         |  51 +++
 pkg/tests/testthat/test-alignMatrices.R      |  69 ++++
 pkg/tests/testthat/test-computeMu.R          |  35 ++
 pkg/tests/testthat/test-hungarianAlgorithm.R |  24 ++
 pkg/tests/testthat/test-jointDiag.R          |  50 +++
 pkg/tests/testthat/test-optimParams.R        |  87 ++++
 to-cran.sh                                   |  22 +
 29 files changed, 2104 insertions(+)
 create mode 100644 .gitattributes
 create mode 100644 .gitfat
 create mode 100644 .gitignore
 create mode 100644 README.md
 create mode 100644 TODO
 create mode 100644 doc/Tensor-Kuleshov2015.pdf
 create mode 100644 patch_Bettina/FLXMRglm.R
 create mode 100644 patch_Bettina/code.R
 create mode 100644 pkg/DESCRIPTION
 create mode 100644 pkg/LICENSE
 create mode 100644 pkg/R/A_NAMESPACE.R
 create mode 100644 pkg/R/computeMu.R
 create mode 100644 pkg/R/multiRun.R
 create mode 100644 pkg/R/optimParams.R
 create mode 100644 pkg/R/plot.R
 create mode 100644 pkg/R/sampleIO.R
 create mode 100644 pkg/R/utils.R
 create mode 100644 pkg/man/morpheus-package.Rd
 create mode 100644 pkg/src/functions.c
 create mode 100644 pkg/src/hungarian.c
 create mode 100644 pkg/tests/testthat.R
 create mode 100644 pkg/tests/testthat/test-Moments_M2.R
 create mode 100644 pkg/tests/testthat/test-Moments_M3.R
 create mode 100644 pkg/tests/testthat/test-alignMatrices.R
 create mode 100644 pkg/tests/testthat/test-computeMu.R
 create mode 100644 pkg/tests/testthat/test-hungarianAlgorithm.R
 create mode 100644 pkg/tests/testthat/test-jointDiag.R
 create mode 100644 pkg/tests/testthat/test-optimParams.R
 create mode 100755 to-cran.sh

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 <Benjamin.Auder@math.u-psud.fr> [aut,cre],
+    Mor-Absa Loum <Mor-Absa.Loum@math.u-psud.fr> [aut]
+Maintainer: Benjamin Auder <Benjamin.Auder@math.u-psud.fr>
+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 <stdlib.h>
+
+//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<d; j++)
+	{
+		for (int i=0; i<n; i++)
+		{
+			M2[mi(j,j,d,d)] -= Y[i] / n;
+			for (int k=0; k<d; k++)
+				M2[mi(j,k,d,d)] += Y[i] * X[mi(i,j,n,d)]*X[mi(i,k,n,d)] / n;
+		}
+	}
+}
+
+// Emprical cross-moment of order 3 between X size nxd and Y size n
+void Moments_M3(double* X, double* Y, int* pn, int* pd, double* M3)
+{
+	int n=*pn, d=*pd;
+	//double* M3 = (double*)calloc(d*d*d,sizeof(double));
+
+	// M3 = E[Y*X^*3] - E[Y*e*X*e] - E[Y*e*e*X] - E[Y*X*e*e]
+	for (int j=0; j<d; j++)
+	{
+		for (int k=0; k<d; k++)
+		{
+			for (int i=0; i<n; i++)
+			{
+				double tensor_elt = Y[i]*X[mi(i,k,n,d)] / n;
+				M3[ti(j,k,j,d,d,d)] -= tensor_elt;
+				M3[ti(j,j,k,d,d,d)] -= tensor_elt;
+				M3[ti(k,j,j,d,d,d)] -= tensor_elt;
+				for (int o=0; o<d; o++)
+					M3[ti(j,k,o,d,d,d)] += Y[i] * X[mi(i,j,n,d)]*X[mi(i,k,n,d)]*X[mi(i,o,n,d)] / n;
+			}
+		}
+	}
+}
diff --git a/pkg/src/hungarian.c b/pkg/src/hungarian.c
new file mode 100644
index 0000000..3a3a5ec
--- /dev/null
+++ b/pkg/src/hungarian.c
@@ -0,0 +1,398 @@
+/********************************************************************
+ ********************************************************************
+ **
+ ** libhungarian by Cyrill Stachniss, 2004
+ ** Modified by Benjamin Auder to work on floating point number,
+ ** and to fit into a R package (2017)
+ **
+ ** Solving the Minimum Assignment Problem using the Hungarian Method.
+ **
+ ** ** This file may be freely copied and distributed! **
+ **
+ ** Parts of the used code was originally provided by the
+ ** "Stanford GraphGase", but I made changes to this code.
+ ** As asked by	the copyright node of the "Stanford GraphGase",
+ ** I hereby proclaim that this file are *NOT* part of the
+ ** "Stanford GraphGase" distrubition!
+ **
+ ** This file is distributed in the hope that it will be useful,
+ ** but WITHOUT ANY WARRANTY; without even the implied
+ ** warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ ** PURPOSE.
+ **
+ ********************************************************************
+ ********************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#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 (a<b)?b:a;
+}
+
+/** This method initialize the hungarian_problem structure and init
+ *  the  cost matrices (missing lines or columns are filled with 0).
+ *  It returns the size of the quadratic(!) assignment matrix. **/
+int hungarian_init(hungarian_problem_t* p, double** cost_matrix, int rows, int cols,
+	int mode)
+{
+	int i,j;
+	double max_cost = 0.;
+	int org_cols = cols,
+	    org_rows = rows;
+
+	// is the number of cols	not equal to number of rows ?
+	// if yes, expand with 0-cols / 0-cols
+	rows = hungarian_imax(cols, rows);
+	cols = rows;
+
+	p->num_rows = rows;
+	p->num_cols = cols;
+
+	p->cost = (double**)calloc(rows,sizeof(double*));
+	p->assignment = (int**)calloc(rows,sizeof(int*));
+
+	for(i=0; i<p->num_rows; i++) {
+		p->cost[i] = (double*)calloc(cols,sizeof(double));
+		p->assignment[i] = (int*)calloc(cols,sizeof(int));
+		for(j=0; j<p->num_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; i<p->num_rows; i++) {
+			for(j=0; j<p->num_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; i<p->num_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; i<p->num_rows; i++) {
+		col_mate[i]=0;
+		unchosen_row[i]=0;
+		row_dec[i]=0.;
+		slack_row[i]=0;
+	}
+	for (j=0; j<p->num_cols; j++) {
+		row_mate[j]=0;
+		parent_row[j] = 0;
+		col_inc[j]=0.;
+		slack[j]=0.;
+	}
+
+	for (i=0; i<p->num_rows; ++i)
+		for (j=0; j<p->num_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; l<n; l++)
+	{
+		s=p->cost[0][l];
+		for (k=1; k<m; k++)
+			if (p->cost[k][l]<s)
+				s=p->cost[k][l];
+		cost+=s;
+		if (s!=0.)
+			for (k=0; k<m; k++)
+				p->cost[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; l<n; l++)
+	{
+		row_mate[l]= -1;
+		parent_row[l]= -1;
+		col_inc[l]=0.;
+		slack[l]=INF;
+	}
+	for (k=0; k<m; k++)
+	{
+		s=p->cost[k][0];
+		for (l=1; l<n; l++)
+			if (p->cost[k][l]<s)
+				s=p->cost[k][l];
+		row_dec[k]=s;
+		for (l=0; l<n; l++)
+			if (s==p->cost[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<t)
+			{
+				// Begin explore node q of the forest 19
+				{
+					k=unchosen_row[q];
+					s=row_dec[k];
+					for (l=0; l<n; l++)
+						if (slack[l] > 0.)
+						{
+							double del=p->cost[k][l]-s+col_inc[l];
+							if (del<slack[l])
+							{
+								if (del==0.)
+								{
+									if (row_mate[l]<0)
+										goto breakthru;
+									slack[l]=0.;
+									parent_row[l]=k;
+									unchosen_row[t++]=row_mate[l];
+								}
+								else
+								{
+									slack[l]=del;
+									slack_row[l]=k;
+								}
+							}
+						}
+				}
+				// End explore node q of the forest 19
+				q++;
+			}
+
+			// Begin introduce a new zero into the matrix 21
+			s=INF;
+			for (l=0; l<n; l++)
+				if (slack[l]>0. && slack[l]<s)
+					s=slack[l];
+			for (q=0; q<t; q++)
+				row_dec[unchosen_row[q]]+=s;
+			for (l=0; l<n; l++)
+				if (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; j<n; j++)
+								if (slack[j]==0.)
+									col_inc[j]+=s;
+							goto breakthru;
+						}
+						else
+						{
+							parent_row[l]=k;
+							unchosen_row[t++]=row_mate[l];
+						}
+						// End look at a new zero 22
+					}
+				}
+				else
+					col_inc[l]+=s;
+			// End introduce a new zero into the matrix 21
+		}
+breakthru:
+		// Begin update the matching 20
+		while (1)
+		{
+			j=col_mate[k];
+			col_mate[k]=l;
+			row_mate[l]=k;
+			if (j<0)
+				break;
+			k=parent_row[j];
+			l=j;
+		}
+		// End update the matching 20
+		if (--unmatched==0)
+			goto done;
+		// Begin get ready for another stage 17
+		t=0;
+		for (l=0; l<n; l++)
+		{
+			parent_row[l]= -1;
+			slack[l]=INF;
+		}
+		for (k=0; k<m; k++)
+			if (col_mate[k]<0)
+				unchosen_row[t++]=k;
+		// End get ready for another stage 17
+	}
+done:
+
+/*	// Begin doublecheck the solution 23
+	for (k=0; k<m; k++)
+		for (l=0; l<n; l++)
+			if (p->cost[k][l]<row_dec[k]-col_inc[l])
+				exit(0);
+	for (k=0; k<m; k++)
+	{
+		l=col_mate[k];
+		if (l<0 || p->cost[k][l]!=row_dec[k]-col_inc[l])
+			exit(0);
+	}
+	k=0;
+	for (l=0; l<n; l++)
+		if (col_inc[l])
+			k++;
+	if (k>m)
+		exit(0);
+	// End doublecheck the solution 23
+*/	// End Hungarian algorithm 18
+
+	for (i=0; i<m; ++i)
+	{
+		p->assignment[i][col_mate[i]]=HUNGARIAN_ASSIGNED;
+		/*TRACE("%d - %d\n", i, col_mate[i]);*/
+	}
+	for (k=0; k<m; ++k)
+	{
+		for (l=0; l<n; ++l)
+		{
+			/*TRACE("%d ",p->cost[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<m; i++)
+		cost+=row_dec[i];
+	for (i=0; i<n; i++)
+		cost-=col_inc[i];
+
+	free(slack);
+	free(col_inc);
+	free(parent_row);
+	free(row_mate);
+	free(slack_row);
+	free(row_dec);
+	free(unchosen_row);
+	free(col_mate);
+}
+
+
+/***************
+ * Usage in R :
+ **************/
+
+// Auxiliary function for hungarianAlgorithm below
+double** array_to_matrix(double* m, int rows, int cols)
+{
+  double** r = (double**)calloc(rows,sizeof(double*));
+  for (int i=0; i<rows; i++)
+  {
+    r[i] = (double*)calloc(cols,sizeof(double));
+    for (int j=0; j<cols; j++)
+      r[i][j] = m[i*cols+j];
+  }
+  return r;
+}
+
+//TODO: re-code this algorithm in a more readable way, based on
+//https://www.topcoder.com/community/data-science/data-science-tutorials/\
+//  assignment-problem-and-hungarian-algorithm/
+// Get the optimal assignment, by calling hungarian_solve above; "distances" in columns
+void hungarianAlgorithm(double* distances, int* pn, int* assignment)
+{
+	int n = *pn;
+	double** mat = array_to_matrix(distances, n, n);
+
+  hungarian_problem_t p;
+	hungarian_init(&p, mat, n, n, HUNGARIAN_MODE_MINIMIZE_COST);
+
+	hungarian_solve(&p);
+
+	// Copy the optimal assignment in a R-friendly structure
+	for (int i=0; i < n; i++) {
+		for (int j=0; j < n; j++) {
+			if (p.assignment[i][j])
+				assignment[i] = j+1; //add 1 since we work with R
+		}
+	}
+
+  hungarian_free(&p);
+	for (int i=0; i<n; i++)
+		free(mat[i]);
+	free(mat);
+}
diff --git a/pkg/tests/testthat.R b/pkg/tests/testthat.R
new file mode 100644
index 0000000..539a565
--- /dev/null
+++ b/pkg/tests/testthat.R
@@ -0,0 +1,5 @@
+library(testthat)
+#library(morpheus)
+load_all()
+
+test_check("morpheus")
diff --git a/pkg/tests/testthat/test-Moments_M2.R b/pkg/tests/testthat/test-Moments_M2.R
new file mode 100644
index 0000000..5201104
--- /dev/null
+++ b/pkg/tests/testthat/test-Moments_M2.R
@@ -0,0 +1,43 @@
+context("Moments_M2")
+
+# (Slower, but trusted) R version of Moments_M2
+.Moments_M2_R = function(X,Y)
+{
+	library(tensor)
+	n = nrow(X)
+	d = ncol(X)
+	v = rep(0,d)
+	e = diag(rep(1,d))
+
+	M21 = M2_1 = tensor(v,v)
+	for (i in 1:n)
+		M21 = M21 + Y[i] * tensor(X[i,],X[i,])
+	M21 = (1/n) * M21
+
+	for (j in 1:d)
+	{
+		L = tensor(v,v)
+		for (i in 1:n)
+			L = L + Y[i]*tensor(e[,j],e[,j])
+		L = (1/n) * L
+		M2_1 = M2_1 + L
+	}
+
+	M2 = M21 - M2_1
+	return (M2)
+}
+
+test_that("both versions of Moments_M2 agree on various inputs",
+{
+	for (n in c(20,200))
+	{
+		for (d in c(2,10,20))
+		{
+			X = matrix( runif(n*d,min=-1,max=1), nrow=n )
+			Y = runif(n,min=-1,max=1)
+			M2 = .Moments_M2(X,Y)
+			M2_R = .Moments_M2_R(X,Y)
+			expect_equal(max(abs(M2 - M2_R)), 0)
+		}
+	}
+})
diff --git a/pkg/tests/testthat/test-Moments_M3.R b/pkg/tests/testthat/test-Moments_M3.R
new file mode 100644
index 0000000..502ceb3
--- /dev/null
+++ b/pkg/tests/testthat/test-Moments_M3.R
@@ -0,0 +1,51 @@
+context("Moments_M3")
+
+# (Slower, but trusted) R version of Moments_M3
+.Moments_M3_R = function(X,Y)
+{
+	library(tensor)
+	n = nrow(X)
+	d = ncol(X)
+	v = rep(0,d)
+	e = diag(rep(1,d))
+
+	M31=M3_1=M3_2=M3_3 = tensor(tensor(v,v),v)
+	for (i in 1:n)
+		M31 = M31 + (Y[i]*tensor(tensor(X[i,],X[i,]),X[i,]))
+	M31 = (1/n)*M31
+
+	for(j in 1:d)
+	{
+		l1=l2=l3 = tensor(tensor(v,v),v)
+		for(i in 1:n)
+		{
+			l1 = l1 + Y[i]*tensor(tensor(e[,j],X[i,]),e[,j])
+			l2 = l2 + Y[i]*tensor(tensor(e[,j],e[,j]),X[i,])
+			l3 = l3 + Y[i]*tensor(tensor(X[i,],e[,j]),e[,j])
+		}
+		l1 = (1/n)*l1
+		l2 = (1/n)*l2
+		l3 = (1/n)*l3
+		M3_1 = M3_1+l1
+		M3_2 = M3_2+l2
+		M3_3 = M3_3+l3
+	}
+
+	M3 = M31-(M3_1+M3_2+M3_3)
+	return (M3)
+}
+
+test_that("both versions of Moments_M3 agree on various inputs",
+{
+	for (n in c(20,200))
+	{
+		for (d in c(2,10,20))
+		{
+			X = matrix( runif(n*d,min=-1,max=1), nrow=n )
+			Y = runif(n,min=-1,max=1)
+			M3 = .Moments_M3(X,Y)
+			M3_R = .Moments_M3_R(X,Y)
+			expect_equal(max(abs(M3 - M3_R)), 0)
+		}
+	}
+})
diff --git a/pkg/tests/testthat/test-alignMatrices.R b/pkg/tests/testthat/test-alignMatrices.R
new file mode 100644
index 0000000..59ad1d1
--- /dev/null
+++ b/pkg/tests/testthat/test-alignMatrices.R
@@ -0,0 +1,69 @@
+context("alignMatrices")
+
+# Helper to generate a random series of matrices to align
+.generateMatrices = function(d, K, N, noise)
+{
+	matrices = list( matrix(runif(d*K, min=-1, max=1),ncol=K) ) #reference
+	for (i in 2:N)
+	{
+		matrices[[i]] <- matrices[[1]][,sample(1:K)]
+		if (noise)
+			matrices[[i]] = matrices[[i]] + matrix(rnorm(d*K, sd=0.05), ncol=K)
+	}
+	matrices
+}
+
+test_that("labelSwitchingAlign correctly aligns de-noised parameters",
+{
+	N = 30 #number of matrices
+	d_K_list = list(c(2,2), c(5,3))
+	for (i in 1:2)
+	{
+		d = d_K_list[[i]][1]
+		K = d_K_list[[i]][2]
+
+		# 1] Generate matrix series
+		matrices_permut = .generateMatrices(d,K,N,noise=FALSE)
+
+		# 2] Call align function with mode=approx1
+		matrices_aligned =
+			alignMatrices(matrices_permut[2:N], ref=matrices_permut[[1]], ls_mode="approx1")
+
+		# 3] Check alignment
+		for (j in 2:N)
+			expect_equal(matrices_aligned[[j-1]], matrices_permut[[1]])
+
+		# 2bis] Call align function with mode=approx2
+		matrices_aligned =
+			alignMatrices(matrices_permut[2:N], ref=matrices_permut[[1]], ls_mode="approx2")
+
+		# 3bis] Check alignment
+		for (j in 2:N)
+			expect_equal(matrices_aligned[[j-1]], matrices_permut[[1]])
+	}
+})
+
+test_that("labelSwitchingAlign correctly aligns noisy parameters",
+{
+	N = 30 #number of matrices
+	d_K_list = list(c(2,2), c(5,3))
+	for (i in 1:2)
+	{
+		d = d_K_list[[i]][1]
+		K = d_K_list[[i]][2]
+		max_error = d * 0.2 #TODO: what value to choose ?
+
+		# 1] Generate matrix series
+		matrices_permut = .generateMatrices(d,K,N,noise=TRUE)
+
+		# 2] Call align function
+		matrices_aligned = alignMatrices(matrices_permut, ref="mean", ls_mode="exact")
+
+		# 3] Check alignment
+		for (j in 2:N)
+		{
+			expect_that( norm(matrices_aligned[[j]] - matrices_permut[[1]]),
+				is_less_than(max_error) )
+		}
+	}
+})
diff --git a/pkg/tests/testthat/test-computeMu.R b/pkg/tests/testthat/test-computeMu.R
new file mode 100644
index 0000000..8234cd7
--- /dev/null
+++ b/pkg/tests/testthat/test-computeMu.R
@@ -0,0 +1,35 @@
+context("computeMu")
+
+test_that("on input of sufficient size, β/||β|| is estimated accurately enough",
+{
+	n = 100000
+	d = 2
+	K = 2
+	p = 1/2
+
+	βs_ref = array( c(1,0,0,1 , 1,-2,3,1), dim=c(d,K,2) )
+	for (i in 1:(dim(βs_ref)[3]))
+	{
+		μ_ref = normalize(βs_ref[,,i])
+		for (model in c("logit","probit"))
+		{
+			cat("\n\n",model," :\n",sep="")
+
+			io = generateSampleIO(n, p, βs_ref[,,i], rep(0,K), model)
+			μ = computeMu(io$X, io$Y, list(K=K))
+			μ_aligned = alignMatrices(list(μ), ref=μ_ref, ls_mode="exact")[[1]]
+
+			#Some traces: 0 is not well estimated, but others are OK
+			cat("Reference normalized matrix:\n")
+			print(μ_ref)
+			cat("Estimated normalized matrix:\n")
+			print(μ_aligned)
+			cat("Difference norm (Matrix norm ||.||_1, max. abs. sum on a column)\n")
+			diff_norm = norm(μ_ref - μ_aligned)
+			cat(diff_norm,"\n")
+
+			#NOTE: 0.5 is loose threshold, but values around 0.3 are expected...
+			expect_that( diff_norm, is_less_than(0.5) )
+		}
+	}
+})
diff --git a/pkg/tests/testthat/test-hungarianAlgorithm.R b/pkg/tests/testthat/test-hungarianAlgorithm.R
new file mode 100644
index 0000000..90cb02f
--- /dev/null
+++ b/pkg/tests/testthat/test-hungarianAlgorithm.R
@@ -0,0 +1,24 @@
+context("hungarianAlgorithm")
+
+test_that("HungarianAlgorithm provides the correct answer on various inputs",
+{
+	for (n in c(2,3,10,50))
+	{
+		for (repetition in 1:3)
+		{
+			# Generate a random vector, and permute it:
+			# we expect the algorithm to retrieve the permutation
+			v = runif(n, min=-1, max=1)
+			permutation = sample(1:n)
+			v_p = v[permutation]
+
+			# v is reference, v_p is v after permutation
+			# distances[i,j] = distance between i-th component of v and j-th component of v_p
+			# "in rows : reference; in columns: after permutation"
+			# "permutation" contains order on v to match v_p as closely as possible
+			distances = sapply(v_p, function(elt) ( abs(elt - v) ) )
+			assignment = .hungarianAlgorithm(distances)
+			expect_equal(assignment, permutation)
+		}
+	}
+})
diff --git a/pkg/tests/testthat/test-jointDiag.R b/pkg/tests/testthat/test-jointDiag.R
new file mode 100644
index 0000000..0eb6ffd
--- /dev/null
+++ b/pkg/tests/testthat/test-jointDiag.R
@@ -0,0 +1,50 @@
+context("jointDiag::ajd")
+
+#auxiliary to test diagonality
+.computeMuCheckDiag = function(X, Y, K, jd_method, β_ref)
+{
+	d = ncol(X)
+	#TODO: redundant code, same as computeMu() main method. Comments are stripped
+	M3 = .Moments_M3(X,Y)
+	M2_t = array(dim=c(d,d,d))
+	for (i in 1:d)
+	{
+		ρ = c(rep(0,i-1),1,rep(0,d-i))
+		M2_t[,,i] = .T_I_I_w(M3,ρ)
+	}
+	jd = jointDiag::ajd(M2_t, method=jd_method)
+	V = if (jd_method=="uwedge") jd$B else solve(jd$A)
+	M2_t = array(dim=c(d,d,K))
+	for (i in 1:K)
+		M2_t[,,i] = .T_I_I_w(M3,V[,i])
+	#END of computeMu() code
+
+	max_error = 0.5 #TODO: tune ?
+	invβ = MASS::ginv(β_ref)
+	for (i in 1:K)
+	{
+		shouldBeDiag = invβ %*% M2_t[,,i] %*% t(invβ)
+		expect_that(
+			mean( abs(shouldBeDiag[upper.tri(shouldBeDiag) | lower.tri(shouldBeDiag)]) ),
+			is_less_than(max_error) )
+	}
+}
+
+test_that("'jedi' and 'uwedge' joint-diagonalization methods return a correct matrix",
+{
+	n = 10000
+	d_K = list( c(2,2), c(5,3), c(20,13) )
+
+	for (dk_index in 1:length(d_K))
+	{
+		d = d_K[[dk_index]][1]
+		K = d_K[[dk_index]][2]
+		#NOTE: sometimes large errors if pr is not balanced enough (e.g. random);
+		#      same note for β. However we could be more random than that...
+		β_ref = rbind(diag(K),matrix(0,nrow=d-K,ncol=K))
+		io = generateSampleIO(n, p=rep(1/K,K-1), β=β_ref, rep(0,K), link="logit")
+		.computeMuCheckDiag(io$X, io$Y, K, jd_method="uwedge", β_ref)
+		#TODO: some issues with jedi method (singular system)
+		#.computeMuCheckDiag(io$X, io$Y, K, jd_method="jedi", β_ref)
+	}
+})
diff --git a/pkg/tests/testthat/test-optimParams.R b/pkg/tests/testthat/test-optimParams.R
new file mode 100644
index 0000000..8d0e7b3
--- /dev/null
+++ b/pkg/tests/testthat/test-optimParams.R
@@ -0,0 +1,87 @@
+context("OptimParams")
+
+naive_f = function(link, M1,M2,M3, p,β,b)
+{
+	d = length(M1)
+	K = length(p)
+	λ <- sqrt(colSums(β^2))
+
+	# Compute β x2,3 (self) tensorial products
+	β2 = array(0, dim=c(d,d,K))
+	β3 = array(0, dim=c(d,d,d,K))
+	for (k in 1:K)
+	{
+		for (i in 1:d)
+		{
+			for (j in 1:d)
+			{
+				β2[i,j,k] = β[i,k]*β[j,k]
+				for (l in 1:d)
+					β3[i,j,l,k] = β[i,k]*β[j,k]*β[l,k]
+			}
+		}
+	}
+
+	res = 0
+	for (i in 1:d)
+	{
+		term = 0
+		for (k in 1:K)
+			term = term + p[k]*.G(link,1,λ[k],b[k])*β[i,k]
+		res = res + (term - M1[i])^2
+		for (j in 1:d)
+		{
+			term = 0
+			for (k in 1:K)
+				term = term + p[k]*.G(link,2,λ[k],b[k])*β2[i,j,k]
+			res = res + (term - M2[i,j])^2
+			for (l in 1:d)
+			{
+				term = 0
+				for (k in 1:K)
+					term = term + p[k]*.G(link,3,λ[k],b[k])*β3[i,j,l,k]
+				res = res + (term - M3[i,j,l])^2
+			}
+		}
+	}
+	res
+}
+
+test_that("naive computation provides the same result as vectorized computations",
+{
+	h <- 1e-7 #for finite-difference tests
+	tol <- .25 * sqrt(h) #about 7.9 e-5
+	for (dK in list( c(2,2), c(5,3)))
+	{
+		d = dK[1]
+		K = dK[2]
+
+		M1 = runif(d, -1, 1)
+		M2 = matrix(runif(d*d,-1,1), ncol=d)
+		M3 = array(runif(d*d*d,-1,1), dim=c(d,d,d))
+
+		for (link in c("logit","probit"))
+		{
+			op = new("OptimParams", "li"=link, "M1"=as.double(M1),
+				"M2"=as.double(M2), "M3"=as.double(M3), "K"=as.integer(K))
+
+			for (var in seq_len((2+d)*K-1))
+			{
+				p = runif(K, 0, 1)
+				p = p / sum(p)
+				β <- matrix(runif(d*K,-5,5),ncol=K)
+				b = runif(K, -5, 5)
+				x <- c(p[1:(K-1)],as.double(β),b)
+
+				# Test functions values
+				expect_equal( op$f(x), naive_f(link,M1,M2,M3, p,β,b) )
+
+				# Test finite differences ~= gradient values
+				dir_h <- rep(0, (2+d)*K-1)
+				dir_h[var] = h
+
+				expect_equal( op$grad_f(x)[var], ( op$f(x+dir_h) - op$f(x) ) / h, tol )
+			}
+		}
+	}
+})
diff --git a/to-cran.sh b/to-cran.sh
new file mode 100755
index 0000000..b7e0ccd
--- /dev/null
+++ b/to-cran.sh
@@ -0,0 +1,22 @@
+#!/bin/bash
+
+# Synchronize package folder
+rsync -a --delete pkg/ pkg-cran/
+cd pkg-cran
+
+# Replace all used greek letters by beta, mu, ...etc
+for file in `find -name *.R`; do
+	if [ -f $file ]; then
+		sed -i 's/μ/mu/g' $file
+		sed -i 's/β/beta/g' $file
+		sed -i 's/λ/lambda/g' $file
+		sed -i 's/Σ/Sigma/g' $file
+	fi
+done
+
+if [ "$1" == "i" ]; then
+	# Install package
+	echo "rp()" | R --slave
+fi
+
+cd ..
-- 
2.44.0