R/classes.R

In DistributionalDDA: Distributional Data Analysis in Bayes Spaces

#' Getting a standard sample from a ReferenceMeasure or distribution
#'
#' This function is here for fast and reproducable Monte Carlo integration,
#' with respect to a measure. Unlike a [simulate] it is allowed to
#' return the same sample in different calls and to provide default
#' size rather than a requested one.
#'
#' @param x the ReferenceMeasure or distribution
#' @param ... further parameters to simulate, which might or might not
#'     be called.
#'
#' @returns a sample of the distribution in the form that would be
#'     expected from a dataset, like e.g. a vector or a data.matrix.
#'
#' @family ReferenceMeasures
#' @export
#'
#' @examples
#' hist(getSample(StdReferenceMeasureUnif01))
#' mean(cos(getSample(StdReferenceMeasureN01)))
#'
getSample = function(x,...) UseMethod("getSample")
#' @describeIn getSample method for discrete measures, returning a weighted sample
#' @export
#' @method getSample ReferenceMeasureDiscrete
getSample.ReferenceMeasureDiscrete = function(x,...) {
    structure(x$X,weights=x$w)
}
#' @describeIn getSample generic method checking whether its weighted sample
#' @export
hasWeights = function(x) UseMethod("hasWeights",x)
#' @describeIn getSample default method checking for a weights attribute.
#' @export
#' @method hasWeights default
hasWeights.default = function(x,...) !is.null(attr(x,"weights"))
#' @describeIn getSample generic method for retrieving weights of a sample
#' @export
getWeights = function(x,...) UseMethod("getWeights",x)
#' @export
#' @method getWeights default
getWeights.default = function(x,...) {
    if(hasWeights(x))
        attr(x,"weights")
    else {
        n=nCases(x)
        rep(1/n,n)
    }
}



# ---- BayesSpace ----
#' Creates the Bayes space for a reference measure
#'
#' A Bayes Space is a Hilbert space of distribution and likelihoods having
#' a non zero density with respect to the reference measure (which needs) to
#' be a distribution itself.
#'
#' @param referenceMeasure a "ReferenceMeasure" object
#' @seealso [ReferenceMeasure]
#' @family BayesSpace
#' @returns An object of class `"BayesSpace"`: a list holding the reference measure.
#' @export
#' @examples
#' bs = BayesSpace(ReferenceMeasureUniform(0,10))
#' bs
#' X = BDataSamples(replicate(10,runif(100,0,10)),basis=BayesSpaceBasis(bs))
#'
BayesSpace <- function(referenceMeasure = ReferenceMeasureUniform(0,1)) {
  structure(
    list(referenceMeasure = referenceMeasure),
    class = "BayesSpace"
  )
}
#' Getter for the ReferenceMeasure of an object
#'
#' BayesSpaces and L2Spaces are defined on underlying reference
#' measures. This function returns this measure.
#'
#' @param x An object that is somehow related to a reference measure,
#'     like e.g. a ReferenceMeasure, a BayesSpace, a BayesSpaceBasis,
#'     or a BData object.
#'
#' @returns The reference-measure object (class `"ReferenceMeasure"`) underlying `x`.
#' @export
#' @seealso [ReferenceMeasure], [BayesSpace], [L2Space]
#' @examples
#' StdReferenceMeasureUnif01
#' getReferenceMeasure(B2Unif01)
#' getReferenceMeasure(B2Unif01Basis5)
#' getReferenceMeasure(L2Space(StdReferenceMeasureUnif01))
#' getReferenceMeasure(BDataSamples(replicate(2,runif(10)),basis=B2Unif01Basis5))
#'
getReferenceMeasure <- function(x) UseMethod("getReferenceMeasure",x)

#' @describeIn getReferenceMeasure method for [BayesSpace]
#' @export
getReferenceMeasure.BayesSpace <- function(x) x$referenceMeasure

#' @describeIn getReferenceMeasure method for [L2Space]
#' @export
getReferenceMeasure.L2Space <- function(x) x$referenceMeasure

#' @describeIn getReferenceMeasure method for ReferenceMeasure, returns the measusre itself
#' @export
getReferenceMeasure.ReferenceMeasure <- function(x) x

#' @describeIn getReferenceMeasure default method tries to get a
#' BayesSpace and the corresponding reference measure from that.
#'
#' @export
#'
getReferenceMeasure.default = function(x) getReferenceMeasure(getBayesSpace(x))

#' Pretty printing of BayesSpaces
#'
#' @param x the Bayes space
#' @param ... further unused parameters for generic use
#' @returns `x`, invisibly (called for the printed side effect).
#' @export
#'
#' @examples
#' print(B2Unif01)
print.BayesSpace <- function(x, ...) {
  cat("BayesSpace with reference measure:")
  print(x$referenceMeasure)
  invisible(x)
  }






# ---- L2Space ----
#' Creates the L2-space for a Reference Measure
#'
#' An L2-space is a Hilbert space of with respect to the reference measure
#' square integrable functions. They are vector addition and multiplication
#' with a scalar is done on the results. The scalar product is defined
#' integrating the product of two functions over the reference measure.
#'
#' @param x the ReferenceMeasure of the L2 space
#' @param ... further arguments for generic use
#'
#' @seealso [ReferenceMeasure]
#' @family BayesSpace
#' @returns An object of class `"L2Space"`: a list holding the reference measure.
#' @export
#' @examples
#' space = L2Space(ReferenceMeasureUniform(0,10))
#' space
#' bs = BayesSpace(ReferenceMeasureUniform(0,10))
#' space2 = L2Space(bs)
#' space2
#'
L2Space = function(x=StdReferenceMeasureUnif01,...) UseMethod("L2Space",x)
#' @describeIn L2Space Method for ReferenceMeasure
#' @export
L2Space.ReferenceMeasure = function(x,...){
  structure(
    .Data=list(referenceMeasure = x),
    class="L2Space"
  )
}


#' @describeIn L2Space Method for BayesSpaces.
#' @export
L2Space.BayesSpace = function(x,...) {
  structure(
    list(referenceMeasure = x$referenceMeasure),
    class="L2Space"
  )
}



#' prints "L2Space with reference measure" and than the reference measure.
#'
#' @param x the L2 space
#' @param ... further unused parameters for generic use
#' @returns `x`, invisibly (called for the printed side effect).
#' @export
#'
#' @examples
#' print(getL2Space(B2Unif01))
print.L2Space <- function(x, ...) {
  cat("L2Space with reference measure:")
  print(x$referenceMeasure)
  invisible(x)
}


# ---- BayesSpaceBasis ----
#' Create a basis for a Bayes space
#'
#' Creates the finite-dimensional basis in which distributional data are
#' represented for computation. Besides the functional spline basis and its
#' [BayesSpace], the returned object stores the transformation matrix `A` and the
#' Gram matrix `S` (the inner products of the basis elements).
#'
#' @details
#' The default basis (constructed from a reference measure or an `fda` B-spline
#' basis) consists of **zero-integral B-splines (ZB-splines)**: spline functions
#' that integrate to zero, which is required to represent centred-log-ratio (clr)
#' transformed densities. The object stores
#' \describe{
#'   \item{`basis`}{a function evaluating the underlying (B-spline) basis,}
#'   \item{`A`}{the transformation matrix mapping those values to the ZB-spline
#'         basis, so the basis functions are `zeta = A \%*\% B`,}
#'   \item{`S`}{the Gram (metric) matrix `S[i,j] = <zeta_i, zeta_j>`,}
#'   \item{`bayesSpace`, `centered`}{the space and a centring flag.}
#' }
#' ZB-splines are **not orthonormal**, so in general `S` is not the identity;
#' inner products, norms, variances and PCA all use `S` (see [getMetricTensor],
#' [princomp.BData]). Pass an explicit `S` only if you have a better value than the
#' internally computed one.
#'
#' @param basis the object describing the basis (a function, an `fda` basis, a
#'     reference measure, or a [BayesSpace]).
#' @param bayesSpace the BayesSpace this basis is designed for
#' @param A the transformation matrix mapping the result of the function to the basis values
#' @param S the Gram (metric) matrix of the given basis, `S[i,j] = <psi_i, psi_j>`. For an
#'     orthonormal basis `S` is the identity; for ZB-splines it is not. Used by all
#'     inner-product based operations.
#' @param ... further generic arguments
#' @param orthonormalize boolean if true creates an orthonormal basis, by a sigualr value decomposition based orthonormalization.

#' @returns An object of class `"BayesSpaceBasis"`: a list with components `basis`,
#'     `bayesSpace`, `A`, `S` and `centered`.
#' @seealso [getMetricTensor], [BayesSpace], [as.BDataSplines], [princomp.BData]
#' @export
#' @examples
#' basis <- BayesSpaceBasis(B2Unif01, nbasis = 7)
#' getMetricTensor(basis)        # Gram matrix S (not the identity for ZB-splines)
BayesSpaceBasis = function(basis,...) UseMethod("BayesSpaceBasis",basis)



#' Creates a BayesSpaceBasis from a function
#'
#' The function needs a vectorized function mapping the support of the
#' Bayes space to real densities.
#'
#' A and S only need to be provided if there is a better result
#' availabel than from the internal computation of the Bayes space
#'
#' @param basis the function representing the log densities
#'     representing the basis.
#' @param bayesSpace the BayesSpace the basis is defined for
#' @param A the matrix mapping the results of f to the log densities.
#' @param S the Gram (metric) matrix of this basis, `S[i,j] = <psi_i, psi_j>` (see
#'     [getMetricTensor]). Stored as the Gram matrix itself, not its inverse.
#' @param ... further unused generic arguments
#' @param centered a boolean, if true the basis is considered centered
#' @param orthonormalize boolean: if TRUE the basis is orthogonalized, using a singluar using a singular value decomposition based matrix square root.
#' @returns A BayesSpaceBasis
#'
#'
#'
#'
#' @seealso BayesSpaceBasis
#'
#' @export
#' @examples
#' BayesSpaceBasis(function(x) outer(x,1:5,function(a,b) cos(a*b*pi)))
#' BayesSpaceBasis(function(x) outer(x,1:5,function(a,b) cos(a*b*pi)),orthonormalize=TRUE)
BayesSpaceBasis.function <- function(basis,
                                     bayesSpace=B2Unif01,
                                     A=diag(rep(1.0,ncol(basis(getSample(bayesSpace))))),
                                     ## S is the Gram (metric) matrix of zeta = A %*% basis
                                     ## (S[i,j] = <zeta_i, zeta_j>), NOT its inverse: every
                                     ## consumer (scalarProduct, mvar, princomp via chol(S))
                                     ## evaluates c^T S d. Same convention as .basisfd; the
                                     ## inverse Gram belongs only in the L2 projection
                                     ## (.SplineCoefsFromLogDensity). Fixes Issue B8.
                                     S=A%*%
                                       scalarProductOfLogDensities(
                                           basis,
                                           space=bayesSpace)%*%
                                       t(A),orthonormalize=FALSE,...,centered=FALSE) {
    if(  orthonormalize ) { #
        MSI = matrixSqrtInv(S)
        A = MSI %*% A
        S = if(length(S)>1)
                diag(rep(1,nrow(S)))
            else
                matrix(1.0,nrow=nrow(S),ncol=nrow(S))
    }
  structure(
    list(
      basis=basis,
      bayesSpace = bayesSpace,
      A = A,
      S = S,
      centered=centered
    ),
    class = "BayesSpaceBasis",
    call=match.call()
  )
}
#' the dimension of the basis
#'
#' @param basis a BayesSpaceBasis
#'
#' @returns a integer providing the dimension of the basis
#' @export
#' @examples
#' getBayesSpaceBasisDimension(B2Unif01Basis5)
#' getBayesSpaceBasisDimension(B2Unif01Basis7)
#'
getBayesSpaceBasisDimension = function(basis) UseMethod("getBayesSpaceBasisDimension",basis)

#' @describeIn getBayesSpaceBasisDimension the method for a standard BayesSpaceBasis
#' @export
getBayesSpaceBasisDimension.BayesSpaceBasis = function(basis) nrow(basis$S)

#' @describeIn BayesSpaceBasis Construct from an basisfd object from functional data
#' @seealso [fda::basisfd]
#' @importFrom fda eval.basis
#' @export
BayesSpaceBasis.basisfd <- function(basis,
                              bayesSpace=BayesSpace(ReferenceMeasureUniform(basis$rangeval[1],basis$rangeval[2])),
                              A=computeL22BayesSpacesBasisMatrix(basis),
                              S=A %*% fda::inprod(basis, basis) %*% t(A) / diff(basis$rangeval),
                              orthonormalize=FALSE,
                              ...) {

                              ## S is the Gram matrix of the ZB-splines zeta = A %*% B with
                              ## respect to the uniform reference MEASURE on rangeval:
                              ## S[i,j] = <zeta_i, zeta_j>_L2(mu) = int zeta_i zeta_j dlambda
                              ## divided by the range width. The division keeps S consistent
                              ## with all inner products computed through CovarianceOperator/
                              ## expectation (which are means w.r.t. mu, not raw Lebesgue
                              ## integrals) and with BayesSpaceBasis.function; without it,
                              ## coefficients from the projection solve(S, <f,zeta>) come out
                              ## too small by the factor diff(rangeval) on domains other than
                              ## the unit interval. fda::inprod() integrates the piecewise-
                              ## polynomial products exactly by quadrature.
                              ## NOT the inverse: all consumers (scalarProduct, mvar, princomp
                              ## via chol(S)) evaluate <c,d>_S = c^T S d. The inverse Gram
                              ## belongs only in the L2 projection coefs = solve(S, <f,zeta>)
                              ## (Issue 3, .SplineCoefsFromLogDensity in helpers.R).
                              ## Verified numerically in tests/testthat/test_GramMatrix.R.

  structure(
    BayesSpaceBasis(
      basis=function(x) fda::eval.basis(x,basis),
      bayesSpace = bayesSpace,
      A = A,
      S = S,
      orthonormalize=orthonormalize,
      centered=TRUE
    ),
    class = "BayesSpaceBasis",
    call=match.call()
  )
}

#' @describeIn BayesSpaceBasis Construct from an fd object from functional data
#' @seealso [fda::fd]
#' @export
BayesSpaceBasis.fd <- function(basis,...) BayesSpaceBasis(basis$basis,...)

#' @describeIn BayesSpaceBasis Calls BayesSpaceBasis on the Reference measure
#' @export
BayesSpaceBasis.BayesSpace = function(basis,...){
   BayesSpaceBasis(basis$referenceMeasure,bayesSpace=basis,...)
}

## NOTE: as.function.BayesSpaceBasis is defined (and exported) in Bayes_spaces.R
## via basisEval(); the earlier duplicate here referenced a free `A` and was a
## dead, shadowed copy -- removed 2026-06-25.

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

#' Gets number of elements in a Basis
#'
#' Returns the dimension of an object.
#'
#'
#' @param x the object
#' @param ... futher unused arguments for generic use
#' @returns an integer giving the number of dimensions of the object
#' @export
#' @examples
#' getDimension(1:3)
#' getDimension(matrix(rnorm(10),ncol=2))
#' getDimension(B2N01Basis5)
#'
getDimension <- function(x,...) UseMethod("getDimension",x)
#' @describeIn getDimension returns number of elements in the basis
#' @export
#' @method getDimension BayesSpaceBasis
getDimension.BayesSpaceBasis = function(x,...) nrow(x$S)
#' @describeIn getDimension returns number of elements in the basis
#' @export
#' @method getDimension BDataSplines
getDimension.BDataSplines = function(x,...) ncol(x)
#' @describeIn getDimension method for matrix returns number of columns
#' @export
#' @method getDimension matrix
getDimension.matrix = function(x,...) ncol(x)
#' @describeIn getDimension method for numeric returns number of elements
#' @export
#' @method getDimension numeric
getDimension.numeric = function(x,...) length(x)



#' Checks that an object is a BayesSpaceBasis
#'
#' @param X the object to be checked
#'
#' @returns boolean
#' @export
#' @examples
#' isValidBayesSpaceBasis(B2Unif01Basis7)
#'
isValidBayesSpaceBasis   = function(X) is(X,"BayesSpaceBasis")

#' Checks whether an object is an element (or dataset of elements) of
#' a BayesSpace
#'
#' May different R objects can be interpreted as elements of a BayesSpace
#' including density functions, histograms, and BData objects.
#'
#' @param X the element to be checked
#' @param basis an optional BayesSpaceBasis in which it should be an
#'     element. If 'basis' is given a TRUE shall only be returned if
#'     can be seen as an element of that BayesSpace. That however is
#'     not fully implemented yet.
#' @returns boolean, FALSE if it is definitely not an element of a
#'     BayesSpace
#' @export
#'
#' @examples
#' isValidBayesSpaceElement(sin)
#' isValidBayesSpaceElement(hist(runif(100), plot=FALSE))
isValidBayesSpaceElement   = function(X,basis=NULL) UseMethod("isValidBayesSpaceElement",X)
#' @export
#' @describeIn isValidBayesSpaceElement functions are valid BayesSpaceElements
#' @method isValidBayesSpaceElement function
isValidBayesSpaceElement.function = function(X,basis=NULL) TRUE
#' @export
#' @describeIn isValidBayesSpaceElement histograms are valid BayesSpaceElements
#' @method isValidBayesSpaceElement histogram
isValidBayesSpaceElement.histogram = function(X,basis=NULL) TRUE
#' @export
#' @describeIn isValidBayesSpaceElement BData objects contain valid BayesSpaceElements if they are for the same BayesSpace
#' @method isValidBayesSpaceElement BData
isValidBayesSpaceElement.BData = function(X,basis=NULL) is.null(basis) || hasSameBayesSpace(X,basis)
#' @export
#' @describeIn isValidBayesSpaceElement method for samples, tries to
#'     check whether the values are in the domain.
#' @method isValidBayesSpaceElement default
isValidBayesSpaceElement.default = function(X,basis=NULL) all(isInDomain(X,getReferenceMeasure(basis)))




#' BayesSpace data in different formats
#'
#' A BDataGeneral object is essentially a list of different
#' representations of BayesSpace elements of the same BayesSpace. They
#' are treated individually as long as possible.
#'
#' @param data the list of different BayesSpace elements
#' @param basis the basis of the common BayesSpace that shall be used,
#'     when a basis is needed.
#' @param close boolean, if true a closed object will be returned
#' @param closed boolean, if true it is assumed that the input is already closed
#' @param centered boolean, if true it is assumed that the input is already centered
#' @param ... further arguments to clo
#'
#' @returns An object of class `c("BDataGeneral","BData")`: a list bundling the
#'     several distributional representations for a common Bayes space.
#' @export
#'
#' @examples
#' X = BDataGeneral(list(sin,hist(runif(100)),BDataSamples(list(runif(100)),basis=B2Unif01Basis7)),
#'              basis=B2Unif01Basis5)
#' Xs= as.BDataSplines(X)
BDataGeneral = function(data,basis,...,close=FALSE,closed=isClosed(data),centered=isCentered(data)) {
    stopifnot(is.list(data))
    stopifnot(isValidBayesSpaceBasis(basis))
    stopifnot(all(sapply(data,isValidBayesSpaceElement,basis=basis)))
    if( close && ! closed )
        data = lapply(data,clo,...)
    structure(data,class=c("BDataGeneral","BData"),closed=close||closed,centered=centered&&!close,basis=basis)
}



# ' # ---- BSSpline ----
# ' #' create a `BSSpline` object: a single distribution represented as a spline
# ' #'
# ' #' BSSpline is a deprecated object, it will probably be removed soon. Use
# ' #' [BDataSplines] instead.
# ' #'
# ' #' @param coefs a vector of spline coefficients
# ' #' @param basis a spline basis - object of class `basisfd` (fda)
# ' #'
# ' #' @returns object of class `BSSpline`
# ' #' @export
# ' #'
# ' #' @examples
# ' #' coefs=1:5
# ' #' bsspline=BSSpline(coefs,B2Unif01Basis5)
# ' BSSpline <- function(coefs, basis) {
# '   if(length(coefs)!=getBayesSpaceBasisDimension(basis)){
# '     stop("`length(coefs)` must be the same as the number of basis functions.")}
# '   structure(
# '     list(coefs = coefs, basis = basis),
# '     class = "BSSpline"
# '   )
# ' }
# '
# #' @export
# print.BSSpline <- function(x, ...) {
#     cat("BSSpline object with", length(x$coefs), "coefficients.\n")
#     invisible(x)
# }


# ---- BDataSplines ----
#' create a `BDataSplines` object: a set of distributions represented as splines;
#' possible representation of individual distribution through a vector of coefficients.
#'
#' @param matrixCoefs a matrix of spline coefficients (nrow: # of distributions;
#'                                                     ncol: # of basis functions)
#' @param basis a spline basis
#' @param x a BDataSplines object to be subsetted
#' @param i index of samples to be selected/filtered out, or logical vector of length equal to `length(x)`
#' @param j index of variables to be selected/filtered out, or logical vector of length equal to `ncol(x[[1]])`
#' @param ... further generic arguments
#'
#' @returns object of class `BDataSplines` describing a set of distributions
#' represented as splines
#' @export
#'
#' @family BDataSplines
#' @seealso BData
#' @examples
#' coefs=matrix(1:15,nrow=3)
#' bssplines=BDataSplines(coefs,B2Unif01Basis5)
#' bssplines
BDataSplines <- function(matrixCoefs, basis) {
  if(is.vector(matrixCoefs)) {
    matrixCoefs=t(as.matrix(matrixCoefs))}
  basis=as.BayesSpaceBasis(basis)
  if(ncol(matrixCoefs)!=getBayesSpaceBasisDimension(basis)){
    stop("`ncol(matrixCoefs)` (",ncol(matrixCoefs),") must be the same as the number of basis functions (",getBayesSpaceBasisDimension(basis),").")}
  structure(
    matrixCoefs,
    class = c("BDataSplines","BData"),
    basis=basis
  )
}

#' Print a BDataSplines object
#'
#' @param x a [BDataSplines] object
#' @param ... further arguments, currently unused.
#' @returns `x`, invisibly (called for the printed side effect).
#' @export
#' @method print BDataSplines
#' @examples
#' X <- BDataSplines(t(1:5), B2Unif01Basis5)
#' print(X)
print.BDataSplines <- function(x, ...) {
  cat("BDataSplines object with", nrow(x), "observations and",
      ncol(x), "coefficients each.\n")
  print(unclass(x))
  invisible(x)
}




# ==== BDataHistograms ====
# ---- |__ BDataFrequencyHistograms ----
#' Data container for distributions as histograms
#'
#' Simple S3 constructor of a data container for distributional data where each
#' distribution is represented as a (absolute) frequency histogram.
#'
#' @param matrixCounts a matrix (or some object convertible to one with `as.matrix()`)
#' containing absolute frequencies, observations by row and histogram bins by column
#' @param breaks defining the histogram bins
#' @param basis the B-spline basis to be used
#' @return An object of class `c("BDataFrequencyHistograms","BDataHistograms","BData")`, with the
#' matrix counts as main value and `breaks` and `basis` as attributres.
#' @export
#' @family BDataHistograms
#'
#' @examples
#' set.seed(47)
#' breaks = seq(-3,3,by=1)
#' inRange = function(x) -3<x & x<3
#' X = t(replicate(10,hist(Filter(inRange,rnorm(1000)),breaks=breaks)$counts))
#' hists = BDataFrequencyHistograms(X,breaks=breaks,basis=B2N01Basis5)
#' plot(hists)
BDataFrequencyHistograms <- function(matrixCounts, breaks, basis) {
    if( as.integer(ncol(matrixCounts)) != as.integer((length(breaks)-1)) ) {
       stop("`length(breaks)-1`= ",length(breaks)-1, " not compatible with `ncol(matrixCounts)`=",ncol(matrixCounts),"->",(ncol(matrixCounts)!=(length(breaks)-1)))
    }
  structure(
    as.matrix(matrixCounts),
    class = c("BDataFrequencyHistograms","BDataHistograms","BData"),
    breaks=breaks,
    basis = basis
  )
}



# ---- |__ BDataDensityHistograms ----
#' Data container for distributions as histograms
#'
#' Simple S3 constructor of a data container for distributional data where each
#' distribution is represented as a density histogram.
#'
#' @param matrixDensities a matrix (or some object convertible to one with `as.matrix()`)
#' containing densities, observations by row and histogram bins by column
#' @param breaks defining the histogram bins
#' @param basis the B-spline basis to be used
#'
#' @return An object of class `c("BDataDensityHistograms","BDataHistograms","BData")`, with the
#' matrix densities as main value and `breaks` and `basis` as attributes.
#' @export
#' @family BData
#' @examples
#' set.seed(47)
#' breaks = seq(-3,3,by=1)
#' inRange = function(x) -3<x & x<3
#' X = t(replicate(10,hist(Filter(inRange,rnorm(1000)),breaks=breaks,plot=FALSE)$density))
#' hists = BDataDensityHistograms(X,breaks=breaks,basis=B2N01Basis5)
#' plot(hists)
BDataDensityHistograms <- function(matrixDensities, breaks, basis) {
  if( as.integer(ncol(matrixDensities)) != as.integer((length(breaks)-1)) ) {
    stop("`length(breaks)-1`= ",length(breaks)-1, " not compatible with `ncol(matrixDensities)`=",ncol(matrixDensities),"->",(ncol(matrixDensities)!=(length(breaks)-1)))
  }
  structure(
    simpleClo(as.matrix(matrixDensities)) %*% diag(1/diff(breaks)),
    class = c("BDataDensityHistograms","BDataHistograms","BData"),
    breaks=breaks,
    basis = basis,
    closed=TRUE
  )
}

# ---- |__ Common utilities for BDataHistograms objects ----
# Suggest to remove, because the default printing is already very clear
# ' Common utilities for BDataHistograms objects
# '
# ' @param x object with a class derived from "BDataHistograms"
# ' @param ... further arguments
# '
# ' @return depeding on the function, this can print the content, extract the breaks
# ' attribute in a convenient fashion or do subsetting.
# ' @export
# ' @family BDataHistograms
# '
# '
# '
# print.BDataHistograms <- function(x, ...) {
#  cat("BDataHistograms with", nrow(x), "observations and",
#      length(getBreaks(x))-1, "bins:\n")
#  print(as.data.frame(unclass(x)))
#  cat("using breaks:\n")
#  print(getBreaks(x))
#  cat(paste("of specific class", class(x)[1]))
#  invisible(x)
#}



#' Subsetting BDataHistograms
#'
#' A subsetting function selecting only a part of the histograms
#' @param x a BDataHistograms object
#' @param i index of rows to be extracted
#' @param j index of the columns, typically missing for the subsetting
#' @param ... further generic parameters such das drop=FALSE
#'
#' @returns A BDataHistograms object with the subset of histograms described by the basis
#'
#' @export
#' @examples
#' set.seed(47)
#' breaks = seq(-3,3,by=1)
#' inRange = function(x) -3<x & x<3
#' X = t(replicate(4,hist(Filter(inRange,rnorm(1000)),breaks=breaks,plot=FALSE)$counts))
#' hists = BDataFrequencyHistograms(X,breaks=breaks,basis=B2N01Basis5)
#' hists[1:2, ]
"[.BDataHistograms" <- function(x, i,j,...) {
  # Subsetting columns is not allowed because it
  # would destroy the meaning of the histogram breaks.
  # Do it manually yourself if you
                                        # know what you are doing.
  if( missing(j) ) {
      breaks = getBreaks(x)
      basis=getBayesSpaceBasis(x)
      structure(unclass(x)[i=i, , drop=FALSE],
                class=class(x),basis=basis, breaks = breaks,closed=attr(x,"closed"))
  } else NextMethod()
}


#' Common utilities for BDataHistograms objects
#'
#' @param x object with a class derived from "BDataHistograms"
#'
#' @return depeding on the function, this can print the content, extract the breaks
#' attribute in a convenient fashion or do subsetting.
#' @export
#' @family BDataHistograms
#'
#' @export
#' @examples
#' set.seed(47)
#' H = hist(rnorm(30))
#' getBreaks(H)
#' X = as.BDataFrequencyHistograms(list(H,H),basis=B2N01Basis5)
#' getBreaks(X)
#' X
#' X[1]
#'
getBreaks = function(x) UseMethod("getBreaks",x)
#' @describeIn getBreaks Method for BDataHistograms
#' @export
getBreaks.BDataHistograms = function(x) attr(x,"breaks")
#' @export
getBreaks.histogram = function(x) x$breaks


### getBayesSpace
#' Generic function returning the Bayes space associated to the given object.
#'
#' @param X the object
#' @param ... further named arguments used by methods e.g. Default
#' @param Default the Bayes Space returned if no Bayes Space is associated
#'
#' @returns The associated [BayesSpace] object (class `"BayesSpace"`).
#' @export
#' @family getBayesSpace
#' @examples
#' B2Unif01Basis5
#' unclass(B2Unif01Basis5)
#' getBayesSpace(B2Unif01Basis5)
#'
#'
getBayesSpace = function(X,...) UseMethod("getBayesSpace",X)
#'
#' @export
#' @describeIn getBayesSpace method for BayesSpace returning itself
getBayesSpace.BayesSpace = function(X,...) ensureType(X,"BayesSpace")
#' @export
#' @describeIn getBayesSpace method for BayesSpaceBasis
getBayesSpace.BayesSpaceBasis = function(X,...) ensureType(X$bayesSpace,"BayesSpace")
#' @export
#' @describeIn getBayesSpace method for NULL, returns default
getBayesSpace.NULL = function(X,...,Default) ensureType(Default,"BayesSpace")
#' @export
#' @describeIn getBayesSpace method for BData
getBayesSpace.BData = function(X,...) getBayesSpace(getBayesSpaceBasis(X))
#' @export
#' @describeIn getBayesSpace creates the BayesSpace for the reference measure
getBayesSpace.ReferenceMeasure = function(X,...) BayesSpace(X,...)
#' @export
#' @describeIn getBayesSpace default method tries to get an attribute
getBayesSpace.default = function(X,...,Default=attr(X,"BayesSpace")) ensureType(Default,"BayesSpace")
#' @export
#' @describeIn getBayesSpace returns BayesSpace with same reference measure
getBayesSpace.L2Space = function(X,...) BayesSpace(X$referenceMeasure,...)

### getL2Space
#' Generic function returning the L2Space associated to the given object.
#'
#' @param X the object
#' @param ... further named arguments used by methods e.g. Default
#' @param Default the Bayes Space returned if no L2 Space is associated
#'
#' @returns The associated [L2Space] object (class `"L2Space"`).
#' @export
#' @family getBayesSpace
#' @examples
#' getL2Space(B2Unif01Basis5)
#' set.seed(47)
#' X = BDataGeneral(list(function(x) x, function(x) x^2),basis=B2Unif01Basis5)
#' getL2Space(X)
#' getL2Space(StdReferenceMeasureUnif01)
#'
getL2Space = function(X,...) UseMethod("getL2Space",X)
#'
#' @export
#' @describeIn getBayesSpace method for BayesSpace returning itself
getL2Space.L2Space = function(X,...) ensureType(X,"L2Space")
#' @export
#' @describeIn getL2Space method for BayesSpaceBasis
getL2Space.BayesSpaceBasis = function(X,...) getL2Space(X$bayesSpace)
#' @export
#' @describeIn getL2Space method for NULL, returns default
getL2Space.NULL = function(X,...,Default) ensureType(Default,"L2Space")
#' @export
#' @describeIn getL2Space method for BData
getL2Space.BData = function(X,...) getL2Space(getBayesSpace(X))
#' @export
#' @describeIn getL2Space creates the BayesSpace for the reference measure
getL2Space.ReferenceMeasure = function(X,...) getL2Space(BayesSpace(X,...))
#' @export
#' @describeIn getL2Space default method tries to get an attribute
getL2Space.default = function(X,...,Default=getL2Space(attr(X,"BayesSpace"))) ensureType(Default,"L2Space")
#' @export
#' @describeIn getL2Space returns L2Space with same reference measure
getL2Space.BayesSpace = function(X,...) L2Space(X$referenceMeasure,...)


#' Generic function returning the BayesSpaceBasis associated to the given object.
#'
#' @param X the object
#' @param ... further named arguments used by methods e.g. Default
#'
#' @returns The associated [BayesSpaceBasis] object (class `"BayesSpaceBasis"`).
#' @export
#' @seealso [BayesSpaceBasis], [getBayesSpace], [BayesSpace]
#' @examples
#' getBayesSpaceBasis(B2Unif01)
#' set.seed(47)
#' X = BDataGeneral(list(function(x) x, function(x) x^2),basis=B2Unif01Basis5)
#' getBayesSpaceBasis(X)
#' getBayesSpace(X)
#' getReferenceMeasure(X)
#' getBayesSpaceBasis(B2Unif01Basis5)
#'
#'
getBayesSpaceBasis = function(X,...) UseMethod("getBayesSpaceBasis",X)
#' @describeIn getBayesSpaceBasis returns itself
#' @export
getBayesSpaceBasis.BayesSpaceBasis = function(X,...) X
#' @describeIn getBayesSpaceBasis method for Datasets
#' @export
getBayesSpaceBasis.BData = function(X,...) attr(X,"basis")
#' @describeIn getBayesSpaceBasis constructs a default vasis from the reference measure. Use ... as further parameters to generic constructor.
#' @export
getBayesSpaceBasis.BayesSpace = function(X,...) BayesSpaceBasis(X,...)



# ---- BDataDensities and BDataLogDensities ----

#' Constructor of BDataDensities objects representing a Dataset of Densities in a BayesSpace
#'
#' The dataset is represented by a list of functions, associated with
#' a Bayes space. The basis is used for functions relying on
#' projections, e.g. when the data is converted to a BDataSpline object.
#'
#' BDataDensities represent densities (aka as Radon-Nikodym
#' derivatives) with repect to the reference measure of the
#' BayesSpace.
#'
#' BDataPDFs for (Probablity Density function) represents densities
#' with respect to the Lebesgue measure (Haar measure) of the
#' support. I.e. the measure that is invariant under shift or
#' permutation of the space. BDataPDFs are only reasonable if the
#' reference measure is absolutely continuous, i.e. itself has a finite
#' Lebesgue density.
#'
#' For the special case of uniform reference measures BDataDensities
#' and BDataPDFs of the same function, represent the same measure.
#'
#' @param f a dataset of densities in some form
#' @param basis the BayesSpaceBasis (or an object that can be converted to one) to be used for projections.
#'
#' @param ... further generic parameters
#' @param close boolean, if true a closed object will be returned
#' @param closed boolean, if true it is assumed that the input is already closed
#' @param centered boolean, if true it is assumed that the input is already centered
#' @seealso [BayesSpaceBasis]
#' @family BData
#'
#' @returns An object of class `c("BDataDensities","BData")` representing the
#'     densities with respect to the reference measure.
#' @export
#' @examples
#' BDataDensities(list(function(x) x, function(x) 1-x),B2Unif01Basis5)
#'
#'
BDataDensities    = function(f,basis,...) UseMethod("BDataDensities",f)
#' @describeIn BDataDensities Object representing Lebesgue densities
#' @export
BDataPDFs         = function(f,basis,...) UseMethod("BDataPDFs",f)

#' Constructor of BDataLogDensities objects representing a Dataset of log Densities in a BayesSpace
#'
#' The dataset is represented by a list of functions, associated with
#' a Bayes space. The basis is used for functions relying on
#' projections, e.g. when the data is converted to a BDataSpline object
#'
#' @param f a dataset of densities in some form
#' @param basis the BayesSpaceBasis (or an object that can be converted to one) to be used for projections.
#' @param ... further generic parameters
#' @param close boolean, if true a closed object will be returned
#' @param closed boolean, if true it is assumed that the input is already closed
#' @param centered boolean, if true the input is assumed to be centered
#'
#' @seealso [BayesSpaceBasis]
#' @family BData
#'
#' @returns An object of class `c("BDataLogDensities","BData")` representing the
#'     log-densities with respect to the reference measure.
#' @export
#' @examples
#' BDataDensities(list(function(x) x, function(x) 1-x),B2Unif01Basis5)
#'
#'
BDataLogDensities = function(f,basis,...) UseMethod("BDataLogDensities",f)
#' @describeIn BDataLogDensities Object representing Lebesgue densities
#' @export
BDataLogPDFs         = function(f,basis,...) UseMethod("BDataLogPDFs",f)

#' @describeIn BDataDensities Method for list of functions interpreted
#'     as log densities.
#' @export
#' @method BDataLogDensities list
BDataLogDensities.list <- function(f, basis,...,close=FALSE,closed=isClosed(f),centered=isCentered(f)) {
  stopifnot(all(sapply(f,is.function)))
  clo(structure(
    f,
    class = c("BDataLogDensities","BData"),
    basis = basis,
    closed=closed,
    centered=centered
  ),...,close=close)
}
#' @describeIn BDataDensities Method for list of functions interpreted
#'     as log probability density functions (PDFs).
#' @export
#' @method BDataLogPDFs list
BDataLogPDFs.list <- function(f, basis,...,close=FALSE,closed=isClosed(f),centered=isCentered(f)) {
  stopifnot(all(sapply(f,is.function)))
  clo(structure(
    f,
    class = c("BDataLogPDFs","BData"),
    basis = basis,
    closed=closed,centered=centered
  ),...,close=close)
}

#' Print a BDataLogDensities object
#'
#' @param x a [BDataLogDensities] object
#' @param ... passed to the next print method.
#' @returns `x`, invisibly.
#' @export
#' @method print BDataLogDensities
print.BDataLogDensities <- function(x, ...) {
     M = getReferenceMeasure(x)
     cat("BDataLogDensities with", length(x$functions), "log-")
     print(M)
     cat("-density functions.\n")
     NextMethod()
}
#' Print a BDataLogPDFs object
#'
#' @param x a [BDataLogPDFs] object
#' @param ... passed to the next print method.
#' @returns `x`, invisibly.
#' @export
#' @method print BDataLogPDFs
print.BDataLogPDFs <- function(x, ...) {
    cat("BDataLogPDF with", length(x), "log-Lebesgue-density functions.\n")
    NextMethod()
}


#' @describeIn BDataDensities Method for a single function
#' @export
BDataDensities.function = function(f,basis,...,close=FALSE,closed=isClosed(f),centered=isCentered(f)) {
    clo(BDataDensities(list(f),basis,...,closed=closed,centered=centered),close=close)
}
#' @describeIn BDataDensities Method for a single function
#' @export
BDataPDFs.function = function(f,basis,...,close=FALSE,closed=isClosed(f),centered=isCentered(f)) {
    clo(BDataPDFs(list(f),basis,...,closed=closed,centered=centered),close=close)
}

#' @describeIn BDataDensities Method for a list of functions
#' @export
BDataDensities.list = function(f,basis,...,close=FALSE,closed=isClosed(f),centered=isCentered(f)) {
    clo(structure(lapply(f,as.function),
              class=c("BDataDensities","BData"),basis=as.BayesSpaceBasis(basis),
              ...,closed=closed,centered=centered
              ),close=close)
}

#' @describeIn BDataDensities Method for a list of functions
#' @export
BDataPDFs.list = function(f,basis,...,close=FALSE,closed=FALSE,centered=FALSE) {
    clo(structure(lapply(f,as.function),
              class=c("BDataPDFs","BData"),basis=as.BayesSpaceBasis(basis),
              ...,closed=closed,centered=centered)
              ,close=close)
}


#' @describeIn BDataLogDensities Method for a list of functions
#' @export
BDataDensities.BDataLogDensities = function(f,basis,...,close=FALSE,closed=isClosed(f),centered=isCentered(f)) {
    clo(structure(lapply(f,as.function),
              class=c("BDataLogDensities","BData"),
              basis=as.BayesSpaceBasis(basis),
              ...,closed=closed,centered=centered),close=close)
    }


#' @describeIn BDataLogDensities conversion method
#' @export
BDataDensities.BDataLogDensities = function(f,basis=getBayesSpaceBasis(f),...,close=FALSE,closed=isClosed(f),centered=isCentered(f)) {
    clo(structure(lapply(f,Exp),
              class=c("BDataLogDensities","BData"),
              basis=as.BayesSpaceBasis(basis),
              ...,closed=closed),close=close,centered=centered)
    }

#' @describeIn BDataLogDensities conversion method
#' @export
BDataLogDensities.BDataDensities = function(f,basis=getBayesSpaceBasis(f),...,close=FALSE,closed=isClosed(f),centered=isCentered(f)) {
    clo(structure(lapply(f,Ln),
              class=c("BDataLogDensities","BData"),
              basis=basis,
              ...,closed=closed,centered=centered),close=close)
    }
#' @describeIn BDataLogDensities conversion method
#' @export
BDataPDFs.BDataLogPDFs = function(f,basis=getBayesSpaceBasis(f),...,close=FALSE,closed=isClosed(f),centered=isCentered(f)) {
    clo(structure(lapply(f,Exp),
              class=c("BDataPDFs","BData"),
              basis=as.BayesSpaceBasis(basis),
              ...,closed=closed,centered=centered),close=close)
    }

#' @describeIn BDataLogDensities conversion method
#' @export
BDataLogPDFs.BDataPDFs = function(f,basis=getBayesSpaceBasis(f),...,close=FALSE,closed=isClosed(f),centered=isCentered(f)) {
    clo(structure(lapply(f,Ln),
              class=c("BDataLogPDFs","BData"),
              basis=basis,
              ...,closed=FALSE,centered=centered),close=close)
    }


#' @export
getBayesSpaceBasis.BDataLogDensities = function(X,...) attr(X,"basis")
#' @export
getBayesSpaceBasis.BDataDensities =    function(X,...) attr(X,"basis")

#' @export
getBayesSpaceBasis.BDataLogPDFs = function(X,...) attr(X,"basis")
#' @export
getBayesSpaceBasis.BDataPDFs =    function(X,...) attr(X,"basis")

#' Returns whether its know that the object is closed or centered
#'
#' @param x the object to be checked
#' @returns a boolean, if true the system believes that the object is
#'     already closed
#'
#' @export
#'
#' @examples
#' isClosed(dnorm) # FALSE, because we have not way of knowing
#' isClosed(BDataDensities(list(dnorm),basis=B2N01Basis5,closed=TRUE)) # TRUE by assertation
#' isClosed(BDataDensities(list(Mul(2,dnorm)),basis=B2N01Basis5,close=TRUE))  # TRUE by making it
#' isCentered(dnorm) # FALSE, because we have not way of knowing
#' isCentered(BDataDensities(list(dnorm),basis=B2N01Basis5,closed=TRUE)) # FALSE
#' isCentered(cen(BDataDensities(list(Mul(2,dnorm)),basis=B2N01Basis5)))  # TRUE by making it
isClosed = function(x) UseMethod("isClosed",x)
#' @describeIn isClosed default methods always returns FALSE
#' @export
#' @method isClosed default
isClosed.default = function(x) FALSE
#' @describeIn isClosed BData returns true if it was closed or declared closed
#' @export
#' @method isClosed BData
isClosed.BData = function(x) isTRUE(attr(x,"closed"))
#' @describeIn isClosed Checks for centering
#' @export
isCentered = function(x) UseMethod("isCentered",x)
#' @describeIn isClosed default methods always returns FALSE
#' @export
#' @method isCentered default
isCentered.default = function(x) FALSE
#' @describeIn isClosed BData returns true if it was closed or declared closed
#' @export
#' @method isCentered BData
isCentered.BData = function(x) isTRUE(attr(x,"centered"))
#' @describeIn isClosed BDataSplines returns true if the basis is centered
#' @export
#' @method isCentered BDataSplines
isCentered.BDataSplines = function(x) isCentered(attr(x,"basis"))
#' @describeIn isClosed BayeSpaceBasis requests whether the basis is considered centered
#' @export
#' @method isCentered BayesSpaceBasis
isCentered.BayesSpaceBasis = function(x) isTRUE(attr(x,"centered"))



#' @export
"[.BDataLogDensities" <- function(x,...) {
    structure(unclass(x)[...],
             class=c("BDataLogDensities","BData"),basis=getBayesSpaceBasis(x),closed=isClosed(x))
    }

#' @export
"[.BDataLogPDFs" <- function(x,...) {
    structure(unclass(x)[...],
             class=c("BDataLogPDFs","BData"),basis=getBayesSpaceBasis(x),closed=isClosed(x),centered=isCentered(x))
    }



#' @export
"[.BDataDensities" <- function(x,...) {
    structure(unclass(x)[...],
             class=c("BDataDensities","BData"),basis=getBayesSpaceBasis(x),closed=isClosed(x),centered=isCentered(x))
    }
#' @export
"[.BDataPDFs" <- function(x,...) {
    structure(unclass(x)[...],
             class=c("BDataPDFs","BData"),basis=getBayesSpaceBasis(x),closed=isClosed(x),centered=isCentered(x))
    }



# ---- BDataSamples and BDataUnivariateSamples ----


#' BData container with distributions as samples
#'
#' Constructs a BData object where each distribution is represented as a sample.
#' Such a sample can be a vector of values (for univariate cases) or some data
#' structure convertible to a "data.frame" (for multivariate cases).
#'
#' @param datalist list of samples
#' @param basis the B-spline basis to be used
#' @param i index of samples to be selected/filtered out, or logical vector of length equal to `length(x)`
#' @param j index of variables to be selected/filtered out, or logical vector of length equal to `ncol(x[[1]])`
#' @return In multivariate cases, the outcome is an object of class
#' `c("BDataSamples","BData")` encoding a list of data.frames of the same number of
#' columns. In univariate cases, the outcome is an object of class
#' `c("BDataUnivariateSamples","BDataSamples","BData")` encoding a list of
#' vectors of samples of each individual distribution.
#'
#' @export
#' @family BDataSamples
#' @examples
#' x = lapply(rep(1000,4), runif) |> BDataSamples(basis=B2Unif01Basis5)
#' summary(x)
BDataSamples <- function(datalist, basis) {
  # compatibility checks:
  n = ncol(datalist[[1]])
  if(is.null(n)){
    if( sapply(datalist, ncol) |> sapply(is.null) |> all() ){
      thisclass = c("BDataUnivariateSamples","BDataSamples","BData")
    }else{
      stop("dimensions of datalist elements do not fit")
    }
  }else{
    if( all(sapply(datalist, ncol) == n) ){
      thisclass = c("BDataSamples","BData")    # TODO: extend this to a subclass of data.frame?
      datalist = lapply(datalist, as.data.frame)
    }else{
      stop("dimensions of datalist elements do not fit")
    }
  }
  # return object:
  structure(
    datalist,
    class = thisclass,
    basis = basis
  )
}

#' @describeIn BDataSamples summary method for univariate cases
#' @export
#' @param object object to be summarized
#' @param ... further arguments passed to summary.default for the
#' individual summaries of each sample/variable, or to vector subsetting
summary.BDataUnivariateSamples = function(object,...){
  cat("BDataUnivariateSamples with", length(object), "observations:\n")
  for(n in 1:length(object)){
    cat("sample", n, "with", length(unclass(object)[[n]]), "individuals,\n")
    print(summary(unclass(object)[[n]], ...))
  }
}

#' @describeIn BDataSamples summary method for multivariate cases
#' @export
summary.BDataSamples = function(object,...){
  cat("BDataSamples with", length(object), "observations:\n")
  for(n in 1:length(object)){
    cat("sample", n, "with", nrow(unclass(object)[n]), "individuals,\n")
    print(sapply(unclass(object)[n], summary, ...))
  }
}


#' @describeIn BDataSamples subsetting method for univariate samples
#' @export
#' @param x object to be subsetted
"[.BDataUnivariateSamples" <- function(x, ...) {
  structure(unclass(x)[...], class=class(x), basis=getBayesSpaceBasis(x))
}

#' @describeIn BDataSamples subsetting method for multivariate samples
#' @export
"[.BDataSamples" <- function(x, i,j) {
                                        # case no j
  basis=getBayesSpaceBasis(x)
  if(missing(j)){
    return(structure(unclass(x)[i],
                     class=class(x), basis=basis))
  }
  j = (1:ncol(x[[1]]))[j]
  # define output class
  if(length(j)>1){
    outclass = c("BDataSamples","BData")
  }else{
    outclass = c("BDataUnivariateSamples","BDataSamples","BData")
  }
  out = lapply(unclass(x)[i], function(xx)xx[,j])
  return(structure(out, class=outclass, basis=basis))
}

#' @describeIn BDataSplines subsetting method, keeps structure if only i is given.
#' @export
#'
"[.BDataSplines" <- function(x,i,j,...) {
    basis=getBayesSpaceBasis(x)
    if( missing(j) ) {
        structure(unclass(x)[i,,drop=FALSE],
                     class=class(x), basis=basis)
    } else {
        unclass(x)[i,j,...]
    }
}