R/Bayes_spaces.R

In DistributionalDDA: Distributional Data Analysis in Bayes Spaces

#' Computes expectations from reference measures
#'
#' @param measure the (probability) measure to compute the expectation
#'     from
#' @param f a real valued random variable, represented as a vectorized
#'     function, for which the expectation shall be computed.
#' @param ... further arguments to f or the method used
#' @param method a method to be selected. Typical methods are
#'     "sample", "MonteCarlo" and "integrate", but only sample is
#'     availabel for all reference measures. "sample" uses the
#'     prestored sample, "MonteCarlo" uses the simulated function of
#'     the reference measure, and "integrate" uses numerical
#'     integration. Typically the prestored integration sample
#'     performs well with derivatives measures derived
#'     ReferenceMeasureUniform, ReferenceMeasusureTransformed, where
#'     it proves a dense grid in the area of measure, and
#'     ReferenceMeasure discrete where it uses precise computation
#'     using a weighted sample. For Cartesian products involving more
#'     than one continues measure due to the curse of dimensionality
#'     typically only MonteCarlo is expected to work reasonably well.
#' @param tol the absolute or relative tolerance (whichever is larger).
#'     accepted for methods with tolerance control.
#'
#' @returns the numerically computed expected value of f with respect
#'     to the measure. The method
#'
#' @examples
#' m <- ReferenceMeasureUniform(0,10)
#' expectation(m, function(x) x^2)
#'
#' @family ReferenceMeasures
#' @export
#'
expectation <- function(measure,f,...,method) UseMethod("expectation",measure)

#' @describeIn expectation general method supporting only the "sample"
#'     method.
#' @export
#' @method expectation ReferenceMeasure
expectation.ReferenceMeasure <-  function(measure,f,...,method=c("default","sample","MonteCarlo","integration"),tol=1E-2) {
    method <- selectIntegrationMethod(measure,match.arg(method))
    if( method=="MonteCarlo" ) {
        .expectationUsingSim(measure=measure,
                             f=as.function(f),...,tol=tol)
    } else {
        if( method!="sample" )
            warning("Unsupported integration method ",method)
        .expectationUsingSample(getSample(measure),as.function(f),...)
    }
}

#' Identifies the integration method to be choosen for a measure
#' @param measure the ReferenceMeasure for which the method should be selected
#' @param method a user provide method hint
#'
#' @return a string identifiying the method
#' @export
#' @examples
#' selectIntegrationMethod(StdReferenceMeasureUnif01,"default")
#'
selectIntegrationMethod <- function(measure,method=c("default","sample","MonteCarlo","integrate")) {
    method = match.arg(method)
    if( method=="default" )
        UseMethod("selectIntegrationMethod",measure)
    else
        method
}

#' @describeIn selectIntegrationMethod standard method using "sample"
#' @export
selectIntegrationMethod.ReferenceMeasure = function(measure,method) "sample"


#' Checks whether values are in the domain of the ReferenceMeasure
#'
#' The function acts vectorized, i.e. in case of numbers or other
#' atomic objects it expects a vector of the given type. In case of
#' values that can not be stored in a vector it expects a list of the elements.
#'
#' @param x values to be checked to be in the reference measure
#' @param referenceMeasure the referenceMeasure to be checked against. It can
#'     be provided indirectly.
#'
#' @returns a logical vector of the same length as `x`, `TRUE` where
#'     the corresponding value lies in the domain (support) of the
#'     reference measure and `FALSE` otherwise.
#'
#' @family ReferenceMeasures
#' @export
#' @seealso isInDomain.ReferenceMeasureUniform
#' @examples
#' isInDomain(1:10,B2Unif01)
#'
isInDomain=function(x,referenceMeasure) UseMethod("isInDomain",referenceMeasure)
#' @describeIn isInDomain method for BayesSpaces
#' @export
isInDomain.BayesSpace=function(x,referenceMeasure) isInDomain(x,referenceMeasure=getReferenceMeasure(referenceMeasure))
#' @describeIn isInDomain method for L2Spaces
#' @export
isInDomain.L2Space=function(x,referenceMeasure) isInDomain(x,referenceMeasure=getReferenceMeasure(referenceMeasure))
#' @describeIn isInDomain method for BayesSpaces
#' @export
isInDomain.BayesSpaceBasis=function(x,referenceMeasure) isInDomain(x,referenceMeasure=getReferenceMeasure(referenceMeasure))


#' create the ReferenceMeasureUniform class
#'
#' @param a the minimum of the domain interval
#' @param b the maximum of the domain interval
#' @param sample the sample to be used for internal computations (use with care)
#' @param seed optional integer. If given, the Monte Carlo sample is
#'     drawn under this seed and the state of the random number
#'     generator is restored afterwards (the same convention as
#'     \code{simulate()}), so that the constructed object is
#'     reproducible independently of the global RNG state and of any
#'     code executed before or after the call. The default sample
#'     however is currently deterministic.
#'
#' @returns a 'ReferenceMeasureUniform class'
#' @family ReferenceMeasures
#' @export
#'
#' @family BDataBayesSpaceSplines
#' @seealso [BData]
#'
#' @examples
#' min <- -5
#' max <- 5
#' RMU <- ReferenceMeasureUniform(min,max)
#' RMU
#' ## reproducible construction, RNG state left untouched:
#' RMU1 <- ReferenceMeasureUniform(0, 1, seed = 20260702)
#' RMU2 <- ReferenceMeasureUniform(0, 1, seed = 20260702)
#' identical(RMU1$sample, RMU2$sample)
#' @importFrom stats runif
ReferenceMeasureUniform = function(a=0,b=1,sample=makeUnifGrid(2^12,a,b),seed=NULL){
    .call = match.call()
    stopifnot(b>a)
    if (!is.null(seed)) {
        if (!exists(".Random.seed", envir = globalenv(), inherits = FALSE))
            stats::runif(1)
        oldSeed <- get(".Random.seed", envir = globalenv())
        # on.exit(assign(".Random.seed", oldSeed, envir = globalenv()))
        set.seed(seed)
    }
    force(sample)
    structure(list(a=a,b=b,sample=sample),
              class=c("ReferenceMeasureUniform","ReferenceMeasure"),
              call=.call)
}

# ' !!! the name for simulations in R is simulate (not random) Please change to simulate including the simulated interface!!!
# '
# ' simulate from the distribution given by the ReferenceMeasureUniform
# '
# '
# '
# ' @param x the ReferenceMeasureUniform object
# ' @param n sample size
# ' @param drop ignored (TODO: remove?)
# ' @param ... ignored
# '
# '
# '
# '
# ' @returns a vector of length n
# ' @export
# '
# ' @examples
# ' RMU <- ReferenceMeasureUniform(0,2)
# ' simulated <- random.ReferenceMeasureUniform(RMU,n=50)
# ' simulated
# ' @importFrom stats runif
# random.ReferenceMeasureUniform=function(x, n = 1L, drop = TRUE, ...){
#  runif(n,x$a,x$b)
# }



#' integrate a function over the uniform reference measure
#'
#' @param measure the ReferenceMeasureUniform object
#' @param f function to be integrated
#' @param ... further arguments to [stats::integrate()]
#' @param method a method for computing the expectation given as a
#'     string. The default argument shows possible choices. Default
#'     chooses the default for reference measure (typically "sample"),
#'     "sample" uses the reference sample of the measure, "MonteCarlo"
#'     uses iterative random simulation until a certain tolerance is
#'     reached, "integrate" is only possible in a few measures and
#'     tries a numerical integration.
#' @returns the integral of f over c(x[1],x[2])
#' @export
#' @examples
#' RMU <- ReferenceMeasureUniform(a=-5,b=5)
#' expectationRMU <- expectation(measure=RMU,f=cos,method="integrate")
#' expectationRMU
#' @importFrom stats integrate
expectation.ReferenceMeasureUniform = function(measure,f,...,method=c("sample","integrate")) {
    method = match.arg(method)
    switch( method,
           integrate=integrate(f,measure$a,measure$b,...),
           sample=NextMethod()
           )
}


#' assess whether a vector of points belong to the domain of the uniform reference measure
#'
#' unlike typical measures the domain is closed, i.e. the extreme
#' values are considered part of the support.
#'
#' @param x a dataset of points to be checked
#' @param referenceMeasure the ReferenceMeasureUniform object,
#'     providing the domain
#'
#' @returns a logical vector of TRUE/FALSE with a length length(x)
#' @export
#'
#' @examples
#' isInDomain(c(1:10),ReferenceMeasureUniform(0,1))
#'
isInDomain.ReferenceMeasureUniform= function(x,referenceMeasure){
  (x>=referenceMeasure$a)&(referenceMeasure$b>=x)
}

#' Printing ReferenceMeasures
#'
#' The measure is printed by the call that created it.
#'
#' @param x the reference measure to be printed
#' @param ... further generic arguments, which will be ignored.
#'
#' @returns an invisible version of the reference measure
#' @export
#'
#' @examples
#' print(StdReferenceMeasureUnif01)
#' print(StdReferenceMeasureN01)
print.ReferenceMeasure = function(x,...){
    .call=attr(x,"call")
    if( is.null(.call) )
        NextMethod()
    else print(.call)
    invisible(x)
}

#' create the `BayesSpaceBasis` object for `ReferenceMeasureUniform`
#'
#' @param basis object of class `ReferenceMeasureUniform`
#' @param bayesSpace object of class `BayesSpace`, of which the result
#'     is a basis.
#' @param nbasis number of basis functions for the B-spline basis
#' @param order order of the polynomial of the basis functions
#' @param breaks ascending nonrepetetive vector of knots
#' @param ...  arguments to be passed to [fda::create.bspline.basis()]
#' @param orthonormalize boolean if true creates and orthonormal basis
#' @returns object of class `BayesSpaceBasis`
#' @export
#'
#' @importFrom utils head
#' @importFrom utils tail
#' @examples
#' basis <- BayesSpaceBasis(ReferenceMeasureUniform(0,1))
#' basis
#' #basis(seq(0,1,length.out=10))
BayesSpaceBasis.ReferenceMeasureUniform = function(basis,
                                                   bayesSpace=getBayesSpace(basis),
                                                   orthonormalize=FALSE,
                                                   nbasis=7,
                                                   order=4,
                                                   breaks=seq(basis$a,
                                                              basis$b,
                                                              length.out=(nbasis-order+3)),...){
    # Remark: to produce a BayesSpaceBasis with nbasis elements
    # we require a bspline basis with one element more.
    structure(
    BayesSpaceBasis(
        fda::create.bspline.basis(rangeval = c(basis$a,basis$b),
                                  nbasis=nbasis+1,
                                  norder=order,
                                  breaks=breaks,
                                  ...),
        bayesSpace=bayesSpace,
        orthonormalize=orthonormalize),
    call=match.call())
}


#' @export
#' @importFrom stats runif
getSample.ReferenceMeasureUniform = function(x,nsim=10000,...) {
  if(is.null(x$sample)) return(runif(nsim,x$a,x$b))
  x$sample
}

#' Evaluates a BayesSpaceBasis
#'
#' The basis is represented by log-densities. The function is
#' vectorized, i.e. multiple x values may be provided.
#'
#' @param basis a description of a basis, e.g. through a
#'     BayesSpaceBasis object
#' @param x the data values on which the basis shall be evaluated
#' @param ... further parameters that a generic can use,
#'     e.g. describing method and precision of evaluation.
#'
#' @returns a matrix with the evaluation for one each of the x in a row.
#'
#' @export
#'
#' @examples
#' x = seq(0,1,by=0.1)
#' basisEval(B2Unif01Basis5,x)
basisEval = function(basis,x,...) UseMethod("basisEval",basis)

#' @describeIn basisEval the standard method for BayesSpaceBasis
#' @export
basisEval.BayesSpaceBasis= function(basis,x,...) .operateAx2v(basis$basis(x),basis$A)

#' @describeIn basisEval a method implementation for basisfd
#' @export
#' @importFrom fda eval.basis
basisEval.basisfd = function(basis,x,...) fda::eval.basis(x,basis,...)

#' @describeIn basisEval the basis can also be converted to a function
#'     or called directly
#' @export
as.function.BayesSpaceBasis = function(x,...){
  f=x
  function(x,...) basisEval(f,x,...)
}

#' Metric tensor of a space or basis
#'
#' The metric tensor is the Gram matrix `S`, used in `|v|^2 = v' S v` to compute the
#' metric of the space when `v` is the vector of basis coefficients. `S` contains the
#' scalar products of the basis elements, `S[i,j] = <psi_i, psi_j>`.
#'
#' @details
#' For an orthonormal basis `S` is the identity matrix. The default zero-integral
#' B-spline (ZB-spline) basis is **not** orthonormal, so `S` is a banded, non-identity
#' matrix; inner products ([scalarProduct]), norms ([norm.BData]), metric variance
#' ([mvar][mvar.BData]) and PCA ([princomp.BData]) all rely on it. You can check whether
#' a basis is orthonormal by comparing the result to `diag(nrow(S))`.
#'
#' @param x the space or basis for which the metric tensor shall be returned.
#' @param ... further generic parameters for the computation
#'
#' @returns A symmetric positive-definite matrix `S` (the Gram matrix / metric tensor),
#'     with `S[i,j]` the inner product of basis elements `i` and `j`.
#'
#' @seealso [BayesSpaceBasis], [scalarProduct], [princomp.BData]
#' @export
#'
#' @examples
#' getMetricTensor(B2Unif01Basis5)   # not the identity: ZB-splines are not orthonormal
#'
getMetricTensor = function(x,...) UseMethod("getMetricTensor",x)

#' @export
getMetricTensor.BayesSpaceBasis = function(x,...) x$S

# ReferenceMeasureTransformed
#' A general continuous reference measure on the real line
#'
#' The reference measure is created by transforming a ReferenceMeasure
#' Uniform to a given distributions based using a distribution
#' function and a quantile function.
#'
#' @param cdf the cumulative distribution function of the reference
#'     measure. For the standard distributions listed as default
#'     argument, it can be provided as a string. In this case the
#'     corresponding quantile and density functions are selected
#'     automatically. The distribution parameters should be provided
#'     in ...
#' @param qf the corresponding quantile function, mapping a value of
#'     the reference measure to its quantiles.
#' @param reference the uniform ReferenceMeasure to be used to do the
#'     computations
#' @param ... further arguments for cdf, qf, and lebesgue_density
#' @param pdf a function representing the Lebesgue Density (aka as
#'     probability density function) of the measure, for conversion to
#'     and from pdfs.
#' @param .args the list of further arguments to cdf and qf
#' @param reference the underlying untransformed reference measure to
#'     be used
#' @param x a ReferenceMeasureTransformed object
#' @param .eps the precision to check the roundtrip between
#'
#' @returns a reference measure with representing the reference
#'     measure described by the distribution and quantile function.
#' @family ReferenceMeasures
#' @export
#' @importFrom stats dbeta
#' @importFrom stats dcauchy
#' @importFrom stats dchisq
#' @importFrom stats dexp
#' @importFrom stats df
#' @importFrom stats dgamma
#' @importFrom stats dlnorm
#' @importFrom stats dlogis
#' @importFrom stats dunif
#' @importFrom stats dweibull
#' @importFrom stats dwilcox
#' @importFrom stats pbeta
#' @importFrom stats pcauchy
#' @importFrom stats pchisq
#' @importFrom stats pexp
#' @importFrom stats pf
#' @importFrom stats pgamma
#' @importFrom stats plnorm
#' @importFrom stats plogis
#' @importFrom stats pnorm
#' @importFrom stats pt
#' @importFrom stats punif
#' @importFrom stats pweibull
#' @importFrom stats pwilcox
#' @importFrom stats qbeta
#' @importFrom stats qcauchy
#' @importFrom stats qchisq
#' @importFrom stats qexp
#' @importFrom stats qgamma
#' @importFrom stats qlnorm
#' @importFrom stats qlogis
#' @importFrom stats qnorm
#' @importFrom stats qt
#' @importFrom stats qunif
#' @importFrom stats qweibull
#' @importFrom stats qwilcox
#' @examples
#' rmt = ReferenceMeasureTransformed(pnorm,qnorm,mean=2,sd=2)
#' isInDomain(0:17,rmt)
#' expectation(rmt,\(x) x^2)
#' hist(getSample(rmt))
ReferenceMeasureTransformed = function(cdf=c("pnorm","pbeta","pcauchy","pchisq","pexp","pf","pgamma","plnorm","plogis","pt","punif","pweibull","pwilcox"),
                                       qf=cdf2qf(cdf,...),
                                       ...,
                                       pdf=cdf2density(cdf,...),
                                       reference=StdReferenceMeasureUnif01,
                                       .args=list(...),
                                       .eps=1E-7) {
    .call=match.call()
    s = getSample(reference)
    if( is.character(cdf) ) {
        useCdf = match.arg(cdf)
        cdf = c("pnorm"=stats::pnorm,"pbeta"=stats::pbeta,"pcauchy"=stats::pcauchy,
                "pchisq"=stats::pchisq,"pexp"=stats::pexp,
                "pf"=stats::pf,"pgamma"=stats::pgamma,"plnorm"=stats::plnorm,
                "plogis"=stats::plogis,"pt"=stats::pt,"punif"=stats::punif,
                "pweibull"=stats::pweibull,"pwilcox"=stats::pwilcox)[[useCdf]]
        qf  =c("pnorm"=stats::qnorm,"pbeta"=stats::qbeta,"pchauchy"=stats::qcauchy,
                "pchisq"=stats::qchisq,"pexp"=stats::qexp,
                "pf"=stats::qf,"pgamma"=stats::qgamma,"plnorm"=stats::qlnorm,
                "plogis"=stats::qlogis,"pt"=stats::qt,"punif"=stats::qunif,
                "pweibull"=stats::qweibull,"pwilcox"=stats::qwilcox)[[useCdf]]
        pdf  =c("pnorm"=stats::dnorm,"pbeta"=stats::dbeta,"pchauchy"=stats::dcauchy,
                "pchisq"=stats::dchisq,"pexp"=stats::dexp,
                "pf"=stats::df,"pgamma"=stats::dgamma,"plnorm"=stats::dlnorm,
                "plogis"=stats::dlogis,"pt"=stats::qt,"punif"=stats::dunif,
                "pweibull"=stats::dweibull,"pwilcox"=stats::dwilcox)[[useCdf]]
    }

    if( any(abs(evalFArgs(cdf,evalFArgs(qf,s,.args=.args),.args=.args)-s) > .eps) )
        warning("cdf and qf do not match in ReferenceMeasureTransformed")
    structure(list(cdf=setFArgs(cdf,.args=.args),qf=setFArgs(qf,.args=.args),
                   reference=reference,
                   sample=evalFArgs(qf,getSample(reference),.args=.args),
                   pdf=setFArgs(pdf,.args=.args)),
              class=c("ReferenceMeasureTransformed","ReferenceMeasure"),
              call=.call)
}

#' @describeIn isInDomain method for ReferenceMeasureTransformed. The
#'     domain is always an open interval, where 0<cdf(x)<1.
#' @export
#'
isInDomain.ReferenceMeasureTransformed = function(x,referenceMeasure) {
    p = referenceMeasure$cdf(x)
    p>0 & p<1
}

## ' @describeIn expectation method for ReferenceMeasureTransformed.
## ' @export
#expectation.ReferenceMeasureTransformed = function(measure,f,...){
# expectation(measure$reference,function(x,...) f(measure$qf(x),...),...)
#}


#' @describeIn ReferenceMeasureTransformed returns transformed sample from the underlying reference
#' @export
getSample.ReferenceMeasureTransformed = function(x,...) {
  x$sample
}


#' Generates a BayesSpaceBasis for transformed Bayes Spaces
#'
#' The basis is a transformation of the Basis that would be computed for the
#' underlying untransformed BayesSpace.
#'
#' @param basis a ReferenceMeasureTransformed, for which the basis
#'     shall be computed.
#' @param bayesSpace the BayesSpace of this ReferenceMeasure
#' @param ... further parameters for the BayesSpaceBasis method of the
#'     underlying untransformed ReferenceMeasure, typically
#'     ReferenceMeasureUniform. The breaks are than to be understood
#'     as quantiles.
#' @param orthonormalize boolean: if TRUE the basis is orthonoramlized, using a singluar using a singular value decomposition based matrix square root.
#'
#' @returns an object of class `"BayesSpaceBasis"`, transforming the
#'     basis of the underlying untransformed BayesSpace.
#' @export
#'
#' @examples
#' rmt <- ReferenceMeasureTransformed("pnorm", mean=2, sd=2)
#' basis <- BayesSpaceBasis(rmt)
#' basis
#'
BayesSpaceBasis.ReferenceMeasureTransformed = function(basis,bayesSpace=getBayesSpace(basis),orthonormalize=FALSE,...) {
    .call = match.call()
    innerBasis = BayesSpaceBasis(basis$reference,
                                 orthonormalize=orthonormalize,...)
    structure(BayesSpaceBasis(
        basis=function(x) innerBasis$basis(basis$cdf(x)),
        bayesSpace=bayesSpace,
        A=innerBasis$A,
        S=innerBasis$S,
        centered=isCentered(innerBasis),
        ),
        call=.call)
}


#' The support of the measures
#'
#' For this package the support tells what is needed to know
#' for plotting
#'
#' @param x the object whos support shall be returned
#' @param X the discrete values of the support, possibly an empty
#'     vector.
#' @param a lower limit of the continuous support interval (a,b), may
#'     be -Inf
#' @param b upper limit of the continuous support interval (a,b), may
#'     be Inf, a==b (which should be a value in x) means no continuous
#'     support)
#' @param lim the effective xlim to use for plots or continuous
#'     evaluations on the support.
#' @param ... further unused arguments for generic use
#' @param alpha a small real number given the portion of the reference
#'     support to be ignored to avoid infinite supports.
#' @param extrems a vector of two numerical values used to compute
#'     values near -Inf and Inf, to look whether the function is a
#'     density (approach zero at both sides) or a cdf (go from 0 to 1)
#' @param starts starting values to find a point where the function is not zero.
#'
#' @returns An object of class "Support" with the given parameters
#' @family Support
#' @export
#' @examples
#' Support(0,1)
#' Support(X=1:4)
#' getSupport(StdReferenceMeasureUnif01)
#' getSupport(StdReferenceMeasureN01)
#'
Support <- function(a,b,X=c(),lim) {
    if( is.matrix(X) ) {
        if( ncol(X)== 1 )
            X=c(X)
    }
    if( missing(a) ) {

        if( missing(X) ) { # take continues from lim
            if( !missing(lim) ) {
                a=lim[1]
                if( missing(b) )
                    b=lim[2]
            } else {
                stop("Support missing parameters: a,X,lim: please give at least one")
            }
        } else { # no contiues area
            if( is.numeric(X) && !is.matrix(X) ) {
                a=min(X)
                b=a
                if( missing(lim) )
                    lim=range(X)
            } else {
                ## Nothing on Real line
                a=0
                b=0
                if( missing(lim) ) {
                    if( is.numeric(X) ) {
                        lim=c(0,nrow(X)+1)
                    } else {
                        lim=c(0,length(X)+1)
                    }
                }
                .call = call("Support",X=X,lim=lim)
                return(structure(list(a=a,b=b,X=X,lim=lim),
                                 class=c("DiscreteSupport","Support"),
                                 call=.call
                       ))
            }
        }
    }
    if( is(a,"Support") )
        return(a)
    if( a==b ) {
        a=X[1]
        b=a
    }
    if( missing(lim) )
        lim=range(c(a,b,X))
    if( !all(is.finite(lim)) ) {
        stop("Infinite plotting limits in Support. Give effective limits.")
    }
    if( length(X)==0 ) {
        X<-NULL
        .call=call("Support",a=a,b=b,lim=lim)
    } else {
        .call=call("Support",a=a,b=b,X=X,lim=lim)
    }
    structure(list(a=a,b=b,X=X,lim=lim),
              class=c("RealSupport","Support"),
              call=.call)
}

#' @export
print.Support <- function(x,...) {
    .call = attr(x,"call")
    if( is.null(.call) )
        NextMethod()
    print(.call,...)
}

#' @describeIn Support method for getting plotting limits of supports
#' @export
getSupport1DPlottingLimits = function(x) UseMethod("getSupport1DPlottingLimits",getSupport(x))
#' @export
#' @method getSupport1DPlottingLimits RealSupport
getSupport1DPlottingLimits.RealSupport = function(x) getSupport(x)$lim
#' @export
#' @method getSupport1DPlottingLimits DiscreteSupport
getSupport1DPlottingLimits.DiscreteSupport = function(x) getSupport(x)$lim
#' @export
#' @method getSupport1DPlottingLimits Support
getSupport1DPlottingLimits.Support = function(x) stop("1D Plotting limits not defined for ",class(x)[1])

#' Gets coordinates on an axis for possibly non real data
#'
#' @param X the data
#' @param support a [Support] object giving the continuous, discrete or mixed support of the data.
#' @param d the dimension of the space, that may be used
#' @return the value(s) of `X` mapped onto the `support`, i.e. possibly
#' the same values of `X` (i.e. if the support is a subset of the real line)
#' or the indices corresponding to the atoms in the support equal to the
#' values of `X` (if the support is a discrete set of points or of categories).
#' @export
mapToRealLineWithSupport = function(X,support,d=1) UseMethod("mapToRealLineWithSupport",getSupport(support))
#' @export
#' @method mapToRealLineWithSupport RealSupport
mapToRealLineWithSupport.RealSupport = function(X,support,d=1) X
#' @export
#' @method mapToRealLineWithSupport DiscreteSupport
mapToRealLineWithSupport.DiscreteSupport = function(X,support,d=1) {
    atoms = getSupportAtoms(getSupport(support))
    if( is.matrix(atoms) ) {
        sapply(1:nrow(X),function(i) which(apply(c(X[i,])==t(atoms),2,all)))
    } else {
        match(X,atoms)
    }
}





#' @describeIn Support method for getting values with mass
#' @export
getSupportAtoms = function(x) UseMethod("getSupportAtoms",getSupport(x))
#' @export
#' @method getSupportAtoms RealSupport
getSupportAtoms.RealSupport = function(x) getSupport(x)$X


#' @describeIn Support get the support associated to an object
#' @export
getSupport <- function(x,...) UseMethod("getSupport",x)
#' @describeIn Support method for Support returns itself
#' @export
#' @method getSupport Support
getSupport.Support = function(x,...) x
#' @describeIn Support method for BayesSpaces
#' @export
#' @method getSupport BayesSpace
getSupport.BayesSpace = function(x,...) getSupport(getReferenceMeasure(x),...)
#' @describeIn Support method for BayesSpaces
#' @export
#' @method getSupport BayesSpaceBasis
getSupport.BayesSpaceBasis = function(x,...) getSupport(getReferenceMeasure(x),...)
#' @describeIn Support method for BData
#' @export
#' @method getSupport BData
getSupport.BData = function(x,...) getSupport(getReferenceMeasure(x),...)
#' @describeIn Support method for ReferenceMeasureUniform
#' @export
#' @method getSupport ReferenceMeasureUniform
getSupport.ReferenceMeasureUniform = function(x,...) Support(x$a,x$b)
#' @describeIn Support method for ReferenceMeasureDiscrete
#' @export
#' @method getSupport ReferenceMeasureDiscrete
getSupport.ReferenceMeasureDiscrete = function(x,...) Support(X=x$X)
#' @describeIn Support method for ReferenceMeasureTransformed
#' @export
#' @method getSupport ReferenceMeasureTransformed
getSupport.ReferenceMeasureTransformed = function(x,...,alpha=0.01) {
    S=getSupport(x$reference)
    a=x$qf(S$a)
    b=x$qf(S$b)
    if( length(x$X) == 0 )
        X=a[c()]
    else
        X=x$qf(x$X)
    if( all(is.finite(a) & is.finite(b)) )
        return(Support(a,b,X=X))
    else {
        aEff = S$a*(1-alpha/2)+S$b*alpha/2
        bEff = S$a*alpha/2+S$b*(1-alpha/2)
        lim = x$qf(c(aEff,bEff))
        if( !all(is.finite(lim) ))
            warning("Non finite values for lim in Support for ReferenceMeasureTransformed.")
        return(Support(a,b,X=X,lim=lim))
    }
}

#' @describeIn Support combines the supports of all Elements in the list
#' @export
#' @method getSupport list
getSupport.list = function(x,...) {
    supports = lapply(x,getSupport)
    Support(a=min(sapply(supports,function(x) x$a)),
            b=max(sapply(supports,function(x) x$b)),
            X=unique(unlist(lapply(supports,function(x) x$X))),
            lim=range(unlist(lapply(supports,function(x) x$lim)))
    )
}

#' @describeIn Support tries to identify pdfs and cdfs and useds
#'     pdf2Support or cdf2Support, a simple 0,1 otherwise.
#' @export
#' @method getSupport function
getSupport.function = function(x,...,extrems=c(-1E20,1E20),
                               starts=c(exp(seq(-20,20,by=0.2)),-exp(seq(-20,20,by=0.2)))
                               ) {
    e = x(extrems)
    y = x(sort(starts))
    if( all(is.finite(y)) && all(is.finite(e)) ) {
        if( isTRUE( all(e>=0) && all(e <= 1E-5*max(y)) ) ) {
                                        # assume its a density
            return(pdf2Support(x,...,starts=starts))
        } else if(isTRUE( all(all(e>=0) && 1-e[2] <= 1E-5 && all(diff(y)>=0) ) )) {
                                        # assume its a cdf
            return(cdf2Support(x),...)
        } else {
                                        # TODO
            Support(0,1)
        }
    } else {
                                        # TODO
        Support(0,1)
    }
}

#' @describeIn Support guesses the support through sample2Support
#' @export
#' @method getSupport numeric
getSupport.numeric <- function(x,...) {
    sample2Support(x,...)
}

#' @describeIn Support just calls NextMethod to hand over to the other classes
#' @export
#' @method getSupport ListData
getSupport.ListData <- function(x,...)  NextMethod("getSupport")

#' Computes Discretization points for nonequidistant discretization
#'
#' This is a helper function to create nonequidistant discretizations
#'
#' @param x typically a ReferenceMeasure
#' @param nEval the number of discretization points
#' @param ... further parameters with details for more complex references
#'
#' @returns an object that can be use together with its evaluation on object$all
#'          to create a DiscreteFunction
#' @export
#'
#' @examples
#' getDiscretizer(B2Unif01Basis5)
#' getDiscretizer(B2N01Basis5)
#'
getDiscretizer = function(x,nEval=513,...) UseMethod("getDiscretizer",x)
#' @describeIn getDiscretizer method for BayesSpaceBasis
#' @export
#' @method getDiscretizer BayesSpaceBasis
#'
getDiscretizer.BayesSpaceBasis = function(x,nEval=513,...) getDiscretizer(getReferenceMeasure(x),nEval=nEval,...)


#' @describeIn getDiscretizer method for ReferenceMeasure discrete
#' @export
#' @method getDiscretizer ReferenceMeasureDiscrete
#'
getDiscretizer.ReferenceMeasureDiscrete = function(x,nEval,...) structure(list(all=x$X,support=getSupport(x),transferCdf=NULL),class="Discretizer")
#' @describeIn getDiscretizer method for ReferenceMeasureUniform
#' @export
#' @method getDiscretizer ReferenceMeasureUniform
#'
getDiscretizer.ReferenceMeasureUniform = function(x,nEval=513,...) structure(list(all=seq(x$a,x$b,length.out=nEval),support=getSupport(x),transferCdf=NULL),class="Discretizer")
#' @describeIn getDiscretizer method for ReferenceMeasureTransformed
#' @export
#' @method getDiscretizer ReferenceMeasureTransformed
getDiscretizer.ReferenceMeasureTransformed = function(x,nEval=513,...) {
    ref = x$reference
    xI  = seq(ref$a,ref$b,length.out=nEval)
    deinfinitize=function(x) {
        if(!is.finite(x[1]))
           x[1]=x[2]-10
        if(!is.finite(x[length(x)]))
           x[length(x)]=x[length(x)-1]+10
        x
    }
    structure(list(all=deinfinitize(x$qf(xI)),
                   support=getSupport(x),
                   transferCdf=x$cdf),class="Discretizer")
}





    #' @export
#' @describeIn B2Unif01 the underlying standard reference measure
StdReferenceMeasureUnif01 = ReferenceMeasureUniform(0,1)
#' @export
#' @describeIn B2Unif01 the underlying standard reference measure
StdReferenceMeasureN01 = ReferenceMeasureTransformed("pnorm",mean=0,sd=1)
#' The standard B2(Unif(0,1)) Bayes space.
#'
#' B2Unif01 is the Bayes space B2(Unif(0,1)), and B2Unif01BasisN its B-Spline basis with N elements of order 2
#' @seealso [BayesSpaceBasis]
#' @family ReferenceMeasures
#' @export
#' @name B2Unif01
#' @examples
#' getReferenceMeasure(B2Unif01)
#' getMetricTensor(B2Unif01Basis5)
B2Unif01 = BayesSpace(referenceMeasure = StdReferenceMeasureUnif01)
#' @export
#' @describeIn B2Unif01 a basis for B2Unif01 with 5 elements
B2Unif01Basis5 = BayesSpaceBasis(B2Unif01,nbasis=5,order=4)
#' @export
#' @describeIn B2Unif01 a basis for B2Unif01 with 7 elements
B2Unif01Basis7 = BayesSpaceBasis(B2Unif01,nbasis=7,order=4)
#' @export
#' @describeIn B2Unif01 a basis for B2Unif01 with 20 elements
B2Unif01Basis20 = BayesSpaceBasis(B2Unif01,nbasis=20,order=4)

#' The standard B2(N(0,1)) Bayes space.
#'
#' B2N01 is the Bayes space B2(N(0,1)), and B2N01BasisN its B-Spline basis with N elements of order 2
#' @seealso [BayesSpaceBasis]
#' @family ReferenceMeasures
#' @export
#' @name B2N01
#' @examples
#' getReferenceMeasure(B2N01)
#' getMetricTensor(B2N01Basis5)
B2N01 = BayesSpace(referenceMeasure = StdReferenceMeasureN01)
#' @export
#' @describeIn B2N01 a basis for B2N01 with 5 elements
B2N01Basis5 = BayesSpaceBasis(B2N01,nbasis=5,order=4)
#' @export
#' @describeIn B2N01 a basis for B2N01 with 7 elements
B2N01Basis7 = BayesSpaceBasis(B2N01,nbasis=7,order=4)
#' @export
#' @describeIn B2N01 a basis for B2N01 with 20 elements
B2N01Basis20 = BayesSpaceBasis(B2N01,nbasis=20,order=4)


#' Discrete Reference Measure for BayesSpaces
#'
#' A discrete reference measure is a measure on a discrete support.
#'
#' @param X either a complete vector or list of all possible values or
#'     a subset that can be used as a sample. In this case the support
#'     shall be given.
#' @param p is either vector of probability values, or a function
#'     mapping a value to its probability. Currently only uniform
#'     distributions are supported
#' @param support the support to be used for the ReferenceMeasure
#' @param sampler a function of a number of random realizations,
#'          that samples ffrom this measure.
#'
#' @return An object of type ReferenceMeasureDiscrete encapsulating this measure
#'
#' @family ReferenceMeasures
#' @export
#' @examples
#' X <- c("A","B","C")
#' ( rmd <- ReferenceMeasureDiscrete(X) )
ReferenceMeasureDiscrete <- function(X,p=rep(1/length(X),length(X)),
                                    sampler=NULL,support=Support(X=X)) {
    .call = match.call()
    if( is.function(p) ) {
        pFunc=p
        p    =p(X)
    } else {
        stopifnot(length(p)==length(X))
        pFunc=function(x) {
            where = match(x,X,0)
            found = where != 0
            ans = numeric(length(x))
            if( any(found) )
                ans[found]=p[where[found]]
            ans
        }
    }
    w = pFunc(X)
    w = w/sum(w)
    return(structure(list(p=p,
                          pFunc=pFunc,
                          support=support,
                          w=pFunc(X),
                          X=X,
                          sampler=sampler),
                     call=.call,
                     class=c("ReferenceMeasureDiscrete","ReferenceMeasure")))
}

#' @describeIn expectation a computing the results explicitly
#' @export
#' @method expectation ReferenceMeasureDiscrete
#'
expectation.ReferenceMeasureDiscrete = function(measure,f,...,method=c("sample")) {
    Y = .sumFirstDimension(f(measure$X)*measure$w)
}
## Function summing over the first dimension
.sumFirstDimension <- function(X) {
    D = length(dim(X))
    if( D<=1 ) sum(X) else apply(X,2:D,sum)
}



#' @export
#' @describeIn CovarianceOperator method for ReferenceMeasuresDiscrete
#'     explicit computation
#' @importFrom stats cov
#' @method CovarianceOperator ReferenceMeasureDiscrete
CovarianceOperator.ReferenceMeasureDiscrete = function(measure,X,Y,...) {
    if( is.list(Y) ) {
        sapply(Y,function(Y1) CovarianceOperator(measure,X,as.function(Y1)))
    } else if( is.list(X) ) {
        sapply(X,function(X1) CovarianceOperator(measure,as.function(X1),Y))
    } else  {
        omega = measure$X
        .wCov(X(omega),Y(omega),measure$w)
    }
}


#' @describeIn weightedStats method for ReferenceMeasureDiscrete exact computation based on X and w of the measure.
#' @export
#'
weightedStats.ReferenceMeasureDiscrete = function(measure,X,f,method,tol) {
    x = measure$X
    XX = as.matrix(X(x))
    wt=c(f(x))*measure$w
    wt=wt/sum(wt)
    reportIfNot(nrow(XX)==length(wt),dimXX=dim(XX),lengthwt=length(wt),dimwt=dim(wt))
    val = stats::cov.wt(XX,wt=wt,method="ML")
    list(E=val$center,V=val$cov)
}
#' @export
#' @describeIn isInDomain method for discrete Reference Measures
isInDomain.ReferenceMeasureDiscrete = function(x,referenceMeasure) {
     x %in% referenceMeasure$X
}

#' Coded Basis for ReferenceMeasureDiscrete
#'
#' @param basis the ReferenceMeasureDiscrete the basis is designed for
#' @param ... further arguments for generic use only.
#' @param bayesSpace the BayesSpace to be used for the BayesSpaceBasis
#' @param orthonormalize boolean: if TRUE the basis is orthonormalized, using a singluar using a singular value decomposition based matrix square root.
#'
#' @returns A BayesSpaceBasis for that discrete reference measure. If
#'     a nonstandard A is used it might be nonorthogonal.
#' @importFrom compositions ilrBase
#' @export
#' @method BayesSpaceBasis ReferenceMeasureDiscrete
#' @examples
#' X=c("A","B","C")
#' rmd = ReferenceMeasureDiscrete(X)
#' basis = BayesSpaceBasis(rmd)
#' as.function(basis)(X)
#' weightedStats(rmd,as.function(basis),function(x) rep(1,length(x)))
# # not yet possible:
# #rmd = ReferenceMeasureDiscrete(X,c(0.2,0.3,0.5))
# #basis = BayesSpaceBasis(rmd)
# #as.function(basis)(X)
# #weightedStats(rmd,as.function(basis),function(x) rep(1,length(x)))
BayesSpaceBasis.ReferenceMeasureDiscrete =
    function(basis,
             bayesSpace=BayesSpace(basis),orthonormalize=FALSE,
             ...) {
        equal=get("==")
        D = length(basis$X)
        if( D == 1 ) {
            return(structure(BayesSpaceBasis(
                function(x) matrix(nrow=length(x),ncol=0),
                bayesSpace=bayesSpace,
                A=matrix(nrow=0,ncol=0),
                S=matrix(nrow=0,ncol=0),
                orthonormalize=orthonormalize,
                ...),
                call=match.call()))
        } else {
#            if( all(basis$w==basis$w[1])) {
#                A=t(compositions::ilrBase(D=length(basis$X)))
#                structure(
#                    BayesSpaceBasis(
#                        function(x) outer(x,basis$X,equal),
#                        bayesSpace=bayesSpace,
#                        A=A,
#                        S=A%*%t(A),
#                        orthonormalize=orthonormalize,
#                        ...),
#                    call=match.call()
#                )
#            } else {
                w = basis$w
                W = cumsum(w)
                A=outer(1:(D-1),1:D,function(j,t) ifelse(t>j+1,0,ifelse(t==j+1,+1/w[t],-1/W[j])*sqrt(W[j]*w[j+1]/W[j+1])))
                if(D>=3) {
                    S=diag(D-1)
                } else {
                    S=matrix(1,nrow=1,ncol=1)
                }
                structure(
                    BayesSpaceBasis(
                        function(x) outer(x,basis$X,equal),
                        bayesSpace=bayesSpace,
                        A=A,
                        S=S,
                        orthonormalize=orthonormalize,
                        ...),
                    call=match.call()
                )
#            }

        }
}


## ReferenceMeasureMixed