## Copyright (C) 2010 Marius Hofert and Martin Maechler
##
## This program is free software; you can redistribute it and/or modify it under
## the terms of the GNU General Public License as published by the Free Software
## Foundation; either version 3 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
## FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
## details.
##
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.

#### Implementation of function evaluations and random number generation for
#### nested Archimedean copulas

##' Returns the copula value at a certain vector u
##' @param x nacopula
##' @param u argument of the copula x
##' @return x(u)
##' @author Marius Hofert, Martin Maechler
## FIXME: maybe make this applicable to a matrix of u's?
pnacopula <- function(x,u) {
    stopifnot(is.numeric(u), 0 <= u, u <= 1,
	      length(u) >= dim(x))	# can be larger
    C <- x@copula
    th <- C@theta
    ## use u[j] for the direct components 'comp':
    C@psi(sum(unlist(lapply(u[x@comp], C@psiInv, theta=th)),
	      C@psiInv(unlist(lapply(x@childCops, pnacopula, u = u)),
		       theta=th)),
	  theta=th)
}

##' Compute the probability P[l < U <= u]  where U ~ copula x.
##' @param x outer nested archimedean copula
##' @param l d-vector of lower "integration" limits
##' @param u d-vector of upper "integration" limits
##' @return the probability that a random vector following the given copula
##'  falls in the hypercube with lower and upper corner l and u, respectively.
##' @author Marius Hofert, Martin Maechler
setGeneric("prob", function(x, l, u) standardGeneric("prob"))

setMethod("prob", signature(x ="outer_nacopula"),
          function(x, l,u) {
              d <- dim(x)
              ## TODO: maybe allow  l & u to be  k x d matrices
              ##        --> return vector of probabilities of length k
              stopifnot(is.numeric(l), is.numeric(u),
                        length(u) == d, d == length(l),
                        0 <= l, l <= u, u <= 1)
              if(d > 30)
                  stop("prob() for copula dimensions > 30 are not supported (yet)")
              D <- 2^d
              m <- 0:(D - 1)
              ## digitsBase() from package 'sfsmisc' {slightly simplified} :
              ## Purpose: Use binary representation of 0:N
              ## Author: Martin Maechler, Date:  Wed Dec  4 14:10:27 1991
              II <- matrix(0, nrow = D, ncol = d)
              for (i in d:1L) {
                  II[,i] <- m %% 2L + 1L
                  if (i > 1) m <- m %/% 2L
              }
              ## Sign: the ("u","u",...,"u") case has +1; = c(2,2,...,2)
              Sign <- c(1,-1)[1L + (- rowSums(II)) %% 2]
              U <- array(cbind(l,u)[cbind(c(col(II)), c(II))], dim = dim(II))
              sum(Sign * apply(U, 1, pnacopula, x=x))
          })

##' Returns (n x d)-matrix of random variates
##' @param n number of vectors of random variates to generate
##' @param x outer_nacopula
##' @return matrix of random variates
##' @author Marius Hofert, Martin Maechler
rnacopula <- function(n, x, ...)
{
    Cp <- x@copula			# outer copula
    theta <- Cp@theta			# theta for outer copula
    V0 <- Cp@V0(n,theta)		# generate V0's
    childL <- lapply(x@childCops, rnchild, # <-- start recursion
		     theta0=theta,V0=V0,...)
    dns <- length(x@comp)	 # dimension of the non-sectorial part
    r <- matrix(rexp(n*dns), n, dns) # generate the non-sectorial part
    ## put pieces together
    mat <- Cp@psi(r/V0, theta=theta)	# transform
    mat <- cbind(mat, do.call(cbind,lapply(childL, `[[`, "U")))
    ## get correct sorting order:
    j <- c(x@comp, unlist(lapply(childL, `[[`, "indCol")))
    ## extra checks TODO: comment
    stopifnot(length(j) == ncol(mat))
    m <- mat[,order(j)]			# permute data and return
    ## extra checks TODO: comment
    stopifnot(length(dm <- dim(m)) == 2, dm == dim(mat))
    m
}

##' Returns a list with an (n x d)-matrix of random variates and a vector of
##' indices.
##' @param x nacopula
##' @param n number of vectors of random variates to generate
##' @param theta0 parameter theta0
##' @param V0 vector of V0's
##' @return list(U = matrix(*,n,d), indCol = vector of length d)
##' @author Marius Hofert, Martin Maechler
rnchild <- function(x, theta0, V0,...)
{
    n <- length(V0)
    Cp <- x@copula # inner copula
    ## Consistency checks -- for now {comment later} :
    stopifnot(is(Cp, "acopula"), is.numeric(n), n == as.integer(n),
              is.numeric(V0), length(V0) == n, is.numeric(theta0))
    theta1 <- Cp@theta # theta_1 for inner copula
    ## generate V01's (only for one sector since the
    ## recursion in rnacopula() takes care of all sectors):
    V01 <- Cp@V01(V0, theta0,theta1,...)
    childL <- lapply(x@childCops, rnchild, # <-- recursion
                     theta0=theta1, V0=V01,...)
    dns <- length(x@comp)	# dimension of the non-sectorial part
    r <- matrix(rexp(n*dns), n, dns) # generate the non-sectorial part
    ## put pieces together: first own comp.s, then the children's :
    mat <- Cp@psi(r/V01, theta1) # transform
    mat <- cbind(mat, do.call(cbind, lapply(childL, `[[`, "U")))
    ## get correct sorting order:
    j <- c(x@comp, unlist(lapply(childL, `[[`, "indCol")))
    list(U = mat, indCol = j) # get list and return
}

if(FALSE) { # evaluate the following into your R session if you need debugging:
    trace(rnacopula, browser, exit=browser, signature=signature(x ="outer_nacopula"))

    debug(rnchild)
}

##' Constructor for outer_nacopula
##' @param family either character string (short or longer form of
##'	 copula family name) or an "acopula" family object
##' @param nacStructure a "formula" of the form C(th, comp, list(C(..), C(..)))
##' @return a valid outer_nacopula object
##' @author Martin Maechler
onacopula <- function(family, nacStructure) {
## , envir = ... , enclos=parent.frame() or
## , envir=environment()
##
### FIXME: base this on  onacopulaL() -- replacing nacStructure by nacList
### ----- : (1) replacing  C() with list() *AND* by
###         (2) wrapping the 3rd argument with list(.) if it's not already
### Use a *recursive* function like this one (but add the "list(.)" wrapping:
###    repC <- function(e) { sC <- as.symbol("C"); if(identical(e[[1]], sC)) e[[1]] <- as.symbol("list"); if(length(e) == 4) e[[4]] <- repC(e[[4]]); e}
    nacl <- substitute(nacStructure)
    stopifnot(identical(nacl[[1]], as.symbol("C")))
    COP <- getAcop(family)
    nacl[[1]] <- as.symbol("oC")
### does not work ..>>>>>>>>>>> we should use onacopulaL() inside functions! <<<<
    ## needed, e.g., when onacopula() is called from inside a user function:
    ## for(j in 2:3) if(is.language(nacl[[j]]))
    ##     nacl[[j]] <- eval(nacl[[j]], parent.frame())
    ## pframe <- parent.frame()
    mkC <- function(cClass, a,b,c) {
	if(missing(b) || length(b) == 0) b <- integer()
	if(missing(c) || length(c) == 0) c <- list()
	else if(length(c) == 1 && !is.list(c)) c <- list(c)
### does not work ...
	## else if(!is.list(c)) {
        ##     c <- eval(c, pframe) ; stopifnot(is.list(c))
        ## }
	## if(!is.numeric(a)) {
        ##     a <- eval(a, pframe) ; stopifnot(is.numeric(a), length(a) == 1)
        ## }
	## if(!is.numeric(b)) {
        ##     b <- eval(b, pframe) ; stopifnot(is.numeric(b))
        ## }
	else stopifnot(is.list(c))
	stopifnot(is.numeric(a), length(a) == 1, is.numeric(b))
	if(any(sapply(c, class) != "nacopula"))
	    stop("third entry of 'nacStructure' must be NULL or list of 'C(..)' terms")
	new(cClass, copula = setTheta(COP, a),
	    comp = as.integer(b), childCops = c)
    }
    C <- function(a,b,c) mkC("nacopula", a,b,c)
    oC <- function(a,b,c) mkC("outer_nacopula", a,b,c)
    eval(nacl)
}

##' @title Constructor for outer_nacopula - with list() input and *recursively*
##' @param family
##' @param nacList
##' @return
##' @author Martin Maechler
onacopulaL <- function(family, nacList) {
    COP <- getAcop(family)
    mkC <- function(cClass, abc) {
        stopifnot(is.list(abc), 2 <= (nL <- length(abc)))
        if(nL == 2) abc[[3]] <- list() else
        if(nL > 3) stop("nacLists must be of length 3 (or 2)")
	new(cClass, copula = setTheta(COP, abc[[1]]),
	    comp = as.integer(abc[[2]]),
            childCops = lapply(abc[[3]], function(.) mkC("nacopula", .)))
    }
    mkC("outer_nacopula", nacList)
}