R/helpers.R

In DistributionalDDA: Distributional Data Analysis in Bayes Spaces

### functionality
# ' @export
# remove later to find missing TODO examples when checking package
#TODO = "TODO"

## a quick function generation tool
#' Shortcut for creating functions of x
#'
#' 'fx(body)' is essentially equal to 'function(x) body'. This can be used
#' to quickly create functions, with minimal syntactic burden.
#'
#' @param .body the expression that will form the body of the function, it
#'     is not evaluated.
#' @param template a template for the function, providing its formal
#'     parameters.
#' @param .env the environment in which the body of the function is
#'     evaluated.
#'
#' @returns a function with .body as body
#' @export
#' @importFrom splines splineDesign
#' @examples
#' f = fx(sin(x)+cos(x))
#' f
#' f(pi)
#'
#'
"fx" <- function(.body,template=function(x) {},.env=parent.frame()) {
   body(template,envir=.env)=substitute(.body)
   template
}

warp <- function(x,i) {
    x[[1+((i-1)%%length(x))]]
}

## equality comparisons
#' @export
`==.ReferenceMeasure` <- function(e1, e2){
   identical(e1,e2)
}


#' @export
`!=.ReferenceMeasure` <- function(e1, e2) !(e1==e2)

#' @export
"==.BayesSpace" = function(x,y){
  class(x)==class(y) && getReferenceMeasure(x)==getReferenceMeasure(y)
}
#' @export
"==.L2Space" = function(x,y){
  class(x)==class(y) && getReferenceMeasure(x)==getReferenceMeasure(y)
}

#' @export
"!=.BayesSpace" = function(x,y) !(x==y)

#' @export
"!=.L2Space" = function(x,y) !(x==y)

# small function to compare quantities that could be equal but for numerical imprecision
equiv = function(x,y, tol=1e-15) abs(x-y)<tol


#computeIntegral and computeScalarProduct seem to duplicate expectation and CovarianceOperator from Bayes_spaces
#computeIntegral
# ' @export
#computeIntegral = function(referenceMeasure,f,...) UseMethod("computeIntegral")
# ' @export
#computeIntegral.ReferenceMeasure = function(referenceMeasure,f,...){
#  mean(f(getSample(referenceMeasure)))
#}

#computeScalarProduct
# ' @export
#computeScalarProduct = function(space,f1,f2,...) UseMethod("computeScalarProduct",space)

# ' @export
#computeScalarProduct.BayesSpace = function(space,f1,f2,...) {
#  referenceMeasure = getReferenceMeasure(space)
#  f1mean = computeIntegral(referenceMeasure,function(x) log(f1(x)),...)
#  f2mean = computeIntegral(referenceMeasure,function(x) log(f2(x)),...)
#  computeIntegral(referenceMeasure,function(x) log(f1(x)-f1mean)*log(f2(x)-f2mean),...)
#}

# ' @export
#computeScalarProduct.ReferenceMeasureUniform =function(space,f1,f2,...){
#  referenceMeasure = getReferenceMeasure(space)
#  computeIntegral(referenceMeasure,function(x) f1(x)*f2(x),...)
#}


# ' @export
#evalBasis = function(basis,X) UseMethod("evalBasis",basis)
# ' @export
#evalBasis.BayesSpaceBasis = function(basis,X) {
#  with(basis,func(X) %*% t(A))
#}

# The evalBasis.BData functions seam not be clearly defined
# ' @export
##evalBasis.BDataSample = function(){}
# ' @export
##evalBasis.BDataPDF = function(){}
# ' @export
##evalBasis.BDataCDF = function(){}
StupidStuff <- function(x) x

#' returns the probability density function
#'
#' The probability density function, also known as the Lebesgue density
#' of the measure is returned. In case of datasets a dataset of such
#' functions is returned. (To be encoded as BDataPDF)
#'
#' @param X the distribution hows Lebesgue density shall be returned
#' @param ... potential futher arguments for the computational method
#'
#' @returns a function (or a list of functions if a dataset of
#'     distributions is given) that can compute the probably density
#'     function.
#'
#' @export
#'
#' @examples
#' x = seq(-1,2,by=0.1)
#' f1 = getPDF(StdReferenceMeasureUnif01)
#' plot(x,f1(x))
#' f2 = getPDF(StdReferenceMeasureN01)
#' plot(x,f2(x),type="l")
#' f3 = getPDF(BDataLogDensities(list(sin),basis=B2Unif01Basis5))
#' plot(x,f3(x),type="l")
#' getPDF(BDataLogDensities(list(sin,cos,\(x)(1)),basis=B2Unif01Basis5))
#'
getPDF = function(X,...) UseMethod("getPDF",X)

#' @export
#' @describeIn getPDF method for ReferenceMeasureUniform
#' @method getPDF ReferenceMeasureUniform
getPDF.ReferenceMeasureUniform = function(X,...) {
  f = function(x) ifelse(x<X$a | x>X$b, 0,1.0/(X$b-X$a))
  return(f)
}
#' @export
#' @describeIn getPDF method for ReferenceMeasureUniform
#' @method getPDF ReferenceMeasureTransformed
getPDF.ReferenceMeasureTransformed = function(X,...) X$pdf

#' @export
#' @describeIn getPDF method for BData objects, returns BDataPDFs (TODO)
#' @method getPDF BData
getPDF.BData = function(X,...) {
    pdfs = as.BDataPDFs(X,...)
    function(x) {
        erg = sapply(pdfs,function(pdf) pdf(x))
        if( ! is.matrix(erg) ) {
            if(numberOfCases(pdfs)==1 )
                erg = matrix(erg,ncol=1)
            else
                erg = matrix(erg,nrow=1)
        }
        erg
    }
}


#' Get a sample from the reference measure
#'
#' This provides a sample of data form the support of the reference
#' measure, e.g. to test functions or to be used as integration
#' points. The sample can optionally come as a weighted sample to
#' ensure integration. In order to achieve good performance and
#' reproducability the sample is typically fixed and not random,
#' except in cases where no good algorithm for a deterministic grid
#' could be implemented.
#'
#'
#' @param x the Reference measure, space for which the support is defined
#' @param ... further generic parameters, typically unused
#'
#' @family ReferenceMeasures
#' @export
#' @examples
#' getSample(B2Unif01)
#'
getSample = function(x,...) UseMethod("getSample")

#' @describeIn getSample method for referencemeasure, please overwrite
#' @export
#' @method getSample ReferenceMeasure
getSample.ReferenceMeasure = function(x,...) x$sample
#' @describeIn getSample method for ReferenceMeasureCartesianProduct
#' @export
#' @method getSample ReferenceMeasureCartesianProduct
getSample.ReferenceMeasureCartesianProduct = function(x,...) attr(x,"sample")

#' @describeIn getSample method for BayesSpace
#' @export
getSample.BayesSpace = function(x,...) getSample(getReferenceMeasure(x,...))
#' @describeIn getSample method for L2Space
#' @export
getSample.L2Space = function(x,...) getSample(getReferenceMeasure(x,...))



#' Checks whether two objects have the same BayesSpace
#'
#' The function fails if no BayesSpace is associated to any of the
#' objects.
#'
#' @param x some object with a BayesSpace
#' @param y some object with a BayesSpace
#'
#' @returns a boolean
#' @export
#'
#' @examples
#' hasSameBayesSpace(B2Unif01Basis5,B2Unif01Basis7) # TRUE
#' hasSameBayesSpace(B2Unif01Basis5,B2N01Basis5) # FALSE
hasSameBayesSpace = function(x,y) {
    b1=getBayesSpace(x)
    b2=getBayesSpace(y)
    stopifnot(!is.null(b1))
    stopifnot(!is.null(b2))
    #if( (b1$referenceMeasure != b2.referenceMeasure) ) return(FALSE)
                                        #TRUE
    if(identical(b1,b2))
        return(TRUE)
    return(b1$referenceMeasure == b2$referenceMeasure)
}

#' Checks whether the associated BayesSpaceBasis of to objects is the same
#'
#' @param x an object with an associated BayesSpaceBasis
#' @param y an object with an associated BayesSpaceBasis
#'
#' @returns a boolean
#'
#' @export
#'
#' @examples
#' x = BDataSamples(list(runif(100)), basis=B2Unif01Basis5)
#' hasSameBayesSpaceBasis(B2Unif01Basis5,x) # TRUE
#' hasSameBayesSpaceBasis(B2Unif01Basis5,B2Unif01Basis7) # FALSE
#' hasSameBayesSpace(B2Unif01Basis5,B2N01Basis5) # FALSE
#'
hasSameBayesSpaceBasis = function(x,y) {
    bX = getBayesSpaceBasis(x)
    bY = getBayesSpaceBasis(y)
    if( identical(bX,bY) )
        return(TRUE)
    if( ! hasSameBayesSpace(bX,bY) )
        return(FALSE)
    # TODO allow for a more loose equality
    return(FALSE)
}


# take a series of numeric vectors and generate the vector of sorted unique values
sortUnique = function(...) sort(unique(c(...)))

### Functions for projection of functions
.SplineCoefsFromLogDensity = function(f,basis,measure=getReferenceMeasure(basis)) {
    ## L2 projection onto the (non-orthonormal) ZB-spline basis: the coefficients
    ## solve the normal equations S c = <f, zeta> with S = getMetricTensor(basis)
    ## the Gram matrix, i.e. c = solve(S, <f, zeta>) (Issue 3). CovarianceOperator
    ## returns the inner products as an n x p matrix -- one row per density when f is
    ## vector-valued (the multi-observation / basis-change reproject path), a single
    ## 1 x p row otherwise. Solve per row via t(solve(S, t(.))), which preserves that
    ## n x p shape (a plain solve(S, as.vector(.)) would flatten and break n > 1).
    inner_products = CovarianceOperator(measure,f,as.function(basis))
    t(solve(getMetricTensor(basis), t(inner_products)))
  }

.SplineCoefsFromDensity = function(f,basis) {
    .SplineCoefsFromLogDensity(Ln(f),basis)
}

#' Compute the density relativ to a measure from histogram weights.
#'
#' The function does not normalize the density.
#'
#' @param weights a value proportional to the density in a histogram
#' @param breaks the breaks of the histogram in accending order. One
#'     more than weights.
#' @param background the value for the density to be chosen outside
#'     of break range.
#' @param dmu a vectorized function representing the Lebesgue density of
#'     the measure relativ to which the density should be valid. The
#'     default computes Lebesgue densities (probability density
#'     functions, PDFs)
#'
#' @returns a function that can be evaluated on the support of the
#'     underlying measure (given by dmu) and provides the a density at
#'     every point.
#'
#' @family .DensitFromHist
#'
#' @export
#'
#' @examples
#' f = .DensityFromHist(1:5,0:5)
#' x = seq(-1,6,by=0.01)
#' plot(x,f(x),type="b")
#'
.DensityFromHist = function(weights,breaks,background=0.0,
                            dmu=function(x) rep(1,numberOfCases(x))) {
    # weights[j] belongs to the j-th bin (breaks[j], breaks[j+1]];
    # findInterval() maps x to its bin index (0 below, length(breaks) above,
    # absorbed by the zero padding).  Versions <= 0.1.0 counted the breaks
    # from the wrong side, which applied the weights in reversed bin order.
    function(x) c((c(0,weights,0)+background)[findInterval(x,breaks,rightmost.closed=TRUE)+1]/dmu(x))
}

#' Compute the log density relativ to a measure from histogram weights.
#'
#' The function does not normalize the density.
#'
#' @param weights a value proportional to the density in a histogram
#' @param breaks the breaks of the histogram in accending order. One
#'     more than weights.
#' @param background a value to be added to ensure that the density is
#'     nowhere 0 and thus the log exists.
#' @param dmu a vectorized function representing the Lebesgue density of
#'     the measure relativ to which the density should be valid. The
#'     default computes Lebesgue densities (probability density
#'     functions, PDFs)
#'
#' @returns a function that can be evaluated on the support of the
#'     underlying measure (given by dmu) and provides the a density at
#'     every point.
#' @family .DensitFromHist
#' @export
#'
#' @examples
#' f = .LogDensityFromHist(1:5,0:5)
#' x = seq(-1,6,by=0.01)
#' plot(x,f(x),type="b")
#'
.LogDensityFromHist = function(weights,breaks,
                               background=mean(weights)/(length(weights)*10),
                               dmu=function(x) rep(1,numberOfCases(x))) {
    # same bin convention as .DensityFromHist (fixed in 0.1.1)
    function(x) c(log(c(0,weights,0)+background)[findInterval(x,breaks,rightmost.closed=TRUE)+1]-log(dmu(x)))
}


#' Coefficients of a single BayesSpace element.
#'
#' The coefficients are computed with respect to the given basis.
#' Only a single density data element my be provide, but in many
#' different froms including single case BData* objects, density
#' functions, samples, frequency and density histograms. A density
#' histogram is a histogram structure without counts. Samples need to
#' provide at least minSampleSize elements, to avoid simply
#' misinterpreting inadmissible object types as samples.
#'
#' @param x a single object representing a distribution. A function
#'     will be interpreted as density unless.
#' @param ... generic arguments supporting the converter. Either
#'     minSampleSize or relevant argumens of as.BDataSpline.
#' @param minSampleSize the minimum length admissible for samples
#' @param basis the BayesSpaceBasis to be used.
#'
#' @returns a numerical vector of coefficients in the given basis
#' @export
#'
#' @examples
#' getSplineCoefs(function(x) x^2,basis=B2Unif01Basis5)
#' getSplineCoefs(rbeta(1000,2,2),basis=B2Unif01Basis5)
#'
getSplineCoefs = function(x,...,basis) UseMethod("getSplineCoefs",x)

#' @describeIn getSplineCoefs the method for BDataSpline objects.
#' @export
getSplineCoefs.BDataSplines = function(x,...,basis=getBayesSpaceBasis(x)) {
    reportIfNot(numberOfCases(x)==1, dimx=dim(x), numberOfCases=numberOfCases(x),
                "getSplineCoefs only works on single observations")
    c(unclass(x))
}

#' @describeIn getSplineCoefs method for density functions. Convert log densities to densities by applying 'Exp()'.
#' @export
getSplineCoefs.function = function(x,...,basis=getBayesSpaceBasis(x)) {
    .SplineCoefsFromDensity(x,basis=basis)
}

#' @describeIn getSplineCoefs method for histograms
#' @export
getSplineCoefs.hist = function(x,...,basis=getBayesSpaceBasis(x)) {
    getSplineCoefs(
    as.BDataSplines(
        if( is.null(x$counts) )
            BDataDensityHistograms(x,basis=basis,...)
        else
            BDataFrequencyHistograms(x,basis=basis,...)
    ))
}

#' @describeIn getSplineCoefs method for distributions represented as
#'     a single element in a BData structure.
#' @export
getSplineCoefs.BData = function(x,...,basis=getBayesSpaceBasis(x)) {
    getSplineCoefs.BDataSplines(as.BDataSplines(x,basis=basis))
}

#' @describeIn getSplineCoefs interpreting it as a sample of data. It
#'     must have at least 10 observations.
#' @export
getSplineCoefs.default = function(x,...,minSampleSize=10,basis=getBayesSpaceBasis(x)) {
    reportIfNot(length(x)>=minSampleSize,
                lengthx=length(x),
                minSampleSize=minSampleSize,
                "At least minSampleSize observations to be condsidered a dataset.")
    getSplineCoefs.BData(BDataSamples(list(x),basis=basis))
}


reproject = function(coefs,basis,origbasis=getBayesSpaceBasis(coefs)) {
    .SplineCoefsFromLogDensity(coef2function(coefs,basis=origbasis),basis=basis)
}




#' A generic function computing the number of cases.
#'
#' Different representations of datasets provide its elements in a very
#' different way. This function provides a uniform interface to figure get
#' the number of cases in a dataset.
#'
#' @param x a dataset, of which the number of cases shall be
#'     counted. It can be a numeric vector, a matrix, a data.frame, a
#'     compositional dataset, or a BData object.
#'
#' @returns an integer
#'
#' @export
#'
#' @examples
#' data(iris)
#' numberOfCases(iris$Sepal.Length)
#' numberOfCases(iris)
#' numberOfCases(lapply(1:numberOfCases(iris),function(i) iris[i,]))
#'
numberOfCases = function(x) UseMethod("numberOfCases",x)

#' @describeIn numberOfCases method for dataframes
#' @export
numberOfCases.data.frame = function(x) nrow(x)
#' @describeIn numberOfCases method for lists and numeric, BDataDensities,...
#' @export
numberOfCases.default    = function(x)  if(is.matrix(x)) nrow(x) else length(x)


FisherScoringDensity <- function(X,basis,...) {
    theta = FisherScoringEstimation(X,basis,...)
    Exp(coef2function(theta,basis))
}

StableSolve <- function(A,b) {
    if( length(A) == 1 ) {
        stopifnot(length(b==1))
        return(b/A)
    } else {
        Asvd=svd(A)
        dmax=max(Asvd$d)
        di = with(Asvd,ifelse(abs(d)<max(abs(d))/1000,0,1/d))
        drop(with(Asvd,v%*%(diag(di)%*%(t(u)%*%b))))
    }
}

getDiag = function(A) {
    if( length(A)>1 ) return(diag(A)) else return(A)
}

makeDiag=function(A) {
    if( length(A)==1 ) matrix(A,nrow=1,ncol=1) else diag(A)
}

FisherScoringEstimation <- function(X,basis,
                                    ...,
                                    T=as.function(basis),
                                    barT = colMeans(T(X)),
                                    g = function(x) rep(1,numberOfCases(x)),
                                    niter=100,
                                    theta0=0*barT,
                                    measure=getReferenceMeasure(basis),
                                    smax=1,
                                    ballMax=1E3,
                                    minStep=smax*1E-16,
                                    verbose=function(...) {}
                                    ) {
    theta = theta0
    for(i in 1:niter) {
        verbose(i,"\n")
        Z = weightedStats(measure,T,Mul(g,Exp(coef2function(theta,basis))))
        Dstats=barT-Z$E
        D = StableSolve(Z$V,Dstats)
        lengthD=sqrt(sum(D^2))
        if( lengthD< minStep ) { # check for covergence
            relDelta = sqrt(sum((Dstats^2)/getDiag(Z$V)))
            if( relDelta>1E-6 )
                warning("Subspace convergence in Fisher Scoring",Dstats)
            verbose("Converged: lengthD=",lengthD,", relDelta=",relDelta,"\n")
            return(theta)
        }
        if( lengthD>smax ) { # don't go to far
            if( lengthD>10*smax) { # Warn if the step is really off
                warning("Large Fisher Scoring Steps",i,", lenD=",lengthD,", theta=",theta)
            }
            D=smax*D/lengthD
        }
        if( sum(theta^2)>ballMax^2 ) { # return early on divergence
            warning("Divergenze in Fisher Scoring",lengthD,theta,D)
            return(theta)
        }
        theta = theta + D
        verbose("theta=",theta,"\n")
    }
    warning("Non-convergence in Fisher Scoring: Theta=",theta,", barT=",barT,", C=",Z$E)
    theta
}

## inverts a symmetric matrix
invertVariance = function(V) {
    # TODO replace with svd based inversion
    # chol2inv(chol(V))
  .powM(V, -1)
}

.powM = function (V, alpha = 1/2){
  if (length(V) == 1) {
    aa = V^alpha
    dim(aa) = dim(V)
    return(aa)
  }
  with(svd(V), u %*% diag(d^alpha) %*% t(v))
}

## Selects the round value of two very similar values. If non is round
## the first is used. This is in internal function for numerical
## utilities trying to infere the boundaries of supports.
.theRoundValue = function(x1,x2,digits=10) {
    r1 = round(x1,10)
    r2 = round(x2,10)
    if( x2==r1 )
        x2
    else if( r1==r2 &&
             ((x1 <= r1 && r1 <= x2) ||
              (x1 >= r1 && r1 >= x2))
            )
        r1
    else
        x1
}

# ' Finds the point where the function gets 0 (or lower) in an interval.
# ' Entry conditions f(x)>0, f(x0)<=0
# .findZeroInInterval(function(x) 4-x,0,7)
# ' @export
.findZeroInInterval = function(f,x,x0,zero=0.0,jumpLimit=1E-10) {
    stopifnot(f(x)>zero)
    stopifnot(f(x0)<=zero)
    while( TRUE ) {
        xs = seq(x,x0,length.out=11)
        if( sum(diff(xs)!=0)<=1 ) {
                    return(.theRoundValue(x0,x))
        }
        y = f(xs)
        g = y>zero
        x0 = xs[!g][1]
        x  = xs[g][sum(g)]
    }
}


## Internal function
##
## Finds the point where the function gets 0
##
## Entry conditions f(x)>0
##
## @param f the function
## @param x the starting point, here f needs to be larger than zero
## @param d the direction given as 1 or -1
## @param limit the value after which infinity is reported
## @param zero the number considered 0
## @returns either the exact point where the function first reaches zero
##          or d*Inf if zero is not reached before limit.
##
## @examples
## .findZero(function(x) 17-x,0,1)
# ' @export
.findZero = function(f,x=0,d=1,limit=1E12,zero=0.0) {
    stopifnot(f(x)>zero)
    h=d
    x0=x+h
    while( abs(x0) < limit ) {
        if( f(x0) <= zero )
            return(.findZeroInInterval(f,x,x0,zero=zero))
        x0=x0+h
        h=2*h
    }
    return(d*Inf)
}

## Internal function
##
## Finds the point where an increasing function first reaches v
##
## @param cdf the function
## @param cdf v value to search
## @param x0 in search direction before the point
## @param x1 in search direction after the point
## @param goInf the absolute which we consider as infinitely far, if we do
##              not get the value in search direction
## @examples
##  .findSupportBoundaryUp(\(x)(ppois(x,2)))
##  .findSupportBoundaryUp(\(x)(pnorm(x,3)))
##  .findSupportBoundaryUp(\(x)(punif(x,3,7)))
# ' @export
.findSupportBoundaryUp = function(cdf,v=1,x0=-1,x1=1,goInf=1E20) {
   ## get an interval containing v
    ## find a value below <=v
    h = 1
    ## find a x0: cdf(x1)<v
    while( cdf(x0)>=v ) {
        x0=x0-h
        h=2*h
        if(abs(x0)>=goInf) {
            stop(".findSupportBoundaryUp: cdf does not go below ",v)
        }
    }
    ## find a value v<=cdf(x1)
    h = 1
    while( v > cdf(x1)  ) {
        x1=x1+h
        h=2*h
        if(abs(x1)>=goInf) {
            return(Inf)
        }
    }
    notHere = TRUE
    while(notHere) {
        ## here cdf(x0)<v<=cdf(x1)
        xm= (x0+x1)/2
        notHere = (xm!=x0)&&(xm!=x1)
        if( cdf(xm) < v )
            x0=xm
        else
            x1=xm
    }
    return(x1)
}

## Internal function
##
## Finds the point where an increasing function first leaves values v
##
## @param cdf the function
## @param cdf v value to search
## @param x0 in search direction before the point
## @param x1 in search direction after the point
## @param goInf the absolute which we consider as infinitely far, if we do
##              not get the value in search direction
## @examples
##  .findSupportBoundaryDown(\(x)(ppois(x,2)))
##  .findSupportBoundaryDown(\(x)(pnorm(x,3)))
##  .findSupportBoundaryDown(\(x)(punif(x,3,7)))
##  .findSupportBoundaryDown(\(x)(punif(x,pi,2*pi)))==pi
# ' @export
.findSupportBoundaryDown = function(cdf,v=0,x0=-1,x1=1,goInf=1E20,minJump=1E-4) {
   ## get an interval containing v
    ## find a value below <=v
    h = 1
    ## find a x0: cdf(x1)<=v
    while( cdf(x0)>v ) {
        x0=x0-h
        h=2*h
        if(abs(x0)>=goInf) {
            return(-Inf)
        }
    }
    ## find a value v<cdf(x1)
    h = 1
    while( v >= cdf(x1)  ) {
        x1=x1+h
        h=2*h
        if(abs(x1)>=goInf) {
            stop(".findSupportBoundaryDown: cdf never goes above ",v)
        }
    }
    notHere = TRUE
    while(notHere) {
        ## here cdf(x0)<=v<cdf(x1)
        xm= (x0+x1)/2
        notHere = (xm!=x0)&&(xm!=x1)
        if( cdf(xm) <= v )
            x0=xm
        else
            x1=xm
        #print(c(x0,x1,cdf(x0),cdf(x1)))
    }
    ## If the lowwer bound has a jump the upper value is probably correct
    if( cdf(x1)-cdf(x0)>minJump )
        return(x1)
    else
        return(x0)
}



#' Try to guess the support from cdf
#'
#' This function essentially tries to infere the support from a cdf.
#'
#' @param cdf the cummulative density function from which the support
#'     shall be guessed.
#' @param alpha the support cutoff probability for effective supports
#'     if the support is infinite
#' @param eps the computational precision for computing the boundaries
#' @param digits the number of digits the boundaries and jumps shall
#'     be rounded to.
#'
#' @returns a `Support` object describing the support inferred from the
#'     `cdf`. Its `a` and `b` elements give the (possibly infinite)
#'     lower and upper support boundaries, its `lim` element the
#'     effective lower and upper boundaries (a numeric vector of length
#'     two, cut at probability `alpha/2` and `1 - alpha/2` when the
#'     support is unbounded), and `X` any detected points of positive
#'     probability mass.
#'
#' @export
#' @family Support
#' @examples
#' cdf2Support(pnorm)
#' cdf2Support(punif)
#' cdf2Support(
#'   function(x){
#'   ifelse(x< -pi,0,((pmax(x,0)+pi)/pi)*0.5) +
#'       ifelse(x<0,0,0.25+(pmax(x,pi)/pi)*0.25)
#'   }
#' )
#'
cdf2Support = function(cdf,alpha=0.01,eps=1E-7,digits=5) {
    Min =.findSupportBoundaryDown(cdf,v=0.0)
    Max =.findSupportBoundaryUp(cdf,v=1.0)
    eMin=.findSupportBoundaryDown(cdf,v=alpha/2)
    eMax=.findSupportBoundaryUp(cdf,v=1-alpha/2)
    DeltaEff = eMax-eMin
    if( is.finite(Min) && (eMin-Min)<DeltaEff*0.5 )
        eMin=Min
    if( is.finite(Max) && (Max-eMax)<DeltaEff*0.5 )
        eMax=Max
    X = findJumps(cdf,eMin,eMax)
    ## TODO: What if the continues support is is smaller than the inner
    ## of the convect hull of the discrete.
    Support(Min,Max,X,c(eMin,eMax))
}

#' Finds jumps of a function of a certain minimum height
#'
#' This is a helper tool to infere supports from cummulative density
#' functions.
#'
#' @param f a vectorized function, whos jumps need to be found. The
#'     function is expected to be deterministic, nondecreasing (or
#'     nonincreasing) and continuous from the right.
#'
#' @param a a numeric value, the lower bound of the range to look for
#'     jumps
#' @param b a numeric value, the upper bound of the range to look for
#'     jumps
#' @param minJump a numeric value giving the minimum jump height, that
#'     is considered a jump.
#' @param eps a numeric value given the distance of evaluations. If
#'     the distance of two evaluations drops to epsilon, and the jump
#'     height is achieved the higher value is rounded to digits and
#'     returned. if eps=0, the search goes until machine precision.
#' @param digits an integer providing the digits of rounding, if
#'     search has finished accordint to the epsilon criterion.
#' @param n the numer discretization intervals considered in each
#'     iteration step.
#'
#' @returns a vector of numerical values at which the function f jumps
#'     by at least minJump in an eps distance, possibly rounded to
#'     digits.
#'
#' @export
#'
#' @examples
#' findJumps(function(x) pbinom(x,4,0.2),-pi,17)
#' findJumps(pnorm,-3,3)
#' p <- findJumps(function(x) ifelse(x<pi,x,x+pi),0,6,eps=0)
#' p==pi # TRUE
#'
#'
findJumps =function(f,a,b,minJump=1/100,
                     eps=1E-7,digits=5,n=100) {
    x  = seq(a,b,length.out=n+1)
    if( sum(diff(x)!=0.0)<2 ) {
        if( abs(f(a)-f(b))>=minJump ) {
            return(b)
        } else {
            return(c())
        }
    }
    dy = diff(f(x))
    possible = which(abs(dy)>=minJump)
    if( length(possible)==0 )
        return(numeric())
    if( length(possible)==1 && abs(a-b) <= eps )
        return(round(x[possible+1],digits))
    unlist(lapply(possible,function(i) findJumps(f,x[i],x[i+1],minJump=minJump,
                                                  eps=eps,digits=digits,n=n)))
}

#' Guesses the support from a sample
#'
#' A support is a range plus possible points with positive
#' probability. The function detects them by bindings.
#'
#' @param X a numerical vector with sample data
#' @param digits the number of digits infered boundaries should be
#'     rounded to.
#' @param minBindings the minimum number of equal values required to
#'     consider them as a value with positive measure.
#'
#' @returns a Support object
#'
#' @family Support
#' @export
#' @importFrom stats quantile
#'
#' @examples
#' sample2Support(rnorm(400))
#' sample2Support(runif(1000))
#' sample2Support(rlnorm(400))
#' sample2Support(rpois(400,2))
#' sample2Support(c(rlnorm(200),rep(0,30)))
#'
#'
sample2Support = function(X,digits=1,minBindings=max(length(X)/50,3)) {
    finits = is.finite(X)
    nf=unique( X[!finits] ) # not finits
    xf = X[finits]          # finits
    xd = sort(unique(xf[rowSums(outer(xf,xf,"=="))>minBindings])) # discrete
    xc = unique(xf[ !( xf %in% xd) ]) # continues
    if( length(xc) == 0 ) {
        return(Support(a=0,b=0,X=sort(c(xd,nf),lim=range(xd))))
    } else if( length(xc)==1 ) {
        return(Support(a=xc,b=xc,X=sort(c(xd,xc,nf)),
                       lim=range(c(xd,xc),na.rm=TRUE)))
    } else {
        ab = range(xc)
        f = 10^digits
        a = floor(ab[1]*f)/f
        b = ceiling(ab[2]*f)/f
        eab = c(stats::quantile(xc,c(0.025,0.975)))
        ea = floor(eab[1]*f)/f
        eb = ceiling(eab[2]*f)/f
        delta = eb-ab
        if( a > 0 && a <= ea / 2 ) {
            a=0
            ea=0
        }
        if( b < 0 && b >= ea / 2 ) {
            b=0
            eb=0
        }
        if( (b-eb)<0.3*(eb-ea) ) {
            eb=b
        }
        if( (ea-a)<0.3*(eb-ea) ) {
            eb=b
        }
        if( length(xd) > 0 ) {
            if( ea > min(xd) )
                ea=min(xd)
            if( eb < max(xd) )
                eb=max(xd)
        }
        Support(a,b,sort(c(xd,nf)),lim=c(ea,eb,use.names=FALSE))
    }
}

#' Guesses the support from a pdf
#'
#' The function is designed to allow a quick start into distributional
#' data analysis and plotting.
#'
#' @param pdf a probability density function of a real random variable
#' @param starts a vector of values, where the algorithm should look
#'     for a nonzero pdf. If the pdf is nowhere positive the
#'     computation will fail.
#' @param minDensityFactor a factor of the maximum density used to
#'     compute an effective range. should be a numeric in (0,0.5).
#' @param alpha a way to specify minDensityFactor by providing an
#'     alpha level for the normal distribution.
#'
#' @returns a Support object with the guessed support
#' @family Support
#' @importFrom stats nlm
#' @export
#' @examples
#' pdf2Support(dnorm)
#' pdf2Support(dunif)
#'
#'
pdf2Support = function(pdf,starts=exp(seq(-20,20,by=0.2)),minDensityFactor=dnorm(qnorm(alpha/2))/dnorm(0),alpha=0.01) {
    y = pdf(starts)
    if( all(y<=0) ) {
        starts = -starts
        y = pdf(starts)
    }
    if( all(y<=0) ) {
        stop("pdf2Support: please provide a staring value where the function is positive")
    }
    x = starts[which.max(y)]
    stopifnot( pdf(x)>0 )
    a = .findZero(pdf,x,-1)
    b = .findZero(pdf,x,1)
    erg = stats::nlm(function(x) -pdf(x),x)
    pdfMax = -erg$minimum
    x=erg$estimate

    ea = .findZero(pdf,x,-1,zero=pdfMax*minDensityFactor)
    eb = .findZero(pdf,x,1,zero=pdfMax*minDensityFactor)
    Support(a,b,lim=c(ea,eb))
}

cdf2qf = function(cdf,...,.args=list(...)) {
    min =.findSupportBoundaryDown(cdf,v=0)
    max =.findSupportBoundaryUp(cdf,v=1)
    x0 = if(min==-Inf) -1 else min
    x1 = if(max==Inf)  1 else max
    F  = function(x) do.call(cdf,c(list(x),.args))
    function(p) {
        K=1:length(p)
        while( any( Q<- (F(x0[K])>p[K]) )) {
            x0[K]=x0[K]*2
            K=K[Q]
        }
        K=1:length(p)
        while( any(Q <- (F(x1[K])<=p[K]))) {
            x1[K]=x1[K]*2
            K=K[Q]
        }
        K=1:length(p)
        p0 = F(x0)
        p1 = F(x1)
        pm = p0
        while( TRUE ) {
            xm[K] = x0 + (x1[K]-x0[K])*(p[K]-p0[K])/(p1[K]-p0[K])
            K = K[is.finite(xm[K]) & x0[K]<xm[K] & xm[K]<x1[K] & p0[K]<p1[K] ]
            if( length(K) == 0 )
                return( x0 )
            ym[K] = F(xm[K])
            lo = (pm[K] <= p[K])
            Klo = K[lo]
            Kup = K[!lo]
            if( length(Klo)>0 ) {
                x0[Klo]=xm[Klo]
                y0[Klo]=ym[Klo]
            }
            if( length(Kup)>0 ) {
                x1[Kup]=xm[Kup]
                y1[Kup]=ym[Kup]
            }
        }
    }
}


## Computer outer difference of datasets
outerDiff2 = function(x,y,Metric) {  # TODO: convert this to a C-function
  ii = expand.grid(ix=1:nrow(x),iy=1:nrow(y))
  ii = x[ii$ix,] - y[ii$iy,]
  ii = apply(ii,1,function(h) sqrt(t(h)%*%Metric%*%h) )
  dim(ii) = c(nrow(x), nrow(y))
  return(ii)
}
# outerDiff2 = function(x,y,Metric) {  # old version
#     outer(1:nrow(x),
#           1:nrow(y),
#           function(i,j) {
#               h = x[i,]-y[j,]
#               sqrt(t(h)%*%Metric%*%h)
#           }
#           )
# }





#' Computes an evaluable kernel density estimate on a dataset
#'
#' The distance is measured in the metric \ifelse{html}{d(x,y)= sqrt((x-y)^t M (x-y))}{\eqn{d(x,y)= \sqrt{(x-y)^t M (x-y)}}{d(x,y)= sqrt((x-y)^t M (x-y))}}
#' given by the 'Metric' tensor M. The kernel density is estimated as
#'
#' \ifelse{html}{f(x) = 1/N sum_i=1^n scale*kernel(d(x,X_i))}{\eqn{f(x) = \frac1N \sum_{i=1}^n scale*kernel(d(x,X_i))}{f(x) = 1/N sum_i=1^n scale*kernel(d(x,X_i))}}
#'
#' The computation is done explicity in all calls to the resulting function.
#'
#' @param X the sample to compute the kernel density estimated from
#' @param ... futher generic parameter
#' @param bandwidth the bandwidth to be used, if null the metric
#'     tensor is computed in a different way.
#' @param support the support on which to compute the kernel density
#'     estimate. Finite supports lead to a mirroring at the support limits. I future also an independent estimate of the density at discrete points is planned.
#' @param Metric the metric tensor used to scale the distances before
#'     passed to the kernel.
#' @param kernel a kernel function evaluated for a scaled distance. In
#'     case of univariate data the kernel is also evaluated with just
#'     the scaled difference (i.e. signed distances) and can thus be
#'     asymetric.
#' @param K the factor to scale the kernel density to get an integral
#'     of 1. The factor can be ignored in base of BData objects.
#' @param Var the variance (matrix) of the dataset to base a Metric
#'     calculation on
#' @param sdFactor the factor to multiply the directional standard
#'     deviation width to the the directional bandwidth.
#' @param minus the function to be used to compute the distance of
#'     elements
#'
#' @returns a function that can be evaluated with a dataset similar to X
#'
#' @export
#'
#' @examples
#' X = rnorm(10)
#' x=seq(-4,4,by=0.1)
#' f = kdeDensity(X,sdFactor=0.5)
#' plot(x,f(x),type="b")
#'
kdeDensity <- function(X,...,bandwidth,Metric,kernel,K,Var,sdFactor) UseMethod("kdeDensity",X)



#' @describeIn kdeDensity method for multivariate datasets provided as matrix
#' @importFrom compositions var
#' @importFrom stats dnorm
#'
#'
#' @export
kdeDensity.matrix =function(X,...,
                            bandwidth=NULL,
                            Metric=if( is.null(bandwidth) )
                                       invertVariance(Var*sdFactor^2)
                                   else
                                       diag(1/bandwidth^2,nrow=ncol(X)),
                            K=1/sqrt(2*pi)^(NCOL(X)-1)*sqrt(det(Metric)),   # corrected from K=1/sqrt(det(Metric))
                            kernel=dnorm,
                            Var=var(X),
                            sdFactor=(0.2)^(1/NCOL(X)) # corrected from sdFactor=0.2
                            ) {
  # warning("function not yet providing correct results")
    N = nrow(X)
    kernel # promises need to be evaluated to avoid unpredictability
    Metric
    K
    function(x) c(kernel(outerDiff2(x,X,Metric))%*%rep(K/N,N))
}

#' @describeIn kdeDensity method for multivariate datasets provided as data.frame
#' @export
#' @examples
#' # data generation and kde function construction
#' xx = cbind(rnorm(100, sd=2.5), rnorm(100, sd=0.5))
#' kdxy = kdeDensity(xx)
#' # check that the integral is approx 1:
#' x0 = seq(from=-2.5, to=2.5, length.out=25)
#' xys = expand.grid(x0*2.5, x0*0.5)
#' zs =  kdxy(xys)
#' (dx0 = x0[2]-x0[1])
#' (A = sum(zs)*dx0^2*2.5*0.5) # bad approx but within expectations given low precision
#' # compare with kde from stats::density for one marginal:
#' dim(zs) = c(25,25)
#' dx = density(xx[,1], kernel="gaussian")
#' dx |> plot()
#' lines(x0*2.5, apply(zs,1,sum)*dx0*0.5/A, col=2)
#' rug(xx[,1])
kdeDensity.data.frame = kdeDensity.matrix

#' @describeIn kdeDensity method for datasets provided as list. The
#'     minus operator must be defined on the elements.
#' @export
kdeDensity.list =function(X,...,
                          bandwidth=NULL,
                          Metric=if( is.null(bandwidth) )
                                     1/(Var*sdFactor^2)
                                 else
                                     1/bandwidth^2,
                          K=sqrt(Metric),   # corrected from K=1/sqrt(Metric)
                          kernel=dnorm,
                          Var=var(X),
                          sdFactor=0.2,
                          minus="-") {
    N = length(X)
    kernel # used promises need to be evaluated to avoid unpredictability
    Metric
    K
    function(x) c(kernel(outer(x,X,minus)*sqrt(Metric))%*%rep(K/N,N))
}


#' @export
#' @describeIn kdeDensity method for datasets provided as a numeric vector.
#' @examples
#' xx = rnorm(1000, sd=2.5) # simulate some data
#' # generate and plot the kde function
#' kdx = kdeDensity(xx)
#' plot(kdx, from=-5, to=5)
#' # ensure that the integral is approx 1
#' xs = seq(from=-5, to=5, by=0.01)
#' ys = kdx(xs)
#' sum(ys)*0.01             #  approx 1.0
#' # graphically compare with kde from stats::density
#' dx = density(xx, kernel="gaussian")
#' with(dx, lines(x, y, col=2))
#' plot(dx)
#' curve(kdx, from=-4, to=4, col=2, add=TRUE)
kdeDensity.numeric = function(X,...,
                          support=Support(a=-Inf,b=Inf,lim=range(X)),
                          bandwidth=NULL,
                          Metric=if( is.null(bandwidth) )
                                     1/(Var*sdFactor^2)
                                 else
                                     1/bandwidth^2,
                          K=sqrt(Metric),   # corrected from K=1/sqrt(Metric)
                          kernel=dnorm,
                          Var=var(X),
                          sdFactor=0.2,
                          minus="-") {
    N = length(X)
    kernel # used promises need to be evaluated to avoid unpredictability
    Metric
    minus
    # TBD: Handle discrete support
    kf = function(x) c(kernel(outer(x,X,minus)*sqrt(Metric))%*%rep(K/N,N))
    a = support$a
    b = support$b
    if( is.finite(a) ) {
        if( is.finite(b) ) {
            KF = function(x) ifelse(x>=a&x<=b,kf(x)+kf(2*a-x)+kf(2*b-x),0)
        } else {
            KF = function(x) ifelse(x>=a,kf(x)+kf(2*a-x),0)
        }
    } else if(is.finite(support$b) ) {
            KF = function(x) ifelse(x<=b,kf(x)+kf(2*b-x),0)
    } else {
            KF = kf
    }
    KF
}


### End Functions for projection of functions.

reportIfNot = function(.x,...,form=substitute(.x),info=list(...)) {
    if( !isTRUE(.x) ) {
      message("Failed check: The expression is not true:")
      message(form)
      message("Reporting")
      message(substitute(list(...)))
      message(info)
         stop("Check failed.")
    }
}

## applies a matrix A to a vector x or if x is a matrix, interpreted as rows of data on each of the
.operateAx2v = function(x,A) {
    if( is.matrix(x) ) {
        reportIfNot(ncol(x)==ncol(A),dimx=dim(x),dimA=dim(A))
        x%*%t(A)
    } else {
        reportIfNot(length(x)==ncol(A),dimx=dim(x),lengthx=length(x),dimA=dim(A))
        A%*%x
    }
}

.operatelx2v = function(x,l) {
    if( is.matrix(x) ) {

        c(x%*%l)
        }
    else
        c(l %*%x )
}



## Converts spline coefficients to a function
coef2function = function(f,basis=getBayesSpaceBasis(f)) {
    if( is.matrix(f) )
            function(x) .operateAx2v(basis$basis(x),basis$A)%*%t(f)
    else
            function(x) .operateAx2v(basis$basis(x),basis$A)%*%f
}


### computing transform matrix A
## renamed from computeA, because it needs to be public as a generic function
#' Computes the matrix mapping an L2-space basis to a Bayes-Space-Basis.
#'
#' This computes a (n-1)x1 matrix mapping the basis results of the L2 basis
#' to L^2_0 basis for a subspace of the Log-Space of a BayesSpace.
#'
#' @param x the L2-basis
#' @param ... further generic arguments potentially parametrizing the
#'     method of computation.
#'
#' @returns the transformation matrix `A` mapping the L2 (B-spline)
#'     basis to the ZB-spline basis of the Bayes space. For a B-spline
#'     basis with `nbasis` functions this is an `(nbasis - 1) x nbasis`
#'     matrix, i.e. one row per ZB-spline basis function.
#'
#' @export
computeL22BayesSpacesBasisMatrix = function(x,...) UseMethod("computeL22BayesSpacesBasisMatrix",x)

#' create the transform matrix A from a B-spline basis (fda::create.bspline.basis())
#'
#' @param x a 'basisfd' object
#' @param ... further unused generic arguments
#'
#' @returns matrix
#' @export
#'
#' @examples
#' computeL22BayesSpacesBasisMatrix(
#'          fda::create.bspline.basis(rangeval = c(0,1),
#'                                   nbasis=9,
#'                                   norder=4,
#'                                   breaks=seq(0,1,length.out=7)
#'                                   )
#'                                 )
#' # a cubic B-spline basis with 8 functions yields a 7 x 8 matrix
#' A = computeL22BayesSpacesBasisMatrix(
#'         fda::create.bspline.basis(c(0,1), nbasis=8, norder=4)
#'     )
#' dim(A)
computeL22BayesSpacesBasisMatrix.basisfd=function(x,...) {
    basis=x
                                        #basis=fda::create.bspline.basis(rangeval = c(x$a,x$b),nbasis=nbasis,norder=order,breaks=knots,...)
    nbasis=basis$nbasis
    breaks=c(basis$rangeval[1],basis$params,basis$rangeval[2])
    nbreaks=length(breaks)
    norder=nbasis-nbreaks+2
    length_inner=length(breaks)-2
    coincident_knots=c(rep(head(breaks,1),norder-1),breaks,rep(tail(breaks,1),norder-1))

    D=norder*diag(1/(tail(coincident_knots,length_inner+norder)-head(coincident_knots,length_inner+norder)))
    matK=diag(length_inner+norder-1)
    matK[cbind(2:(length_inner+norder-1), 1:(length_inner+norder-2))] = -1
    K=rbind(matK,c(rep(0,length_inner+norder-2),-1))
    A=t(D%*%K)
                                        #ZB=fda::fd(A,basis) -> functional dataobject for the ZB-spline basis
    reportIfNot(nrow(A)+1==ncol(A),dimA=dim(A))
    return(A)
}

# ' create the transform matrix A from a UniformReferenceMeasure object
# '
# ' @param x 'UniformReferenceMeasure' object
# ' @param nbasis number of basis functions
# ' @param order  order of basis functions
# ' @param knots  ascending nonrepetetive vector of knots
# ' @param ... further arguments to [fda::create.bspline.basis]
# '
# ' @returns matrix
# ' @export
# '
# ' @examples
# '
                                        # '
# I don't think this function should be needed
#computeL22BayesSpacesBasisMatrix.ReferenceMeasureUniform=function(x,nbasis,order=4,
#                                          knots=seq(x$a,x$b,length.out=(nbasis-order+2)),...){
#  basis=fda::create.bspline.basis(rangeval = c(x$a,x$b),nbasis=nbasis,norder=order,breaks=knots,...)

#computeL22BayesSpacesBasisMatrix(basis)
#}

#' Returns the argument, but ensures that it has a certain class.
#'
#' If the x has not the given type it, the function prints information
#' about the problem and stops with an error.
#'
#' @param x the objects that shall be checked.
#' @param type a string representing the type to be checked with 'is(x,type)'
#'
#' @returns x
#' @export
#' @examples
#' ensureType(3,"numeric")
#' ensureType(diag(2),"matrix")
#'
ensureType = function(x,type) {
    if( !is(x,type) ) {
      message("EnsureType error")
      message(x)
      message(class(x))
         stop("Object should have class", type," but has ", class(x) ,"in ensureTye")
    }
    x
}


#' Computes the scalar product of log densities in a BayesSpace
#'
#' This uses the Covariance Operator capability of the
#' ReferenceMeasure, x and y are converted to functions.
#'
#' @param x a log density function (centering is not required)
#' @param y a log density function (centering is not required)
#' @param space a BayesSpace (or anything that has a ReferenceMeasure)
#'
#' @returns the scalar product of the (centred) log-densities `x` and
#'     `y` with respect to the reference measure of `space`, computed
#'     via the covariance operator of that measure. A single numeric
#'     value for scalar functions, or a matrix of pairwise scalar
#'     products when `x` and `y` are vector-valued (e.g. basis
#'     functions).
#'
#' @export
#' @examples
#' scalarProductOfLogDensities(function(x) rep(1,length(x)),space=B2Unif01)
#'
scalarProductOfLogDensities = function(x,y=x,space) {
    x = as.function(x)
    y = as.function(y)
    measure = getReferenceMeasure(space)
    CovarianceOperator(measure,x,y)
}


#' numerically compute a density from a cdf
#'
#' The function tries to keep the (open) support, by using a one sided
#' derivative at limits of the support, while using a two sided in the
#' middle of the support.
#'
#' @param cdf the cummulative distribution function to be derived
#' @param ... futher named parameters for cdf
#' @param .args the further named parameters given as a list
#' @param eps the distance to be used in the numerical derivative
#' @param cdf2 a function representing the cdf of the reference
#'     measure. NULL uses a simplification for Lebesque reference.
#'
#' @returns a vectorized function giving the density obtained by
#'     numerically differentiating `cdf` (a one-sided difference at the
#'     boundaries of the support, a two-sided difference in the
#'     interior). If `cdf2` is supplied the result is the density
#'     relative to that reference measure instead of the Lebesgue
#'     density.
#'
#' @export
#'
#' @examples
#' dUnif = cdf2density(punif,min=-0.2,0.4)
#' x=seq(-1,2,by=0.01)
#' plot(x,dUnif(x))
#'
cdf2density = function(cdf,...,.args=list(...),
                       eps=1E-7,
                       cdf2=NULL) {
    if( is.null(cdf2) )
        function(x) {
            cdfX = evalFArgs(cdf,x,.args=.args)
            cdfXp= evalFArgs(cdf,x+eps,.args=.args)
            cdfXm= evalFArgs(cdf,x-eps,.args=.args)
            ifelse(cdfXm<=0.0,
            ifelse(cdfX<=0.0,
                   0,
                   (cdfXp-cdfX)/eps
                   ),
            ifelse(cdfXp>=1.0,
            (cdfX-cdfXm)/eps,
            (cdfXp-cdfXm)/(2*eps)
            )
            )
        }
    else
        function(x) {
            cdfX = evalFArgs(cdf,x,.args=.args)
            cdfXp= evalFArgs(cdf,x+eps,.args=.args)
            cdfXm= evalFArgs(cdf,x-eps,.args=.args)
            cdf2x = cdf2(x)
            cdf2p= cdf2(x+eps)
            cdf2m= cdf2(x-eps)
            ifelse(cdfXm<=0.0,
            ifelse(cdfX<=0.0,
                   0,
                   (cdfXp-cdfX)/(cdf2p-cdf2x)
                   ),
            ifelse(cdfXp>=1.0,
            (cdfX-cdfXm)/(cdf2x-cdf2m),
            (cdfXp-cdfXm)/(cdf2p-cdf2m)
            )
            )
        }

}

#' Closure methods for BData objects
#'
#' Closure functions to enforce a distributional representation to integrate to
#' one, if needed.
#'
#' @param X distribution representation, typically an object of a derived class of "BData"
#' @param ... further arguments, typically ignores
#'
#' @return A fitting object with elements now integrating to one, depending on
#' the specific method.
#'
#' @importFrom compositions clo
#' @export
#' @method clo BDataHistograms
#' @examples
#' # close frequency histograms to density histograms integrating to one
#' counts = rbind(c(10,20,30,40), c(25,25,25,25))
#' fh = BDataFrequencyHistograms(counts, breaks = seq(0,1,by=0.25),
#'                               basis = B2Unif01Basis5)
#' clo(fh)
#'
clo.BDataHistograms <- function(X,...){
  dd = clo(unclass(X))
  BDataDensityHistograms(dd, breaks=getBreaks(X), basis=getBayesSpaceBasis(X))
}


evalFArgs = function(.f,...,.args=NULL) do.call(.f,c(list(...),.args))

setFArgs = function(.f,...,.args=NULL) {
    if( length(.args) > 0 ) {
        reportIfNot(all(names(args) %in% names(formals(.f))), names(args)[!(names(args) %in% names(formals(.f)))] )
        formals(.f)[names(.args)]=.args
    }
    .f
}


simpleClo = function(x) {
    if( is.matrix(x) )
        x/rowSums(x)
    else
        t(c(x/sum(x)))
}

deNanFunction = function(f,NaNValue=0) {
    if( is.list(f) )
        return(lapply(f,deNanFunction,NaNValue=0))
    function(x) {
        y = f(x)
        good=is.finite(y)
        if(!all(good))
            y[!good]=NaNValue
        y
    }
}

#' Finds the closing constant(s)
#'
#' A closing constant is a real value, by which the function
#' (or its Exp in case of log functions) needs to be divided
#' to get the representative with measure 1.
#'
#' @param f a representation of a measure
#' @param ... further arguments for generic use
#' @param measure the ReferenceMeasure, needed especially for functions
#'
#' @returns a (vector) of numeric closing constants. Might contain Inf
#'     in case of infinite measures and 0 in case of functions that
#'     are constant zeros
#'
#' @export
#'
#' @examples
#' # closing constant of f(x) = x^2 on Uniform(0,1)
#' closingConstant(function(x) x^2, measure = ReferenceMeasureUniform(0,1))
#'
closingConstant = function(f,...) UseMethod("closingConstant",f)
#' @describeIn closingConstant method for function
#' @export
#' @method closingConstant function
closingConstant.function = function(f,measure=getReferenceMeasure(f),...) {
    expectation(measure,deNanFunction(f),...)
}
#' @describeIn closingConstant method for BayesSpaceBasis
#' @export
#' @method closingConstant BayesSpaceBasis
closingConstant.BayesSpaceBasis = function(f,measure=getReferenceMeasure(f),...) {
    expectation(measure,deNanFunction(Exp(as.function(f))),...)
}
#' @describeIn closingConstant method for BDataDensities
#' @export
#' @method closingConstant BDataDensities
closingConstant.BDataDensities = function(f,measure=getReferenceMeasure(f),...) {
    sapply(deNanFunction(f),closingConstant,measure=measure,...)
}
#' @describeIn closingConstant method for BDataLogDensities
#' @export
#' @method closingConstant BDataLogDensities
closingConstant.BDataLogDensities = function(f,measure=getReferenceMeasure(f),...) {
    sapply(deNanFunction(Exp(unclass(f))),closingConstant,measure=measure,...)
}
#' @describeIn closingConstant method for BDataLogPDFs
#' @export
#' @method closingConstant BDataLogPDFs
closingConstant.BDataLogPDFs = function(f,measure=getReferenceMeasure(f),...) {
    lapply(deNanFunction(Div(Exp(unclass(f)),getPDF(measure))),closingConstant,measure=measure,...)
}
#' @describeIn closingConstant method for BDataPDFs
#' @export
#' @method closingConstant BDataPDFs
closingConstant.BDataPDFs = function(f,measure=getReferenceMeasure(f),...) {
    sapply(deNanFunction(Div(unclass(f),getPDF(measure))),closingConstant,measure=measure,...)
}
#' @describeIn closingConstant method for BDataSplines
#' @export
#' @method closingConstant BDataSplines
closingConstant.BDataSplines = function(f,measure=getReferenceMeasure(f),...) {
    sapply(deNanFunction(as.BDataDensities(f)),closingConstant,measure=measure,...)
}

## TODO CenteringConstant DONE

#' Finds the centering constant(s)
#'
#' A centering constant is a real value, by which needs to be
#' substracted from a function (of its log in case of densities)
#' get an integral of 0.
#'
#' @param f a representation of a measure
#' @param ... further arguments for generic use
#' @param measure the ReferenceMeasure, needed especially for functions
#'
#' @returns a (vector) of numeric centering constants. Might contain Inf
#'     in case of infinite measures and 0 in case of functions that
#'     are constant zeros
#'
#' @export
#'
#' @examples
#' # centering constant of f(x) = x^2 on Uniform(0,1)
#' centeringConstant(function(x) x^2, measure = ReferenceMeasureUniform(0,1))
#'
centeringConstant = function(f,...) UseMethod("centeringConstant",f)
#' @describeIn centeringConstant method for function
#' @export
#' @method centeringConstant function
centeringConstant.function = function(f,measure=getReferenceMeasure(f),...) {
    expectation(measure,deNanFunction(f),...)
}
#' @describeIn centeringConstant method for BayesSpaceBasis
#' @export
#' @method centeringConstant BayesSpaceBasis
centeringConstant.BayesSpaceBasis = function(f,measure=getReferenceMeasure(f),...) {
    expectation(measure,deNanFunction(as.function(f)),...)
}
#' @describeIn centeringConstant method for BDataDensities
#' @export
#' @method centeringConstant BDataDensities
centeringConstant.BDataDensities = function(f,measure=getReferenceMeasure(f),...) {
    sapply(Ln(deNanFunction(f)),centeringConstant,measure=measure,...)
}
#' @describeIn centeringConstant method for BDataLogDensities
#' @export
#' @method centeringConstant BDataLogDensities
centeringConstant.BDataLogDensities = function(f,measure=getReferenceMeasure(f),...) {
    sapply(deNanFunction(unclass(f)),centeringConstant,measure=measure,...)
}
#' @describeIn centeringConstant method for BDataLogPDFs
#' @export
#' @method centeringConstant BDataLogPDFs
centeringConstant.BDataLogPDFs = function(f,measure=getReferenceMeasure(f),...) {
    lapply(deNanFunction(Sub(unclass(f),Ln(getPDF(measure)))),centeringConstant,measure=measure,...)
}
#' @describeIn centeringConstant method for BDataPDFs
#' @export
#' @method centeringConstant BDataPDFs
centeringConstant.BDataPDFs = function(f,measure=getReferenceMeasure(f),...) {
    sapply(deNanFunction(Ln(Div(unclass(f),getPDF(measure)))),centeringConstant,measure=measure,...)
}
#' @describeIn centeringConstant method for BDataSplines
#' @export
#' @method centeringConstant BDataSplines
centeringConstant.BDataSplines = function(f,measure=getReferenceMeasure(f),...) {
    rep(0,nrow(f))
}



## TODO cen
#' Centers the representation
#'
#' This centers an representation in a way that a direct converstion to a log Density has integral 0.
#'
#' @param X the object to be centered
#' @param ... further arguments for generic use
#' @param measure the reference measure for which the closing constant
#'     shall be computed. This is especially required for functions
#' @param center a boolean, if FALSE the object is returned directly,
#'     if true uncentered objects will be centered.
#'
#' @returns and object of the same type as the incoming, but closed.
#'
#' @details
#' Centring subtracts the integral mean (with respect to the reference measure) so
#' that the result integrates to zero; equivalently it makes the clr representative
#' of the distribution centred.
#'
#' @seealso [clo]
#' @export
#'
#' @examples
#' cen(function(x) x^2,measure=ReferenceMeasureUniform(0,1))
#'
cen = function(X,...,center=TRUE) if(center&&!isCentered(X)) UseMethod("cen",X) else X
#' @describeIn cen method for functions
#' @export
#' @method cen function
cen.function = function(X,measure=getReferenceMeasure(X),...,center) {
    Sub(unclass(X),centeringConstant(X,measure=measure))
}
#' @describeIn cen method for BDataDensities
#' @export
#' @method cen BDataDensities
cen.BDataDensities = function(X,...,center) {
    BDataDensities(Div(unclass(X),exp(centeringConstant(X,...))),
                   basis=getBayesSpaceBasis(X),centered=TRUE)
}
#' @describeIn cen method for BDataPDFs
#' @export
#' @method cen BDataPDFs
cen.BDataPDFs = function(X,...,center) {
    BDataPDFs(Div(unclass(X),exp(centeringConstant(X,...))),
                   basis=getBayesSpaceBasis(X),centered=TRUE)
}
#' @describeIn cen method for BDataLogDensities
#' @export
#' @method cen BDataLogDensities
cen.BDataLogDensities = function(X,...,center) {
    BDataLogDensities(Sub(unclass(X),centeringConstant(X,...)),
                      basis=getBayesSpaceBasis(X),centered=TRUE)
}
#' @describeIn cen method for BDataLogPDFs
#' @export
#' @method cen BDataLogPDFs
cen.BDataLogPDFs = function(X,...,center) {
    BDataLogPDFs(Sub(unclass(X),centeringConstant(X,...)),
                      basis=getBayesSpaceBasis(X),centered=TRUE)
}

#' @describeIn cen method for BayesSpaceBasis returns the centered logDensities
#' @export
#' @method cen BayesSpaceBasis
cen.BayesSpaceBasis = function(X,...,center) as.BDataLogDensities(X,...)


## TODO remember that something is centered or closed? Probably Done, needs checks

#' compute the clr transform
#'
#' This is a generic version of clr overwriting the one in compositions
#'
#' @param x the object to be clr transformed
#' @param ... further arguments for generic use
#'
#' @returns the centred log-ratio (clr) representation of `x`. For
#'     `BData` objects this is the centred log-density representation
#'     (an object of class `BDataLogDensities` whose log-density
#'     integrates to zero); the default method returns the clr
#'     transform computed by [compositions::clr].
#'
#' @details
#' The centred log-ratio map sends a distribution to its log-density centred to
#' integrate to zero, i.e. to its coordinates in the L2 space underlying the Bayes
#' space. It is an isometry from the Bayes space onto the centred subspace of that
#' L2 space, so Bayes-space distances and inner products become the ordinary L2 ones.
#'
#' @export
#'
#' @examples
#' X1 = BDataPDFs(basis=B2Unif01Basis5,
#'    list(
#'     \(x)(dbeta(x,1,2)),
#'     \(x)(dbeta(x,1,1)),
#'     \(x)(dbeta(x,1.5,0.7))
#'    ) )
#' plot(X1,main="Beta Densities on Unif(0,1)")
#' plot(clr(X1),main="CLRs of Beta Densities on Unif(0,1)")
#'
clr = function(x,...) UseMethod("clr",x)
#' @describeIn clr the default is [compositions::clr]
#' @export
#' @method clr default

clr.default = function(x,...) compositions::clr(x,...)
#' @describeIn clr returns a centered representation as BDataLogDensities
#' @export
#' @method clr BData
clr.BData = function(x,...) cen(as.BDataLogDensities(x,...))

#' Closes to a given total measure
#'
#' This is a generic version of clo overwriting the one in compositions
#'
#' @param X the object to be closed
#' @param ... further arguments for generic use
#' @param measure the reference measure for which the closing constant
#'     shall be computed. This is especially required for functions
#' @param close a boolean, if FALSE the object is returned directly,
#'     if true unclosed objects will be closed.
#' @param total the integral that is desired
#'
#' @returns and object of the same type as the incoming, but closed.
#'
#' @details
#' Closure rescales a (positive) function so that it integrates to `total` (one by
#' default) with respect to the reference measure, i.e. it projects the function to a
#' proper density and thereby picks the representative of its Bayes-space class.
#'
#' @seealso [compositions::clo] for parameters of the default, [cen]
#' @export
#'
#' @examples
#' clo(function(x) x^2,measure=ReferenceMeasureUniform(0,1))
#'
clo = function(X,...,close=TRUE) if(close&&!isClosed(X)) UseMethod("clo",X) else X
#' @describeIn clo the default is [compositions::clo]
#' @export
#' @method clo default
clo.default = function(X,...,total=1,close) compositions::clo(X,...,total=total)
#' @describeIn clo method for functions
#' @export
#' @method clo function
clo.function = function(X,measure=getReferenceMeasure(X),...,total=1,close) {
    Div(unclass(X),Div(closingConstant(X,measure=measure),total))
}
#' @describeIn clo method for BDataDensities
#' @export
#' @method clo BDataDensities
clo.BDataDensities = function(X,...,total=1,close) {
    BDataDensities(Div(unclass(X),closingConstant(X,...)/total),
                   basis=getBayesSpaceBasis(X),closed=TRUE)
}
#' @describeIn clo method for BDataPDFs
#' @export
#' @method clo BDataPDFs
clo.BDataPDFs = function(X,...,total=1,close) {
    BDataPDFs(Div(unclass(X),Div(closingConstant(X,...),total)),
                   basis=getBayesSpaceBasis(X),closed=TRUE)
}
#' @describeIn clo method for BDataLogDensities
#' @export
#' @method clo BDataLogDensities
clo.BDataLogDensities = function(X,...,total=1,close) {
    BDataLogDensities(Sub(unclass(X),Ln(Div(closingConstant(X,...),total))),
                      basis=getBayesSpaceBasis(X),closed=TRUE)
}
#' @describeIn clo method for BDataLogPDFs
#' @export
#' @method clo BDataLogPDFs
clo.BDataLogPDFs = function(X,...,total=1,close) {
    BDataLogPDFs(Sub(unclass(X),Ln(Div(closingConstant(X,...),total))),
                      basis=getBayesSpaceBasis(X),closed=TRUE)
}


rowVecSub = function(x,y) {
    if( is.matrix(x) ) {
        if( is.matrix(y)) {
            x-y
        } else {
            t(t(x)-y)
        }
    } else if(is.matrix(y)) {
        t(x-t(y))
    } else {
        x-y
    }
}

#' @describeIn clo method for BayesSpaceBasis returns the closed logDensities
#' @export
#' @method clo BayesSpaceBasis
clo.BayesSpaceBasis = function(X,...,total=1) {
    C=log(closingConstant(X)/total)
    f = as.function(X)
    function(x) rowVecSub(f(x),C)
}

#' Creates a function based on a regular discretization
#'
#' A simple fast tool for defining a discretized function defined
#' continuously on [a,b] and discrete values in X, with equidistant
#' data values y in between.  The resulting function is vectorized and
#' interpolates linearily with errors or NA outside the region. It
#' handles infinite supports, but not infinite results on the
#' continuous region.
#'
#' @param y the function to be approximated, or evaluation values,
#'     first those for X and than those equidistantly spaced xmin and
#'     xmax (in the image of transferCdf), or ReferenceMeasure on which to define
#'     the function (in this case the function is reported is given in f)
#' @param Y replaces y, if y is a Discretizer
#' @param f the function to be approximated (only if y is a ReferenceMeasure)
#' @param a the lower limit of the domain
#' @param b the upper limit of the domain
#' @param n the number of discretization points on [a,b]
#' @param X the discrete elements of the support (i.e. positions,
#'     where the function may jump.
#' @param support the [Support] on which the function shall be
#'     defined.
#' @param ... further arguments for Support (a,b,X)
#' @param transferCdf a transfer function mapping the interval a,b to [0,1]
#'     (or another interval), in which the nots are place
#'     equidistantly and are interpolated linearily.
#' @param transferQf (qunantile function), must be the inverse of transferCdf and is
#'     used to determine the evaluation points for the function.
#'
#'
#'
#' @returns a function interpolating linearly between the given
#'     evaluation points and working in a vectorized form.
#'
#' @export
#'
#' @examples
#' x=seq(0,1,length.out=100)
#' f = DiscretizedFunction(sin(x),0,1)
#' x=seq(0,1,length.out=1000)
#' max(abs(f(x)-sin(x)))
DiscretizedFunction = function(y,...) UseMethod("DiscretizedFunction",y)

#' @describeIn DiscretizedFunction method for precomputed values
#' @export
#' @method DiscretizedFunction numeric
DiscretizedFunction.numeric = function(y,...,support=Support(...),transferCdf=NULL) {
    # TODO correct for discontinuity in the support
             X = support$X
             if( length(X)>0 ) {
                 if(support$a==support$b) { # no continues part
                     return(InnerDiscretizedFunctionDiscrete(y,X))
                 } else { # continues part
                     finner=InnerDiscretizedFunction(y[-(1:length(X))],support$a,support$b,transferCdf=transferCdf)
                     y=y[1:length(X)] # remove unnecaary data from namespace
                     return(function(x) {
                         what = match(x,X,nomatch=0)
                         cont = what==0
                         if(any(cont)) {
                             yr=rep(0,length(x))
                             yr[cont] = finner(x[cont])
                             yr[!cont] = y[what]
                         } else {
                             y[what]
                         }
                     }
                     )
                 }
             } else
                 return(InnerDiscretizedFunction(y,
                                                 support$a,
                                                 support$b,
                                                 transferCdf=transferCdf))
 }

## This function returns a lookup returning the corresponding value
## of y for each value in X. It is here to keep the namespace
## stored with the lookup clean.
InnerDiscretizedFunctionDiscrete = function(y,X) {
    ## TODO generalize for arbitrary objects
    function(x) y[match(x,X)]
}

## This function returns function returning a vectorized linear
## interpoliation from equispaced lookup table, possibly after a
## transform. NULL stands for no transform and is implemented more
## efficiently.
InnerDiscretizedFunction = function(y,a,b,transferCdf=NULL) {
    if( is.null(transferCdf) ) {
        function(x) {
            k = 1+(length(y)-1)*((x-a)*(b-a))
            l = floor(k)
            u = ceiling(k)
            q=k-l
            y[l]*(1-q)+y[u]*q
        }
    } else {
        a = transferCdf(a)
        b = transferCdf(b)
        function(x) {
            k = 1+(length(y)-1)*((transferCdf(x)-a)*(b-a))
            l = floor(k)
            u = ceiling(k)
            q=k-l
            y[l]*(1-q)+y[u]*q
        }
    }
}



#' computes the limit of a function
#'
#' @param f the function whos limit shall be computed
#' @param x the location, where the limit shall be computed
#' @param d the direction of the limit
#'
#' @returns the limit it will still return
#'
#' @export
#'
#' @examples
#' computeLimit(\(x)(2*(x-pi)/(x-pi)),x=pi)
#' np = computeLimit(\(x)(pi*(x-pi)/(x-pi)),x=pi)
#' np-pi
#'
#'
computeLimit = function(f,x,d=1) {
    ## The multiplyer to look for a place where things start to be
    ## look like a function
    d=sign(d)
    if( !is.finite(d) || d==0 )
        d=1
    m = 1000
    if( x==Inf ) {
        fo=f
        f = function(x) fo(1/x)
        d=-1
        x=0
    } else if( x==-Inf ) {
        fo=f
        f = function(x) fo(1/x)
        d=1
        x=0
    }
    while(TRUE) {
        ## Look from here to x
        if( abs(x)>=1 )
            xa = x*(1+d*m*.Machine$double.eps)
        else
            xa = x+d*m*.Machine$double.eps
        if( !is.finite(xa) )
            stop("computeLimit failed identifying a point with a value")
        xs = seq(x,xa,length.out=11)[-1]
        ys = f(xs)
        finite =is.finite(xs) & xs!=x
        ## We would like to see a
        ## a gradient, to be able to extend it.
        if( all(finite) && (ys[1]!=ys[length(ys)] || m*.Machine$double.eps > 1000 )) {
            yo=ys[10]
            xo=xs[10]
            dydx = (ys[1]-ys[10])/(xs[1]-xs[10])

            ## Loop contracting the interval again to find the nearest
            ## computable value.
            while(TRUE) {
                finite =is.finite(ys)
                xa=xs[finite][1]
                ya=ys[finite][1]
                if( sum(finite) < 10 ) {
                    if( ya > .Machine$double.xmax/2 )
                        return(Inf)
                    if( ya < -(.Machine$double.xmax/2) )
                        return(-Inf)
                    ## the possible value trespasing a boundary is computed
                    ## based on the active gradient
                    return(.theRoundValue(ya,ya+2*abs(xa-xo)*dydx,digits=10))
                }
                if(abs((ya-yo)/yo)>.Machine$double.eps*2) {
                    ndydx = (ya-yo)/(xa-xo)
                    yo=ya
                    xo=xa
                    if( abs(ndydx)>abs(dydx) || sign(ndydx)<sign(dydx) )
                        dydx=ndydx
                    if( abs(ndydx)>0 && abs(ndydx)<abs(dydx)/2 )
                        dydx=ndydx
                }
                xs=seq(x,xa,length.out=11)[-1]
                ys=f(xs)
            }
        }
        m=m*10
    }

}

#' Makes a real function returning limits at undefined points
#'
#' The function identifies Punctures in the definition of a function,
#' by getting a non finite result (Inf,-Inf,Nan) and than asks
#' computeLimit to compute the corresponding directional limit (for
#' finite values given by d). The value can still be infinite.
#'
#' @param f a function, to be make stable against problems with limits
#' @param x0 the direction from which limits shall be taken.
#'
#' @returns a numerically stabilised version of `f`: a function that
#'     evaluates `f` but replaces any non-finite result (Inf, -Inf,
#'     NaN) by the corresponding directional limit obtained from
#'     [computeLimit] (which may itself still be infinite).
#'
#'
makeLimitStable = function(f,x0=0.01) {
    function(x) {
        y = f(x)
        bad = !is.finite(y)
        if( any(bad) )
            y[bad] = sapply(x[bad],function(x) computeLimit(f=f,d=sign(x-x0)))
        y
    }
}

## evaluates the function y vectorized at x, but if x in X it will replace
## it with a limit.
deXEval = function(y,x,X) {
    yC=y(x)
    if( length(yC) == length(x) ) {
        # function is real
        nC= x%in%X | !is.finite(yC)
        if( any(nC) ) {
            d=c(-1,rep(1,length(x)-1))
            if( sum(nC)>1 )
                yC[nC] = mapply(computeLimit,x=x[nC],d=d[nC],MoreArgs=list(f=y))
            else
                yC[nC] = computeLimit(x=x[nC],d=d[nC],f=y)

        }
        yC
    } else {
        # function is itself returning multiple values
        nC= x%in%X | !is.finite(yC)
        if( any(nC) ) {
            d=c(-1,rep(1,length(x)-1))
            for( i in 1:nrow(yC) ) {
                nCi=nC[i,]
                yi =function(x) {if(length(x)==1) y(x)[i] else y(x)[i,]}
                if( sum(nCi)>1 )
                    yC[i,][nCi] = mapply(computeLimit,x=x[nCi],d=d[nCi],MoreArgs=list(f=yi))
                else
                    yC[i,][nCi] = computeLimit(x=x[nCi],d=d[nCi],f=yi)
            }
        }
        yC
    }
}

#' @describeIn DiscretizedFunction method for functions that need to
#'     be computed first.
#' @export
#' @method DiscretizedFunction function
DiscretizedFunction.function = function(y,a=0,b=a,n=501,X=c(),transferCdf=NULL,
                                        transferQf=cdf2qf(transferCdf,...),...) {
    if( length(X)>0 ) {
        yX=y(X)
        if( a==b ) {
            return(function(x) yX[match(x,X)])
        }
    } else {
        yX=c()
    }
    if( is.null(transferCdf) ) {
        x = seq(a,b,length.out=n)
        yC= deXEval(y,x,X)
        return(DiscretizedFunction(c(yX,yC),a=a,b=b,X=X))
    } else {
        cdfS = makeLimitStable(transferCdf)
        qfS  = makeLimitStable(transferQf)
        yS   = makeLimitStable(yS)
        yC = deXEval(yS,qfS(seq(cdfS(a),cdfS(b),length.out=n)),X)
        return(DiscretizedFunction(c(yX,yC),a=a,b=b,X=X,transferCdf=transferCdf))
    }
}

#' @describeIn DiscretizedFunction method for Discretizer
#' @export
#' @method DiscretizedFunction Discretizer

DiscretizedFunction.Discretizer = function(y,Y,...) {
    DiscretizedFunction(Y,support=y$support,transferCdf=y$transferCdf)
}

#' @describeIn DiscretizedFunction method for functions on Transformed
#'     Reference measures.
#' @export
#' @method DiscretizedFunction ReferenceMeasureTransformed
DiscretizedFunction.ReferenceMeasureTransformed = function(y,
                                                           f,
                                                           n=513,
                                                           support=getSupport(y),
                                                           ...) {
    DiscretizedFunction(f,
                        a=support$a,
                        b=support$b,
                        X=support$X,
                        transferCdf=y$cdf,
                        transferQf=y$qf,...)
}

#' @describeIn DiscretizedFunction method for functions on Transformed
#'     Reference measures.
#' @export
#' @method DiscretizedFunction ReferenceMeasureUniform
DiscretizedFunction.ReferenceMeasureUniform = function(y,
                                                           f,
                                                           n=513,
                                                           support=getSupport(y),
                                                           ...) {
    DiscretizedFunction(f,
                        a=support$a,
                        b=support$b,
                        X=support$X,...)
}


#' Operator to act on a list
#'
#' The operator acts on all elements of a list and
#' returns a list of the results. The right hand side
#' can be a name in which case the '$' operator is used to
#' e
#'
#' @param X the list on which to operate
#' @param expr an expression, which is either a single name, in which
#'     case a the correponding '$'-operator is applied to the list, or a
#'     full expression, in which '.' takes the role of the selected element.
#' @returns a list, with results of the operations on the individual
#'     list elements.
#' @export
#' @examples
#' data(DTD)
#' DTD$data %.% shape
#' DTD$data %.% summary(.$shape)
#' DTD$data %..% summary(.$shape)
"%.%" <- function(X,expr) {
    listDataElementAccess(X,substitute(expr),parent.frame())
}
#' @describeIn grapes-.-grapes same as '%.%', but tries to simplify the result to an array
#' @export
"%..%" <- function(X,expr) simplify2array(listDataElementAccess(X,substitute(expr),parent.frame()))


#' Accesses list elements with an expression
#'
#' This is the generic behind '%.%'. You need to overwrite this
#' to modify the behaviour of '%.%'.
#'
#' @param X the list on which to operate
#' @param .body an expression, which is either a single name, in which
#'     case a the correponding '$'-operator is applied to the list, or
#'     a full expression, in which '.' takes the role of the selected
#'     element.
#' @param .env the environment the .body expression shall be evaluated in
#' @returns a list from the results of body evaluated in .env, with . standing in for the elements of X, or if .body is a name with .$name evaluted.
#' @export
#'
#' @examples
#' L = list(list(name = "a", val = 1), list(name = "b", val = 2))
#' # a bare name applies the corresponding `$` accessor to each element
#' listDataElementAccess(L, quote(name), environment())
#' # a full expression uses `.` for each list element
#' listDataElementAccess(L, quote(.$val * 10), environment())
#'
listDataElementAccess<- function(X,.body,.env) UseMethod("listDataElementAccess",X)
# ' @describeIn listDataElementAccess default method implementing the behaviour for lists, vectors, BDataObjects,...
#' @export
#' @method listDataElementAccess default
listDataElementAccess.default <- function(X,.body,.env) {
    .f <- function(.) .$name
    environment(.f)=.env
    if( is.name(.body) ) {
        body(.f)[[3]] = .body
    } else {
        body(.f) = .body
    }
    lapply(X,.f)
}




boundedRange = function(...) {
    X = unlist(c(...),use.names=FALSE)
    r = range(X[is.finite(X)])
    if( all(is.finite(r)) )
        return(r)
    else
        return(c(0,10))
}

## This gets to datasets, or single vectors
## A single vector is added to each row of the dataset
## The vector might be without frame or a row
addDatasets = function(x,y) {
    if( is.matrix(x) ) {
        if(is.matrix(y)) {
            if( ncol(x)!=ncol(y)) {
                stop("different number of columns in addDataset:",ncol(x),"!=",ncol(y))
            }
            if( nrow(x)==1 ) {
                if( nrow(y) == 1 ) {
                    return(x+y)
                } else {
                    return(t(c(x)+t(y)))
                }
            } else if( nrow(y)==1 ) {
                    return(t(c(y)+t(x)))
            } else {
                if( nrow(x)!=nrow(y) ) {
                    stop("different number of rows in addDatasets:",nrow(x),"!=",nrow(y))
                }
                return(x+y)
            }
        } else {
            if( ncol(x)!=length(y) ) {
                stop("different number of columns do not match vector legnth:",ncol(x),"!=",length(y))
            }
            return(t(t(x)+c(y)))
        }
    } else {
        if(is.matrix(y)) {
            if( length(x)!=ncol(y) ) {
                stop("vectorlength does not match ncol in addDatasets:",length(x),"!=",ncol(y))
            }
            return( t(t(y)+c(x)))
        } else {
            if( length(x)!=length(y) ) {
                stop("Different lengths in addDatasets: ", length(x),"!=",length(y))
            }
            return(matrix(c(x)+c(y),nrow=1))
        }
    }
}


## computeFor
##
## @param test a logical vector of the length of the result
## @param values computed values for the positions of the test
## @param Default the
##
## @returns returns a vector of length test. The untrue values
##      are filled with default, the true are filled with values in the sequence
##      of the results
computedFor = function(test,values,Default=0) {
    test <- as.logical(test)
    ans <- rep(Default,length.out=length(test))
    len <- length(ans)
    ypos <- which(test)
    if (length(ypos) > 0L)
        ans[ypos] <- values
    ans
}


#' Makes a vectorized 2-ary function into one that returns a matrix
#'
#' A vectorized function of two arguments is made returning the result
#' for each pair structured as a matrix, or if it already returns an
#' n-ary array, as an n+1-ary array with the last dimension splitted
#'
#' @param f the 2-ary working on vectors, lists or rows of matrices or
#'     datasets
#' @returns a function like f but returning the result for each pair of
#'     input structured with one more dimension.
#' @export
#' @examples
#' makeFunctionOuter(\(x,y) exp(-abs(x-y)))(1:4,1:4)
#'
makeFunctionOuter <- function(f) {
    f=f
    function(x,y) {
        nDimX=length(dim(x))
        nDimY=length(dim(y))
        stopifnot(nDimX==nDimY)
        nx = if(nDimX==2) nrow(x) else length(x)
        ny = if(nDimY==2) nrow(y) else length(y)
        iX = rep(1:nx,ny)
        iY = rep(1:ny,each=nx)
        if( nDimX==2 ) {
            Z = f(x[,iX,drop=FALSE],y[,iY,drop=FALSE])
        } else {
            Z = f(x[iX],y[iY])
        }
        if( length(dim(Z)) < 2 )
            return(matrix(Z,nrow=nx,ncol=ny))
        else
            return(structure(Z,dim=c(dim(Z)[-length(dim(Z))],nx,ny)))
    }
}


#' Computes the discretization of a mixed support
#'
#' This provides discretization points for functions with a individual
#' jumps for the purpose of graphical display, i.e. with equispaced
#' points in lim
#'
#' @param support the [Support] to be discretized (or an object which
#'     has a support)
#' @param nEval the number of gridpoints to discretize the continuous
#'     part (defaults to 513)
#' @param ... further generic arguments
#'
#' @returns a list with the elements
#'  \item{x}{the continuous discretization values}
#'  \item{X}{the discrete values}
#'  \item{Xplot}{the location to plot the discrete values}
#'  \item{Xlab}{labels to plot the discrete values}
#'  \item{all}{X and x combined}
#'  \item{iX}{the indices to access X in all}
#'  \item{ix}{the indices to access x in all}
#'
#' @export
#'
#' @examples
#' supportDiscretization(Support(1,3,X=seq(1,2,by=0.1)))
#'
supportDiscretization = function(support,...) UseMethod("supportDiscretization",getSupport(support))
#' @describeIn supportDiscretization method for real supports
#' @export
supportDiscretization.RealSupport = function(support,nEval=513,...) {
    support=getSupport(support)
    X      = getSupportAtoms(support)
    if( support$a==support$b ) {
        xs=X[0]
    } else {
        x     = seq(support$lim[1],support$lim[2],length.out=nEval)
        if( length(X) > 0 )
            x     = x[! (x %in% X) & x>=support$a & x<= support$b]
    }
    all=c(X,x)
    structure(c(support,list(x=x,all=all,Xplot=X,Xlab=X,
         iX=if(length(X)>0) 1:length(X) else c(),
         ix=if(length(x)>0) (length(X)+1):length(all) else c(),
         nEval=nEval
         )),class=c("SupportDiscretization",class(support)))
}

#' @describeIn supportDiscretization method for real supports
#' @method supportDiscretization DiscreteSupport
#' @export
supportDiscretization.DiscreteSupport = function(support,...) {
    support=getSupport(support)
    X      = getSupportAtoms(support)
    if( is.matrix(X) )
        Xlab=paste("X",1:nrow(X),sep="")
    else
        Xlab=X
    structure(c(support,list(x=numeric(0),all=X,Xplot=1:length(X),Xlab=Xlab,
         iX=1:length(X),
         ix=c()
         )),class=c("SupportDiscretization",class(support)))
}




## #' #' SmoothSplineFromHist
## #' #'
## #' #' @param x vector of x-values to be evaluated over the defined B-spline basis
## #' #' @param f vector of clr functional values/evaluations of 'x' over the defined B-spline basis
## #' #' @param knots vector of knots (unique, increasing) defining the B-spline basis
## #' #' @param k degree of the polynomial (order = k+1)
## #' #' @param l degree of a smooth derivation
## #' #' @param resolution argument for the the grid for the numerical integration
## #' #' @param method method used for smoothing: (1) penalization ; (2) convex combination
## #' #' @param alpha mixture parameter
## #' #' @param w vector of weights
## #' #' @param x_eval vector for evaluation of the final smoothed spline
## #' #'
## #' #' @returns  a list with zb-spline coefficients, evaluation of the spline in x_eval and
## #' #'           value of the J functional (see Machalova, Talska, Hron, Gaba:
## #' # Compositional splines for representation of density functions
## #' # Computational Statistics (2021))
## #' #' @export
## #' #'
## #' #' @examples
## #' #' # x=c(0,0.1,0.2,0.35,0.65,0.7,0.8,0.95,0.99)
## #' # f=c(-1,0,0.2,1,1,-0.5,-0.2,-1.5,0)
## #' # knots=seq(0,1,length=6)
## #' # k=3
## #' # l=2
## #' # resolution=1000
## #' # method=2
## #' # alpha=0.99
## #' # w=rep(1,length(x))
## #' # x_eval=seq(0,1,length=100)
## #' # ZB_coefs=SmoothSplineFromHist(x,f,knots,k,l,resolution,method,alpha,w,x_eval=x)
## #' #
## #' # plot(x,f,pch=4,lwd=2,col=2)
## #' # points(x_eval,ZB_coefs$zb_spline,ylim = range(c(f,ZB_coefs$zb_spline)),
## #' #      xlab="x",ylab="f(x)",main="Smoothing spline")
## #' # abline(h=0,col="darkgrey")
## #' # abline(v=knots,col="lightgrey",lty="dashed")
## #'
## #' SmoothSplineFromHist=function(x,f,knots,k,l,resolution,method,alpha,w,x_eval){
## #'   if (!requireNamespace("splines", quietly = TRUE)) {
## #'     stop("Package 'splines' is required but not availabel.")
## #'   }
## #'
## #'   #Bbasis - to be replaced from existing functions
## #'   Bbasis=function(x,knots,k,...){
## #'
## #'   coincident_knots=c(rep(head(knots,1),k),knots,rep(tail(knots,1),k))
## #'   splines::splineDesign(knots=coincident_knots,
## #'                         x=x,
## #'                         ord=k+1,
## #'                         outer.ok = TRUE)
## #' }
## #'
## #'   #L20basis
## #'   L20basis=function(x,knots,k,...){
## #'     g=length(knots)-2
## #'     coincident_knots=c(rep(head(knots,1),k),knots,rep(tail(knots,1),k))
## #'
## #'     D=(k+1)*diag(1/(tail(coincident_knots,g+k+1)-head(coincident_knots,g+k+1)))
## #'
## #'     {matK=diag(g+k)
## #'       matK[cbind(2:(g+k), 1:((g+k)-1))] <- -1
## #'       K=rbind(matK,c(rep(0,g+k-1),-1))}
## #'
## #'     Bbasis(x,knots,k)%*%D%*%K
## #'   }
## #'
## #'   #smoothing
## #'   inner_knots=knots[-c(1,length(knots))]
## #'   g=length(inner_knots)
## #'   coincident_knots=c(rep(head(knots,1),k),knots,rep(tail(knots,1),k))
## #'   N = g+k+1
## #'
## #'   # Matrix S
## #'   S = array(0)
## #'   if (l == 0){
## #'     S = diag(1, c(N,N))
## #'   }
## #'   if (l > 0){
## #'     j=l
## #'     while (j>0){
## #'       Dj = array(0)
## #'       Dj=(k+1-j)*diag(1/(coincident_knots[(k+1)+(-k+j+k+1-j):(g+k+1-j)]-coincident_knots[(k+1)+(-k+j):(g)]))
## #'
## #'       Lj = array(0, c(N-j,N-j+1))
## #'       for (jj in (1:(N-j))){
## #'         Lj[jj,jj] = (-1)
## #'         Lj[jj,jj+1] = 1
## #'       }
## #'       if (j==l){
## #'         S = Dj%*%Lj
## #'       } else{
## #'         S = S%*%Dj%*%Lj
## #'       }
## #'       j=j-1
## #'     }
## #'   }
## #'
## #'   # Matrix M
## #'   augment_knots=coincident_knots[(1+l):(g+2*k+2-l)]
## #'   grid = seq(min(augment_knots), max(augment_knots),  length = resolution)
## #'
## #'   M_dim = g+abs(-k+l)+1
## #'   CC = splineDesign(augment_knots, grid, (k+1-l), outer.ok=TRUE)
## #'   step=diff(grid[1:2])
## #'
## #'   M=array(0, c(M_dim,M_dim))
## #'   for (i in 1:M_dim){
## #'     for (j in 1:M_dim){
## #'       M[i,j]=trapzc_w(step, CC[,i]*CC[,j])
## #'     }
## #'   }
## #'
## #'   # Matrix D
## #'   D=(k+1)*diag(1/(tail(coincident_knots,g+k+1)-head(coincident_knots,g+k+1)))
## #'
## #'   # Matrix K
## #'   {matK=diag(g+k)
## #'     matK[cbind(2:(g+k), 1:((g+k)-1))] <- -1
## #'     K=rbind(matK,c(rep(0,g+k-1),-1))
## #'   }
## #'
## #'   #Martix U
## #'   U = S%*%D%*%K
## #'
## #'   # test
## #'   if (method==2&&(alpha>=1|alpha<=0))stop ('For the selected method, alpha must be from (0,1).')
## #'   if (length(x) <= N) stop ('Length of t must be higher then dimension of the space of splines.')
## #'   if (qr(Bbasis(x,knots,k))$rank != N) stop ('Collocaton matrix B does not have full column rank.')
## #'   if(any(w<=0)) stop ('Weights must be positive.')
## #'
## #'   W = diag(w)
## #'   Bbasis_eval=Bbasis(x,knots,k)
## #'
## #'   if (method==1){G = alpha*t(U)%*%M%*%U+ t(K)%*%t(D)%*%t(Bbasis_eval)%*%W%*%Bbasis_eval%*%D%*%K
## #'   g_bold = t(K)%*%t(D)%*%t(Bbasis_eval)%*%W%*%f}
## #'
## #'   if (method==2){G = t(U)%*%((1-alpha)*M)%*%U + alpha * t(K)%*%t(D)%*%t(Bbasis_eval)%*%W%*%Bbasis_eval%*%D%*%K
## #'   g_bold = alpha*t(K)%*%t(D)%*%t(Bbasis_eval)%*%W%*%f}
## #'
## #'
## #'   # z - vector of z coefficients
## #'   z = solve(G)%*%g_bold
## #'
## #'   # functional J
## #'   if (method==1) {J =     alpha*t(z)%*%t(U)%*%M%*%U%*%z +       t(f-Bbasis_eval%*%D%*%K%*%z)%*%W%*%(f-Bbasis_eval%*%D%*%K%*%z)}
## #'   if (method==2) {J = (1-alpha)*t(z)%*%t(U)%*%M%*%U%*%z + alpha*t(f-Bbasis_eval%*%D%*%K%*%z)%*%W%*%(f-Bbasis_eval%*%D%*%K%*%z)}
## #'
## #'   #
## #'   zb_spline=L20basis(x_eval,knots,k)%*%z
## #'
## #'   return(list(zb_coefficients=z,zb_spline=zb_spline,J_functional=J))
## #' }

matrixSqrtInv = function(S) {
    if( length(S)==0 ) {
        return(S)
    } else if(length(S)==1) {
        return(1/sqrt(S))
    } else {
        SVD =svd(S)
        with(SVD,v%*%diag(ifelse(d>01E-16*min(d[d>0]),1/sqrt(d),0))%*%t(u))
    }
}


#' Accessing cases regardless how a dataset is stored.
#'
#' @param X the dataset, may be a vector a matrix, a dataframe, or
#'          if implmenented something more advanced.
#' @param I the Index vector indexing the cases
#'
#' @return returns the subset of cases.
#' @export
#' @examples
#' nCases(1:5)
#' nCases
#'
nCases = function(X) UseMethod("nCases",X)
#' @describeIn nCases default method working for most know types of datasets, by using the first array dimension (if applicable) or the length as number of cases.
#' @export
#' @method  nCases default
nCases.default = function(X) {
    D = dim(X)
    if(length(D)>1) D[1] else length(X)
}
#' @describeIn nCases selects the cases given by the index I
#' @export
selectCases = function(X,I) UseMethod("selectCases",X)
#' @describeIn nCases selects the cases given by the index I
#' @export
#' @method  selectCases default
selectCases.default = function(X,I) {
    D = dim(X)
    if(length(D)>1) X[I,,drop=FALSE] else X[I]
}