R/Cartesian.R

In DistributionalDDA: Distributional Data Analysis in Bayes Spaces

#' Compute the cartesian product
#'
#' @param ... objects of the same kind (space, measure, support, basis, ...) to be used for a cartesian product, and named parameters of the method functions. The first will determine the generic behaviour. if exactly one argument is given and is a pure list it is interpreted as the list of objects given. The individual spaces can be named.
#' @family Cartesian
#' @return the Cartesian product of the objects given. Check [SupportCartesianProduct()]
#' for details about this structure.
#'
#' The cartesian product of a space (like a vector space or support) is the space of tuples, where each tuple element comes from the space of the same position or name.
#' The cartesian product of measures is the product measure, which gives a cartesian product of sets the product of the measures.
#' The cartesian product of functions(like basis) is a function on the cartesian product of supports and gives the product of the function values for each of the components.
#' @export
#'
#'
#'
#'
cartesianProduct <- function(...) {
    .LIST=list(...)
    l = length(.LIST)
    if( l==0 )
        return(NULL)
    x=.LIST[[1]]
    if( l==1 && is(x,"list") )
        return(do.call("cartesianProduct",x))
    UseMethod("cartesianProduct",x)
}
#' @describeIn cartesianProduct method for Supports
#' @export
cartesianProduct.Support <- function(...) SupportCartesianProduct(...)
#' @describeIn cartesianProduct method for ReferenceMeasure
#' @export
cartesianProduct.ReferenceMeasure <- function(...) ReferenceMeasureCartesianProduct(...)
#' @describeIn cartesianProduct method for BayesSpace
#' @export
cartesianProduct.BayesSpace <- function(...) BayesSpaceCartesianProduct(...)
#' @describeIn cartesianProduct method for BayesSpaceBasis
#' @export
cartesianProduct.BayesSpaceBasis <- function(...) BayesSpaceBasisCartesianProduct(...)

#' Decomposes a cartesian product into a list of its factors.
#'
#' @param ... the cartesian product or products to be decomposed
#'
#' @return If x is not a cartesian product the function should return
#' list(x). if x is a cartesian product it should provide its factors.
#' Check [SupportCartesianProduct()] for details about this latter
#' structure.
#'
#' @family Cartesian
#'
#' @export
getCartesianFactors <- function(...) {
    L = list(...)
    if(length(L)==1)
        UseMethod("getCartesianFactors",L[[1]])
    else
        unlist(lapply(L,getCartesianFactors))
}
#' @export
getCartesianFactors.default <- function(...) list(...)
#' @export
getCartesianFactors.SupportProduct <-
    function(...) structure(...,class=NULL)
#' @export
getCartesianFactors.ReferenceMeasureCartesianProduct <-
    function(...) structure(...,class=NULL)
#' @export
getCartesianFactors.BayesSpaceCartesianProduct <-
    function(...) structure(...,class=NULL)
#' @export
getCartesianFactors.BayesSpaceBasisCartesianProduct <-
    function(...) attr(...,which="marginalBases")

.listCatcher <- function(...) {
    L = list(...)
    if( length(L) == 1 && is(L[[1]],"list") )
        L[[1]]
    else
        L
}

#' Creates a cartesian product of supports
#'
#' @param ... Give the supports, either directly or as a simple list. If they correspond to named columns in a dataframe, name them
#' @return A Support of type SupportCartesianProduct that provides the Cartesian product of the supports
#'
#' Supports are essentially encoded sets or spaces of possible values.
#' The cartesian product of such spaces is the space containing the tuples of
#' with one entry from each space. In R it will be encoded by the columns of
#' a data.frame.
#'
#' @export
#'
#' @seealso [Support], [data.frame]
#' @family Cartesian
#'
#' @examples
#'  s1 = Support(-Inf,Inf,lim=c(-2,2))
#'  s2 = Support(0,2)
#'  s3 = Support(X=c(1,2,3))
#'  s4 = Support(X=c("A","B","C","D"))
#'  s  = SupportCartesianProduct(s1,s2,s3,s4)
#'
SupportCartesianProduct = function(...) {
    .LIST = .listCatcher(...)
    .call=as.call(c(list(as.name("SupportCartesianProduct")),
                    .LIST %..% call))
    structure(.LIST,
              class=c("SupportCartesianProduct","Support"),
              call = .call)
}

## ReferenceMeasureCartesianProduct
#' Compute a cartesian Product of measures
#'
#' @param ... the measures, possibly named
#' @param nSamples the (maximum) number of samples to be used to represent
#'        the measure in the fixed sample for Monte Carlo integrations.
#'
#' @return A ReferenceMeasureCartesianProduct, which is essentially a
#'     list of the measures, in an object of that class, that acts like a
#'     ReferenceMeasure for the Cartesian Product of the supports of the
#'     individual measures.
#'
#' The product measure corresponds gives a any product set the measure
#' of the that corresponds to the product of the individual
#' measures. In case of a probablity distribution, that would be the
#' measure representing an independent distribution, with the initial
#' measures as marginal distributions. This concept is used to
#' establish a reference measure for multivariate BayesSpaces which
#' can also contain non independent measures.
#'
#' @export
#' @family Cartesian
#'
#' @examples
#'  RMCP3 = ReferenceMeasureCartesianProduct(A=StdReferenceMeasureN01,
#'                                           B=StdReferenceMeasureUnif01,
#'                                           C=ReferenceMeasureDiscrete(c("A","B","C"),c(0.2,0.3,0.5))
#'                                            )
#'
ReferenceMeasureCartesianProduct = function(...,nSamples=100000) {
    .LIST = .listCatcher(...)
    structure(.LIST,
              class=c("ReferenceMeasureCartesianProduct","ReferenceMeasure"),
              sample=.crossSampler(.LIST=lapply(.LIST,getSample),
                                   nSamples=nSamples),
              call = as.call(c(list(as.name("ReferenceMeasureCartesianProduct")),
                               lapply(.LIST,attr,which="call"))
                             )
              )
}


#' The Cartesian Product of Bayes Spaces
#'
#' @param ... the Bayes Spaces to be used in the product either given directly or as a simple list. Names will to used to name the variable they are operating on.
#'
#' @return the Product space of the BayesSpaces, a BayesSpace itself.
#' @family Cartesian
#'
#' @export
BayesSpaceCartesianProduct = function(...) {
    structure(.listCatcher(...),
              class=c("BayesSpaceCartesianProduct","BayesSpace"),
              call = match.call())
}
#' @describeIn BayesSpaceCartesianProduct Product of L2 spaces
#' @export
L2SpaceCartesianProduct = function(...) {
    structure(.listCatcher(...),
              class=c("L2SpaceCartesianProduct","BayesSpace"),
              call = match.call())
}


#' @export
getReferenceMeasure.BayesSpaceCartesianProduct = function(x) {
    . <- NULL # Trick to avoid warning
   ReferenceMeasureCartesianProduct(x %.% getReferenceMeasure(.))
}

#' Construct the product Bayes Space for a product reference measure
#'
#' @param X a ReferenceMeasureCartesianProduct to construct the
#'     BayesSpace for
#' @param ... further arguments to the getBayesSpace of the individual
#'     spaces
#' @return A BayesSpaceCartesianProduct corresponding to the cartesian
#'     Product of Bayes measures, itself also a Bayes measure.
#'
#' @export
#' @method getBayesSpace ReferenceMeasureCartesianProduct
#' @family Cartesian
#'
getBayesSpace.ReferenceMeasureCartesianProduct = function(X,...) {
    structure(
        BayesSpaceCartesianProduct(lapply(X,getBayesSpace,...)),
        call = match.call())
}



#' Computes the product basis of BayesSpaceBasis for the product BayesSpace
#'
#' @param ... the BayesSpacesBases to operate, either directly or as a simple list.
#' @param .LIST the list of these bases, use either ... or .LIST
#' @param bayesSpace the Product BayesSpace for the basis. The product
#'     space is computed automatically, if not given. Try to create
#'     this space only once.
#'
#'
#'
#' @return A Bayes Space Basis for the product space as defined the DDA Book.
#' @family Cartesian
#'
#' @seealso [BayesSpaceBasis]
#'
#'
BayesSpaceBasisCartesianProduct = function(...,
                                           .LIST=.listCatcher(...),
                                           bayesSpace=BayesSpaceCartesianProduct(.LIST=lapply(.LIST,getBayesSpace))) {
    f = function(x) {
        A1=.cartesianDimensionIndices(unclass(.LIST))
        A2=lapply(A1,\(i) cbind(1,basisEval(unclass(.LIST)[[i]],x[[i]])))
        .outerBasisProduct(A2)[,-1,drop=FALSE]
        #.LIST |>
        #    .cartesianDimensionIndices()  |>
        #    lapply(\(i) cbind(1,basisEval(.LIST[[i]],x[[i]]))) |>
        #    .outerBasisProduct()
    }
    S <- NULL # Trick to avoid warning
    basisDims  = sapply(.LIST,function(x) nrow(x$S))
    structure(BayesSpaceBasis(
        f,
        bayesSpace=bayesSpace,
        A=diag(prod(basisDims+1)-1),
        S=.LIST %.% S |> .SProduct()),
        marginalBases=.LIST,
        call=match.call()
        )

}

#' Create a BayesSpaceBasis for a Cartesian Product of Bayes Spaces
#'
#' @param basis A cartesian product of Bayes Spaces (BayesSpaceCartesianProduct)
#' @param ... further parameter controling the BayesSpaceBasis generation, used by the BayesSpaces basis generation tools of the factors of the cartesian product.
#' @param useBases a list of marginal BayesSpaceBases corresponding to the BayesSpaces.
#'
#' @return A BayesSpaceBasisCartesianProduct of the bases in useBases, for the BayesSpace basis.
#' @family Cartesian
#'
#' @seealso BayesSpaceBasis
#' @method BayesSpaceBasis BayesSpaceCartesianProduct
#' @export
#' @examples
#' BayesSpaceBasis(BayesSpaceCartesianProduct(
#'     A=BayesSpace(StdReferenceMeasureUnif01),
#'     B=BayesSpace(ReferenceMeasureDiscrete(c("A","B","C"),c(0.2,0.3,0.5))),
#'     C=BayesSpace(StdReferenceMeasureN01)
#'     )
#' )
#'
BayesSpaceBasis.BayesSpaceCartesianProduct <-
    function(basis,...,
             useBases=lapply(basis,BayesSpaceBasis,...)) {
        structure(
            BayesSpaceBasisCartesianProduct(useBases,bayesSpace=basis),
            call=match.call())
    }
#' @describeIn BayesSpaceBasis.BayesSpaceCartesianProduct creates the BayesSpaceCartesian Product for this reference measure and the corresponding product basis.
#' @param bayesSpace the BayesSpace to be used for this CartesianProduct. It has a good default, but should be given in case you want to create the basis for an already existing ProductSpace.
#' @export
#' @method BayesSpaceBasis ReferenceMeasureCartesianProduct
BayesSpaceBasis.ReferenceMeasureCartesianProduct <-
    function(basis,...,
             useBases=lapply(bayesSpace,BayesSpaceBasis,...),
             bayesSpace=BayesSpaceCartesianProduct(lapply(basis,BayesSpace))) {
        structure(
            BayesSpaceBasisCartesianProduct(useBases,bayesSpace=bayesSpace),
            call=match.call()
            )
    }




## Indices to the right columns of a dataset.
.cartesianDimensionIndices <- function(X) {
   NAMES <- names(X)
   if(is.null(NAMES)) {
      1:length(X)
   } else {
      NAMES
    }
}

## Takes a list of boolean vectors can makes a rowwise logical and.
.parallelAnd <- function(...) {
    .LIST=.listCatcher(...)
    if( length(.LIST) > 0 )
        .LIST[[1]] & .parallelAnd(.LIST[-1])
    else
        TRUE
}

#' @describeIn isInDomain method for product measures
#' @export
#' @method isInDomain ReferenceMeasureCartesianProduct
isInDomain.ReferenceMeasureCartesianProduct <-
  function(x,referenceMeasure) {
    .parallelAnd(.LIST=lapply(.cartesianDimensionIndices(referenceMeasure),
                              \(j) isInDomain(x[,j],referenceMeasure[j])))
  }

## Create cross samples from marginal samples.
.fullCrossSampler = function(...,.LIST=list(...)) {
    if(length(.LIST)==0)
        return(data.frame())
    if( is.null(names(.LIST)) )
        names(.LIST)=paste("t",1:length(.LIST),sep="")
    Y <- .fullCrossSampler(.LIST=.LIST[-1])
    nY=nrow(Y)
    if(hasWeights(Y) || hasWeights(.LIST[[1]])) {
        wX = getWeights(.LIST[1])
        X <- as.data.frame(structure(.LIST[[1]],weights=NULL))
        if(length(.LIST)==1)
            return(structure(X,weights=wX))
        nX=nrow(X)
        wY=getWeights(Y)
        i = rep(1:nX,nY)
        j = rep(1:nY,each=nX)
        structure(cbind(structure(X[i,],weights=NULL),
                        structure(Y[j,],weights=NULL)),
                  weights=wX[i]*wY[j])
    } else {
        X <- as.data.frame(.LIST[1])
        if(length(.LIST)==1)
            return(X)
        nX=nrow(X)
        i = rep(1:nX,nY)
        j = rep(1:nY,each=nX)
        cbind(structure(X[i,]),structure(Y[j,]))
    }
}

.randomCrossSampler = function(...,.LIST=list(...),nSamples=1000000) {
    if( is.null(names(.LIST)) )
        names(.LIST) <- paste("t",1:length(.LIST),sep="")
    as.data.frame(lapply(.LIST,function(X) {
        i = sample(1:nCases(X),nSamples,replace=TRUE,getWeights(X))
        selectCases(X,i)
    }))
}

.subSampleTo <- function(s,n) {
    N = length(s)

}

# In development
#.systematicCrossSampler1 <- function(...,.LIST=list(...),nSamples=1000000) {
#    ns = sapply(.LIST,nCases)
#    reduce = ns<30
#    nns = length(nns)
#    undersamplingFactor=ceiling((prod(ns)/nSamples)^(1/nns))
#    newNs = floor(ns / undersamplingFactor)
#    .crossSampler(.LIST=mapply(.subSampleTo,.LIST,newNs),nSamples=nSamples)
#}

.crossSampler <- function(...,.LIST=list(...),nSamples=1000000) {
    sizeOfFullCrossSample = prod(sapply(.LIST,nCases))
    if(  sizeOfFullCrossSample <= nSamples ) {
        .fullCrossSampler(.LIST=.LIST)
    } else {
        .randomCrossSampler(.LIST=.LIST,nSamples=nSamples)
    }
}

#' @describeIn selectIntegrationMethod select "MonteCarlo" if the getSample is an exhausted random default sample. Otherwise use "sample"
#' @export
selectIntegrationMethod.CartesianProductMeasure = function(measure,method) {
    ifelse(nCases(getSample(measure))==1000000,"MonteCarlo","sample")
}


## OuterBasisProducts

.outerBasisProduct <- function(bases) {
    if( length(bases) == 1 )
        return(bases[[1]])
    Y = .outerBasisProduct(bases[-1])
    X = bases[[1]]
    mX=ncol(X)
    mY=ncol(Y)
    (X[,rep(1:mX,mY),drop=FALSE]*Y[,rep(1:mY,each=mX),drop=FALSE])
}
#' Computes the indices of a certain marginal basis
#'
#' @param basis the full cartesian BayesSpaceBasis
#' @param indices giving which the marginal subspaces of interest, given as booelan vector or vector of names of positive or negative integers .
#' @param basisDims the dimensions of the basises for the various subspaces
#'
#' @return an integer vector with the indices of the interaction described
#'         by the indices.
#' @export
#' @family Cartesian
subBasisIndices <- function(basis,indices,
                            basisDims  = sapply(attr(basis,".LIST"),
                                                function(x) nrow(x$S))) {
    if( !is.logical(indices) ) {
        boolIndices = sapply(basisDims,\(x) FALSE)
        boolIndices[indices]=TRUE
        indices=boolIndices
    }
    stopifnot(length(indices)==length(basisDims))
    .subBasisIndicesGetter(basisDims,indices)[-1] |> which()
}

## Helplerfunktion for computing the indices through indicators.
.subBasisIndicesGetter = function(basisDims,indices) {
    stopifnot(length(basisDims)==length(indices))
    if(length(indices)>1) {
        kronecker(.subBasisIndicesGetter(basisDims[-1],indices[-1]),
                  c(!indices[1],rep(indices[1],basisDims[1])))
    } else {
        c(!indices,rep(indices,basisDims))
    }
}

## Helper functions to compute the Scalar products of the tensor basis
### Compute the block matrix [[1 0] [0 S]]
.SExtend <- function(S) rbind(c(1,rep(0,NCOL(S))),cbind(0,S))

### Compute the joint scalar product for the tensor product basis by
### computing the unconstrained and removing the constant
.SProduct <- function(L) .SProductIntern(L)[-1,-1]


### Compute the unconstraint product scalar product through
### Recursive Kroneker products of the extended matrices
.SProductIntern <- function(L) {
    if( length(L) > 1  ) {
        kronecker(.SProductIntern(L[-1]),.SExtend(L[[1]]))
    } else {
        .SExtend(L[[1]])
    }
}

#' @describeIn Support returns the product support
#' @export
#' @seealso [ReferenceMeasureCartesianProduct], [SupportCartesianProduct]
#' @method getSupport ReferenceMeasureCartesianProduct
getSupport.ReferenceMeasureCartesianProduct <- function(x,...) {
    SupportCartesianProduct(.LIST=lapply(x,getSupport,...))
}

.listTypeChecker <- function(class2,...) {
    L = list(...)
    if(length(L)==1 && is(L[[1]],"list") && !is(L[[1]],class2))
        L = L[[1]]
    fits = sapply(L,is,class2)
    if( !all(fits) ) {
        if( is.null(names(L)) )
            n=(1:length(L))[fits]
        else
            n=1:names(L)[fits]
        stop("Expected ",class2," but got wrong classes at ",paste(n)," with classes ",deparse(L[!fits] %..% class))
    }
    L
}