R/operations.R

In DistributionalDDA: Distributional Data Analysis in Bayes Spaces

## Operations on BData objects

# ====== algebra overwritting ==========
#' Operations for BData objects
#'
#' BData objects are considered a elements of a BayesSpace and
#' operated accordingly to Bayes space rules, using the functions
#' [Add], [Sub], [Mul], [Div]. For operating with a scalar, please use
#' the BData object first.
#'
#' @param x a BData value to be operated
#' @param y another value, like e.g. another BData value or a scalar
#' @returns A BData object of the same class as `x`, resulting from the
#'   corresponding Bayes-space operation (addition, subtraction, scalar
#'   multiplication, or scalar division).
#' @export
#' @rdname BDataOperations
#'
#' @examples
#' X <- BDataSplines(t(1:5),B2Unif01Basis5)
#' X+X
"+.BData" <- function(x,y) Add(x,y)
#' @export
#' @describeIn BDataOperations Bayes space subtraction divides densities
"-.BData" <- function(x,y) Sub(x,y)
#' @export
#' @describeIn BDataOperations Bayes space multiplication works only with scalars
"*.BData" <- function(x,y) Mul(x,y)
#' @export
#' @method / BData
#' @describeIn BDataOperations Bayes space division only works with scalars
#'     and leads to do the y-th root
"/.BData" <- function(x,y) Div(x,y)

#' Computes the scalar product of (datasets of) vectors from various spaces
#'
#' The scalar product might depend on the space in which the operation
#' takes place.
#'
#' @param x an element of the space (or a dataset of such elements)
#' @param y an element of the space (or a dataset of such elements)
#' @param space the space in which the scalar product shall be computed
#' @param ... further generic arguments, like e.g. method for the compution
#'
#' @returns A numeric scalar, vector, or matrix of inner products, depending on
#'   the inputs.
#' @export
#'
#' @examples
#' X <- BDataSplines(matrix(c(1,2,3,4,5, 2,1,0,1,2), nrow=2, byrow=TRUE),
#'                   B2Unif01Basis5)
#' scalarProduct(X, X)
scalarProduct = function(x,y,space=NULL,...) UseMethod("scalarProduct",x)


#' Pool absolute frequencies
#'
#' Sum together the absolute frequencies of two representations with this information
#'
#' @param x a compatible BData representation
#' @param y another compatible BData representation
#' @param ... additional arguments (mostly internal)
#' @export
#'
#' @return An object of the least informative kind from the two pooled together,
#' which contains the added absolute frequencies of both input objects. The Bayes Space
#' basis is controlled by the first argument.
addFrequency = function(x,y, ...){
  UseMethod("addFrequency",x)
}


# ====== Exp-Ln mappings ==========


## Operations on functions
### Ln
#' @describeIn Exp Inverse of Exp
#' @export
Ln = function(f) UseMethod("Ln",f)
#' @describeIn Exp method for numerics just applies the natural logarithm
#' @export
Ln.numeric = function(f) log(f)
#' @describeIn Exp Inverse functional mapping
#' @export
Ln.function = function(f) {
    b = as.list(body(f))
    if( length(b) > 0 && b[[1]]==as.name("Exp") ) {
        body(f)=b[[2]]
        f
    } else function(...) Ln(f(...))
}
#' @describeIn Exp list method operates in a vectorized way over the list
#' @export
Ln.list = function(f) lapply(f,Ln)


### Exp
#' Functional mapping a function to its exp version
#'
#' Exp maps f(x) to Exp(f(x)). Ln is its inverse operation.
#'
#' @param f the function (or a constant)
#'
#' @return same type as f but exponentiated
#' @family Exp
#' @export
#'
#' @examples
#' Exp(1)
#' f = function(x) x^2
#' expf = Exp(f)
#' expf
#' lnf = Ln(f)
#' Exp(lnf)
#' Ln(expf)
#'
#'
Exp = function(f) UseMethod("Exp",f)
#' @describeIn Exp numeric method just evaluates exp
#' @export
Exp.numeric = function(f) exp(f)

#' @describeIn Exp Function method maps the function
#' @export
Exp.function = function(f) {
    b = as.list(body(f))
    if( length(b) > 0 && b[[1]]==as.name("Ln") ) {
        body(f)=b[[2]]
        f
    } else function(...) Exp(f(...))
}

#' @describeIn Exp list method iterates in a vectorized manner over the list
#' @export
Exp.list = function(f) {
    lapply(f,Exp)
}

# ====== basic-R Add definitions ==========

# Add
#' Adds functional objects
#'
#' This additional is also able to add functions and BayesSpace
#' objects.  The sum of two functions is a function that results in
#' the sum of the two results. The sum of two BayesSpace objects (in
#' BData), is given by the BayesSpace operations, which is given by
#' the product of the densities with respect to the reference measure.
#'
#' @param x a functional object
#' @param y a functional object
#'
#' @returns A BData object of the same subclass as the first argument (the
#'   Bayes-space sum / difference / scalar multiple / quotient).
#' @export
#' @family Add
#'
#' @examples
#' Add(sin,cos)
#' Add(list(sin,cos),sin)
#' Add(cos,list(sin,cos))
#' Add(list(sin,cos),list(sin,cos))
#' Add(list(sin,cos),10)
#'
Add = function(x,y) UseMethod("Add",x)
#' @export
#' @describeIn Add method for functions
#' @method Add function
Add.function = function(x,y) UseMethod("Add.function",y)
#' @export
#' @describeIn Add method for lists, iterating over them
#' @method Add list
Add.list = function(x,y) UseMethod("Add.list",y)
#' @export
#' @describeIn Add method for numerics, considering them as constant functions
#' @method Add numeric
Add.numeric = function(x,y) UseMethod("Add.numeric",y)
#' @export
#' @describeIn Add method for adding two functions, doing it by creating a function adding their values.
#' @method Add.function function
Add.function.function = function(x,y) function(...) x(...)+y(...)
#' @export
#' @describeIn Add method adding a function to a list, creating a list with entries for each list
#' @method Add.function list
Add.function.list = function(x,y) lapply(y,function(y) Add(x,y))
#' @export
#' @describeIn Add method adding a constant to a function
#' @method Add.function numeric
Add.function.numeric = function(x,y) {
        function(...) x(...)+y
}
#' @export
#' @describeIn Add method for adding a constant to a function
#' @method Add.numeric function
Add.numeric.function = function(x,y) {
    if( length(x)==1 ) {
        function(...) x+y(...)
    } else {
        mapply(Add,x,y)
    }
}
#' @export
#' @describeIn Add method for two numerics, simply adding the values
#' @method Add.numeric numeric
Add.numeric.numeric = function(x,y) x+y
#' @export
#' @describeIn Add method for adding a numeric to every value of the list
#' @method Add.list numeric
Add.list.numeric = function(x,y) mapply(Add,x,y)
#' @export
#' @describeIn Add method for adding a numeric to every value of the list
#' @method Add.numeric list
Add.numeric.list = function(x,y) mapply(Add,x,y)
#' @export
#' @describeIn Add method for vectorized addition of list elements
#' @method Add.list list
Add.list.list = function(x,y) mapply(Add,x,y)
#' @export
#' @describeIn Add method for adding a function to every element of a list
#' @method Add.list function
Add.list.function = function(x,y) lapply(x,function(x) Add(x,y))



# ====== basic-R Sub definitions ==========

# Sub
#' Substract functional objects
#'
#' This subtraction is also able to add functions and BayesSpace
#' objects.  The sum of two functions is a function that results in
#' the difference of the two results. The difference of two BayesSpace
#' objects (in BData), is given by the BayesSpace operations, which is
#' given by the quotient of the densities with respect to the reference
#' measure.
#'
#' @param x a functional object
#' @param y a functional object
#'
#' @returns A BData object of the same subclass as the first argument (the
#'   Bayes-space sum / difference / scalar multiple / quotient).
#' @export
#' @family Add
#' @examples
#' Sub(sin,cos)
#' Sub(list(sin,cos),sin)
#' Sub(cos,list(sin,cos))
#' Sub(list(sin,cos),list(sin,cos))
#' Sub(list(sin,cos),10)
#'
Sub = function(x,y) UseMethod("Sub",x)
#' @export
#' @describeIn Sub method for functions
#' @method Sub function
Sub.function = function(x,y) UseMethod("Sub.function",y)
#' @export
#' @describeIn Sub method for lists, iterating over the list
#' @method Sub list
Sub.list = function(x,y) UseMethod("Sub.list",y)
#' @export
#' @describeIn Sub method for numerics, treating them as constant
#'     functions
#' @method Sub numeric
Sub.numeric = function(x,y) UseMethod("Sub.numeric",y)
#' @export
#' @describeIn Sub method for two functions, returning a function
#'     giving the difference of the two as a result
#' @method Sub.function function
Sub.function.function = function(x,y)
function(...) x(...)-y(...)
#' @export
#' @describeIn Sub method for a function and list, iterating over the
#'     list
#' @method Sub.list function
Sub.list.function = function(x,y) lapply(x,function(x) Sub(x,y))
#' @export
#' @describeIn Sub method for a function and list, iterating over the
#'     list
#' @method Sub.function list
Sub.function.list = function(x,y) lapply(y,function(y) Sub(x,y))
#' @export
#' @describeIn Sub method substracting a constant value from a
#'     function
#' @method Sub.function numeric
Sub.function.numeric = function(x,y) {
   function(...) x(...)-y
}
#' @export
#' @describeIn Sub method substracting a function from a constant
#' @method Sub.numeric function
Sub.numeric.function = function(x,y) {
    if( length(x)==1 ) { function(...) x-y(...)  } else { mapply(Sub,x,y)
} }
#' @export
#' @describeIn Sub method substracting two values
#' @method Sub.numeric numeric
Sub.numeric.numeric = function(x,y) x-y
#' @export
#' @describeIn Sub method for substracting a value from each element
#'     of a list, returning a list.
#' @method Sub.list numeric
Sub.list.numeric = function(x,y)
mapply(Sub,x,y)
#' @export
#' @describeIn Sub method each element of a list from a constant
#' @method Sub.numeric list
Sub.numeric.list = function(x,y)
mapply(Sub,x,y)
#' @export
#' @describeIn Sub method for vectorized application of sub on the
#'     elements of lists
#' @method Sub.list list
Sub.list.list = function(x,y) mapply(Sub,x,y)


# ====== basic-R Mul definitions ==========

### Mul
#' Multiplication of functional objects
#'
#' This multiplication is also able to multiply functions and
#' BayesSpace objects with constants.  The product of two functions is a
#' function that results in the product of the two results. Product
#' of a BayesSpace objects (in BData), is given by the
#' BayesSpace operations, which is given powering the
#' density with respect to the reference measure, with the given scalar.
#'
#' @param x a functional object
#' @param y a functional object
#'
#' @returns A BData object of the same subclass as the first argument (the
#'   Bayes-space sum / difference / scalar multiple / quotient).
#' @export
#' @family Add
#' @examples
#' Mul(sin,cos)
#' Mul(list(sin,cos),sin)
#' Mul(cos,list(sin,cos))
#' Mul(list(sin,cos),list(sin,cos))
#' Mul(list(sin,cos),10)
Mul = function(x,y) UseMethod("Mul",x)
#' @export
#' @describeIn Mul method for multiplying functions
#' @method Mul function
Mul.function = function(x,y) UseMethod("Mul.function",y)
#' @export
#' @describeIn Mul method for vectorized operation over lists
#' @method Mul list
Mul.list = function(x,y) UseMethod("Mul.list",y)
#' @export
#' @describeIn Mul method for numerical values treated as constant functions
#' @method Mul numeric
Mul.numeric = function(x,y) UseMethod("Mul.numeric",y)
#' @export
#' @describeIn Mul method for operating two functions
#' @method Mul.function function
Mul.function.function = function(x,y) function(...) x(...)*y(...)
#' @export
#' @describeIn Mul method for multiplying a function with each element of a list
#' @method Mul.list function
Mul.list.function = function(x,y) lapply(x,function(x) Mul(x,y))
#' @export
#' @describeIn Mul method for multiplying a function with each element of a list
#' @method Mul.function list
Mul.function.list = function(x,y) lapply(y,function(y) Mul(x,y))
#' @export
#' @describeIn Mul method for multipying a function with a numeric value
#' @method Mul.function numeric
Mul.function.numeric = function(x,y) {
    function(...) x(...)*y
}
#' @export
#' @describeIn Mul method for multiplying numeric constant with a function
#' @method Mul.numeric function
Mul.numeric.function = function(x,y) {
    if( length(x)==1 ) {
        function(...) x*y(...)
    } else {
        mapply(Mul,x,y)
    }
}
#' @export
#' @describeIn Mul method for multiplying two constant values
#' @method Mul.numeric numeric
Mul.numeric.numeric = function(x,y) x*y

#' @export
#' @describeIn Mul method for multiplying a numeric object by a BData object
#' @method Mul.numeric BData
Mul.numeric.BData = function(x,y) y*x

#' @export
#' @describeIn Mul method for multiplying each element of a list with a numeric value
#' @method Mul.list numeric
Mul.list.numeric = function(x,y) mapply(Mul,x,y)
#' @export
#' @describeIn Mul method for multiplying a numeric value with each element of a list
#' @method Mul.numeric list
Mul.numeric.list = function(x,y) mapply(Mul,x,y)
#' @export
#' @describeIn Mul method for vectorized multiplication of the elements of a list.
#' @method Mul.list list
Mul.list.list = function(x,y) mapply(Mul,x,y)



# ====== basic-R Div definitions ==========


### Div
#' Division of functional objects
#'
#' This division is also able to divide functions and
#' BayesSpace objects with constants.  The quotient of two functions is a
#' function that results in the quotient of the two results. The quotient of
#' a BayesSpace objects (in BData) with a constant, is given by the
#' BayesSpace operations of powering with the inverse value 1/y.
#'
#' @param x a functional object
#' @param y a functional object
#'
#' @returns A BData object of the same subclass as the first argument (the
#'   Bayes-space sum / difference / scalar multiple / quotient).
#' @export
#' @family Add
#' @examples
#' Div(sin,cos)
#' Div(list(sin,cos),sin)
#' Div(cos,list(sin,cos))
#' Div(list(sin,cos),list(sin,cos))
#' Div(list(sin,cos),10)
#'
Div = function(x,y) UseMethod("Div",x)
#' @export
#' @describeIn Div method for functions
#' @method Div function
Div.function = function(x,y) UseMethod("Div.function",y)
#' @export
#' @describeIn Div method for vectorizing over lists
#' @method Div list
Div.list = function(x,y) UseMethod("Div.list",y)
#' @export
#' @describeIn Div method for numerical values seen as constant functions
#' @method Div numeric
Div.numeric = function(x,y) UseMethod("Div.numeric",y)
#' @export
#' @describeIn Div method for operating any BData object (argument `x`)
#' by a numeric value (argument `y`); defers calculations to `Mul()`
#' @method Div BData
Div.BData = function(x,y) Mul(x,1/y)
#' @export
#' @describeIn Div method for deviding two functions
#' @method Div.function function
Div.function.function = function(x,y) function(...) x(...)/y(...)
#' @export
#' @describeIn Div method for multiplying a function with each element of a list
#' @method Div.list function
Div.list.function = function(x,y) lapply(x,function(x) Div(x,y))
#' @export
#' @describeIn Div method for multiplying a function with each element of a list
#' @method Div.function list
Div.function.list = function(x,y) lapply(y,function(y) Div(x,y))
#' @export
#' @describeIn Div method for dividing a function with a constant
#' @method Div.function numeric
Div.function.numeric = function(x,y) {
    function(...) x(...)/y
}
#' @export
#' @describeIn Div method for deviding a constant with a function
#' @method Div.numeric function
Div.numeric.function = function(x,y) {
    if( length(x)==1 ) {
        function(...) x/y(...)
    } else {
        mapply(Div,x,y)
    }
}
#' @export
#' @describeIn Div method for deviding two numeric values
#' @method Div.numeric numeric
Div.numeric.numeric = function(x,y) x/y
#' @export
#' @describeIn Div method for dividing each element of a list with a constant
#' @method Div.list numeric
Div.list.numeric = function(x,y) mapply(Div,x,y)
#' @export
#' @describeIn Div method for deviding a constant with each element of a list
#' @method Div.numeric list
Div.numeric.list = function(x,y) mapply(Div,x,y)
#' @export
#' @describeIn Div method for vectorized division of the elements of lists
#' @method Div.list list
Div.list.list = function(x,y) mapply(Div,x,y)





# ====== Add for BData Objects ==========

#' @export
#' @describeIn Add method for BDataDensities (effectively multiplying the densities as addition in the BayesSpace)
#' @method Add BDataDensities
Add.BDataDensities = function(x,y) UseMethod("Add.BDataDensities",y)
#' @export
#' @describeIn Add method for BDataLogDensities
#' @method Add BDataLogDensities
Add.BDataLogDensities = function(x,y) UseMethod("Add.BDataLogDensities",y)

#' @export
#' @describeIn Add method for BDataSplines
#' @method Add BDataSplines
Add.BDataSplines = function(x,y) UseMethod("Add.BDataSplines",y)

#' @export
#' @describeIn Add method for BDataHistograms
#' @method Add BDataHistograms
Add.BDataHistograms = function(x,y) UseMethod("Add.BDataHistograms",y)


#' @export
#' @describeIn Add method for two BDataDensities
#' @method Add.BDataDensities BDataDensities
Add.BDataDensities.BDataDensities = function(x,y) {
    stopifnot(hasSameBayesSpace(x,y))
    BDataDensities(mapply(Mul,unclass(x),unclass(y)),
                   basis=getBayesSpaceBasis(x))
}

#' @export
#' @describeIn Add method for adding two BDataLogDensities
#' @method Add.BDataLogDensities BDataLogDensities
Add.BDataLogDensities.BDataLogDensities = function(x,y) {
    stopifnot(hasSameBayesSpace(x,y))
    BDataLogDensities(mapply(Add,x,y),basis=getBayesSpaceBasis(x))
}

#' @export
#' @describeIn Add method converting the second argument to a BDataLogDensity and then adding
#' @method Add.BDataLogDensities default
Add.BDataLogDensities.default = function(x,y) {
    Add(x,as.BDataLogDensities(y))
}

#' @export
#' @describeIn Add method converting the second argument to a BDataDensity object and then adding.
#' @method Add.BDataDensities default
Add.BDataDensities.default = function(x,y) {
    Add(x,as.BDataDensities(y))
}

#' @export
Add.BDataSplines = function(x,y) UseMethod("Add.BDataSplines",y)

#' @export
#' @describeIn Add method for two BDataSplines
#' @method Add.BDataSplines BDataSplines
Add.BDataSplines.BDataSplines = function(x,y) {
    stopifnot(hasSameBayesSpaceBasis(x,y))
    stopifnot(dim(x)==dim(y)||(ncol(x)==ncol(y)&&(nrow(x)==1||nrow(y)==1)))
    BDataSplines(t(unclass(t(x))+unclass(t(y))),basis=getBayesSpaceBasis(x))
}



#' @export
#' @describeIn Add method for BDataSplines default
#' @method Add.BDataSplines default
Add.BDataSplines.default = function(x,y) {
    stopifnot(hasSameBayesSpaceBasis(x,y))
    Add(x,as.BDataSplines(y,basis=getBayesSpaceBasis(x)))
}


#' @export
#' @describeIn Add method for BDataHistograms
#' @method Add.BDataHistograms BDataHistograms
Add.BDataHistograms.BDataHistograms = function(x,y) {
  b1 = getBreaks(x)
  b2 = getBreaks(y)
  stopifnot(equiv(b1[1], b2[1]), equiv(b1[length(b1)],b2[length(b2)]), hasSameBayesSpace(x,y))
  # operate on transposed data, to ensure that row-observations can operate with matrix-datasets
  if(all( equiv(b1,b2) )){
     return(BDataDensityHistograms( .apply2OneOrMany(unclass(x), unclass(y), FUN="*"), b1, getBayesSpaceBasis(x)))
  }
  # sum logs, then do exp
  NewLogDens = .apply2OneOrMany(FUN="+", log(unclass(x)) %*% sign( .distributingMatrix(b1, b2)), log(unclass(y)) %*% sign( .distributingMatrix(b2, b1))  )
  return(BDataDensityHistograms(exp(NewLogDens), sortUnique(b1,b2), getBayesSpaceBasis(x)))
}

# should not be exported, it is an internal thing
# @export
 .distributingMatrix = function(b1,b2){
  b0 = sortUnique(b1,b2)
  sapply(1:(length(b1)-1), function(i1){
    i0 = which.min((b1[i1]-b0)^2)  # == has a precision problem
    rs = (b0-b0[i0])/(b1[i1+1]-b1[i1])
    rs = ifelse(rs>0, ifelse(rs<=1,rs,1), 0)
    return(rs)
  }) |> apply(2, diff) |> t()
}


.apply2OneOrMany = function(x,y,FUN){
  # assumes x and y both to be matrices, eventually one of them 1-row matrix
  if(nrow(x)==nrow(y)){
    return( do.call(FUN, list(x, y)) )
  }else if(nrow(x)==1){
    return( sweep(y, 2, c(x), FUN) )
  }else if(nrow(y)==1){
    return( .apply2OneOrMany(y,x,FUN=FUN) )
  }else{
    stop("Arguments 'x' and 'y' must either have the same row.number/sample.length, or one of them must have one row/length")
  }
}

# ====== addFrequency for BData Objects ==========


#' @export
#' @method addFrequency BDataFrequencyHistograms
addFrequency.BDataFrequencyHistograms = function(x,y,...) UseMethod("addFrequency.BDataFrequencyHistograms",y)

#' @export
#' @method addFrequency BDataUnivariateSamples
addFrequency.BDataUnivariateSamples = function(x,y,...) UseMethod("addFrequency.BDataUnivariateSamples",y)


#' @describeIn addFrequency method for two "BDataFrequencyHistograms" object. Returns
#' a "BDataFrequencyHistograms" where the frequencies are summed (if breaks are
#' equivalent), or where they are distributed uniformly along the interval bins and
#' rebinned by summing with a break set resulting from the combination of both
#' input bin sets.
#' @export
#' @method addFrequency.BDataFrequencyHistograms BDataFrequencyHistograms
#' @examples
#' xh = lapply(rep(1000,4), runif) |> as.BDataFrequencyHistograms(basis=B2Unif01Basis5)
#' yh = lapply(rep(1000,4), runif) |> as.BDataFrequencyHistograms(basis=B2Unif01Basis5)
#' addFrequency(xh,yh)
#' xh1 = list(runif(300)) |> as.BDataFrequencyHistograms(basis=B2Unif01Basis5)
#' addFrequency(xh,xh1)
#' xh1 = list(runif(300)) |> as.BDataFrequencyHistograms(basis=B2Unif01Basis5, breaks=c(0:4)/4)
#' addFrequency(xh,xh1)
addFrequency.BDataFrequencyHistograms.BDataFrequencyHistograms = function(x,y,...) {
  b1 = getBreaks(x)
  b2 = getBreaks(y)
  stopifnot( equiv(b1[1],b2[1]), equiv(b1[length(b1)],b2[length(b2)]), hasSameBayesSpace(x,y) )
  # operate on transposed data, to ensure that row-observations can operate with matrix-datasets
  if(length(b1)==length(b2) && all(equiv(b1,b2))){
    return(BDataFrequencyHistograms(.apply2OneOrMany(unclass(x),unclass(y),"+"), b1, getBayesSpaceBasis(x)))
  }
  NewCounts = .apply2OneOrMany(FUN="+", unclass(x) %*% .distributingMatrix(b1, b2), unclass(y) %*% .distributingMatrix(b2, b1))
  return(BDataFrequencyHistograms(NewCounts, sortUnique(b1,b2), getBayesSpaceBasis(x)))
}


#' @describeIn addFrequency method for two "BDataUnivariateSamples" object. Returns
#' a "BDataUnivariateSamples" where the observations are pooled together (in the order
#' they are given), pair by pair. It can recycle a "BDataUnivariateSamples" of one
#' single population with a multi-population "BDataUnivariateSamples", but will
#' not work with any other combination of "BDataUnivariateSamples" of different
#' lengths.
#' @export
#' @method addFrequency.BDataUnivariateSamples BDataUnivariateSamples
#' @examples
#' xs = lapply(rep(200,4), runif) |> BDataSamples(B2Unif01Basis5)
#' ys = lapply(rep(10,4), runif) |> BDataSamples(B2Unif01Basis5)
#' addFrequency(xs,ys) |> summary()
#' x1 = BDataSamples(list(runif(100)), B2Unif01Basis5)
#' addFrequency(xs, x1)  |> summary()
addFrequency.BDataUnivariateSamples.BDataUnivariateSamples = function(x,y,...){
  stopifnot(hasSameBayesSpace(x,y))
  if( length(x)==1 ){
    res = lapply(unclass(y), function(yy) c(x[[1]],yy) )
  }else if( length(y)==1 ){
    res = lapply(unclass(x), function(xx) c(xx,y[[1]]) )
  }else if( length(x)==length(y)  ){
    res = lapply(1:length(x), function(i) c(x[[i]], y[[i]])) # purrr::map2(unclass(x),unclass(y),c)
  }else{
    stop("lengths of 'x' and 'y' provided do not match")
  }
  BDataSamples(res, basis=getBayesSpaceBasis(x))
}



#' @describeIn addFrequency method for a "BDataFrequencyHistograms" and a "BDataUnivariateSamples"
#' objects; classifies the data sets from "BDataUnivariateSamples" into the bins
#' of "BDataFrequencyHistograms". As with other methods of this function, it only works for
#' `x` and `y` of the same length, or else either `x` or `y` must have length 1.
#' @param .baseSpaceBasis which base space basis to use, defaults to the one specified by the first argument `x`
#' @export
#' @method addFrequency.BDataFrequencyHistograms BDataUnivariateSamples
#' @examples
#' xs = lapply(rep(200,4), runif) |> BDataSamples(B2Unif01Basis5)
#' yh = lapply(rep(1000,4), runif) |> as.BDataFrequencyHistograms(basis=B2Unif01Basis5)
#' addFrequency(xs,yh)
addFrequency.BDataFrequencyHistograms.BDataUnivariateSamples = function(x,y, .baseSpaceBasis = getBayesSpaceBasis(x),...) {
  stopifnot(hasSameBayesSpace(x,y))
  addFrequency(x, as.BDataFrequencyHistograms(y, breaks =getBreaks(x), basis=.baseSpaceBasis) )
}


#' @describeIn addFrequency method for a "BDataUnivariateSamples" and a "BDataFrequencyHistograms" objects; same as in the reverse order
#' @export
#' @method addFrequency.BDataUnivariateSamples BDataFrequencyHistograms
addFrequency.BDataUnivariateSamples.BDataFrequencyHistograms = function(x,y, .baseSpaceBasis = getBayesSpaceBasis(x),...) addFrequency(y,x, .baseSpaceBasis)


# ====== Sub for BData Objects ==========


#' @export
#' @describeIn Sub method for BDataHistograms (operating in BayesSpaces)
#' @method Sub BDataHistograms
Sub.BDataHistograms = function(x,y) UseMethod("Sub.BDataHistograms",y)


#' @export
#' @describeIn Sub method for two BDataHistograms
#' @method Sub.BDataHistograms BDataHistograms
Sub.BDataHistograms.BDataHistograms = function(x,y) {
  b1 = getBreaks(x)
  b2 = getBreaks(y)
  stopifnot(equiv(b1[1], b2[1]), equiv(b1[length(b1)],b2[length(b2)]), hasSameBayesSpace(x,y))
  # operate on transposed data, to ensure that row-observations can operate with matrix-datasets
  if(all( equiv(b1,b2) )){
    return(BDataDensityHistograms( .apply2OneOrMany(unclass(x), unclass(y), FUN="/"), b1, getBayesSpaceBasis(x)))
  }
  # subtract logs, then do exp
  NewLogDens = .apply2OneOrMany(FUN="-", log(unclass(x)) %*% sign( .distributingMatrix(b1, b2)), log(unclass(y)) %*% sign( .distributingMatrix(b2, b1))  )
  return(BDataDensityHistograms(exp(NewLogDens), sortUnique(b1,b2), getBayesSpaceBasis(x)))
}


# ====== Scalar Product for BData Objects ==========

#' @export
#' @describeIn scalarProduct of vectors (or rows of a matrix) potentially over the metric tensor of the space.
#' @method scalarProduct numeric
scalarProduct.numeric = function(x,y=x,space,...) UseMethod("scalarProduct.numeric",y)
#' @export
scalarProduct.numeric.numeric = function(x,y=x,space=NULL,...) {
    if( is.matrix(x) ) x=t(x)
    if( is.matrix(y) ) y=t(y)
    if( is.null(space) )
        t(x) %*% y
    else {
        V = getMetricTensor(space)
        c(t(x) %*% V %*% y)
    }
}

#' @export
#' @describeIn scalarProduct in BayesSpaces for BDataSplines objects
scalarProduct.BDataSplines = function(x,y,space=getBayesSpace(x),...) {
        stopifnot(hasSameBayesSpace(x,space))
        basis=getBayesSpaceBasis(x)
        V = getMetricTensor(basis)
        unclass(x) %*% V %*% t(unclass(as.BDataSplines(y,basis=basis)))
}

#' @export
#' @describeIn scalarProduct in BayesSpaces for BDataDensities objects
scalarProduct.BDataDensities = function(x,y=x,space=getBayesSpace(x),...) {
    x = as.BDataLogDensities(x)
    y = as.BDataDensities(y)
    scalarProduct(x,y,space=space)
}


#' @export
#' @describeIn scalarProduct in BayesSpaces for BDataLogDensities objects
scalarProduct.BDataLogDensities = function(x,y=x,space=getBayesSpace(x),...) {
    y = as.BDataLogDensities(y)
    stopifnot(hasSameBayesSpace(x,y))
    stopifnot(hasSameBayesSpace(x,space))
    measure = getReferenceMeasure(space)
    CovarianceOperator(measure,unclass(x),unclass(y))
}



## mean
#' Computes the distribution mean for BDataSplines.
#'
#' This function simply averages the spline coefficients and returns a
#' correponding spline.
#'
#' @param x the BDataSplines object over which to average
#' @param ... unused arguments required for generic use
#' @param method the method to be chosen (as in [cov.rob]), functions
#'     does operates according to BayesSpace formulas directly on the
#'     densities, classical uses the spline representation, mcd and
#'     mve apply these robust methods on the spline coefficients.
#'
#' @returns a BDataSplines with a single observation representing the mean
#' @export
#' @importFrom MASS cov.mcd
#' @importFrom MASS cov.mve
#' @method mean BDataSplines
#' @examples
#' data(DTD)
#' Xpdf    = BDataPDFs(DTD$shapePDF,basis=B2Unif01Basis5)
#' Xspline = as.BDataSplines(Xpdf)
#' plot(Xspline,col=2:12)
#' meanXspline = mean(Xspline)
#' plot(meanXspline,col=1,lwd=3,add=TRUE)
#'
#' plot(Xspline,as="PDF")
#' plot(meanXspline,col="red",lwd=3,add=TRUE,as="PDF")
#'
#' Xdensities = as.BDataDensities(Xpdf)
#'
#' plot(Xdensities)
#' meanXdensities = clo(mean(Xdensities))
#' plot(meanXdensities,col="green",lwd=7,add=TRUE)
#' plot(meanXspline,col="red",lwd=3,add=TRUE,as="density")
#'
#' plot(Xpdf)
#' meanXpdf = mean(Xpdf)
#' plot(clo(meanXpdf),col="green",lwd=7,add=TRUE)
#' plot(meanXspline,col="red",lwd=3,add=TRUE,as="PDF",close=TRUE)
#' plot(as.BDataPDFs(meanXdensities),col="blue",lwd=3,add=TRUE,close=TRUE)
#' plot(as.BDataPDFs(mean(as.BDataLogPDFs(Xpdf))),col="orange",lwd=3,add=TRUE,close=TRUE)
#' plot(as.BDataPDFs(mean(as.BDataLogDensities(Xpdf))),col="magenta",lwd=3,add=TRUE,close=TRUE)
#'
mean.BDataSplines = function(x,...,method=c("classical","mcd","mve")) {
    switch(match.arg(method),
           classical=BDataSplines(t(colMeans(x)),basis=getBayesSpaceBasis(x)),
           mcd=BDataSplines(cov.mcd(unclass(x),...)$center,basis=getBayesSpaceBasis(x)),
           mve=BDataSplines(cov.mve(unclass(x),...)$center,basis=getBayesSpaceBasis(x))
           )
}

#' @describeIn mean.BDataSplines method for full densities
#' @export
#' @method mean BDataDensities
#' @importFrom compositions geometricmean
mean.BDataDensities = function(x,...,method=c("functions","classical","mcd","mve")) {
    method=match.arg(method)
    if(method=="functions") {
        densities=x
        BDataDensities(list(function(x) {
            if( numberOfCases(x)==1 ) {
                compositions::geometricmean(sapply(densities,function(f) f(x)),...)
            } else {
                apply(sapply(densities,function(f) f(x)),1,compositions::geometricmean,...)
            }
        }),basis=getBayesSpaceBasis(x))
    }  else {
        mean(as.BDataSplines(x),...,method=method)
    }
}

#' @describeIn mean.BDataSplines method for full PDFs
#' @export
#' @method mean BDataPDFs
#' @importFrom compositions geometricmean
mean.BDataPDFs = function(x,...,method=c("functions","classical","mcd","mve")) {
    method=match.arg(method)
    if(method=="functions") {
        densities=x
        pdf = getPDF(getReferenceMeasure(getBayesSpaceBasis(x)))
        BDataPDFs(list(function(x) {
            if( numberOfCases(x)==1 ) {
                geometricmean(sapply(densities,function(f) f(x)),...)
            } else {
                apply(sapply(densities,function(f) f(x)),1,geometricmean,...)
            }
        }),basis=getBayesSpaceBasis(x))
    }  else {
        mean(as.BDataSplines(x),...,method=method)
    }
}

#' @describeIn mean.BDataSplines method for full logPDFs
#' @export
#' @method mean BDataLogPDFs
mean.BDataLogPDFs = function(x,...,method=c("functions","classical","mcd","mve")) {
    method=match.arg(method)
    if(method=="functions") {
        densities=x
        pdf = getPDF(getReferenceMeasure(getBayesSpaceBasis(x)))
        BDataLogPDFs(list(function(x) {
            if( numberOfCases(x)==1 ) {
                mean(sapply(densities,function(f) f(x)),...)
            } else {
                apply(sapply(densities,function(f) f(x)),1,mean,...)
            }
        }),basis=getBayesSpaceBasis(x))
    }  else {
        mean(as.BDataSplines(x),...,method=method)
    }
}

#' @describeIn mean.BDataSplines method for full logDensities
#' @export
#' @method mean BDataLogDensities
mean.BDataLogDensities = function(x,...,method=c("functions","classical","mcd","mve")) {
    method=match.arg(method)
    if(method=="functions") {
        densities=x
        pdf = getPDF(getReferenceMeasure(getBayesSpaceBasis(x)))
        BDataLogDensities(list(function(x) {
            if( numberOfCases(x)==1 ) {
                mean(sapply(densities,function(f) f(x)),...)
            } else {
                apply(sapply(densities,function(f) f(x)),1,mean,...)
            }
        }),basis=getBayesSpaceBasis(x))
    }  else {
        mean(as.BDataSplines(x),...,method=method)
    }
}


## end mean

## varKernel
#' Compute the kernel of the Variance Kovariance operator
#'
#' Variance in Bayes spaces can be expressed in various forms. The
#' variance kernel is essentially the covariance function of the
#' clr-transform of a random variable:
#'
#'     c(x,y) = cov( clr(X)(x),clr(X)(y))
#'
#' @param X a dataset of BData objects
#' @param ... further parameters detailing computation
#' @param basis the BayesSpaceBasis for which this covariance kernel is valid
#'
#' @return a function of two arguments with a basis attribute representing the BayesSpaceBasis, and varKernel class.
#' @export
#' @seealso plot.varKernel, variationKernel
#'
#' @examples
#' data(DTD)
#' X = BDataPDFs(DTD$shapePDF,basis=B2Unif01Basis5)
#' C = varKernel(X)
#' C
#' plot(C)
#'
#' K = BDataSplines(matrix(rnorm(10*5),ncol=5),basis=B2Unif01Basis5)
#' var(K)
#' C2 = varKernel(K)
#' C2
#' plot(C2)
#'
varKernel = function(X,...) UseMethod("varKernel",X)
#' @describeIn varKernel method for BDataSplines uses the covariance of the spline coefficients to compute
#' @export
# ' method varKernel BDataSplines
varKernel.BDataSplines = function(X,...) {
    basis=getBayesSpaceBasis(X)
    C=var(unclass(X))
    b=as.function(basis)
    structure(function(x,y) {
        b(x) %*% C %*% t(b(y))
    },C=C,
    basis=basis,
    class="varKernel")
}

#' @describeIn varKernel default method for BData computes implicitly using clr transforms.
#' @export
# ' @method varKernel BData
varKernel.BData = function(X,...) {
    basis=getBayesSpaceBasis(X)
    . = clr(X)
    CLR = function(x) t(. %..% c(.(x)))
    structure(function(x,y) {
        cov(CLR(x),CLR(y))
    },
    basis=basis,
    class="varKernel"
    )
}

#' @describeIn varKernel constructor
#' @export
#' @method varKernel function
varKernel.function = function(X,basis=getBayesSpaceBasis(X),...) structure(X,basis=basis,class="varKernel")


#' @describeIn varKernel method for varKernel, keeps the kernel, possibly overwrites the basis.
#' @export
#' @method varKernel varKernel
#'
varKernel.varKernel= function(X,basis=getBayesSpaceBasis(X),...) structure(X,basis=basis)





#' @describeIn getBayesSpaceBasis method for varKernel
#' @method getBayesSpaceBasis varKernel
#' @export
getBayesSpaceBasis.varKernel = function(X,...) attr(X,"basis")
#' @describeIn Support getter method for varKernel
#' @method getSupport varKernel
#' @export
getSupport.varKernel = function(x,...) getSupport(getBayesSpaceBasis(x))

###################### variation

## variationKernel
#' Compute the variation kernel
#'
#' Variance in Bayes spaces can be expressed in various forms. The
#' variation kernel is essentially the variance of the difference of
#' the log-densities or clr-transforms at evaluated at different values
#'
#'     v(x) = (var( clr(X)(x_i) - clr(X)(y_j)))_ij
#'
#' It relates to the covariance like the variogram of the residual
#' mean corrected field corresponds to the covariance. The advantage
#' of the variation kernel is its subcompositional coherence, i.e. it
#' is independent of support and reference measure.
#'
#' @param X a dataset of BData objects
#' @param ... further parameters detailing computation
#' @param basis the BayesSpaceBasis for which this covariance kernel is valid
#'
#' @return a function of two arguments with a basis attribute representing the BayesSpaceBasis, and variationKernel class.
#' @export
#' @seealso plot.variationKernel, varKernel
#'
#' @examples
#' data(DTD)
#' X = BDataPDFs(DTD$shapePDF,basis=B2Unif01Basis5)
#' V = variationKernel(X)
#' V
#' plot(V)
#'
#' K = BDataSplines(matrix(rnorm(10*5),ncol=5),basis=B2Unif01Basis5)
#' V2 = variationKernel(K)
#' V2
#' plot(V2)
#'
variationKernel = function(X,...) UseMethod("variationKernel",X)
#' @describeIn variationKernel method for BDataSplines uses the covariance of the spline coefficients to compute
#' @export
# ' method variationKernel BDataSplines
variationKernel.BDataSplines = function(X,...) {
    basis=getBayesSpaceBasis(X)
    C=var(unclass(X))
    b=as.function(basis)
    structure(function(x) {
        B = b(x)
        CC = B %*% C %*% t(B)
        if( length(CC) == 1 )
            return( 0 )
        v = diag(CC)
        one = rep(1,ncol(CC))
        outer(v,one)+outer(one,v)-2*CC
    },C=C,
    basis=basis,
    class="variationKernel")
}

#' @describeIn variationKernel default method for BData computes implicitly using clr transforms.
#' @export
# ' @method variationKernel BData
variationKernel.BData = function(X,...) {
    basis=getBayesSpaceBasis(X)
    . = as.BDataLogDensities(X)
    LogDensities = function(x) t(. %..% c(.(x)))
    structure(function(x) {
        C = var(LogDensities(x))
        if( length(C) == 1 )
            return( 0 )
        v = diag(C)
        one = rep(1,ncol(C))
        outer(v,one)+outer(one,v)-2*C
    },
    basis=basis,
    class="variationKernel"
    )
}

#' @describeIn variationKernel constructor
#' @export
#' @method variationKernel function
variationKernel.function = function(X,basis=getBayesSpaceBasis(X),...) structure(X,basis=basis,class="variationKernel")


#' @describeIn variationKernel method for varKernel, keeps the kernel, possibly overwrites the basis.
#' @export
#' @method variationKernel variationKernel
#'
variationKernel.variationKernel= function(X,basis=getBayesSpaceBasis(X),...) structure(X,basis=basis)





#' @describeIn getBayesSpaceBasis method for variationKernel
#' @method getBayesSpaceBasis variationKernel
#' @export
getBayesSpaceBasis.variationKernel = function(X,...) attr(X,"basis")
#' @describeIn Support getter method for variationKernel
#' @method getSupport variationKernel
#' @export
getSupport.variationKernel = function(x,...) getSupport(getBayesSpaceBasis(x))







## mvar - Metric Variance
#' Compute the metric variance
#'
#' The metric variance is the expected squared distance from the mean
#' in the metric of the enclosing space. By construction the metric
#' variance computed from a spline representation should be smaller than
#' metrix variances computed on functional data, because the spline captures
#' only a part of the variance.
#'
#' @param x some BData dataset.
#' @param ... further parameters for the methods to compute it (expectation)
#' @param meanx an external mean estimate
#' @param method the method of variance estimation
#'
#' @returns a positive real number representing the metric variance
#' @export
#' @method mvar BData
#' @importFrom compositions mvar
#' @examples
#' data(DTD)
#' Xpdf    = BDataPDFs(DTD$shapePDF,basis=B2Unif01Basis5)
#' Xsplines= as.BDataSplines(Xpdf)
#' mvar(Xsplines)
#' \donttest{ mvar(Xsplines,method="mcd")
#' mvar(as.BDataDensities(Xpdf))
#' mvar(as.BDataLogDensities(Xpdf))
#' mvar(as.BDataLogPDFs(Xpdf))
#' mvar(Xpdf)}
#'
mvar.BData = function(x,...,meanx=mean(x),method="classical") {
    RM = getReferenceMeasure(x)
    method=match.arg(method)
    sum(sapply(clr(x-meanx),
                    function(f) expectation(RM,function(x) f(x)^2,...)))/(numberOfCases(x)-1)
}

#' @describeIn mvar.BData a specialized method for BDataSpline, which
#'     does not require integration.
#' @export
#' @importFrom MASS cov.mcd
#' @importFrom MASS cov.mve
#' @method mvar BDataSplines
mvar.BDataSplines = function(x,...,method=c("classical","mcd","mve"),
                             meanx=mean(x,method)) {
    method = match.arg(method)
    basis = getBayesSpaceBasis(x)
    metric=getMetricTensor(basis)
    switch(method,
           classical={
               xc = (t(unclass(x))-c(unclass(as.BDataSplines(meanx,basis=basis))))
               sum(xc * (metric %*% xc))/(numberOfCases(x)-1)
           },
           #mcd=sum(diag(metric %*% MASS::cov.mcd(unclass(x),...)$cov)),
           #mve=sum(diag(metric %*% MASS::cov.mve(unclass(x),...)$cov))
           mcd=sum(diag(metric %*% covRobCollinear(unclass(x),...,method=method))),
           mve=sum(diag(metric %*% covRobCollinear(unclass(x),...,method=method)))
           )
}


#' @describeIn mvar.BData compute the metric standard deviation, proportional to the
#' square root of the metric variance
#' @export
#' @method msd BData
#' @importFrom compositions msd
msd.BData = function(x, ...){
  sqrt(mvar(x,...))
}



#' Scale and/or center a BData set
#'
#' Apply scaling and/or centering to a dataset of BData objects.
#' Diverts all calculations to `mean.BData`-instances as well as
#' to `msd.BData()` and the arithmetic operations.
#'
#' @param x  data set to transform
#' @param center either logical (default: TRUE) or the mean BData object to subtract
#' from the elements of `x`
#' @param scale either logical (default: TRUE) or the scalar to divide
#'              the (centered?) elements of `x` with
#' @param ... further generic arguments
#'
#' @returns a data set of the same kind and size of `x`
#' @export
#' @importFrom compositions scale
#'
#' @examples
#' bdl = BDataDensities(list(\(x)(dnorm(x,1,1)),\(x)(dnorm(x,-1,1)) ),basis=B2N01Basis5,closed=TRUE)
#' bdlsc = scale(bdl)
#' msd(bdlsc) # should now be 1
#' mean(bdlsc)[[1]](-5:5) # should be 1 everywhere, i.e. the neutral element
scale.BData = function(x,center=TRUE, scale=TRUE,...){
  if(is.logical(scale)){
    if(scale){
      scale = msd(x)
    }else{
      scale = 1
    }
  }
  if(is.logical(center)){
    if(center){
      return((x-mean(x))/scale)
    }else{
      return(x/msd)
    }
  }
}


#' @importFrom MASS cov.rob
covRobCollinear = function(x,...,method=c("mcd","mve","classical"),condition=1E-5) {
    SVD=svd(t(t(x)-colMeans(x)))
    lambdaMax=max(abs(SVD$d))
    take=(abs(SVD$d)>=lambdaMax*condition)
    if( all(take) )
        return(MASS::cov.rob(x,...,method=method))
    if( sum(take) == 0 )
        return( matrix(0,ncol=ncol(x),nrow=nrow(x)))
    xr = SVD$u[,take]
    C=MASS::cov.rob(xr,...,method=method)$cov
    D=SVD$d[take]
    SVD$v[,take]%*%(C * outer(D,D))%*%t(SVD$v[,take])
}
## Add BData

#' @describeIn Add BayesSpaceAddition for BDataDensities
#' @export
#' @method Add BDataDensities
Add.BDataDensities = function(x,y) {
    x=as.BDataDensities(x)
    basis=getBayesSpaceBasis(x)
    y=as.BDataDensities(y,basis=basis)
    BDataDensities(Mul(unclass(x),unclass(y)),basis=basis)
}

#' @describeIn Add BayesSpaceAddition for BDataSplines
#' @export
#' @method Add BDataSplines
Add.BDataSplines = function(x,y) {
    x=as.BDataSplines(x)
    basis=getBayesSpaceBasis(x)
    y=as.BDataSplines(y,basis=basis)
    BDataSplines(addDatasets(unclass(x),unclass(y)),basis=basis)
}



#' @describeIn Add BayesSpaceAddition for BDataLogDensities
#' @export
#' @method Add BDataLogDensities
Add.BDataLogDensities = function(x,y) {
    x=as.BDataLogDensities(x)
    basis=getBayesSpaceBasis(x)
    y=as.BDataLogDensities(y,basis=basis)
    BDataLogDensities(Add(unclass(x),unclass(y)),basis=basis)
}

#' @describeIn Add BayesSpaceAddition for BDataPDFs
#' @export
#' @method Add BDataPDFs
Add.BDataPDFs = function(x,y) {
    x=as.BDataPDFs(x)
    basis=getBayesSpaceBasis(x)
    y=as.BDataDensities(y,basis=basis) # sic! Densities
    BDataPDFs(Mul(unclass(x),unclass(y)),basis=basis)
}

#' @describeIn Add BayesSpaceAddition for BDataLogPDFs
#' @export
#' @method Add BDataLogPDFs
Add.BDataLogPDFs = function(x,y) {
    x=as.BDataLogPDFs(x)
    basis=getBayesSpaceBasis(x)
    y=as.BDataLogDensities(y,basis=basis) # sic! LogDensities
    BDataLogPDFs(Add(unclass(x),unclass(y)),basis=basis)
}

## Substract
#' @describeIn Sub BayesSpace Difference for BDataDensities
#' @export
#' @method Sub BDataDensities
Sub.BDataDensities <- function(x,y) {
    x=as.BDataDensities(x)
    basis=getBayesSpaceBasis(x)
    y=as.BDataDensities(y,basis=basis)
    BDataDensities(Div(unclass(x),unclass(y)),basis=basis)
}

#' @describeIn Sub BayesSpace Difference for BDataSplines
#' @export
#' @method Sub BDataSplines
Sub.BDataSplines = function(x,y) {
    x=as.BDataSplines(x)
    basis=getBayesSpaceBasis(x)
    y=as.BDataSplines(y,basis=basis)
    BDataSplines(addDatasets(unclass(x),-unclass(y)),basis=basis)
}



#' @describeIn Sub BayesSpace Difference for BDataLogDensities
#' @export
#' @method Sub BDataLogDensities
Sub.BDataLogDensities = function(x,y) {
    x=as.BDataLogDensities(x)
    basis=getBayesSpaceBasis(x)
    y=as.BDataLogDensities(y,basis=basis)
    BDataLogDensities(Sub(unclass(x),unclass(y)),basis=basis)
}

#' @describeIn Sub BayesSpace Difference for BDataPDFs
#' @export
#' @method Sub BDataPDFs
Sub.BDataPDFs = function(x,y) {
    x=as.BDataPDFs(x)
    basis=getBayesSpaceBasis(x)
    y=as.BDataDensities(y,basis=basis)
    BDataPDFs(Div(unclass(x),unclass(y)),basis=basis)
}

#' @describeIn Sub BayesSpaceAddition for BDataLogPDFs
#' @export
#' @method Sub BDataLogPDFs
Sub.BDataLogPDFs = function(x,y) {
    x=as.BDataLogPDFs(x)
    basis=getBayesSpaceBasis(x)
    y=as.BDataLogDensities(y,basis=basis)
    BDataLogPDFs(Sub(unclass(x),unclass(y)),basis=basis)
}


## Scale

#' @describeIn Mul BayesSpace scaling multiplication
#' @export
#' @method Mul BDataDensities
Mul.BDataDensities = function(x,y) {
    if( !is.numeric(y) ) stop("Mul.BDataDensities works only with numeric `y`") # return(Mul(y,x))
    x=as.BDataDensities(x)
    basis=getBayesSpaceBasis(x)

    BDataDensities(mapply(function(f,y) function(x) f(x)^y,x,y),
                   basis=basis, closed = isClosed(x)*all(y==1))
}

#' @describeIn Mul BayesSpace scaling multiplication
#' @export
#' @method Mul BDataLogDensities
Mul.BDataLogDensities = function(x,y) {
  if( !is.numeric(y) ) stop("Mul.BDataLogDensities works only with numeric `y`") # return(Mul(y,x))
  x=as.BDataLogDensities(x)
    basis=getBayesSpaceBasis(x)
    BDataLogDensities(mapply(function(f,y) function(x) f(x)*y,x,y),
                   basis=basis)
}

#' @describeIn Mul BayesSpace scaling multiplication
#' @export
#' @method Mul BDataPDFs
Mul.BDataPDFs = function(x,y) {
    #if( is.numeric(x) )
    #    return(Mul(y,x))
  if( !is.numeric(y) ) stop("Mul.BDataPDFs works only with numeric `y`") # return(Mul(y,x))
    as.BDataPDFs(Mul(as.BDataDensities(x),y))
}

#' @describeIn Mul BayesSpace scaling multiplication for log-PDFs
#' @export
#' @method Mul BDataLogPDFs
Mul.BDataLogPDFs = function(x,y) {
#    if( is.numeric(x) )
#        return(Mul(y,x))
  if( !is.numeric(y) ) stop("Mul.BDataLogPDFs works only with numeric `y`") # return(Mul(y,x))
  as.BDataLogPDFs(Mul(as.BDataLogDensities(x),y))
}

#' @describeIn Mul BayesSpace scaling multiplication
#' @export
#' @method Mul BDataSplines
Mul.BDataSplines = function(x,y) {
    # NOTE: do not test is.numeric(x) here -- a BDataSplines object is a
    # numeric matrix, so that guard would always fire and bounce between
    # Mul.BDataSplines and Mul.numeric.BData in infinite mutual recursion
    if( !is.numeric(y) ) stop("Mul.BDataSplines works only with numeric `y`") # return(Mul(y,x))
    lenY=length(y)
    if( lenY != nrow(x) && lenY!=1)
        stop("length of Y does not match nrow(x) in Mul.BDataSplines",lenY,"!=",nrow(x))
    BDataSplines(unclass(x)*y,basis=getBayesSpaceBasis(x))
}



#' @describeIn norm.BData generic functionality computing Euclidean norms of
#' vectors (default method) and redirecting to `base::norm()` for matrices.
#' @export
norm <- function(x,...) UseMethod("norm",x)

#' @method norm default
#' @export
norm.default <- function(x,...) {
  sqrt( sum(x^2) )
}

#' @method norm matrix
#' @export
norm.matrix <- function(x,...){
  base::norm(x=x, ...)
}

#' Computes the Bayes-Space norm of BData objects
#' @param x a BData object
#' @param measure the reference measure to on which to compute the L2 norm of the clr
#' @param ... further arguments to expectation.
#' @returns A numeric vector of the L2 (Bayes-space) norms, one per observation.
#' @importFrom compositions norm
#' @export
#' @method norm BData
#' @examples
#' lbdl = BDataLogDensities(basis=B2N01Basis5, closed=TRUE,
#'   list(
#'     \(x)(dnorm(x,2,1,log=TRUE)),
#'     \(x)(dnorm(x,1,1,log=TRUE)),
#'     \(x)(dnorm(x,-1,1,log=TRUE)),
#'     \(x)(dnorm(x,-2,1,log=TRUE))
#'     ))
#' norm(lbdl)
norm.BData = function(x,...,measure=getReferenceMeasure(x)) {
    sqrt(sapply(clr(x), function(f,...) expectation(measure,function(x) f(x)^2,...),...))
}

#' @describeIn norm.BData fast exact norm for BDataSplines via the Gram
#' matrix of the basis: the squared norm of each observation is the
#' quadratic form of its coefficient vector with the metric tensor S.
#' @export
#' @method norm BDataSplines
norm.BDataSplines = function(x,...) {
    S = getMetricTensor(getBayesSpaceBasis(x))
    coefs = unclass(x)
    if( !is.matrix(coefs) )
        coefs = matrix(coefs,nrow=1)
    sqrt(pmax(rowSums((coefs %*% S) * coefs),0))
}



### distance ===========



#' Distance method for BData objects
#'
#' Computes inter-distances between all objects of `x`
#'
#' @param x a data container of class `BData` or derived
#' @param ... further parameters to `stats::dist()`, i.e. among other `method`
#' allows to use Manhattan or other Minkowski norms.
#'
#' @returns a distance object, as `stats::dist()`. The implementations
#' rely on transformation to BDataSplines; the spline coefficients are
#' first mapped to orthonormal coordinates via the Gram matrix of the
#' basis, so that the default `method = "euclidean"` yields the exact
#' Bayes-space (clr-L2) distance also for non-orthonormal bases such as
#' ZB-splines. Other `method`s (Manhattan, Minkowski, ...) operate on
#' these orthonormal coordinates and depend on the chosen basis.
#' @export
#' @method dist BData
#' @importFrom compositions dist
#' @examples
#' lbdl = BDataLogDensities(basis=B2N01Basis5, closed=TRUE,
#'   list(
#'     \(x)(dnorm(x,2,1,log=TRUE)),
#'     \(x)(dnorm(x,1,1,log=TRUE)),
#'     \(x)(dnorm(x,-1,1,log=TRUE)),
#'     \(x)(dnorm(x,-2,1,log=TRUE))
#'     ))
#' dist(lbdl)
dist.BData = function(x, ...){
  dist(as.BDataLogDensities(x),...)
}


#' @describeIn dist.BData calculations for BDataLogDensities objects
#' @export
#' @method dist BDataLogDensities
dist.BDataLogDensities = function(x, ...){
  dist(as.BDataSplines(x),...)
}


#' @describeIn dist.BData calculations for BDataSplines objects
#' @export
#' @method dist BDataSplines
dist.BDataSplines = function(x, ...){
  coefs <- unclass(idt(x))
  S <- getMetricTensor(getBayesSpaceBasis(x))
  if (!isTRUE(all.equal(S, diag(nrow(S)), tolerance = .Machine$double.eps^0.5))) {
    # map coefficients to orthonormal coordinates: with S = M %*% t(M),
    # (c1-c2) %*% S %*% t(c1-c2) equals the squared Euclidean norm of (c1-c2) %*% M
    M <- tryCatch(t(chol(S)), error = function(e) {
      eig <- eigen(S, symmetric = TRUE)
      if (any(eig$values < -sqrt(.Machine$double.eps)))
        stop("Gram matrix of the basis is not positive semi-definite.")
      eig$vectors %*% diag(sqrt(pmax(eig$values, 0)))
    })
    coefs <- coefs %*% M
  }
  dist(coefs, ...)
}