From 6e9cf699337288482f1d6aa531e3e55229aa307e Mon Sep 17 00:00:00 2001
From: Raimon Tolosana-Delgado <tolosa53@fwg206.ad.fz-rossendorf.de>
Date: Fri, 18 Sep 2020 14:47:07 +0200
Subject: [PATCH] minor changes

---
 DESCRIPTION  |   2 +-
 NAMESPACE    |   3 +
 NEWS.md      |   8 +
 R/genDiag.R  | 964 +++++++++++++++++++++++++--------------------------
 R/geostats.R |  34 +-
 R/grids.R    |  15 +-
 6 files changed, 520 insertions(+), 506 deletions(-)
 create mode 100644 NEWS.md

diff --git a/DESCRIPTION b/DESCRIPTION
index 50a8413..ded559b 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,5 +1,5 @@
 Package: gmGeostats
-Version: 0.10-6
+Version: 0.10-6.9000
 Date: 2020-09-07
 Title: Geostatistics for Compositional Analysis
 Authors@R: c(person(given = "Raimon", 
diff --git a/NAMESPACE b/NAMESPACE
index 10cf645..77c6ff8 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -59,7 +59,10 @@ S3method(getStackElement,default)
 S3method(getStackElement,list)
 S3method(gmApply,DataFrameStack)
 S3method(gmApply,default)
+S3method(gsi.gstatCokriging2compo,data.frame)
 S3method(gsi.gstatCokriging2compo,default)
+S3method(gsi.gstatCokriging2rmult,data.frame)
+S3method(gsi.gstatCokriging2rmult,default)
 S3method(has.missings,data.frame)
 S3method(image,logratioVariogramAnisotropy)
 S3method(image,mask)
diff --git a/NEWS.md b/NEWS.md
new file mode 100644
index 0000000..e4b75b3
--- /dev/null
+++ b/NEWS.md
@@ -0,0 +1,8 @@
+# gmGeostats 0.10.x
+
+* (2020-09-18) conversion of variogram models of class CompLinModCoReg -> LMCAnisCompo -> variogramModel(List)
+* (2020-09-18) `gsi.gstatCokriging2compo()` and `gsi.gstatCokriging2rmult()` now work with systems of 3 components resp. 2 variables (systems of 2 components or 1 variable still fail)
+
+# gmGeostats 0.10.6
+
+* (2020-09-17) we are in CRAN! 
\ No newline at end of file
diff --git a/R/genDiag.R b/R/genDiag.R
index f6472d9..ef21837 100644
--- a/R/genDiag.R
+++ b/R/genDiag.R
@@ -1,483 +1,481 @@
-#### MAF calculations (2 versions) ----------------
-# compute the MAF basis from two matrices of (alr/ilr/clr)-covariance 
-gsi.getMAFbasis = function(C0, CT, tol=1e-12){
-  svdT = svd(CT)
-  A1 = with(svdT, u %*% diag(ifelse(d>tol, 1/sqrt(d), 0)) )
-  B = t(A1) %*% C0 %*% A1
-  svdB = svd(B)
-  o = nrow(svdB$u):1  
-  A = A1[,o] %*% svdB$u[o,]
-  return(A)
-}
-
-  
-gsi.computeExplainedVariance <- function(X, Winv){
-  rs = sapply(1:ncol(X), function(i) mvar(rmult(outer(X[,i], Winv[i,]))))
-  return(rs)
-}
-
-
-#' @describeIn Maf.data.frame  Generalised diagonalisations
-#' @param ... generic functionality arguments
-#' @export
-Maf = function(x,...) UseMethod("Maf",x)
-
-#' Generalised diagonalisations
-#' Calculate several generalized diagonalisations out of a data set and its empirical 
-#' variogram
-#' @param x a data set, typically of class "data.frame" or of a compositional class 
-#' @param vg empirical variogram, of a kind fitting to the data provided
-#' @param i a slicer for the variogram, typically this will be one or more indices of 
-#' the lag distance to take. %For other options see code{getEmpVariogramSlice}.
-#'
-#' @return An object extending \code{c("princomp.CLASSOF(x)",}"\code{\link{princomp}}") 
-#' with classes "\code{genDiag}", and an extra class argument depending on the 
-#' diagonalisation method chosen.
-#' 
-#' Function \code{Maf} results carry the extra class "\code{maf}", and they correspond
-#' to minimum/maximum autocorrelation factors (MAF) as proposed by Switzer and Green 
-#' (XXXX) or XXXXX (XXXX). In this case, the
-#' slicer is typically just the index of one lag distance (defaults to i=2). MAF 
-#' provides the analytical solution to the joint diagonalisation of two matrices,
-#' the covariance of increments provided by the slicing and the conventional covariance
-#' matrix (of the idt transformed values, if appropriate). Resulting factors are ordered
-#' in decreasing order of spatial continuity, which might produce surprising 
-#' scree-plots for those who are used to see a screeplot of a principal component analysis.
-#' 
-#' Function \code{UWEDGE} (Uniformly Weighted Exhaustive Diagonalization with Gauss 
-#' iterations) results carry the extra class "\code{uwedge}". The function 
-#' is a wrapper on \code{jointDiag::uwedge} from package \code{jointDiag}. In this case the 
-#' slicer is typically just a subset of indices of lag distances to consider 
-#' (defaults to the nearest indexes to mininum, maximum and mean lag distances of the 
-#' variogram). UWEDGE iteratively seeks for a pair of matrices (a mixing and a
-#' demixing matrices) diagonalises the set of matrices \eqn{M_1, M_2, \ldots, M_K} 
-#' given, by minimizing the target quantity 
-#' \deqn{Q_{uwedge}(A, W) = \sum_{k=1}^K Tr[E_k^t\cdot E_k],}
-#' where \eqn{E_k = (W^t\cdot M_k \cdot W- A\cdot \Lambda_k\cdot A^t)} and 
-#' \eqn{\Lambda_k = diag(W^t\cdot M_k \cdot W)} is the resulting diagonalized version of 
-#' each matrix. Obtained factors are ordered
-#' in decreasing order of spatial continuity, which might produce surprising 
-#' scree-plots for those who are used to see a screeplot of a principal component analysis.
-#' 
-#' Function \code{RJD} results carry the extra class "\code{rjd}". The function 
-#' is a wrapper on \code{JADE::rjd}, implementing the Rotational joint diagonalisation method . In this case the 
-#' slicer is typically just a subset of indices of lag distances to consider 
-#' (defaults to the nearest indexes to mininum, maximum and mean lag distances).
-#' This algorithm also served for quasi-diagonalising a set of matrices as in UWEDGE,
-#' just that in this case the quantity to minimise is the sum of sequares of all off-diagonal
-#' elements of \eqn{A^t\cdot M_k\cdot A} for all \eqn{k=1, 2, \ldots K}.
-#' 
-#' All these functions produce output mimicking \code{\link{princomp}}, i.e. with 
-#' elements
-#' \describe{
-#'   \item{sdev}{contrary to the output in PCA, this contains the square root of the
-#'   metric variance of the predictions obtained for each individual factor; this is the
-#'   quantity needed for \code{\link{screeplot}} to create plots of explained variance
-#'   by factor}
-#'   \item{loadings}{matrix of contributions of each (cdt-transformed) original variable to the new factors}
-#'   \item{center}{center of the data set (eventually, represented through \code{\link{cdt}}), 
-#'   in compositional methods}
-#'   \item{scale}{the scalings applied to each original variable}
-#'   \item{n.obs}{number of observations}
-#'   \item{scores}{the scores of the supplied data on the new factors}
-#'   \item{call}{the call to the function (attention: it actually may come much later)}
-#' } 
-#' and additionally some of the following arguments, in different order
-#' \describe{
-#'   \item{invLoadings}{matrix of contributions of each factor onto each original variable}
-#'   \item{Center}{compositional methods return here the cdt-backtransformed center}
-#'   \item{InvLoadings}{compositional methods return here the clr-backtransformed inverse loadings, so that
-#'   each column of this matrix can be understood as a composition on itself}
-#'   \item{DownInvLoadings}{compositional methods return here the clr-backtransformed "minus inverse loadings", so that
-#'   each column of this matrix can be understood as a composition on itself; details in 
-#'   \code{\link{princomp.acomp}} }
-#'   \item{C1, C2}{Maf returns the two matrices that were diagonalised}
-#'   \item{eigenvalues}{Maf returns the generalized eigenvalues of the diagonalisation of C1 and C2}
-#'   \item{gof}{UWEDGE returns the values of the goodness of fit criterion across sweeps}
-#'   \item{diagonalized}{RJD returns the diagonalized matrices, in an array of (K,D,D)-dimensions, being
-#'   D the number of variables in \code{x}}
-#'   \item{type}{a string describing which package and which function was used as a workhorse for
-#'   the calculation}
-#' }
-#'  
-#' 
-#' NOTE: if the arguments you provide to RJD and UWEDGE are not of the appropriate type
-#' (i.e. data.frames or equivalent) the default method of these functions will just attempt
-#' to call the basis functions JADE:rjd and JointDiag:uwedge respectively. 
-#' This will be the case if you provide \code{x} as a "\code{matrix}", or as
-#' an "\code{array}". In those cases, the output will NOT be structured as an extension
-#' to princomp results; instead they will be native output from those functions.
-#' 
-#'  
-#' @export
-#' @method Maf data.frame
-#' @aliases genDiag
-#' 
-#' @family generalised Diagonalisations
-#'
-#' @examples
-#' require("magrittr")
-#' require("gstat")
-#' require("compositions")
-#' data("jura", package="gstat")
-#' gs = gstat(id="Cd", formula=log(Cd)~1, locations=~Xloc+Yloc, data=jura.pred) %>% 
-#' gstat(id="Co", formula=log(Cd)~1, locations=~Xloc+Yloc, data=jura.pred) %>% 
-#' gstat(id="Cr", formula=log(Cr)~1, locations=~Xloc+Yloc, data=jura.pred) %>% 
-#' gstat(id="Cu", formula=log(Cu)~1, locations=~Xloc+Yloc, data=jura.pred) %>% 
-#' gstat(id="Ni", formula=log(Ni)~1, locations=~Xloc+Yloc, data=jura.pred) %>% 
-#' gstat(id="Pb", formula=log(Pb)~1, locations=~Xloc+Yloc, data=jura.pred) %>% 
-#' gstat(id="Zn", formula=log(Zn)~1, locations=~Xloc+Yloc, data=jura.pred)
-#' vg = variogram(gs)
-#' mf = Maf(aplus(jura.pred[, -(1:6)]), vg=vg)
-#' mf
-#' mf$loadings
-#' biplot(mf)
-Maf.data.frame <- function(x, vg, i=2,...){
-    cl <- match.call()
-  C1 = unclass(var(idt(x)))
-  mn = colMeans(cdt(x))
-  C2 = 0*C1
-  ## BEGIN: to change in the future to use gmEVario 
-  vg = as.gstatVariogram(vg)
-  #dists = unique(vg$dist)
-  #vv = vg[abs(dists[i]-vg$dist)<1e-9*dists[i],]
-  vv = sapply(split(vg$gamma, vg$id), function(x)x[i])
-  vv = vv[order(names(vv))]
-  C2[lower.tri(C2,diag=TRUE)]=vv
-  dd = diag(diag(C2))
-  C2 = unclass(C2 + t(C2) - dd)
-  ## END: to change in the future to use gmEVario 
-  
-  
-  dt.centred = scale(x, center = TRUE, scale=FALSE)
-  eT = eigen(C1)
-  A1 = with(eT, vectors %*% diag(ifelse(values>1e-12, 1/sqrt(values), 0)) ) ## C1^{-0.5}
-  B = t(A1) %*% C2 %*% A1 
-  eB = eigen(B)
-  A2 = eB$vectors
-  
-  eigenvalues = eB$values
-  ord = order(eigenvalues)
-  A2 = A2[, ord]
-  W = unclass(A1) %*% unclass(A2)
-  rownames(W) = colnames(x)
-  # do not normalize
-  # W = unclass(t(normalize(rmult(t(W))))) 
-  W = unclass(W)
-  Winv = solve(W)
-   colnames(Winv) = colnames(W)
-  facscores = unclass(dt.centred) %*% W
-  colnames(facscores) = paste("maf", 1:ncol(facscores), sep="")
-  sdevs = sqrt(gsi.computeExplainedVariance(facscores, Winv))
-  res = list(sdev=sdevs, loadings = W, center=mn, scale=rep(1, ncol(x)),
-             n.obs=nrow(x), scores = facscores, C1=C1, C2=C2, eigenvalues = eigenvalues[ord], 
-             invLoadings=Winv)
-  if(class(x)!="data.frame"){
-    res$Center = cdtInv(mn)
-    res$InvLoadings = cdtInv(Winv, orig=x)
-    res$InvDownLoadings = cdtInv(-Winv, orig=x)
-  }
-  res$call = cl
-  res$type="gmGeostats::MAF"
-  class(res)= c("maf","genDiag",paste("princomp", class(x), sep="."), "princomp")
-  return(res)
-}
-
-#' @describeIn Maf.data.frame  Generalised diagonalisations
-#' @method Maf rmult
-#' @export
-Maf.rmult <- Maf.data.frame
-
-#' @describeIn Maf.data.frame  Generalised diagonalisations
-#' @method Maf aplus
-#' @export
-Maf.aplus <- Maf.data.frame
-
-#' @describeIn Maf.data.frame  Generalised diagonalisations
-#' @method Maf rplus
-#' @export
-Maf.rplus <- Maf.data.frame
-
-#' @describeIn Maf.data.frame  Generalised diagonalisations
-#' @method Maf ccomp
-#' @export
-Maf.ccomp <- Maf.data.frame
-
-
-#' @describeIn Maf.data.frame  Generalised diagonalisations
-#' @method Maf rcomp
-#' @export
-Maf.rcomp  <- function(x, vg, i=2,...){
-  requireNamespace("compositions", quietly=TRUE)
-  cl <- match.call()
-  C1 = unclass(clrvar2ilr(cov(x)))
-  C2 = unclass(clrvar2ilr(variation2clrvar(vg$vg[i,,]))) 
-  clr.centred = scale(cdt(x), center = TRUE, scale=FALSE)
-  
-  eT = eigen(C1)
-  A1 = with(eT, vectors %*% diag(ifelse(values>1e-12, 1/sqrt(values), 0)) ) ## C1^{-0.5}
-  B = t(A1) %*% C2 %*% A1 
-  eB = eigen(B)
-  A2 = eB$vectors
-  
-  eigenvalues = eB$values
-  ord = order(eigenvalues)
-  A2 = A2[, ord]
-  W = unclass(A1) %*% unclass(A2)
-  # do not normalize
-  #W = unclass(t(normalize(rmult(t(W))))) 
-  Wclr = unclass(t(ilr2clr(t(W))))
-  rownames(Wclr) = colnames(x)
-  Winv = solve(W)
-  Wclr.inv = unclass(ilr2clr(Winv))
-    colnames(Wclr.inv) = colnames(x)
-  facscores = unclass(clr.centred) %*% Wclr
-  colnames(facscores) = paste("maf", 1:ncol(facscores), sep="")
-  sdevs = sqrt(gsi.computeExplainedVariance(facscores, Wclr.inv))
-  res = list(sdev=sdevs, loadings = Wclr, center=mean(cdt(x),...), scale=rep(1, ncol(x)),
-             n.obs=nrow(x), scores = facscores, 
-             C1=C1, C2=C2, eigenvalues = eigenvalues[ord], invLoadings=Wclr.inv,
-             Center = mean(x, ...),
-             InvLoadings = cdtInv(Wclr.inv, orig=x),
-             InvDownLoadings = cdtInv(-Wclr.inv, orig=x),
-             call = cl, type="gmGeostats::MAF")
-  kkk = paste("princomp", class(x), sep=".")
-  class(res)= c("maf","genDiag", kkk,  "princomp")
-  return(res)
-} 
-
-#' @describeIn Maf.data.frame  Generalised diagonalisations
-#' @method Maf acomp
-#' @export
-Maf.acomp <- Maf.rcomp
-
-
-
-#### UWEDGE calculations ----------
-
-#' @describeIn Maf.data.frame  Generalised diagonalisations
-#' @export
-UWEDGE = function(x,...){
-  if(!requireNamespace("jointDiag", quietly=FALSE)) stop("UWEDGE requires package 'jointDiag' installed")
-  UseMethod("UWEDGE",x)
-}
-
-#' @describeIn Maf.data.frame  Generalised diagonalisations
-#' @method UWEDGE default
-#' @export
-UWEDGE.default <- function(x, ...) jointDiag::uwedge(M=x, ...)
-
-#' @describeIn Maf.data.frame  Generalised diagonalisations
-#' @method UWEDGE acomp
-#' @export
-UWEDGE.acomp <-  function(x, vg, i=NULL, ...){
-  requireNamespace("compositions", quietly=TRUE)
-  cl <- match.call()
-  vg = as.logratioVariogram(vg)
-  if(is.null(i)){
-    i = c(1, floor(c(0.5,1)*dim(vg$vg)[1]))
-  }
-  V = compositions::ilrBase(D=ncol(x)) 
-  # alternative: V = loadings(princomp(x))
-  C1 = -clrvar2ilr(cov(x), V = V)
-  d = ncol(x)-1
-  M = -0.5*apply(vg$vg[c(1,i),,], 1, clrvar2ilr, V=V )
-  dim(M)=c(d,d,length(i)+1)
-  M[,,1]=C1
-  res = jointDiag::uwedge(M)
-  
-  ord = order(diag(res$B %*% M[,,length(i)] %*% t(res$B)))
-  res$B = res$B[ord,]
-  W = t(res$B)
-  # do not normalize
-  # W = unclass(t(normalize(rmult(t(W))))) 
-  W = unclass(W)
-  ilrx = ilr(x, V=V)
-  mns = mean(ilrx)
-  ilrxc = ilrx-mns
-  facscores = unclass(ilrxc) %*% W
-  Wclr = V %*% W
-  rownames(Wclr) = colnames(x)
-  colnames(Wclr) <- paste("uwedge", 1:d, sep="")
-  Wclr.inv = gsiInv(Wclr)
-  colnames(Wclr.inv) = colnames(x)
-  colnames(facscores) <- paste("uwedge", 1:d, sep="")
-  sdevs = sqrt(gsi.computeExplainedVariance(facscores, Wclr.inv))
-  out = list(sdev=sdevs, loadings = Wclr, center=mean(cdt(x)), 
-             scale=rep(1, ncol(x)),
-             n.obs=nrow(x), scores = facscores, 
-             gof=res$criter, invLoadings=Wclr.inv,
-             Center = mean(x),
-             InvLoadings = cdtInv(Wclr.inv, orig=x),
-             InvDownLoadings = cdtInv(-Wclr.inv, orig=x),
-             call = cl, type="jointDiag::uwedge")
-  kkk = paste("princomp", class(x), sep=".")
-  class(out)= c("uwedge","genDiag", kkk,  "princomp")
-  return(out)
-}  
-
-#' @describeIn Maf.data.frame  Generalised diagonalisations
-#' @method UWEDGE rcomp
-#' @export
-UWEDGE.rcomp <- UWEDGE.acomp
-
-
-
-
-### RJD calculation -----------------------
-
-#' @describeIn Maf.data.frame  Generalised diagonalisations
-#' @export
-RJD = function(x,...){
-  if(!requireNamespace("JADE", quietly=FALSE)) stop("RJD requires package 'JADE' installed")
-  UseMethod("RJD",x)
-}
-
-#' @describeIn Maf.data.frame  Generalised diagonalisations
-#' @method RJD default
-#' @export
-RJD.default <- function(x, ...){
-  JADE::rjd(X=x, ...)
-} 
-
-
-#' @describeIn Maf.data.frame  Generalised diagonalisations
-#' @method RJD acomp
-#' @export
-RJD.acomp <- function(x, vg, i=NULL,  ...){
-  requireNamespace("compositions", quietly=TRUE)
-  cl <- match.call()
-  vg = as.logratioVariogram(vg)
-  if(is.null(i))
-    i = c(1, floor(c(0.5,1)*dim(vg$vg)[1]))
-  d = ncol(x)-1
-  i = sort(i,decreasing=TRUE)
-  V = compositions::ilrBase(D=d+1)
-  # alternative: V = loadings(princomp(x))
-  M = -0.5*apply(vg$vg[i ,,], 1, clrvar2ilr, V=V )
-  dim(M)=c(d,d,length(i))
-  res = JADE::rjd(M)
-  res$B = t(res$V) 
-  
-  ord = order(diag(res$B %*% M[,,2] %*% t(res$B)))
-  res$B = res$B[ord,]
-  W = t(res$B)
-  # do not normalize
-  # W = unclass(t(normalize(rmult(t(W))))) 
-  W = unclass(W)
-  ilrx = ilr(x, V=V)
-  mns = mean(ilrx)
-  ilrxc = ilrx-mns
-  facscores = unclass(ilrxc) %*% W
-  Wclr = V %*% W
-  rownames(Wclr) = colnames(x)
-  Wclr.inv = t(Wclr)
-  colnames(Wclr.inv) = colnames(x)
-  sdevs = sqrt(gsi.computeExplainedVariance(facscores, Wclr.inv))
-  o = order(sdevs, decreasing = TRUE)
-  Wclr = Wclr[,o]
-  facscores = facscores[,o]
-  Wclr.inv = Wclr.inv[o,]
-  colnames(Wclr) <- rownames(Wclr.inv) <- paste("rjd", 1:d, sep="")
-  colnames(facscores) <- paste("rjd", 1:d, sep="")
-  out = list(sdev=sdevs[o], loadings = Wclr, center=mean(cdt(x)), 
-             scale=rep(1, ncol(x)),
-             n.obs=nrow(x), scores = facscores, 
-             diagonalized = res$D[o,o,], invLoadings=Wclr.inv,
-             Center = mean(x),
-             InvLoadings = cdtInv(Wclr.inv, orig=x),
-             InvDownLoadings = cdtInv(-Wclr.inv, orig=x),
-             call = cl, type="JADE::rjd")
-  kkk = paste("princomp", class(x), sep=".")
-  class(out)= c("rjd","genDiag", kkk,  "princomp")
-  return(out)
-}  
-
-#' @describeIn Maf.data.frame  Generalised diagonalisations
-#' @method RJD rcomp
-#' @export
-RJD.rcomp <- RJD.acomp
-
-
-### genDiag methods ----
-
-#' Predict method for generalised diagonalisation objects
-#'
-#' @param object a generalized diagonalisation object, as obtained from a call to
-#' \code{\link{Maf}}, and on the same page, information on the other diagonalisation
-#' methods \code{UWEDGE} or \code{RJD}
-#' @param newdata a matrix or data.frame of factor scores to convert back to the original
-#' scale (default: the scores element from `object`)
-#' @param ... not used, kept for generic compatibility
-#'
-#' @return A data set or compositional object of the nature of the original data
-#' used for creating the genDiag object.
-#' @export
-#' @family generalised Diagonalisations
-#' @examples
-#' data("jura", package="gstat")
-#' juracomp = compositions::acomp(jura.pred[, -(1:6)]) 
-#' lrvg = logratioVariogram(data=juracomp, loc=jura.pred[,1:2])
-#' mf = Maf(juracomp, vg=lrvg)
-#' mf
-#' biplot(mf)
-#' predict(mf) 
-#' unclass(predict(mf)) - unclass(juracomp) # predict recovers the original composition
-predict.genDiag = function (object, newdata=NULL, ...) {
-  requireNamespace("compositions", quietly=TRUE)
-  if(is.null(newdata))
-    newdata = object$scores
-  if("data.frame" %in% class(newdata))
-    newdata = as.matrix(newdata)
-  Z = newdata %*% unclass(object$invLoadings)
-  if("Center" %in% names(object)){
-    Z = cdtInv(Z, orig = object$Center) + object$Center
-  }else{
-    Z = Z + outer(rep(1, nrow(Z)), object$center)
-    Z = data.frame(Z)
-    colnames(Z) = names(object$center)
-  }
-  return(Z)
-}
-
-
-
-#' Colored biplot for gemeralised diagonalisations
-#' Colored biplot method for objects of class genDiag
-#' @param x a generalized diagonalisation object, as obtained from a call to
-#' \code{\link{Maf}} (or to \code{UWEDGE} or \code{RJD}, on the help page of \code{\link{Maf}}).
-#' @param choices which factors should be represented? vector of 2 indices; defaults to
-#' c(1,2)
-#' @param scale deprecated, kept for coherence with \code{link{biplot.princomp}}
-#' @param pc.biplot same as the homonimous argument from \code{link{biplot.princomp}}:
-#'  boolean, to scale variables down by sqrt(n) and observations up by the same factor. 
-#' @param ... further arguments to \code{\link{coloredBiplot}}
-#'
-#' @return nothing. Function is called exclusively to produce the plot
-#' @export
-#' @family generalised Diagonalisations
-#' @importFrom compositions coloredBiplot
-#' @method coloredBiplot genDiag
-#'
-#' @examples
-#' data("jura", package="gstat")
-#' juracomp = compositions::acomp(jura.pred[, -(1:6)]) 
-#' lrvg = logratioVariogram(data=juracomp, loc=jura.pred[,1:2])
-#' mf = Maf(juracomp, vg=lrvg)
-#' mf
-#' compositions::coloredBiplot(mf, xlabs.col=as.integer(jura.pred$Rock)+2)
-coloredBiplot.genDiag <- function(x, choices = 1:2, scale=0, pc.biplot=FALSE, ...){
-  if(length(choices) != 2) stop("length of choices must be 2")
-  if(!length(scores <- x$scores))
-    stop(gettextf("object '%s' has no scores", deparse(substitute(x))),
-         domain = NA)
-  lam <- rep(1, along.with=choices)
-  if(is.null(n <- x$n.obs)) n <- 1
-  lam <- lam * sqrt(n)
-  if(scale < 0 || scale > 1) warning("'scale' is outside [0, 1]")
-  if(scale != 0) lam <- lam^scale else lam <- 1
-  if(pc.biplot) lam <- lam / sqrt(n)
-  coloredBiplot.default(t(t(scores[,choices]) / lam),
-                         t(x$invLoadings[choices,] * lam), ...)
-  invisible()
-}
+#### MAF calculations (2 versions) ----------------
+# compute the MAF basis from two matrices of (alr/ilr/clr)-covariance 
+gsi.getMAFbasis = function(C0, CT, tol=1e-12){
+  svdT = svd(CT)
+  A1 = with(svdT, u %*% diag(ifelse(d>tol, 1/sqrt(d), 0)) )
+  B = t(A1) %*% C0 %*% A1
+  svdB = svd(B)
+  o = nrow(svdB$u):1  
+  A = A1[,o] %*% svdB$u[o,]
+  return(A)
+}
+
+  
+gsi.computeExplainedVariance <- function(X, Winv){
+  rs = sapply(1:ncol(X), function(i) mvar(rmult(outer(X[,i], Winv[i,]))))
+  return(rs)
+}
+
+
+#' @describeIn Maf.data.frame  Generalised diagonalisations
+#' @param ... generic functionality arguments
+#' @export
+Maf = function(x,...) UseMethod("Maf",x)
+
+#' Generalised diagonalisations
+#' Calculate several generalized diagonalisations out of a data set and its empirical 
+#' variogram
+#' @param x a data set, typically of class "data.frame" or of a compositional class 
+#' @param vg empirical variogram, of a kind fitting to the data provided
+#' @param i a slicer for the variogram, typically this will be one or more indices of 
+#' the lag distance to take. %For other options see code{getEmpVariogramSlice}.
+#'
+#' @return An object extending \code{c("princomp.CLASSOF(x)",}"\code{\link{princomp}}") 
+#' with classes "\code{genDiag}", and an extra class argument depending on the 
+#' diagonalisation method chosen.
+#' 
+#' Function \code{Maf} results carry the extra class "\code{maf}", and they correspond
+#' to minimum/maximum autocorrelation factors (MAF) as proposed by Switzer and Green 
+#' (XXXX) or XXXXX (XXXX). In this case, the
+#' slicer is typically just the index of one lag distance (defaults to i=2). MAF 
+#' provides the analytical solution to the joint diagonalisation of two matrices,
+#' the covariance of increments provided by the slicing and the conventional covariance
+#' matrix (of the idt transformed values, if appropriate). Resulting factors are ordered
+#' in decreasing order of spatial continuity, which might produce surprising 
+#' scree-plots for those who are used to see a screeplot of a principal component analysis.
+#' 
+#' Function \code{UWEDGE} (Uniformly Weighted Exhaustive Diagonalization with Gauss 
+#' iterations) results carry the extra class "\code{uwedge}". The function 
+#' is a wrapper on \code{jointDiag::uwedge} from package \code{jointDiag}. In this case the 
+#' slicer is typically just a subset of indices of lag distances to consider 
+#' (defaults to the nearest indexes to mininum, maximum and mean lag distances of the 
+#' variogram). UWEDGE iteratively seeks for a pair of matrices (a mixing and a
+#' demixing matrices) diagonalises the set of matrices \eqn{M_1, M_2, \ldots, M_K} 
+#' given, by minimizing the target quantity 
+#' \deqn{Q_{uwedge}(A, W) = \sum_{k=1}^K Tr[E_k^t\cdot E_k],}
+#' where \eqn{E_k = (W^t\cdot M_k \cdot W- A\cdot \Lambda_k\cdot A^t)} and 
+#' \eqn{\Lambda_k = diag(W^t\cdot M_k \cdot W)} is the resulting diagonalized version of 
+#' each matrix. Obtained factors are ordered
+#' in decreasing order of spatial continuity, which might produce surprising 
+#' scree-plots for those who are used to see a screeplot of a principal component analysis.
+#' 
+#' Function \code{RJD} results carry the extra class "\code{rjd}". The function 
+#' is a wrapper on \code{JADE::rjd}, implementing the Rotational joint diagonalisation method . In this case the 
+#' slicer is typically just a subset of indices of lag distances to consider 
+#' (defaults to the nearest indexes to mininum, maximum and mean lag distances).
+#' This algorithm also served for quasi-diagonalising a set of matrices as in UWEDGE,
+#' just that in this case the quantity to minimise is the sum of sequares of all off-diagonal
+#' elements of \eqn{A^t\cdot M_k\cdot A} for all \eqn{k=1, 2, \ldots K}.
+#' 
+#' All these functions produce output mimicking \code{\link{princomp}}, i.e. with 
+#' elements
+#' \describe{
+#'   \item{sdev}{contrary to the output in PCA, this contains the square root of the
+#'   metric variance of the predictions obtained for each individual factor; this is the
+#'   quantity needed for \code{\link{screeplot}} to create plots of explained variance
+#'   by factor}
+#'   \item{loadings}{matrix of contributions of each (cdt-transformed) original variable to the new factors}
+#'   \item{center}{center of the data set (eventually, represented through \code{\link{cdt}}), 
+#'   in compositional methods}
+#'   \item{scale}{the scalings applied to each original variable}
+#'   \item{n.obs}{number of observations}
+#'   \item{scores}{the scores of the supplied data on the new factors}
+#'   \item{call}{the call to the function (attention: it actually may come much later)}
+#' } 
+#' and additionally some of the following arguments, in different order
+#' \describe{
+#'   \item{invLoadings}{matrix of contributions of each factor onto each original variable}
+#'   \item{Center}{compositional methods return here the cdt-backtransformed center}
+#'   \item{InvLoadings}{compositional methods return here the clr-backtransformed inverse loadings, so that
+#'   each column of this matrix can be understood as a composition on itself}
+#'   \item{DownInvLoadings}{compositional methods return here the clr-backtransformed "minus inverse loadings", so that
+#'   each column of this matrix can be understood as a composition on itself; details in 
+#'   \code{\link{princomp.acomp}} }
+#'   \item{C1, C2}{Maf returns the two matrices that were diagonalised}
+#'   \item{eigenvalues}{Maf returns the generalized eigenvalues of the diagonalisation of C1 and C2}
+#'   \item{gof}{UWEDGE returns the values of the goodness of fit criterion across sweeps}
+#'   \item{diagonalized}{RJD returns the diagonalized matrices, in an array of (K,D,D)-dimensions, being
+#'   D the number of variables in \code{x}}
+#'   \item{type}{a string describing which package and which function was used as a workhorse for
+#'   the calculation}
+#' }
+#'  
+#' 
+#' NOTE: if the arguments you provide to RJD and UWEDGE are not of the appropriate type
+#' (i.e. data.frames or equivalent) the default method of these functions will just attempt
+#' to call the basis functions JADE:rjd and JointDiag:uwedge respectively. 
+#' This will be the case if you provide \code{x} as a "\code{matrix}", or as
+#' an "\code{array}". In those cases, the output will NOT be structured as an extension
+#' to princomp results; instead they will be native output from those functions.
+#' 
+#'  
+#' @export
+#' @method Maf data.frame
+#' @aliases genDiag
+#' 
+#' @family generalised Diagonalisations
+#'
+#' @examples
+#' require("magrittr")
+#' require("gstat")
+#' require("compositions")
+#' data("jura", package="gstat")
+#' gs = gstat(id="Cd", formula=log(Cd)~1, locations=~Xloc+Yloc, data=jura.pred) %>% 
+#' gstat(id="Co", formula=log(Cd)~1, locations=~Xloc+Yloc, data=jura.pred) %>% 
+#' gstat(id="Cr", formula=log(Cr)~1, locations=~Xloc+Yloc, data=jura.pred) %>% 
+#' gstat(id="Cu", formula=log(Cu)~1, locations=~Xloc+Yloc, data=jura.pred) %>% 
+#' gstat(id="Ni", formula=log(Ni)~1, locations=~Xloc+Yloc, data=jura.pred) %>% 
+#' gstat(id="Pb", formula=log(Pb)~1, locations=~Xloc+Yloc, data=jura.pred) %>% 
+#' gstat(id="Zn", formula=log(Zn)~1, locations=~Xloc+Yloc, data=jura.pred)
+#' vg = variogram(gs)
+#' mf = Maf(aplus(jura.pred[, -(1:6)]), vg=vg)
+#' mf
+#' mf$loadings
+#' biplot(mf)
+Maf.data.frame <- function(x, vg, i=2,...){
+    cl <- match.call()
+  C1 = unclass(var(idt(x)))
+  mn = colMeans(cdt(x))
+  C2 = 0*C1
+  ## BEGIN: to change in the future to use gmEVario 
+  vg = as.gstatVariogram(vg)
+  #dists = unique(vg$dist)
+  #vv = vg[abs(dists[i]-vg$dist)<1e-9*dists[i],]
+  vv = sapply(split(vg$gamma, vg$id), function(x)x[i])
+  vv = vv[order(names(vv))]
+  C2[lower.tri(C2,diag=TRUE)]=vv
+  dd = diag(diag(C2))
+  C2 = unclass(C2 + t(C2) - dd)
+  ## END: to change in the future to use gmEVario 
+  
+  
+  dt.centred = scale(x, center = TRUE, scale=FALSE)
+  eT = eigen(C1)
+  A1 = with(eT, vectors %*% diag(ifelse(values>1e-12, 1/sqrt(values), 0)) ) ## C1^{-0.5}
+  B = t(A1) %*% C2 %*% A1 
+  eB = eigen(B)
+  A2 = eB$vectors
+  
+  eigenvalues = eB$values
+  ord = order(eigenvalues)
+  A2 = A2[, ord]
+  W = unclass(A1) %*% unclass(A2)
+  rownames(W) = colnames(x)
+  # do not normalize
+  # W = unclass(t(normalize(rmult(t(W))))) 
+  W = unclass(W)
+  Winv = solve(W)
+   colnames(Winv) = colnames(W)
+  facscores = unclass(dt.centred) %*% W
+  colnames(facscores) = paste("maf", 1:ncol(facscores), sep="")
+  sdevs = sqrt(gsi.computeExplainedVariance(facscores, Winv))
+  res = list(sdev=sdevs, loadings = W, center=mn, scale=rep(1, ncol(x)),
+             n.obs=nrow(x), scores = facscores, C1=C1, C2=C2, eigenvalues = eigenvalues[ord], 
+             invLoadings=Winv)
+  if(class(x)!="data.frame"){
+    res$Center = cdtInv(mn)
+    res$InvLoadings = cdtInv(Winv, orig=x)
+    res$InvDownLoadings = cdtInv(-Winv, orig=x)
+  }
+  res$call = cl
+  res$type="gmGeostats::MAF"
+  class(res)= c("maf","genDiag",paste("princomp", class(x), sep="."), "princomp")
+  return(res)
+}
+
+#' @describeIn Maf.data.frame  Generalised diagonalisations
+#' @method Maf rmult
+#' @export
+Maf.rmult <- Maf.data.frame
+
+#' @describeIn Maf.data.frame  Generalised diagonalisations
+#' @method Maf aplus
+#' @export
+Maf.aplus <- Maf.data.frame
+
+#' @describeIn Maf.data.frame  Generalised diagonalisations
+#' @method Maf rplus
+#' @export
+Maf.rplus <- Maf.data.frame
+
+#' @describeIn Maf.data.frame  Generalised diagonalisations
+#' @method Maf ccomp
+#' @export
+Maf.ccomp <- Maf.data.frame
+
+
+#' @describeIn Maf.data.frame  Generalised diagonalisations
+#' @method Maf rcomp
+#' @export
+Maf.rcomp  <- function(x, vg, i=2,...){
+  requireNamespace("compositions", quietly=TRUE)
+  cl <- match.call()
+  C1 = unclass(clrvar2ilr(cov(x)))
+  C2 = unclass(clrvar2ilr(variation2clrvar(vg$vg[i,,]))) 
+  clr.centred = scale(cdt(x), center = TRUE, scale=FALSE)
+  
+  eT = eigen(C1)
+  A1 = with(eT, vectors %*% diag(ifelse(values>1e-12, 1/sqrt(values), 0)) ) ## C1^{-0.5}
+  B = t(A1) %*% C2 %*% A1 
+  eB = eigen(B)
+  A2 = eB$vectors
+  
+  eigenvalues = eB$values
+  ord = order(eigenvalues)
+  A2 = A2[, ord]
+  W = unclass(A1) %*% unclass(A2)
+  # do not normalize
+  #W = unclass(t(normalize(rmult(t(W))))) 
+  Wclr = unclass(t(ilr2clr(t(W))))
+  rownames(Wclr) = colnames(x)
+  Winv = solve(W)
+  Wclr.inv = unclass(ilr2clr(Winv))
+    colnames(Wclr.inv) = colnames(x)
+  facscores = unclass(clr.centred) %*% Wclr
+  colnames(facscores) = paste("maf", 1:ncol(facscores), sep="")
+  sdevs = sqrt(gsi.computeExplainedVariance(facscores, Wclr.inv))
+  res = list(sdev=sdevs, loadings = Wclr, center=mean(cdt(x),...), scale=rep(1, ncol(x)),
+             n.obs=nrow(x), scores = facscores, 
+             C1=C1, C2=C2, eigenvalues = eigenvalues[ord], invLoadings=Wclr.inv,
+             Center = mean(x, ...),
+             InvLoadings = cdtInv(Wclr.inv, orig=x),
+             InvDownLoadings = cdtInv(-Wclr.inv, orig=x),
+             call = cl, type="gmGeostats::MAF")
+  kkk = paste("princomp", class(x), sep=".")
+  class(res)= c("maf","genDiag", kkk,  "princomp")
+  return(res)
+} 
+
+#' @describeIn Maf.data.frame  Generalised diagonalisations
+#' @method Maf acomp
+#' @export
+Maf.acomp <- Maf.rcomp
+
+
+
+#### UWEDGE calculations ----------
+
+#' @describeIn Maf.data.frame  Generalised diagonalisations
+#' @export
+UWEDGE = function(x,...){
+  if(!requireNamespace("jointDiag", quietly=FALSE)) stop("UWEDGE requires package 'jointDiag' installed")
+  UseMethod("UWEDGE",x)
+}
+
+#' @describeIn Maf.data.frame  Generalised diagonalisations
+#' @method UWEDGE default
+#' @export
+UWEDGE.default <- function(x, ...) jointDiag::uwedge(M=x, ...)
+
+#' @describeIn Maf.data.frame  Generalised diagonalisations
+#' @method UWEDGE acomp
+#' @export
+UWEDGE.acomp <-  function(x, vg, i=NULL, ...){
+  requireNamespace("compositions", quietly=TRUE)
+  cl <- match.call()
+  vg = as.logratioVariogram(vg) 
+  if(is.null(i)){
+    i = c(1, floor(c(0.5,1)*dim(vg$vg)[1]))
+  }
+  V = loadings(princomp(x))
+  C1 = clrvar2ilr(cov(x), V = V)
+  d = ncol(x)-1
+  M = -0.5*apply(vg$vg[c(1,i),,], 1, clrvar2ilr, V=V )
+  dim(M)=c(d,d,length(i)+1)
+  M[,,1]=C1
+  res = jointDiag::uwedge(M)
+  
+  ord = order(diag(res$B %*% M[,,length(i)] %*% t(res$B)))
+  res$B = res$B[ord,]
+  W = t(res$B)
+  # do not normalize
+  # W = unclass(t(normalize(rmult(t(W))))) 
+  W = unclass(W)
+  ilrx = ilr(x, V=V)
+  mns = mean(ilrx)
+  ilrxc = ilrx-mns
+  facscores = unclass(ilrxc) %*% W
+  Wclr = V %*% W
+  rownames(Wclr) = colnames(x)
+  colnames(Wclr) <- paste("uwedge", 1:d, sep="")
+  Wclr.inv = gsiInv(Wclr)
+  colnames(Wclr.inv) = colnames(x)
+  colnames(facscores) <- paste("uwedge", 1:d, sep="")
+  sdevs = sqrt(gsi.computeExplainedVariance(facscores, Wclr.inv))
+  out = list(sdev=sdevs, loadings = Wclr, center=mean(cdt(x)), 
+             scale=rep(1, ncol(x)),
+             n.obs=nrow(x), scores = facscores, 
+             gof=res$criter, invLoadings=Wclr.inv,
+             Center = mean(x),
+             InvLoadings = cdtInv(Wclr.inv, orig=x),
+             InvDownLoadings = cdtInv(-Wclr.inv, orig=x),
+             call = cl, type="jointDiag::uwedge")
+  kkk = paste("princomp", class(x), sep=".")
+  class(out)= c("uwedge","genDiag", kkk,  "princomp")
+  return(out)
+}  
+
+#' @describeIn Maf.data.frame  Generalised diagonalisations
+#' @method UWEDGE rcomp
+#' @export
+UWEDGE.rcomp <- UWEDGE.acomp
+
+
+
+
+### RJD calculation -----------------------
+
+#' @describeIn Maf.data.frame  Generalised diagonalisations
+#' @export
+RJD = function(x,...){
+  if(!requireNamespace("JADE", quietly=FALSE)) stop("RJD requires package 'JADE' installed")
+  UseMethod("RJD",x)
+}
+
+#' @describeIn Maf.data.frame  Generalised diagonalisations
+#' @method RJD default
+#' @export
+RJD.default <- function(x, ...){
+  JADE::rjd(X=x, ...)
+} 
+
+
+#' @describeIn Maf.data.frame  Generalised diagonalisations
+#' @method RJD acomp
+#' @export
+RJD.acomp <- function(x, vg, i=NULL,  ...){
+  requireNamespace("compositions", quietly=TRUE)
+  cl <- match.call()
+  vg = as.logratioVariogram(vg)
+  if(is.null(i))
+    i = c(1, floor(c(0.5,1)*dim(vg$vg)[1]))
+  d = ncol(x)-1
+  i = sort(i,decreasing=TRUE)
+  V = loadings(princomp(x))
+  M = -0.5*apply(vg$vg[i ,,], 1, clrvar2ilr, V=V )
+  dim(M)=c(d,d,length(i))
+  res = JADE::rjd(M)
+  res$B = t(res$V) 
+  
+  ord = order(diag(res$B %*% M[,,2] %*% t(res$B)))
+  res$B = res$B[ord,]
+  W = t(res$B)
+  # do not normalize
+  # W = unclass(t(normalize(rmult(t(W))))) 
+  W = unclass(W)
+  ilrx = ilr(x, V=V)
+  mns = mean(ilrx)
+  ilrxc = ilrx-mns
+  facscores = unclass(ilrxc) %*% W
+  Wclr = V %*% W
+  rownames(Wclr) = colnames(x)
+  Wclr.inv = t(Wclr)
+  colnames(Wclr.inv) = colnames(x)
+  sdevs = sqrt(gsi.computeExplainedVariance(facscores, Wclr.inv))
+  o = order(sdevs, decreasing = TRUE)
+  Wclr = Wclr[,o]
+  facscores = facscores[,o]
+  Wclr.inv = Wclr.inv[o,]
+  colnames(Wclr) <- rownames(Wclr.inv) <- paste("rjd", 1:d, sep="")
+  colnames(facscores) <- paste("rjd", 1:d, sep="")
+  out = list(sdev=sdevs[o], loadings = Wclr, center=mean(cdt(x)), 
+             scale=rep(1, ncol(x)),
+             n.obs=nrow(x), scores = facscores, 
+             diagonalized = res$D[o,o,], invLoadings=Wclr.inv,
+             Center = mean(x),
+             InvLoadings = cdtInv(Wclr.inv, orig=x),
+             InvDownLoadings = cdtInv(-Wclr.inv, orig=x),
+             call = cl, type="JADE::rjd")
+  kkk = paste("princomp", class(x), sep=".")
+  class(out)= c("rjd","genDiag", kkk,  "princomp")
+  return(out)
+}  
+
+#' @describeIn Maf.data.frame  Generalised diagonalisations
+#' @method RJD rcomp
+#' @export
+RJD.rcomp <- RJD.acomp
+
+
+### genDiag methods ----
+
+#' Predict method for generalised diagonalisation objects
+#'
+#' @param object a generalized diagonalisation object, as obtained from a call to
+#' \code{\link{Maf}}, and on the same page, information on the other diagonalisation
+#' methods \code{UWEDGE} or \code{RJD}
+#' @param newdata a matrix or data.frame of factor scores to convert back to the original
+#' scale (default: the scores element from `object`)
+#' @param ... not used, kept for generic compatibility
+#'
+#' @return A data set or compositional object of the nature of the original data
+#' used for creating the genDiag object.
+#' @export
+#' @family generalised Diagonalisations
+#' @examples
+#' data("jura", package="gstat")
+#' juracomp = compositions::acomp(jura.pred[, -(1:6)]) 
+#' lrvg = logratioVariogram(data=juracomp, loc=jura.pred[,1:2])
+#' mf = Maf(juracomp, vg=lrvg)
+#' mf
+#' biplot(mf)
+#' predict(mf) 
+#' unclass(predict(mf)) - unclass(juracomp) # predict recovers the original composition
+predict.genDiag = function (object, newdata=NULL, ...) {
+  requireNamespace("compositions", quietly=TRUE)
+  if(is.null(newdata))
+    newdata = object$scores
+  if("data.frame" %in% class(newdata))
+    newdata = as.matrix(newdata)
+  Z = newdata %*% unclass(object$invLoadings)
+  if("Center" %in% names(object)){
+    Z = cdtInv(Z, orig = object$Center) + object$Center
+  }else{
+    Z = Z + outer(rep(1, nrow(Z)), object$center)
+    Z = data.frame(Z)
+    colnames(Z) = names(object$center)
+  }
+  return(Z)
+}
+
+
+
+#' Colored biplot for gemeralised diagonalisations
+#' Colored biplot method for objects of class genDiag
+#' @param x a generalized diagonalisation object, as obtained from a call to
+#' \code{\link{Maf}} (or to \code{UWEDGE} or \code{RJD}, on the help page of \code{\link{Maf}}).
+#' @param choices which factors should be represented? vector of 2 indices; defaults to
+#' c(1,2)
+#' @param scale deprecated, kept for coherence with \code{link{biplot.princomp}}
+#' @param pc.biplot same as the homonimous argument from \code{link{biplot.princomp}}:
+#'  boolean, to scale variables down by sqrt(n) and observations up by the same factor. 
+#' @param ... further arguments to \code{\link{coloredBiplot}}
+#'
+#' @return nothing. Function is called exclusively to produce the plot
+#' @export
+#' @family generalised Diagonalisations
+#' @importFrom compositions coloredBiplot
+#' @method coloredBiplot genDiag
+#'
+#' @examples
+#' data("jura", package="gstat")
+#' juracomp = compositions::acomp(jura.pred[, -(1:6)]) 
+#' lrvg = logratioVariogram(data=juracomp, loc=jura.pred[,1:2])
+#' mf = Maf(juracomp, vg=lrvg)
+#' mf
+#' compositions::coloredBiplot(mf, xlabs.col=as.integer(jura.pred$Rock)+2)
+coloredBiplot.genDiag <- function(x, choices = 1:2, scale=0, pc.biplot=FALSE, ...){
+  if(length(choices) != 2) stop("length of choices must be 2")
+  if(!length(scores <- x$scores))
+    stop(gettextf("object '%s' has no scores", deparse(substitute(x))),
+         domain = NA)
+  lam <- rep(1, along.with=choices)
+  if(is.null(n <- x$n.obs)) n <- 1
+  lam <- lam * sqrt(n)
+  if(scale < 0 || scale > 1) warning("'scale' is outside [0, 1]")
+  if(scale != 0) lam <- lam^scale else lam <- 1
+  if(pc.biplot) lam <- lam / sqrt(n)
+  coloredBiplot.default(t(t(scores[,choices]) / lam),
+                         t(x$invLoadings[choices,] * lam), ...)
+  invisible()
+}
diff --git a/R/geostats.R b/R/geostats.R
index 90f0abb..8e1ae85 100644
--- a/R/geostats.R
+++ b/R/geostats.R
@@ -207,22 +207,24 @@ gsi.CondTurningBands <- function(Xout, Xin, Zin, vgram,
 ### tests --
 # if(!exists("do.test")) do.test=FALSE
 # if(do.test){
-#   library("gstat")
-#   library("magrittr")
-#   data("jura", package = "gstat")
-#   dt = jura.pred %>% dplyr::select(Cd:Zn)
-#   X = jura.pred[,1:2]
-#   spc = SpatialPointsDataFrame(SpatialPoints(X),acomp(dt))
-#   a2 = make.gmCompositionalGaussianSpatialModel(spc)
-#   #attach(gsi)
-#   a3gs = as.gstat(a2)
-#   a3vg = variogram(a3gs) 
-#   a3gs = fit_lmc(a3vg, a3gs, vgm("Sph", psill=1, nugget=0, range=1))
-#   a4 = make.gmCompositionalGaussianSpatialModel(spc, V="alr", model=a3gs$model)
-#   a4gs = as.gstat(a4)
-#   p1 = predict(a4gs, newdata = jura.grid)
-#   p2 = predict(a4gs, newdata = SpatialPoints(jura.grid[,1:2]))
-#   p3 = predict(a4gs, newdata = SpatialPixels(SpatialPoints(jura.grid[,1:2])))
+  # library("compositions")
+  # library("sp")
+  # library("magrittr")
+  # library("gstat")
+  # data("jura", package = "gstat")
+  # dt = jura.pred %>% dplyr::select(Cd:Cr)
+  # X = jura.pred[,1:2]
+  # spc = SpatialPointsDataFrame(SpatialPoints(X),acomp(dt))
+  # a2 = make.gmCompositionalGaussianSpatialModel(spc)
+  # #attach(gsi)
+  # a3gs = as.gstat(a2)
+  # a3vg = variogram(a3gs)
+  # a3gs = fit_lmc(a3vg, a3gs, vgm("Sph", psill=1, nugget=0, range=1))
+  # a4 = make.gmCompositionalGaussianSpatialModel(spc, V="alr", model=a3gs$model)
+  # a4gs = as.gstat(a4)
+  # p1 = predict(a4gs, newdata = jura.grid)
+  # p2 = predict(a4gs, newdata = jura.grid[,1:2])
+  # p3 = predict(a4gs, newdata = SpatialPixels(SpatialPoints(jura.grid[,1:2])))
 #   
 # }
 
diff --git a/R/grids.R b/R/grids.R
index 5abfc1f..b1d396e 100644
--- a/R/grids.R
+++ b/R/grids.R
@@ -152,6 +152,7 @@ gsi.gstatCokriging2compo.default <- function(COKresult, ...){
 #' @param tol for generalized inversion of the matrix (**rarely touched!**)
 #' @param nscore boolean, were the data normal score-transformed? (**deprecated**)
 #' @param gg in the case that normal score transformation was applied, provide the gstat object! (**deprecated**)
+#' @export
 gsi.gstatCokriging2compo.data.frame = function(COKresult, # output of predict.gstat
                                                V=NULL,    # string or matrix describing which logratio was applied
                                                orignames=NULL, # names of the original components (OPTIONAL)
@@ -171,7 +172,7 @@ gsi.gstatCokriging2compo.data.frame = function(COKresult, # output of predict.gs
   if(is.null(orignames)) orignames = paste("v", 1:D, sep="")
   if(length(orignames)!=D) stop("names provided not consistent with number of logratio variables. Did you forget the rest?")
   prefix = "ilr"
-  prediccions = COKresult[,prednames]
+  prediccions = COKresult[,prednames, drop=FALSE]
   if(is.character(V)){
     if(V=="ilr"){
       V = ilrBase(prediccions)
@@ -205,12 +206,12 @@ gsi.gstatCokriging2compo.data.frame = function(COKresult, # output of predict.gs
     # add cokriging variance matrices as an attribute as well
     noms = sub(".pred", "", prednames)
     cvmat = array(0, dim=c(nrow(rg), D-1, D-1), dimnames=list(NULL, noms, noms))
-    vrs = COKresult[,varnames]
+    vrs = COKresult[,varnames, drop=FALSE]
     colnames(vrs) = sub(".var", "", varnames)
     for(ivr in noms){
       cvmat[ ,ivr, ivr] = vrs[,ivr]
     }
-    cvs = COKresult[,covnames]
+    cvs = COKresult[,covnames, drop=FALSE]
     colnames(cvs) = sub("cov.", "", covnames)
     for(ivr in noms){
       for(jvr in noms){
@@ -234,6 +235,7 @@ gsi.gstatCokriging2rmult <- function(COKresult, ...) UseMethod("gsi.gstatCokrigi
 
 #' @describeIn gsi.gstatCokriging2compo Reorganisation of cokriged multivariate data
 #' @method gsi.gstatCokriging2rmult default
+#' @export
 gsi.gstatCokriging2rmult.default <- function(COKresult, ...){
   stopifnot(is(COKresult, "Spatial"), "data" %in% slotNames(COKresult))
   coord = sp::coordinates(COKresult)
@@ -251,6 +253,7 @@ gsi.gstatCokriging2rmult.default <- function(COKresult, ...){
 
 #' @describeIn gsi.gstatCokriging2compo Reorganisation of cokriged multivariate data
 #' @method gsi.gstatCokriging2rmult data.frame
+#' @export
 gsi.gstatCokriging2rmult.data.frame = function(COKresult, # output of predict.gstat
                                     nscore=FALSE, # were the data NS-transformed?
                                     gg=NULL, # in the case that NS was applied, provide the gstat object!
@@ -264,7 +267,7 @@ gsi.gstatCokriging2rmult.data.frame = function(COKresult, # output of predict.gs
   D = length(prednames)
   noms = sub(".pred", "", prednames)
   
-  prediccions = COKresult[,prednames]
+  prediccions = COKresult[,prednames, drop=FALSE]
   if(nscore){
     ## space to back-transform the predictions
     #if(is.null(gg))stop("To apply a nscore backtransformation, the gstat object must be provided!")
@@ -283,12 +286,12 @@ gsi.gstatCokriging2rmult.data.frame = function(COKresult, # output of predict.gs
   if(!nscore){
     # add cokriging variance matrices as an attribute as well
     cvmat = array(0, dim=c(nrow(rg), D, D), dimnames=list(NULL, noms, noms))
-    vrs = COKresult[,varnames]
+    vrs = COKresult[,varnames, drop=FALSE]
     colnames(vrs) = sub(".var", "", varnames)
     for(ivr in noms){
       cvmat[ ,ivr, ivr] = vrs[,ivr]
     }
-    cvs = COKresult[,covnames]
+    cvs = COKresult[,covnames, drop=FALSE]
     colnames(cvs) = sub("cov.", "", covnames)
     for(ivr in noms){
       for(jvr in noms){
-- 
GitLab