From: Benjamin Auder <benjamin.auder@somewhere>
Date: Tue, 10 Dec 2019 11:42:29 +0000 (+0100)
Subject: Replace tabs by 2spaces
X-Git-Url: https://git.auder.net/doc/html/app_dev.php/current/css/rpsls.css?a=commitdiff_plain;h=6dd5c2acccd10635449230faa824b7e8906911bf;p=morpheus.git

Replace tabs by 2spaces
---

diff --git a/pkg/R/computeMu.R b/pkg/R/computeMu.R
index 2a62793..f7f82ad 100644
--- a/pkg/R/computeMu.R
+++ b/pkg/R/computeMu.R
@@ -26,67 +26,67 @@
 #' @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")
+  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 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 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 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])
+  # 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]
-	μ
+  # 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
index 8b70d6d..56c3e19 100644
--- a/pkg/R/multiRun.R
+++ b/pkg/R/multiRun.R
@@ -82,48 +82,48 @@
 #'   res[[i]] <- alignMatrices(res[[i]], ref=β, ls_mode="exact")}
 #' @export
 multiRun <- function(fargs, estimParams,
-	prepareArgs = function(x,i) x, N=10, ncores=3, agg=lapply, verbose=FALSE)
+  prepareArgs = function(x,i) 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")
+  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, index)
-		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
-		})
-	}
+  estimParamAtIndex <- function(index)
+  {
+    fargs <- prepareArgs(fargs, index)
+    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)
+  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]]) )
+  # 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
index a45f71a..505b665 100644
--- a/pkg/R/optimParams.R
+++ b/pkg/R/optimParams.R
@@ -33,17 +33,17 @@
 #' @export
 optimParams <- function(X, Y, K, link=c("logit","probit"))
 {
-	# Check arguments
+  # Check arguments
   if (!is.matrix(X) || any(is.na(X)))
     stop("X: numeric matrix, no NAs")
   if (!is.numeric(Y) || any(is.na(Y)) || any(Y!=0 & Y!=1))
     stop("Y: binary vector with 0 and 1 only")
-	link <- match.arg(link)
+  link <- match.arg(link)
   if (!is.numeric(K) || K!=floor(K) || K < 2)
     stop("K: integer >= 2")
 
-	# Build and return optimization algorithm object
-	methods::new("OptimParams", "li"=link, "X"=X,
+  # Build and return optimization algorithm object
+  methods::new("OptimParams", "li"=link, "X"=X,
     "Y"=as.integer(Y), "K"=as.integer(K))
 }
 
@@ -60,30 +60,30 @@ optimParams <- function(X, Y, K, link=c("logit","probit"))
 #' @field W Weights matrix (iteratively refined)
 #'
 setRefClass(
-	Class = "OptimParams",
+  Class = "OptimParams",
 
-	fields = list(
-		# Inputs
-		li = "character", #link function
-		X = "matrix",
-		Y = "numeric",
+  fields = list(
+    # Inputs
+    li = "character", #link function
+    X = "matrix",
+    Y = "numeric",
     Mhat = "numeric", #vector of empirical moments
-		# Dimensions
-		K = "integer",
+    # Dimensions
+    K = "integer",
     n = "integer",
-		d = "integer",
+    d = "integer",
     # Weights matrix (generalized least square)
     W = "matrix"
-	),
+  ),
 
-	methods = list(
-		initialize = function(...)
-		{
-			"Check args and initialize K, d, W"
+  methods = list(
+    initialize = function(...)
+    {
+      "Check args and initialize K, d, W"
 
       callSuper(...)
-			if (!hasArg("X") || !hasArg("Y") || !hasArg("K") || !hasArg("li"))
-				stop("Missing arguments")
+      if (!hasArg("X") || !hasArg("Y") || !hasArg("K") || !hasArg("li"))
+        stop("Missing arguments")
 
       # Precompute empirical moments
       M <- computeMoments(X, Y)
@@ -92,28 +92,28 @@ setRefClass(
       M3 <- as.double(M[[3]])
       Mhat <<- c(M1, M2, M3)
 
-			n <<- nrow(X)
-			d <<- length(M1)
+      n <<- nrow(X)
+      d <<- length(M1)
       W <<- diag(d+d^2+d^3) #initialize at W = Identity
-		},
+    },
 
-		expArgs = function(v)
-		{
-			"Expand individual arguments from vector v into a list"
+    expArgs = function(v)
+    {
+      "Expand individual arguments from vector v into a list"
 
-			list(
-				# p: dimension K-1, need to be completed
-				"p" = c(v[1:(K-1)], 1-sum(v[1:(K-1)])),
-				"β" = matrix(v[K:(K+d*K-1)], ncol=K),
-				"b" = v[(K+d*K):(K+(d+1)*K-1)])
-		},
+      list(
+        # p: dimension K-1, need to be completed
+        "p" = c(v[1:(K-1)], 1-sum(v[1:(K-1)])),
+        "β" = matrix(v[K:(K+d*K-1)], ncol=K),
+        "b" = v[(K+d*K):(K+(d+1)*K-1)])
+    },
 
-		linArgs = function(L)
-		{
-			"Linearize vectors+matrices from list L into a vector"
+    linArgs = function(L)
+    {
+      "Linearize vectors+matrices from list L into a vector"
 
-			c(L$p[1:(K-1)], as.double(L$β), L$b)
-		},
+      c(L$p[1:(K-1)], as.double(L$β), L$b)
+    },
 
     computeW = function(θ)
     {
@@ -130,32 +130,32 @@ setRefClass(
       "Vector of moments, of size d+d^2+d^3"
 
       p <- θ$p
-			β <- θ$β
-			λ <- sqrt(colSums(β^2))
-			b <- θ$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)
-
-			c(
-				β  %*% (p * .G(li,1,λ,b)),
-				β2 %*% (p * .G(li,2,λ,b)),
-				β3 %*% (p * .G(li,3,λ,b)))
+      β <- θ$β
+      λ <- sqrt(colSums(β^2))
+      b <- θ$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)
+
+      c(
+        β  %*% (p * .G(li,1,λ,b)),
+        β2 %*% (p * .G(li,2,λ,b)),
+        β3 %*% (p * .G(li,3,λ,b)))
     },
 
     f = function(θ)
     {
-			"Product t(hat_Mi - Mi) W (hat_Mi - Mi) with Mi(theta)"
+      "Product t(hat_Mi - Mi) W (hat_Mi - Mi) with Mi(theta)"
 
       L <- expArgs(θ)
-			A <- as.matrix(Mhat - Moments(L))
+      A <- as.matrix(Mhat - Moments(L))
       t(A) %*% W %*% A
     },
 
-		grad_f = function(θ)
-		{
-			"Gradient of f, dimension (K-1) + d*K + K = (d+2)*K - 1"
+    grad_f = function(θ)
+    {
+      "Gradient of f, dimension (K-1) + d*K + K = (d+2)*K - 1"
 
       L <- expArgs(θ)
       -2 * t(grad_M(L)) %*% W %*% as.matrix((Mhat - Moments(L)))
@@ -165,80 +165,80 @@ setRefClass(
     {
       "Gradient of the vector of moments, size (dim=)d+d^2+d^3 x K-1+K+d*K"
 
-			p <- θ$p
-			β <- θ$β
-			λ <- sqrt(colSums(β^2))
-			μ <- sweep(β, 2, λ, '/')
-			b <- θ$b
+      p <- θ$p
+      β <- θ$β
+      λ <- sqrt(colSums(β^2))
+      μ <- sweep(β, 2, λ, '/')
+      b <- θ$b
 
       res <- matrix(nrow=nrow(W), ncol=0)
 
-			# 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)
+      # 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)
+      # 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)
 
       # Gradient on p: K-1 columns, dim rows
-			km1 = 1:(K-1)
-			res <- cbind(res, rbind(
+      km1 = 1:(K-1)
+      res <- cbind(res, rbind(
         sweep(as.matrix(β [,km1]), 2, G1[km1], '*') - G1[K] * β [,K],
         sweep(as.matrix(β2[,km1]), 2, G2[km1], '*') - G2[K] * β2[,K],
         sweep(as.matrix(β3[,km1]), 2, G3[km1], '*') - G3[K] * β3[,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, '*')
-
-				res <- cbind(res, rbind(dβ, dβ2, dβ3))
-			}
+      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, '*')
+
+        res <- cbind(res, rbind(dβ, dβ2, dβ3))
+      }
 
       # Gradient on b
-			res <- cbind(res, rbind(
-				sweep(β,  2, p * G2, '*'),
-				sweep(β2, 2, p * G3, '*'),
-				sweep(β3, 2, p * G4, '*') ))
+      res <- cbind(res, rbind(
+        sweep(β,  2, p * G2, '*'),
+        sweep(β2, 2, p * G3, '*'),
+        sweep(β3, 2, p * G4, '*') ))
 
-			res
-		},
+      res
+    },
 
-		run = function(θ0)
-		{
-			"Run optimization from θ0 with solver..."
+    run = function(θ0)
+    {
+      "Run optimization from θ0 with solver..."
 
-	    if (!is.list(θ0))
-		    stop("θ0: list")
+      if (!is.list(θ0))
+        stop("θ0: list")
       if (is.null(θ0$β))
         stop("At least θ0$β must be provided")
-			if (!is.matrix(θ0$β) || any(is.na(θ0$β)) || ncol(θ0$β) != K)
-				stop("θ0$β: matrix, no NA, ncol == K")
+      if (!is.matrix(θ0$β) || any(is.na(θ0$β)) || ncol(θ0$β) != K)
+        stop("θ0$β: matrix, no NA, ncol == K")
       if (is.null(θ0$p))
         θ0$p = rep(1/K, K-1)
       else if (length(θ0$p) != K-1 || sum(θ0$p) > 1)
@@ -265,9 +265,9 @@ setRefClass(
         print(expArgs(op_res$par))
       }
 
-			expArgs(op_res$par)
-		}
-	)
+      expArgs(op_res$par)
+    }
+  )
 )
 
 # Compute vectorial E[g^{(order)}(<β,x> + b)] with x~N(0,Id) (integral in R^d)
@@ -281,9 +281,9 @@ setRefClass(
 #
 .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...)
+  # 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...)
   sapply( seq_along(λ), function(k) {
     res <- NULL
     tryCatch({
@@ -306,24 +306,24 @@ setRefClass(
 # Derivatives list: g^(k)(x) for links 'logit' and 'probit'
 #
 .deriv <- list(
-	"probit"=list(
-		# 'probit' derivatives list;
-		# NOTE: exact values for the integral E[g^(k)(λz+b)] could be computed
-		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)
-	)
+  "probit"=list(
+    # 'probit' derivatives list;
+    # NOTE: exact values for the integral E[g^(k)(λz+b)] could be computed
+    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)
@@ -332,6 +332,6 @@ setRefClass(
 #
 .zin <- function(x)
 {
-	x[is.nan(x)] <- 0.
-	x
+  x[is.nan(x)] <- 0.
+  x
 }
diff --git a/pkg/R/plot.R b/pkg/R/plot.R
index 29a254e..9727493 100644
--- a/pkg/R/plot.R
+++ b/pkg/R/plot.R
@@ -6,10 +6,10 @@
 #
 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])
-	} )
+  # 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
@@ -31,12 +31,12 @@ extractParam <- function(mr, x=1, y=1)
 #' @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")
+  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
@@ -50,12 +50,12 @@ plotHist <- function(mr, x, y)
 #' @export
 plotBox <- function(mr, x, y, xtitle="")
 {
-	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]], xlab=xtitle, ylab="Parameter value")
+  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]], xlab=xtitle, ylab="Parameter value")
 }
 
 #' plotCoefs
@@ -71,32 +71,32 @@ plotBox <- function(mr, x, y, xtitle="")
 #' @export
 plotCoefs <- function(mr, params, idx, xtitle="Parameter")
 {
-	L <- nrow(mr[[1]][[1]])
-	K <- ncol(mr[[1]][[1]])
+  L <- nrow(mr[[1]][[1]])
+  K <- ncol(mr[[1]][[1]])
 
-	params_hat <- matrix(nrow=L, ncol=K)
-	stdev <- matrix(nrow=L, ncol=K)
-	for (x in 1:L)
-	{
-		for (y in 1:K)
-		{
-			estims <- extractParam(mr, x, y)
-			params_hat[x,y] <- mean(estims[[idx]])
-#			stdev[x,y] <- sqrt( mean( (estims[[idx]] - params[x,y])^2 ) )
-			# HACK remove extreme quantile in estims[[i]] before computing sd()
-			stdev[x,y] <- sd( estims[[idx]] ) #[ estims[[idx]] < max(estims[[idx]]) & estims[[idx]] > min(estims[[idx]]) ] )
-		}
-	}
+  params_hat <- matrix(nrow=L, ncol=K)
+  stdev <- matrix(nrow=L, ncol=K)
+  for (x in 1:L)
+  {
+    for (y in 1:K)
+    {
+      estims <- extractParam(mr, x, y)
+      params_hat[x,y] <- mean(estims[[idx]])
+#      stdev[x,y] <- sqrt( mean( (estims[[idx]] - params[x,y])^2 ) )
+      # HACK remove extreme quantile in estims[[i]] before computing sd()
+      stdev[x,y] <- sd( estims[[idx]] ) #[ estims[[idx]] < max(estims[[idx]]) & estims[[idx]] > min(estims[[idx]]) ] )
+    }
+  }
 
-	par(cex.axis=1.5, cex.lab=1.5, mar=c(4.7,5,1,1))
-	params <- as.double(params)
-	o <- order(params)
-	avg_param <- as.double(params_hat)
-	std_param <- as.double(stdev)
-	matplot(cbind(params[o],avg_param[o],avg_param[o]+std_param[o],avg_param[o]-std_param[o]),
-		col=1, lty=c(1,5,3,3), type="l", lwd=2, xlab=xtitle, ylab="")
+  par(cex.axis=1.5, cex.lab=1.5, mar=c(4.7,5,1,1))
+  params <- as.double(params)
+  o <- order(params)
+  avg_param <- as.double(params_hat)
+  std_param <- as.double(stdev)
+  matplot(cbind(params[o],avg_param[o],avg_param[o]+std_param[o],avg_param[o]-std_param[o]),
+    col=1, lty=c(1,5,3,3), type="l", lwd=2, xlab=xtitle, ylab="")
 
-	#print(o) #not returning o to avoid weird Jupyter issue... (TODO:)
+  #print(o) #not returning o to avoid weird Jupyter issue... (TODO:)
 }
 
 #' plotQn
@@ -113,19 +113,19 @@ plotCoefs <- function(mr, params, idx, xtitle="Parameter")
 #' @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...
+  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
index 6fa38ae..3d46092 100644
--- a/pkg/R/sampleIO.R
+++ b/pkg/R/sampleIO.R
@@ -20,54 +20,54 @@
 #' @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")
+  # 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)
+  #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)
-		#probas = rowSums(p * probas)
-		X = rbind(X, newXblock)
-		#Y = c( Y, vapply(probas, function(p) (ifelse(p >= .5, 1, 0)), 1) )
-		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])
+  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)
+    #probas = rowSums(p * probas)
+    X = rbind(X, newXblock)
+    #Y = c( Y, vapply(probas, function(p) (ifelse(p >= .5, 1, 0)), 1) )
+    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
index 6d1c361..ff988a0 100644
--- a/pkg/R/utils.R
+++ b/pkg/R/utils.R
@@ -9,9 +9,9 @@
 #' @export
 normalize = function(X)
 {
-	X = as.matrix(X)
-	norm2 = sqrt( colSums(X^2) )
-	sweep(X, 2, norm2, '/')
+  X = as.matrix(X)
+  norm2 = sqrt( colSums(X^2) )
+  sweep(X, 2, norm2, '/')
 }
 
 # Computes a tensor-vector product
@@ -23,11 +23,11 @@ normalize = function(X)
 #
 .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
+  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
@@ -39,11 +39,11 @@ normalize = function(X)
 #
 .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)
+  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
@@ -55,11 +55,11 @@ normalize = function(X)
 #
 .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) )
+  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
@@ -72,7 +72,7 @@ normalize = function(X)
 #'
 #' @export
 computeMoments = function(X, Y)
-	list( colMeans(Y * X), .Moments_M2(X,Y), .Moments_M3(X,Y) )
+  list( colMeans(Y * X), .Moments_M2(X,Y), .Moments_M3(X,Y) )
 
 # Find the optimal assignment (permutation) between two sets (minimize cost)
 #
@@ -83,9 +83,9 @@ computeMoments = function(X, Y)
 #
 .hungarianAlgorithm = function(distances)
 {
-	n = nrow(distances)
-	.C("hungarianAlgorithm", distances=as.double(distances), pn=as.integer(n),
-		assignment=integer(n), PACKAGE="morpheus")$assignment
+  n = nrow(distances)
+  .C("hungarianAlgorithm", distances=as.double(distances), pn=as.integer(n),
+    assignment=integer(n), PACKAGE="morpheus")$assignment
 }
 
 #' alignMatrices
@@ -105,65 +105,65 @@ computeMoments = function(X, Y)
 #' @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'}")
+  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
+  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]
-			}
-		}
+    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
+    # 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/src/functions.c b/pkg/src/functions.c
index f311a34..e534c57 100644
--- a/pkg/src/functions.c
+++ b/pkg/src/functions.c
@@ -3,62 +3,62 @@
 // Index matrix (by columns)
 int mi(int i, int j, int d1, int d2)
 {
-	return j*d1 + i;
+  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;
+  return k*d1*d2 + j*d1 + i;
 }
 
 // Empirical 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));
+  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;
-		}
-	}
+  // 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;
+    }
+  }
 }
 
 // Empirical 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));
+  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;
-			}
-		}
-	}
+  // 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;
+      }
+    }
+  }
 }
 
 // W = 1/N sum( t(g(Zi,theta)) g(Zi,theta) )
 // with g(Zi, theta) = i-th contribution to all moments (size dim) - real moments
 void Compute_Omega(double* X, double* Y, double* M, int* pn, int* pd, double* W)
 {
-	int n=*pn, d=*pd;
+  int n=*pn, d=*pd;
   int dim = d + d*d + d*d*d;
   //double* W = (double*)malloc(dim*dim*sizeof(double));
 
@@ -81,7 +81,7 @@ void Compute_Omega(double* X, double* Y, double* M, int* pn, int* pd, double* W)
       g[j] = 0.0;
       if (idx1 == idx2)
         g[j] -= Y[i];
-			g[j] += Y[i] * X[mi(i,idx1,n,d)]*X[mi(i,idx2,n,d)] - M[j];
+      g[j] += Y[i] * X[mi(i,idx1,n,d)]*X[mi(i,idx2,n,d)] - M[j];
     }
     for (int j=d+d*d; j<dim; j++)
     {
diff --git a/pkg/src/hungarian.c b/pkg/src/hungarian.c
index 7850f64..9671a43 100644
--- a/pkg/src/hungarian.c
+++ b/pkg/src/hungarian.c
@@ -11,7 +11,7 @@
  **
  ** 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",
+ ** As asked by  the copyright node of the "Stanford GraphGase",
  ** I hereby proclaim that this file are *NOT* part of the
  ** "Stanford GraphGase" distrubition!
  **
@@ -36,318 +36,318 @@
 #define verbose (0)
 
 typedef struct {
-	int num_rows;
-	int num_cols;
-	double** cost;
-	int** assignment;
+  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;
+  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 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;
+  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;
+  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;
+  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
-		}
+    ;
+  }
+  // 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
-	}
+    // 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);
+/*  // 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);
 }
 
 
@@ -373,24 +373,24 @@ double** array_to_matrix(double* m, int rows, int cols)
 // 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);
+  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_init(&p, mat, n, n, HUNGARIAN_MODE_MINIMIZE_COST);
 
-	hungarian_solve(&p);
+  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
-		}
-	}
+  // 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);
+  for (int i=0; i<n; i++)
+    free(mat[i]);
+  free(mat);
 }
diff --git a/pkg/src/morpheus_init.c b/pkg/src/morpheus_init.c
index 5352ee8..b872488 100644
--- a/pkg/src/morpheus_init.c
+++ b/pkg/src/morpheus_init.c
@@ -12,16 +12,16 @@ extern void Moments_M3(void*, void*, void*, void*, void*);
 extern void Compute_Omega(void*, void*, void*, void*, void*, void*);
 
 static const R_CMethodDef CEntries[] = {
-    {"hungarianAlgorithm", (DL_FUNC) &hungarianAlgorithm, 3},
-    {"Moments_M2",         (DL_FUNC) &Moments_M2,         5},
-    {"Moments_M3",         (DL_FUNC) &Moments_M3,         5},
-    {"Compute_Omega",      (DL_FUNC) &Compute_Omega,      6},
-    {NULL, NULL, 0}
+  {"hungarianAlgorithm", (DL_FUNC) &hungarianAlgorithm, 3},
+  {"Moments_M2",         (DL_FUNC) &Moments_M2,         5},
+  {"Moments_M3",         (DL_FUNC) &Moments_M3,         5},
+  {"Compute_Omega",      (DL_FUNC) &Compute_Omega,      6},
+  {NULL, NULL, 0}
 };
 
 void R_init_morpheus(DllInfo *dll)
 {
-    R_registerRoutines(dll, CEntries, NULL, NULL, NULL);
-    R_useDynamicSymbols(dll, FALSE);
+  R_registerRoutines(dll, CEntries, NULL, NULL, NULL);
+  R_useDynamicSymbols(dll, FALSE);
 }
 
diff --git a/pkg/tests/testthat/test-Moments_M2.R b/pkg/tests/testthat/test-Moments_M2.R
index 5201104..0395d4e 100644
--- a/pkg/tests/testthat/test-Moments_M2.R
+++ b/pkg/tests/testthat/test-Moments_M2.R
@@ -3,41 +3,41 @@ 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))
+  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
+  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
-	}
+  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)
+  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)
-		}
-	}
+  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
index 502ceb3..8ac3e75 100644
--- a/pkg/tests/testthat/test-Moments_M3.R
+++ b/pkg/tests/testthat/test-Moments_M3.R
@@ -3,49 +3,49 @@ 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))
+  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
+  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
-	}
+  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)
+  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)
-		}
-	}
+  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
index 59ad1d1..8ce7abf 100644
--- a/pkg/tests/testthat/test-alignMatrices.R
+++ b/pkg/tests/testthat/test-alignMatrices.R
@@ -3,67 +3,67 @@ 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
+  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]
+  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)
+    # 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")
+    # 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]])
+    # 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")
+    # 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]])
-	}
+    # 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 ?
+  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)
+    # 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")
+    # 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) )
-		}
-	}
+    # 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
index 8234cd7..fccef55 100644
--- a/pkg/tests/testthat/test-computeMu.R
+++ b/pkg/tests/testthat/test-computeMu.R
@@ -2,34 +2,34 @@ context("computeMu")
 
 test_that("on input of sufficient size, β/||β|| is estimated accurately enough",
 {
-	n = 100000
-	d = 2
-	K = 2
-	p = 1/2
+  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="")
+  β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]]
+      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")
+      #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) )
-		}
-	}
+      #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
index 90cb02f..591c3c1 100644
--- a/pkg/tests/testthat/test-hungarianAlgorithm.R
+++ b/pkg/tests/testthat/test-hungarianAlgorithm.R
@@ -2,23 +2,23 @@ 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]
+  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)
-		}
-	}
+      # 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
index 0eb6ffd..bb12a6c 100644
--- a/pkg/tests/testthat/test-jointDiag.R
+++ b/pkg/tests/testthat/test-jointDiag.R
@@ -3,48 +3,48 @@ 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
+  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) )
-	}
+  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) )
+  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)
-	}
+  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
index 6864804..8f65e46 100644
--- a/pkg/tests/testthat/test-optimParams.R
+++ b/pkg/tests/testthat/test-optimParams.R
@@ -2,86 +2,86 @@ context("OptimParams")
 
 naive_f = function(link, M1,M2,M3, p,β,b)
 {
-	d = length(M1)
-	K = length(p)
-	λ <- sqrt(colSums(β^2))
+  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]
-			}
-		}
-	}
+  # 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
+  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 <- 1e-3 #large tolerance, necessary in some cases... (generally 1e-6 is OK)
-	for (dK in list( c(2,2), c(5,3)))
-	{
-		d = dK[1]
-		K = dK[2]
+  h <- 1e-7 #for finite-difference tests
+  tol <- 1e-3 #large tolerance, necessary in some cases... (generally 1e-6 is OK)
+  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))
+    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 (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)
+      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 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
+        # 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 )
-			}
-		}
-	}
+        expect_equal( op$grad_f(x)[var], ( op$f(x+dir_h) - op$f(x) ) / h, tol )
+      }
+    }
+  }
 })