variograms.R

#### theoretical variogram ---------------------

#' Generate D-variate variogram models
#'
#' @description Function to set up D-variate variogram models based on model type, the variogram parameters sill and nugget and a matrix describing the anisotropy of the range.
#'
#' @param type model of correlation function. The function expects a constant, e.g. the internal constants 'vg.Gau' for Gaussian model or 'vg.Exp'. for exponential models. See examples for usage.
#' @param nugget (DxD)-matrix for the nugget effect. Default is a muted nugget (0).
#' @param sill (DxD)-matrix for the partial sills of the correlation function
#' @param anisRanges 2x2 or 3x3 matrix of ranges (see details)
#' @param extraPar for certain correlation functions, extra parameters (smoothness, period, etc)
#'
#' @return an object of class "gmCgram" containing the linear model of corregionalization
#' of the nugget and the structure given.
#' @details The argument `type` must be an integer indicating the model to be used.
#' Some constants are available to make reading code more understandable. That is, you can
#' either write `1`, `vg.sph`, `vg.Sph` or `vg.Spherical`, they will all work and produce
#' a spherical model. The same applies for the following models:
#' `vg.Gauss = vg.Gau = vg.gau = 0`;
#' `vg.Exponential = vg.Exp = vg. exp = 2`.
#' These constants are available after calling `data("variogramModels")`.
#' No other model is currently available, but this data object will be
#' regularly updated.
#' The constant vector `gsi.validModels` contains all currently valid models.
#'
#' Argument `anisRange` expects a matrix $M$ such that
#' \deqn{
#' h^2 = (\mathbf{x}_i-\mathbf{x}_j)\cdot M^{-1}\cdot (\mathbf{x}_i-\mathbf{x}_j)^t
#' }
#' is the (square of) the lag distance to be fed into the correlation function.
#' @family gmCgram
#' @export
#' @aliases vg.Exp vg.exp vg.Exponential vg.Gau vg.gauss
#' vg.Gauss vg.Sph vg.sph vg.Spherical gsi.validModels
#'
#' @examples
#' utils::data("variogramModels") # shortcut for all model constants
#' v1 = setCgram(type=vg.Gau, sill=diag(2), anisRanges = 3*diag(c(3,1)))
#' v2 = setCgram(type=vg.Exp, sill=0.3*diag(2), anisRanges = 0.5*diag(2))
#' vm = v1+v2
#' plot(vm)
setCgram = function(type, nugget=sill*0, sill, anisRanges, extraPar=0){
  utils::data("variogramModels")
  stopifnot(all(dim(sill)==dim(nugget)),
            ncol(sill)==nrow(sill),
            type %in% gsi.validModels)
  dim(sill) = c(1,dim(sill))
  dim(anisRanges) = c(1,dim(anisRanges))
  vgout <- list(type=type,
                data=extraPar,
                nugget=nugget,
                sill=sill, #(nstru, nvar, nvar)
                M=AnisotropyRangeMatrix(anisRanges)     # these are "classical" ranges (i.e. distances)
  )
  class(vgout) = "gmCgram"
  return(vgout)
}



#' Subsetting of gmCgram variogram structures
#'
#' Extraction or combination of nested structures of a gmCgram object
#'
#' @param x `gmCgram` variogram object
#' @param i indices of the structures that are desired to be kept (0=nugget) or removed (see details)
#' @param ... extra arguments for generic functionality
#'
#' @return a `gmCgram` variogram object with the desired structures only.
#' @details This function can be used to: extract the nugget (i=0), extract some
#' structures (i=indices of the structures, possibly including 0 for the nugget),
#' or filter some structures out (i=negative indices of the structures to remove;
#' nugget will always removed in this case). If you want to extract "slots" or
#' "elements" of the variogram, use the $-notation. If you want to extract variables of the
#' variogram matrices, use the `[`-notation. The contrary operation (adding structures together)
#' is obtained by summing (+) two `gmCgram` objects.
#' @export
#' @method [[ gmCgram
#' @family gmCgram functions
#' @examples
#' utils::data("variogramModels")
#' v1 = setCgram(type=vg.Gau, sill=diag(2), anisRanges = 3*diag(c(3,1)))
#' v2 = setCgram(type=vg.Exp, sill=0.3*diag(2), anisRanges = 0.5*diag(2))
#' vm = v1+v2
#' vm[[1]]
"[[.gmCgram"<- function(x,i,...){
  nG = dim(x$M)[3]
  nD = dim(x$nugget)[1]
  nSo = length(i)
  nSi = dim(x$M)[1]
  if(i==0){# case only nugget wanted
    out = with(x, list(type=type[i],
                       data=data[i],
                       nugget=nugget,
                       sill=structure(0*sill[1,,], dim=c(nSo, nD, nD)),
                       M=structure(M[1,,], dim=c(nSo, nG, nG))))
  }else if( (0 %in% i) & !any(i<0)){
    j = i[i>0] # case some structures wanted, including nugget
    out = with(x, list(type=type[j],
                       data=data[j],
                       nugget=nugget,
                       sill=structure(sill[j,,], dim=c(nSo-1, nD, nD)),
                       M=structure(M[j,,], dim=c(nSo-1, nG, nG))))
  }else if(all(i>0) | all(i<0)){# case some structured (un)wanted, nugget surely not wanted
    if(all(i<0)) nSo = nSi-nSo
    out = with(x, list(type=type[i],
                       data=data[i],
                       nugget=nugget*0,
                       sill=structure(sill[i,,], dim=c(nSo, nD, nD)),
                       M=structure(M[i,,], dim=c(nSo, nG, nG))))
  }else{# unsolvable case
    stop("index set i cannot merge 0 and negative numbers")
  }
  class(out) = class(x)
  return(out)
}


#' Combination of gmCgram variogram structures
#'
#' combination of nested structures of a gmCgram object
#'
#' @param x `gmCgram` variogram object
#' @param y `gmCgram` variogram object
#' @export
#' @return The combined nested structures
#' @method + gmCgram
#' @examples
#' utils::data("variogramModels")
#' v1 = setCgram(type=vg.Gau, sill=diag(2), anisRanges = 3*diag(c(3,1)))
#' v2 = setCgram(type=vg.Exp, sill=0.3*diag(2), anisRanges = 0.5*diag(2))
#' vm = v1+v2
"+.gmCgram" <- function(x,y) {
  y = as.gmCgram(y)
  stopifnot(is(y,"gmCgram"),
            dim(x$sill)[-1]==dim(y$sill)[-1],
            dim(x$M)[-1]==dim(y$M)[-1])
  myfun = function(A,B){
    D = dim(A)[2]
    nA = dim(A)[1]
    nB = dim(B)[1]
    dim(A) = c(nA, D^2)
    dim(B) = c(nB, D^2)
    out = rbind(A, B)
    dim(out) = c(nA+nB, D, D)
    return(out)
  }
  x$type = c(x$type, y$type)
  x$data= c(x$data, y$data)
  x$nugget = x$nugget + y$nugget
  x$sill = myfun(x$sill, y$sill)
  x$M = myfun(x$M, y$M)
  return(x)
}




#' Subsetting of gmCgram variogram structures
#'
#' Extraction of some variables of a gmCgram object
#'
#' @param x \code{gmCgram} variogram object
#' @param i row-indices of the variables to be kept/removed
#' @param j column-indices of the variables to be kept/removed (if only \code{i}
#' is specified, \code{j} will be taken as equal to \code{i}!)
#' @param ...  extra arguments for generic functionality
#'
#' @return a \code{gmCgram} variogram object with the desired variables only.
#' @details This function can be used to extract the model for a a subset of variables.
#' If only \code{i} is specified, \code{j} will be taken as equal to \code{i}.
#' If you want to select all \code{i}'s for certain \code{j}'s or vice versa, give
#' \code{i=1:dim(x$nugget)[1]} and \code{j=} your desired indices, respectively
#' \code{j=1:dim(x$nugget)[2]} and \code{i=} your desired indices; replace \code{x} by the
#' object you are giving. If \code{i!=j}, the output will be a \code{c("gmXCgram","gmCgram")}
#' object, otherwise it will be a regular class \code{"gmCgram"} object.
#' If you want to extract "slots" or
#' "elements" of the variogram, use the $-notation. If you want to extract variables of the
#' variogram matrices, use the `[`-notation.
#' @export
#' @method [ gmCgram
#' @family gmCgram functions
#' @examples
#' utils::data("variogramModels")
#' v1 = setCgram(type=vg.Gau, sill=diag(2), anisRanges = 3*diag(c(3,1)))
#' v2 = setCgram(type=vg.Exp, sill=0.3*diag(2), anisRanges = 0.5*diag(2))
#' vm = v1+v2
#' vm[1,1]
"[.gmCgram"<- function(x,i,j=i,...){
  nDi = length(i)
  nDj = length(j)
  nS = dim(x$M)[1]
  out = with(x, list(type=type,
                     data=data,
                     nugget=structure(nugget[i,j, drop=F], dim=c(nDi, nDj)),
                     sill=structure(sill[,i,j, drop=F], dim=c(nS, nDi, nDj)),
                     M=M))
  class(out) = class(x)
  if(!all(i==j)) class(out) = unique(c("gmXCgram", class(out)))
  return(out)
}

#' Length, and number of columns or rows
#'
#' Provide number of structures, and nr of variables of an LMC of class gmCgram
#'
#' @param x gmCgram object
#'
#' @return \code{length} returns the number of structures (nugget not counted), while
#' \code{ncol} and \code{nrow} return these values for the nugget (assuming that they will
#' be also valid for the sill).
#' @export
#' @aliases ncol.gmCgram nrow.gmCgram
#' @family gmCgram functions
#'
#' @method length gmCgram
#' @examples
#' utils::data("variogramModels")
#' v1 = setCgram(type=vg.Gau, sill=diag(3)+0.5, anisRanges = 2*diag(c(3,0.5)))
#' v2 = setCgram(type=vg.Exp, sill=0.3*diag(3), anisRanges = 0.5*diag(2))
#' vm = v1+v2
#' length(vm)
#' ncol(vm)
#' nrow(vm)
length.gmCgram = function(x)  length(x$type)
if(!isGeneric("ncol")){
  ncol <- function(x) UseMethod("ncol",x)
  ncol.default <- base::ncol
}
if(!isGeneric("nrow")){
  nrow <- function(x) UseMethod("nrow",x)
  nrow.default <- base::nrow
}
ncol.gmCgram = function(x) ncol(x$nugget)
nrow.gmCgram = function(x) nrow(x$nugget)


#' Convert a gmCgram object to an (evaluable) function
#'
#' Evaluate a gmCgram on some h values, or convert the gmCgram object into an evaluable function
#'
#' @param x a gmCgram object
#' @param ... extra arguments for generic functionality
#'
#' @return a \code{function} that can be evaluated normally, with an argument \code{X}
#' and possibly another argument \code{Y}; both must have the same number of columns
#' than the geographic dimension of the variogram (i.e. \code{dim(x$M)[3]}).
#' @export
#' @method as.function gmCgram
#' @family gmCgram functions
#' @examples
#' utils::data("variogramModels")
#' v1 = setCgram(type=vg.Gau, sill=diag(2)+0.5, anisRanges = 2*diag(c(3,0.5)))
#' v2 = setCgram(type=vg.Exp, sill=0.3*diag(2), anisRanges = 0.5*diag(2))
#' vm = v1+v2
#' vgf = as.function(vm)
#' (h = rbind(c(0,1), c(0,0), c(1,1)))
#' vgf(h)
#' predict(vm, h)
as.function.gmCgram = function(x,...){
    f <- function(X,Y=X){
      if(is(X,"Spatial")) X = sp::coordinates(X)
      X = as.matrix(X)
      if(is(Y,"Spatial")) Y = sp::coordinates(Y)
      Y = as.matrix(Y)
      stopifnot(ncol(X)==ncol(Y), ncol(X)==dim(x$M)[3])
      ijEqual = ifelse(nrow(X)==nrow(Y), all(X==Y), FALSE)
      o = gsi.calcCgram(X,Y,x,ijEqual)
      return(o)
    }
  return(f)
}

#' @describeIn as.function.gmCgram predict a gmCgram object on some lag vector coordinates
#' @param object gmCgram object
#' @param newdata matrix, data.frame or Spatial object containing coordinates
#' @include gmAnisotropy.R
#' @method predict gmCgram
#' @export
predict.gmCgram = function(object, newdata, ...){
  as.function(object)(X=newdata)
}


#' Convert theoretical structural functions to gmCgram format
#'
#' Convert covariance function or variogram models to the format gmCgram
#' of package gmGeostats
#'
#' @param m model to be converted
#' @param ... further parameters
#'
#' @return the covariance/variogram model, recasted to class \code{gmCgram}.
#' This is a generic function. Methods exist for objects of class
#' \code{LMCAnisCompo} (for compositional data) and \code{variogramModelList}
#' (as provided by package \code{gstat}).
#' @export
#' @family gmCgram functions
as.gmCgram <- function(m, ...)  UseMethod("as.gmCgram",m)

#' @describeIn as.gmCgram Convert theoretical structural functions to gmCgram format
#' @method as.gmCgram default
#' @export
as.gmCgram.default <- function(m,...) m



#' Draw cuves for covariance/variogram models
#'
#' Represent a gmCgram object as a matrix of lines in several plots
#'
#' @param x object to draw, of class gmCgram // curently only valid for symmetric functions
#' @param xlim.up range of lag values to use in plots of the upper triangle
#' @param xlim.lo range of lag values to use in plots of the lower triangle
#' @param vdir.up geograohic directions to represent in the upper triangle
#' @param vdir.lo geograohic directions to represent in the lower triangle
#' @param xlength number of discretization points to use for the curves (defaults to 200)
#' @param varnames string vector, variable names to use in the labelling of axes
#' @param add logical, should a new plot be created or stuff be added to an existing one?
#' @param commonAxis logical, is a common Y axis for all plots in a row desired?
#' @param cov logical, should the covariance function (=TRUE) or the variogram (=FALSE) be plotted?
#' @param closeplot logical, should the plot be left open (FALSE) for further changes, or be frozen (TRUE)?
#' defaults to TRUE
#' @param ... further graphical parameters for the plotting function
#'
#' @return This function is called for its side effect of producing a plot: the plot will be
#' open to further changes if you provide `closeplot=FALSE`. Additionally, the function
#' invisibly returns the graphical parameters that were active before starting the plot. Hence,
#' if you want to freeze a plot and not add anymore to it, you can do \code{par(plot(x, closeplot=FALSE, ...))},
#' or \code{plot(x, closeplot=TRUE, ...)}.
#' If you want to further add stuff to it, better just call \code{plot(x, closeplot=FALSE,...)}. The difference
#' is only relevant when working with the screen graphical device.
#' @export
#' @method plot gmCgram
#' @family gmCgram functions
#' @examples
#' utils::data("variogramModels")
#' v1 = setCgram(type=vg.Gau, sill=diag(3)-0.5, anisRanges = 2*diag(c(3,0.5)))
#' v2 = setCgram(type=vg.Exp, sill=0.3*diag(3), anisRanges = 0.5*diag(2))
#' vm = v1+v2
#' plot(vm)
#' plot(vm, cov=FALSE)
plot.gmCgram = function(x, xlim.up=NULL, xlim.lo=NULL, vdir.up= NULL, vdir.lo= NULL, xlength=200, varnames = colnames(x$nugget),
                        add=FALSE, commonAxis=TRUE, cov =TRUE, closeplot=TRUE, ...){
  Dg = dim(x$M)[2]
  Ns = dim(x$M)[1]
  Dv = dim(x$nugget)[1]
  if(is.null(varnames)) varnames = paste("v", 1:Dv, sep="")
  if(is.null(vdir.up) & is.null(vdir.lo)){
    vdir.lo = rep(0, Dg)
    vdir.lo[1] = 1
    dim(vdir.lo) = c(1,Dg)
    if(!is.isotropic(x)){
      aux = diag(Dg)
      if(Dg==3){
        vdir.lo = aux[3,]
        vdir.up = aux[-3,]
      }else
        vdir.lo = aux
    }
  }
  if(!is.null(vdir.up)) vdir.up = compositions::oneOrDataset(compositions::normalize(vdir.up))
  if(!is.null(vdir.lo)) vdir.lo = compositions::oneOrDataset(compositions::normalize(vdir.lo))
  fk = c(sqrt(3),1,3)[x$type+1]*1.25  # range to effective range factor expansion
  if(is.null(xlim.up)){
    if(add){
      par(mfg=c(1,2))
      xlim.up = par()$usr[1:2]
    }else if(!is.null(vdir.up)){
      maxdist = sapply(1:Ns, function(i) max(sapply(1:nrow(vdir.up), function(j) vdir.up[j,]%*%x$M[i,,]%*%vdir.up[j,])))
      # here we must compute the max M projected on vd.up
      xlim.up = c(0, max(fk*maxdist))
    }
  }
  if(is.null(xlim.lo)){
    if(add){
      par(mfg=c(2,1))
      xlim.lo = par()$usr[1:2]
    }else if(!is.null(vdir.lo)){
      maxdist = sapply(1:Ns, function(i) max(sapply(1:nrow(vdir.lo), function(j) vdir.lo[j,]%*%x$M[i,,]%*%vdir.lo[j,])))
      # here we must compute the max M projected on vd.lo
      xlim.lo = c(0, max(fk*maxdist))
    }
  }
  if(!is.null(xlim.up)){
    xseq.up = seq(from=xlim.up[1], to=xlim.up[2], length.out=xlength)
  }else{ xseq.up=NULL}
  if(!is.null(xlim.lo)){
    xseq.lo = seq(from=xlim.lo[1], to=xlim.lo[2], length.out=xlength)
  }else{ xseq.lo=NULL}

  opar = par()
  opar = par_remove_readonly(opar)

  if(closeplot) on.exit(par(opar))

  getVdens = function(vdir, xseq){
    if(is.null(vdir)|is.null(xseq)) return(NULL)
    Vdens = sapply(1:nrow(vdir), function(k){
      X = outer(xseq, vdir[k,])
      Y = X[1,,drop=F]*0
      gsi.calcCgram(X,Y,x,FALSE)
    })
    dim(Vdens) = c(Dv,xlength,Dv,nrow(vdir))
    if(!cov){
      # convert to variogram if cov=FALSE
      Y = matrix(rep(0,Dg), ncol=Dg)
      C0 = gsi.calcCgram(Y,Y,x,FALSE)
      dim(C0) = c(Dv,Dv)
      Vdens = sweep(-Vdens, c(1,3), C0, "+") ## this must be corrected when we allow non-symmetric covariances
    }
    return(Vdens)
  }
  Vdens.up = getVdens(vdir.up, xseq.up)
  Vdens.lo = getVdens(vdir.lo, xseq.lo)

  myplot = function(...) matplot(type="l",ylab="", xlab="",xaxt="n", ...)
  if(add) myplot = function(...) matlines(...)
  if(!add){
    par(mfrow=c(Dv+1,Dv+1), mar=c(2,3,0,0), oma=c(1,4,1,1), xpd=NA)
    myplot(c(0,0), c(0,0), pch="", ann=FALSE, bty="n", yaxt="n")
  }
  for(i in 1:Dv){
      for(j in 1:Dv){
        if((i>=j)&!(is.null(vdir.lo)|is.null(xlim.lo))){
            par(mfg=c(i+1,j,Dv+1,Dv+1))
            ylim = range(Vdens.lo[,,j,])
            if(commonAxis) ylim=range(Vdens.lo[,,j,])
            myplot(xseq.lo, Vdens.lo[i,,j,], ylim=ylim, ...)
            if(i==j & !add){
              axis(side = 3)
              mtext(text=varnames[i], side = 3, line=3)
              mtext(text=varnames[i], side = 4, line=3)
            }
        }
        if((i<=j)&!(is.null(vdir.up)|is.null(xlim.up))){
            par(mfg=c(i,j+1,Dv+1,Dv+1))
            ylim = range(Vdens.up[,,j,])
            if(commonAxis) ylim=range(Vdens.up[,,j,])
            myplot(xseq.up, Vdens.up[i,,j,], ylim=ylim,...)
            if(i==j & !add){
              axis(side = 1)
              mtext(text=varnames[i], side = 1, line=3)
              mtext(text=varnames[i], side = 2, line=3)
            }
        }
      }}
  mtext(text="lag distance", side=1, outer = TRUE, line=0)
  mtext(text=c("semivariogram","covariance")[cov+1], side=2, outer = TRUE, line=2)
  invisible(opar)
}





#' Check for anisotropy of a theoretical variogram
#'
#' Checks for anisotropy of a theoretical variogram or covariance function model

#' @param x variogram or covariance model object
#' @param tol tolerance for
#' @param ... extra arguments for generic  functionality
#'
#' @return Generic function. Returns of boolean answering the question of the name,
#' or NA if object \code{x} does not contain a known theoretical variogram
#' @export
is.isotropic <- function(x, tol=1e-10, ...){ UseMethod("is.isotropic", x) }

#' @method is.isotropic default
#' @export
is.isotropic.default = function(x, tol=1e-10, ...) NA

#' @method is.isotropic gmCgram
#' @export
is.isotropic.gmCgram = function(x, tol=1e-10, ...){
  all(apply(x$M, 1, function(y){
    ev = eigen(y, only.values=TRUE)[[1]]
    all(abs(ev-ev[1])<tol)
  }))
}

#' @method is.isotropic variogramModel
#' @export
is.isotropic.variogramModel = function(x, tol=1e-10, ...){
  anis = x[,grep("anis", colnames(x))]
  all(apply(anis, 2, function(y) all(abs(y-y[1])<tol) ) )
}

#' @method is.isotropic variogramModelList
#' @export
is.isotropic.variogramModelList = function(x, tol=1e-10, ...) is.isotropic(x[[1]], tol=tol)

#' @method is.isotropic LMCAnisCompo
#' @export
is.isotropic.LMCAnisCompo = function(x, tol=1e-10, ...){
  all(sapply(x["A",], 1, function(y){
    ev = eigen(y$A, only.values=TRUE)[[1]]
    all(abs(ev-ev[1])<tol)
  }))
}





#' Check for any anisotropy class
#'
#' Check that an object contains a valid specification of anisotropy, in any form
#'
#' @param x object to check
#'
#' @return a logical, TRUE if the object is an anisotropy specification; FALSE otherwise
#' @export
#' @family anisotropy
#'
#' @examples
#' a =  anis2D_par2A(0.5, 30)
#' a
#' is.anisotropySpecification(a)
is.anisotropySpecification = function(x){
  length(grep("Anisotropy", class(x))) |  inherits(x, "Anisotropy")
}







#### empirical variogram ---------------------




#' Variogram method for gmSpatialModel objects
#'
#' Compute the empirical variogram of the conditioning data contained in a [gmSpatialModel-class] object
#'
#' @param object a gmSpatialModel object containing spatial data.
#' @param methodPars (currently ignored)
#' @param ... further parameters to [gstat::variogram()]
#'
#' @return Currently the function is just a convenience wrapper on
#' the variogram calculation functionalities of package "gstat",
#' and returns objects of class "\code{gstatVariogram}". Check the
#' help of \code{gstat::variogram} for further information.
#' In the near future, methods will be created, which will depend on
#' the properties of the two arguments provided,  \code{object} and
#' \code{methodPars}.
#' @export
#' @importFrom gstat variogram
variogram_gmSpatialModel <-  function(object, methodPars=NULL, ...){
  if(!is.null(methodPars)) stop("use 'variogram' with named parameters only")
  gstat::variogram(as.gstat(object), ...)
}


# Variogram calculations
#
# Compute empirical variograms out of a spatial data object
#
# @param object spatial data container
# @param ... further parameters for variogram calculation
#
# @return depending on the input data, different kinds of empirical variograms
# will be produced. See appropriate method descriptions.
#
# @importFrom sp variogram
# @export
##variogram <- function(object, ...) UseMethod("variogram", object)

# @describeIn variogram
# @method variogram default
# @export
#variogram.default <- function(object, ...){
#  return(variogram_gmSpatialModel(object, ...))
#}




#' Empirical variogram or covariance function in 2D
#'
#' compute the empirical variogram or covariance function in a 2D case study
#'
#' @param X matrix of Nx2 columns with the geographic coordinates
#' @param Z matrix or data.frame of data with dimension (N,Dv)
#' @param Ff for variogram, matrix of basis functions with nrow(Ff)=N (can be a N-vector of 1s);
#' for covariance function, a (N,Dv)-matrix or a  Dv-vector giving the mean values
#' @param maxdist maximum lag distance to consider
#' @param lagNr number of lags to consider
#' @param lags if maxdist and lagNr are not specified, either: (a) a matrix of 2 columns giving
#' minimal and maximal lag distance defining the lag classes to consider, or (b) a vector of lag breaks
#' @param azimuthNr number of azimuths to consider
#' @param azimuths if azimuthNr is not specified, either: (a) a matrix of 2 columns giving
#' minimal and maximal azimuth defining the azimuth classes to consider, or (b) a vector of azimuth breaks
#' @param maxbreadth maximal breadth (in lag units) orthogonal to the lag direction
#' @param minpairs minimal number of pairs falling in each class to consider the calculation sufficient; defaults to 10
#' @param cov boolean, is covariance (TRUE) or variogram (FALSE) desired? defaults to variogram
#'
#' @return An empirical variogram for the provided data. NOTE: avoid using directly gsi.* functions! They
#' represent either internal functions, or preliminary, not fully-tested functions. Use \code{\link{variogram}} instead.
#' @export
#' @family gmEVario functions
#'
#' @examples
#' library(gstat)
#' data("jura", package = "gstat")
#' X = as.matrix(jura.pred[,1:2])
#' Z = as.matrix(jura.pred[,c("Zn","Cd","Pb")])
#' vge = gsi.EVario2D(X,Z)
#' dim(vge)
#' dimnames(vge)
#' class(vge["gamma",1])
#' dim(vge["gamma",1][[1]])
#' vge["npairs",1]
#' vge["lags",1]
gsi.EVario2D = function(X,Z,Ff=rep(1, nrow(X)),
                      maxdist= max(dist(X[sample(nrow(X),min(nrow(X),1000)),]))/2,
                      lagNr = 15, lags = seq(from=0, to=maxdist, length.out=lagNr+1),
                      azimuthNr=4, azimuths = seq(from=0, to=180, length.out=azimuthNr+1)[1:azimuthNr],
                      maxbreadth=Inf, minpairs=10, cov=FALSE){

  # dimensions
  N = nrow(X)
  Dv = ncol(Z)
  Dg = ncol(X)
  stopifnot(N==nrow(Z))
  if(length(dim(Ff))==0){
    stopifnot(N==length(Ff))
  }else{
    stopifnot(N==nrow(Ff))
  }
  # expand the information given into a set of columns stating conditions
  if(length(dim(lags))==0){
    lags = data.frame(minlag=lags[-length(lags)], maxlag=lags[-1])
    if(maxbreadth!=Inf) lags[,"maxbreadth"]=maxbreadth
  }else if(dim(lags)==2){
    lags = data.frame(lags)
    colnames(lags) = c("minlag","maxlag","maxbreadth")[1:ncol(lags)]
  }else stop("lags can be either a vector of lags or a data.frame, see ?gsi.EVario2D")

  if(length(dim(azimuths))==0){
    tol = (azimuths[2]-azimuths[1])/2
    if(is.na(tol)) tol=180
    azimuths = data.frame(minaz=azimuths-tol, maxaz=azimuths+tol)
  }else if(dim(azimuths)==2){
    azimuths = data.frame(azimuths)
    colnames(azimuths) = c("minaz","maxaz")
  }else stop("azimuths can be either a vector of lags or a data.frame, see ?gsi.EVario2D")

  # compute pairs of locations
  ij = expand.grid(1:nrow(X), 1:nrow(X))# indices
  XX = X[ij[,1],]-X[ij[,2],] # locations
  XXabs = gmApply(XX, 1, function(x) sqrt(sum(x^2)))
  XXaz = gmApply(XX, 1, function(x) pi/2-atan2(x[2],x[1])) +2*pi
  XXaz = XXaz %% pi  # residual to 180??

  # compute residuals and pairs of variables appropriate to the structural function
  if(cov){
    if(all(dim(Z)==dim(Ff))){
      Z = as.matrix(Z-Ff)
    }else if(Dv==length(c(unlist(Ff)))){
      Z = as.matrix(sweep(Z, 2, Ff, "-"))
    }
    ZZ = outer(Z,Z) # variables
    ZZ = aperm(ZZ, c(1,3,2,4))
    dim(ZZ) = c(N*N, Dv, Dv)
  }else{
    Z = lm(as.matrix(Z)~as.matrix(Ff)+0)$residuals ## ideally this should be a GLS fit
    ZZ = Z[ij[,1],]-Z[ij[,2],]
  }


  # output
  ## ATTENTION: needs to be changed to return a structure (3,Na)-matrix of objects,
  #     like logratioVariogramAnisotropy
  Nh = nrow(lags)
  Na = nrow(azimuths)
  vg = array(0, dim=c(Nh, Dv, Dv, Na))
  n = array(0, dim=c(Nh, Na))
  azs = azimuths * pi/180
  res = sapply(1:Na, function(i){
    tk_a = (azs[i,1]<=XXaz) & (azs[i,2]>=XXaz)
    zz = ZZ[tk_a,]
    xxabs = XXabs[tk_a]
    xxaz = XXaz[tk_a]
    tk_h = outer(xxabs, lags[,1],">=") & outer(xxabs, lags[,2],"<=")
    if(ncol(lags)>2){
      tk_b = outer(xxabs * abs(sin((xxaz-(azs[i,2]-azs[i,1])))), lags[,3], "<=")
      tk_h = tk_h & tk_b
    }
    n[,i] = colSums(tk_h)
    for(j in 1:Nh){
      if(!is.na(n[j,i]) && n[j,i]>minpairs){
        if(cov){ # covariance function
          vg[j,,,i] = gmApply((ZZ[tk_a,][tk_h[,j],,]), c(2,3),"sum")/(n[j,i])
        }else{   # semi-variogram
          aux = rowSums(
            gmApply(X=ZZ[tk_a,][tk_h[,j],], 1, function(x)outer(x,x,"*"))
          )/(2*n[j,i])
          vg[j,,,i] = aux
        }
      }else{
        vg[j,,,i]=NA
      }
    }
    return(list(gamma=vg[,,,i], lags=gsi.lagClass(lags), npairs =n[,i]))
  })

  # output
  attr(res, "directions") = gsi.azimuthInterval(azimuths)
  # attr(res, "lags") = gsi.lagClass(lags)
  attr(res, "type") = ifelse(cov, "covariance","semivariogram")
  class(res) = "gmEVario"
  return(res)
}

#' Empirical variogram or covariance function in 3D
#'
#' compute the empirical variogram or covariance function in a 3D case study
#'
#' @param X matrix of Nx3 columns with the geographic coordinates
#' @param Z matrix or data.frame of data with dimension (N,Dv)
#' @param Ff for variogram, matrix of basis functions with nrow(Ff)=N
#' (can be a N-vector of 1s; should include the vector of 1s!!);
#' for covariance function, a (N,Dv)-matrix or a Dv-vector giving the mean values
#' @param maxdist maximum lag distance to consider
#' @param lagNr number of lags to consider
#' @param lags if maxdist and lagNr are not specified, either: (a) a matrix of 2 columns giving
#' minimal and maximal lag distance defining the lag classes to consider, or (b) a vector of lag breaks
#' @param dirvecs matrix which rows are the director vectors along which variograms will be computed (these will be normalized!)
#' @param angtol scalar, angular tolerance applied (in degrees; defaults to 90??, ie. isotropic)
#' @param maxbreadth maximal breadth (in lag units) orthogonal to the lag direction (defaults to `Inf`, ie. not used)
#' @param minpairs minimal number of pairs falling in each class to consider the calculation sufficient; defaults to 10
#' @param cov boolean, is covariance (TRUE) or variogram (FALSE) desired? defaults to variogram
#'
#' @return An empirical variogram for the provided data. NOTE: avoid using directly gsi.* functions! They
#' represent either internal functions, or preliminary, not fully-tested functions. Use \code{\link{variogram}} instead.
#'
#'
#' @export
#' @family gmEVario functions
#'
# @examples
# dt <- readr::read_csv("~/Geochem_sum_imp2.csv")
#
# X = as.matrix(dt[,c("X","Y","Z")])
# Z = as.matrix(compositions::clr(dt[,c("Cu","Zn","Pb", "As", "Cd", "In", "Other")]))
# dirvecs = rbind( c(1,0,0), c(1,1,0), c(0,1,0), c(-1,1,0), c(0,0,1))
# vge = gsi.EVario3D(X, Z, dirvecs=dirvecs, maxdist=50, angtol=10, maxbreadth=5,
# cov = TRUE, Ff = colMeans(Z))
# dim(vge)
# dimnames(vge)
# class(vge["gamma",1])
# dim(vge["gamma",1][[1]])
# vge["npairs",1]
# vge["lags",1]
# plot.gmEVario(vge, varnames = colnames(Z), commonAxis=FALSE,
#    vdir.up = 1:4, vdir.lo=5, xlim.up=c(0,50), xlim.lo=c(0,15)
#    )
## ie. the 4 planar directions on the upper triangle,
## the downhole direction in the lower triangle
gsi.EVario3D = function(X,Z,Ff=rep(1, nrow(X)),
                        maxdist= max(dist(X[sample(nrow(X),min(nrow(X),1000)),]))/2,
                        lagNr = 15, lags = seq(from=0, to=maxdist, length.out=lagNr+1),
                        dirvecs=t(c(1,0,0)), angtol=90,
                        maxbreadth=Inf, minpairs=10, cov=FALSE){
  # dimensions
  N = nrow(X)
  Dv = ncol(Z)
  Dg = ncol(X)
  stopifnot(N==nrow(Z))
  if(length(dim(Ff))==0){
    stopifnot(length(Ff) %in% c(N, Dv))
  }else{
    stopifnot(nrow(Ff) %in% c(N, 1) )
  }
  ## expand the information given into a set of columns stating conditions
  # prepare lags
  if(length(dim(lags))==0){
    lags = data.frame(minlag=lags[-length(lags)], maxlag=lags[-1])
    # if(maxbreadth!=Inf)
    lags[,"maxbreadth"] = maxbreadth
  }else if(dim(lags)==2){
    lags = data.frame(lags)
    colnames(lags) = c("minlag","maxlag","maxbreadth")[1:ncol(lags)]
  }else stop("lags can be either a vector of lags or a data.frame, see ?gsi.EVario3D")
  lagsSq = lags^2

  # prepare directions with tolerance
  if(length(dim(dirvecs))==0){
    if(length(dirvecs)!=3) stop("dirvecs must be a matrix with 3 columns!")
    dirvecs = t(dirvecs)
  }
  dirvecs = dirvecs / sqrt(rowSums(dirvecs^2))
  if(length(angtol) %in% c(1, nrow(dirvecs))){
    dirvecs = cbind(dirvecs, tol=cos(angtol*pi/180) )
  }else stop("angtol can be either a scalar or a vector of length=nrow(dirvecs), see ?gsi.EVario3D")

  # compute pairs of locations
  ij = expand.grid(1:nrow(X), 1:nrow(X))# indices
  XX = X[ij[,1],]-X[ij[,2],] # locations
  XXmodSq = gmApply(XX, 1, function(x) sum(x^2))
  XXmod = sqrt(XXmodSq)

  # compute residuals and pairs of variables appropriate to the structural function
  if(cov){
    if(all(dim(Z)==dim(Ff))){
      Z = as.matrix(Z-Ff)
    }else if(Dv==length(c(unlist(Ff)))){
      Z = as.matrix(sweep(Z, 2, Ff, "-"))
    }
    ZZ = outer(Z,Z) # variables
    ZZ = aperm(ZZ, c(1,3,2,4))
    dim(ZZ) = c(N*N, Dv, Dv)
  }else{
    Z = lm(as.matrix(Z)~as.matrix(Ff)+0)$residuals ## ideally this should be a GLS fit
    ZZ = Z[ij[,1],]-Z[ij[,2],]
  }

  # output
  ## ATTENTION: needs to be changed to return a structure (3,Na)-matrix of objects,
  #     like logratioVariogramAnisotropy
  Nh = nrow(lags)
  Na = nrow(dirvecs)
  vg = array(0, dim=c(Nh, Dv, Dv, Na),
             dimnames=list(rownames(lags), colnames(Z), colnames(Z), rownames(dirvecs))
  )
  n = array(0, dim=c(Nh, Na))
  res = sapply(1:Na, function(i){
    projections = abs(XX %*% dirvecs[i,1:3])
    cosinus = ifelse(XXmod==0,0,projections/XXmod)
    tk_a = cosinus > dirvecs[i,"tol"]
    xxabs = XXmod[tk_a]
    tk_h = outer(xxabs, lags[,1],">=") & outer(xxabs, lags[,2],"<=")
    if(ncol(lags)>2){
      residualsSq = XXmodSq[tk_a]-projections[tk_a]^2
      tk_b = outer(residualsSq, lagsSq[,3], "<=")
      tk_h = tk_h & tk_b
    }
    n[,i] = colSums(tk_h)
    for(j in 1:Nh){
      if(!is.na(n[j,i]) && n[j,i]>minpairs){
        if(cov){ # covariance function
          vg[j,,,i] = gmApply((ZZ[tk_a,,][tk_h[,j],,]), c(2,3),"sum")/(n[j,i])
        }else{   # semi-variogram
          aux = rowSums(
            gmApply(X=ZZ[tk_a,][tk_h[,j],], 1, function(x)outer(x,x,"*"))
          )/(2*n[j,i])
          vg[j,,,i] = aux
        }
      }else{
        vg[j,,,i]=NA
      }
    }
    return(list(gamma=vg[,,,i], lags=gsi.lagClass(lags), npairs =n[,i]))
  })

  # output
  attr(res, "directions") = gsi.directorVector(dirvecs[,1:3])
  # attr(res, "lags") = gsi.lagClass(lags)
  attr(res, "type") = ifelse(cov, "covariance","semivariogram")
  class(res) = "gmEVario"
  return(res)
}





#' Plot empirical variograms
#'
#' Flexible plot of an empirical variogram of class gmEVario
#'
#' @param x object to print, of class gmEVario
#' @param xlim.up range of X values to be used for the diagrams of the upper triangle
#' @param xlim.lo range of X values to be used for the diagrams of the lower triangle
#' @param vdir.up in case of anisotropic variograms, indices of the directions to be plotted
#' on the upper triangle
#' @param vdir.lo ..., indices of the directions to be plotted on the lower triangle
#' @param varnames variable names to be used
#' @param type  string, controlling whether to plot lines, points, etc (see \code{\link{plot}})
#' @param add boolean, add stuff to an existing diagram?
#' @param commonAxis boolean, should vertical axes be shared by all plots in a row?
#' @param cov boolean, is this a covariance? (if FALSE, it is a variogram)
#' @param closeplot logical, should the plot be left open (FALSE) for further changes, or be frozen (TRUE)?
#' defaults to TRUE
#' @param ... further parameters to \code{\link{matplot}}
#'
#' @return invisibly, the graphical parameters active before calling the function.
#' This is useful for freezing the plot if you provided `closeplot=FALSE`.
#'
#' How to use arguments `vdir.lo` and `vdir.up`? Each empirical variogram \code{x} has been
#' computed along certain distances, recorded in its attributes and retrievable with command
#' \code{\link{ndirections}}.
#' @export
#' @family gmEVario functions
#' @method plot gmEVario
#'
#' @examples
#' library(gstat)
#' data("jura", package = "gstat")
#' X = as.matrix(jura.pred[,1:2])
#' Z = as.matrix(jura.pred[,c("Zn","Cd","Pb")])
#' vge = gsi.EVario2D(X,Z)
#' plot(vge)
#' plot(vge, pch=22, lty=1, bg="grey")
plot.gmEVario = function(x, xlim.up=NULL, xlim.lo=NULL, vdir.up= NULL, vdir.lo= NULL,
                         varnames = dimnames(x$gamma)[[2]], type="o",
                         add=FALSE, commonAxis=TRUE, cov =attr(x,"type")=="covariance",
                         closeplot=TRUE, ...){
  Dv = dim(x[1,1][[1]])[2]
  if(is.null(varnames)) varnames = paste("v", 1:Dv, sep="")
  if(is.null(vdir.up)&is.null(vdir.lo)) vdir.lo <- 1:ndirections(x)
  if( any(c(vdir.up, vdir.lo)>ndirections(x))){
    stop("indicated directions (vdir.up or vdir.lo) do not exist in x")
  }
  if(is.null(xlim.up)){
    if(add){
      par(mfg=c(1,2))
      xlim.up = par()$usr[1:2]
    }else if(!is.null(vdir.up)){
      maxdist = max(sapply(x["lags",], gsi.midValues.lagClass ) )
      xlim.up = c(0, maxdist)
    }
  }
  if(is.null(xlim.lo)){
    if(add){
      par(mfg=c(2,1))
      xlim.lo = par()$usr[1:2]
    }else if(!is.null(vdir.lo)){
      maxdist = max(sapply(x["lags",], gsi.midValues.lagClass ) )
      xlim.lo = c(0, maxdist)
    }
  }

  opar = par()
  opar = par_remove_readonly(opar)

  if(closeplot) on.exit(par(opar))

  myplot = function(...) matplot(type=type, ylab="", xlab="",xaxt="n", ...)
  if(add) myplot = function(...) matpoints(type=type, ...)
  if(!add){
    myplot(c(0,0), c(0,0), pch="", ann=FALSE, bty="n", yaxt="n")
    par(mfrow=c(Dv+1,Dv+1), mar=c(2,3,0,0), oma=c(1,4,1,1), xpd=NA)
  }
  for(i in 1:Dv){
    for(j in 1:Dv){
      if((i>=j)&(!is.null(vdir.lo))){
        par(mfg=c(i+1,j,Dv+1,Dv+1))
        ylim = range(sapply(vdir.lo, function(kk) x["gamma",kk][[1]][,i,j]), na.rm = TRUE)
        if(commonAxis) ylim=range(sapply(vdir.lo, function(kk) x["gamma",kk][[1]][,,j]), na.rm=TRUE)
        myplot(
          sapply(vdir.lo, function(kk) gsi.midValues.lagClass(x["lags",kk][[1]])),
          sapply(vdir.lo, function(kk) x["gamma",kk][[1]][,i,j]),  ylim=ylim, ...)
        if(i==j){
          axis(side = 3)
          mtext(text=varnames[i], side = 3, line=3)
          mtext(text=varnames[i], side = 4, line=3)
        }
      }
      if((i<=j)&(!is.null(vdir.up))){
        par(mfg=c(i,j+1,Dv+1,Dv+1))
        ylim = range(sapply(vdir.up, function(kk) x["gamma",kk][[1]][,i,j]), na.rm = TRUE)
        if(commonAxis) ylim=range(sapply(vdir.up, function(kk) x["gamma",kk][[1]][,,j]), na.rm=TRUE)
        myplot(
          sapply(vdir.up, function(kk) gsi.midValues.lagClass(x["lags",kk][[1]])),
          sapply(vdir.up, function(kk) x["gamma",kk][[1]][,i,j]),  ylim=ylim, ...)
        if(i==j){
          axis(side = 1)
          mtext(text=varnames[i], side = 1, line=3)
          mtext(text=varnames[i], side = 2, line=3)
        }
      }
    }}
  mtext(text="lag distance", side=1, outer = TRUE, line=0)
  mtext(text=c("semivariogram","covariance")[cov+1], side=2, outer = TRUE, line=2)
  attr(opar, "vdir.up")=vdir.up
  attr(opar, "vdir.lo")=vdir.lo
  invisible(opar)
}






#' Convert empirical structural function to gmEVario format
#'
#' Convert empirical covariance functions or variograms to the format gmEVario
#' of package gmGeostats
#'
#' @param vgemp variogram/covariance function to be converted
#' @param ... further parameters
#'
#' @return the empirical covariance function or variogram, recasted to class
#' \code{gmEVario}. This is a generic function. Methods exist for objects of