version 1.0-0

DESCRIPTION

@@ -0,0 +1,16 @@
+Package: R2Cuba
+Title: Multidimensional Numerical Integration
+Version: 1.0-0
+Date: 2010-01-01
+Author: The Cuba library has been written by Thomas Hahn
+(http: //wwwth.mppmu.mpg.de/members/hahn); Interface to R was written
+        by Annie Bouvier and Ki�n Ki�u
+Maintainer: Annie Bouvier <Annie.Bouvier@jouy.inra.fr>
+Depends: methods
+Description: R2Cuba implements four general-purpose multidimensional
+        integration algorithms: Vegas, Suave, Divonne and Cuhre
+License: GPL (>= 3)
+Encoding: latin1
+Packaged: Tue Jan 5 15:11:01 2010; abouvier
+Repository: CRAN
+Date/Publication: 2010-01-06 15:31:51

NAMESPACE

@@ -0,0 +1,4 @@
+# Les points d'entree C ou Fortran appel�s par les fonctions R;
+useDynLib(R2Cuba, Rcuhre, Rdivonne, Rsuave, Rvegas)
+# Les fonctions R a disposition de l'utilisateur
+export(cuhre, divonne, suave, vegas, print.cuba)

R/commoncuba.R

@@ -0,0 +1,137 @@
+# ++++++++++++++++++++++++++++++++++++++++++
+# Functions shared by all the R2Cuba functions
+# ++++++++++++++++++++++++++++++++++++++++++
+verif <- function(ndim, lower, upper, rel.tol, abs.tol,
+            flags, min.eval,  max.eval) {
+  # Verification of the input arguments
+  # Issue a warning for each error and return T if none
+  bon <- TRUE
+
+  if (ndim<=0) {
+    bon <- FALSE
+      warning("ndim should be positive")
+  }
+
+  if ((length(lower) != ndim) ||
+      (length(upper) != ndim)) {
+       bon <- FALSE
+      warning("lower and upper should be vectors of length ndim")
+  }
+
+    if (any(lower >= upper)) {
+      bon <- FALSE
+      warning("Lower bounds should be less than upper bounds")
+    }
+  if ( !is.null(flags$verbose) &&
+      ((flags$verbose <0) || (flags$verbose >3))) {
+     bon <- FALSE
+      warning("flags$verbose should be in [0,3]")
+   }
+
+  if ( !is.null(flags$final) &&
+       ((flags$final <0) || (flags$final >1))) {
+     bon <- FALSE
+      warning("flags$final should be in [0,1]")
+   }
+
+      if ( !is.null(flags$pseudo.random) &&
+       ((flags$pseudo.random <0) || (flags$pseudo.random >1))) {
+     bon <- FALSE
+      warning("flags$pseudo.random should be in [0,1]")
+   }
+      if ( !is.null(flags$smooth) &&
+       ((flags$smooth <0) || (flags$smooth >1))) {
+     bon <- FALSE
+      warning("flags$smooth should be in [0,1]")
+   }
+
+  if ( (min.eval<0) || (min.eval>max.eval)) {
+    bon <- FALSE
+      warning("Error in min.eval or max.eval")
+   }
+
+     return(bon)
+  } # End verif
+# -------------------------------------------
+decodflags <- function(flags) {
+# Decode the flags
+    if (!is.null(flags$verbose))
+      lesflags <- flags$verbose
+    else
+      lesflags <- 1 # valeur par défaut
+    if (!is.null(flags$final))
+      lesflags <- lesflags + (flags$final*4)
+    if (!is.null(flags$pseudo.random))
+      lesflags <- lesflags + (flags$pseudo.random*8)
+    if (!is.null(flags$smooth))
+      lesflags <- lesflags + (flags$smooth*16)
+return(lesflags)
+  } # End decodflags
+
+# -------------------------------------------------------
+
+crff <- function(lecall, integrand, nomf, libargs, ...) {
+  # Determine how to call the user function according to
+  # the list of its arguments and the current list of arguments
+
+ #nfarg: number of the R user function formals arguments
+ nfarg <- length(formals(integrand))
+
+ if (nfarg <1)
+   stop("Function integrand should have one argument at least")
+
+ # nargsup:  number of additional arguments in the current call
+ a <- deparse(lecall, width=500)
+ # Mettre les arguments courants dans une liste
+ zl <- NULL
+ eval(parse(text= sub(nomf, "zl<-list", a)))
+# determiner ceux additionnels
+ a <- (names(zl) %in% libargs)
+ nargsup <- length(a[a==FALSE])
+
+ if (nargsup >0) {
+   # call with additional arguments
+   if (nfarg == (nargsup+1))
+     ffintegrand <- function(x, phw=0) integrand(x, ...)
+   else
+     if (nfarg >= (nargsup+2)) {
+       ffintegrand <- function(x, phw=0) integrand(x,phw, ...)
+     }
+   else 
+     stop(paste("Additional argument", names(zl)[a==FALSE], "not expected in the integrand function\n"))
+
+ } # End (nargsup >0)
+ else {
+  if (nfarg == 1)
+     ffintegrand <- function(x, phw=0) integrand(x)
+   else {
+      ffintegrand <- function(x, phw=0) integrand(x,phw)
+   }
+ }
+ return(ffintegrand)
+} # End crff
+
+# ++++++++++++++++++++++++++++++++++++++++++
+# The methods of the class "cuba"
+# ++++++++++++++++++++++++++++++++++++++++++
+
+print.cuba <- 
+          function(x,  ...)
+  {
+                for (i in 1:length(x$value)) {
+       cat("integral: ", format(x$value[i], ...),
+           " (+-",
+            format(x$abs.error[i], digits = 2), ")\n", sep = "")
+            if (!is.null(x$nregions))
+              cat("nregions: ", x$nregions, "; ",sep="")
+     }
+            cat("number of evaluations: ", x$neval)
+             cat("; probability: ", format(x$prob, ...), "\n")
+
+              if (x$message !="OK") 
+         cat("failed with message ", sQuote(x$message), "\n", 
+        sep = "")
+    invisible(x)
+
+          }
+

R/cuhre.R

@@ -0,0 +1,80 @@
+# -------------------------------------------------------
+
+cuhre <- function(ndim, ncomp,
+                 integrand, ...,
+                 lower=rep(0,ndim), upper=rep(1,ndim),
+                  rel.tol= 1e-3,
+                  abs.tol = 0,
+                  flags=list(verbose=1,
+                    final=1, pseudo.random=0,
+                    mersenne.seed=NULL),
+                  min.eval=0,  max.eval=50000,
+                    key=0)
+
+  {
+
+     # Verification
+ if (!verif(ndim, lower, upper, rel.tol, abs.tol,
+            flags, min.eval,  max.eval))
+      stop("Error in input: see the warnings")
+ if (is.null(flags$mersenne.seed))
+   flags$mersenne.seed <- NA
+
+  # Decode the flags
+    lesflags <- decodflags(flags)
+    # Allocate output structures
+    nregions <- 0
+    neval <- 0
+    fail <- 0
+    integral <- rep(0, ncomp)
+    error <- rep(0, ncomp)
+    prob <- rep(0, ncomp)
+
+# Determine how to call the user function according to
+  # the list of its arguments and the current list of arguments
+ libargs <- c("ndim", "ncomp", 
+                 "integrand","lower", "upper",
+                  "rel.tol", "abs.tol", "flags",
+                  "min.eval","max.eval", "key")
+ ffintegrand <- crff(match.call(), integrand, "cuhre", libargs, ...)
+
+  prdbounds <- prod(upper-lower)
+    ret <-  .C("Rcuhre", as.integer(ndim),
+             as.integer(ncomp),
+             ffintegrand, new.env(),
+             as.double(lower), as.double(upper),as.double(prdbounds),
+             as.double(rel.tol), as.double(abs.tol),
+             as.integer(flags$mersenne.seed),
+               as.integer(lesflags),
+             as.integer(min.eval),  as.integer(max.eval),
+             as.integer(key),
+             nregions=as.integer(nregions),
+             neval=as.integer(neval), fail=as.integer(fail),
+             integral=as.double(integral),
+             error=as.double(error), prob=as.double(prob),
+               NAOK=TRUE)
+#Add to finish the last print:
+cat("\n")
+
+# To homogeneize with the R function "integrate", add
+    # message and call into the output,
+    # ifail rather than fail , abs.error rather than  error,
+    # value rather than integral
+
+
+    if (ret$fail ==0)
+      mess ="OK" # OK required to be printed by print.cuba
+    else
+      mess="Dimension out of range"
+
+    retour <- list(neval=ret$neval, nregions=ret$nregions,
+                   ifail=ret$fail, value=ret$integral,
+                abs.error=ret$error,
+                   neval=ret$neval,
+                   prob=ret$prob, message=mess,
+                call=match.call())
+    attr(retour, "class") = c("cuba")
+    return(retour)
+  } # End cuhre
+
+

R/divonne.R

@@ -0,0 +1,113 @@
+divonne <- function(ndim, ncomp,
+                   integrand, ...,
+                 lower=rep(0,ndim), upper=rep(1,ndim), 
+                  rel.tol= 1e-3,
+                  abs.tol = 0,
+                  flags=list(verbose=1,
+                    final=1, pseudo.random=0,
+                    mersenne.seed=NULL),
+                    min.eval=0,  max.eval=50000,
+                    key1=47,
+                    key2=1, key3=1, 
+                    max.pass=5, border=0, max.chisq=10,
+                    min.deviation=0.25,
+                    xgiven=NULL,  nextra=0,
+                    peakfinder=NULL)
+
+  {
+     # Verification
+ if (!verif(ndim, lower, upper, rel.tol, abs.tol,
+            flags, min.eval,  max.eval))
+      stop("Error in input: see the warnings")
+     # Decode the flags
+    lesflags <- decodflags(flags)
+if (is.null(flags$mersenne.seed))
+  flags$mersenne.seed <- NA
+
+ #  xgiven dimensions:
+ if (!is.null(xgiven)) {
+   if (!is.matrix(xgiven))
+     stop("xgiven should be a matrix")
+   ngiven <- ncol(xgiven)
+   ldxgiven <- nrow(xgiven)
+   if (ldxgiven != ndim)
+     stop("Matrix xgiven should have ndim rows")
+   # Rescaler xgiven dans l'hypercube unité
+   for (i in 1:ndim) {
+     xgiven[i,] <- xgiven[i,]/(upper[i]-lower[i])
+   }
+ }
+ else {
+    ngiven=0; ldxgiven=ndim
+  }
+      if (nextra <0)
+       stop("nextra should be positive or zero")
+
+ if ((nextra >0) && is.null(peakfinder))
+   stop("nextra positive but not peakfinder subroutine")
+ if ((nextra==0) && !is.null(peakfinder))
+   warning("peakfinder ignored: nextra is zero")
+
+
+
+       # Allocate outputs
+    nregions <- 0
+    neval <- 0
+    fail <- 0
+    integral <- rep(0, ncomp)
+    error <- rep(0, ncomp)
+    prob <- rep(0, ncomp)
+  libargs <- c("ndim", "ncomp", 
+                 "integrand","lower", "upper",
+                  "rel.tol", "abs.tol", "flags",
+                    "min.eval","max.eval", "key1",
+                    "key2", "key3", 
+                    "max.pass", "border", "max.chisq",
+                    "min.deviation", "xgiven",  "nextra",
+                    "peakfinder")
+
+ ffintegrand <- crff(match.call(), integrand, "divonne", libargs, ...)
+  prdbounds <- prod(upper-lower)
+  ret <-  .C("Rdivonne", as.integer(ndim),
+             as.integer(ncomp),
+             ffintegrand, new.env(),
+             as.double(lower), as.double(upper),as.double(prdbounds),
+             as.double(rel.tol), as.double(abs.tol),
+             as.integer(flags$mersenne.seed),
+             as.integer(lesflags),
+             as.integer(min.eval),  as.integer(max.eval),
+             as.integer(key1), as.integer(key2), as.integer(key3),
+             as.integer(max.pass), as.double(border),
+             as.double(max.chisq), as.double(min.deviation),
+             as.integer(ngiven), as.integer(ldxgiven),
+             as.double(xgiven), as.integer(nextra),
+             peakfinder,
+             nregions=as.integer(nregions),
+             neval=as.integer(neval), fail=as.integer(fail),
+             integral=as.double(integral),
+             error=as.double(error), prob=as.double(prob),
+             NAOK=TRUE)
+#Add to finish the last print:
+cat("\n")
+
+# To homogeneize with the R function "integrate", add
+    # message and call into the output,
+    # ifail rather than fail , abs.error rather than  error,
+    # value rather than integral    
+
+    if (ret$fail ==0)
+      mess ="OK" # OK required to be printed by print.cuba
+    else
+      mess="Dimension out of range"
+
+    retour <- list(neval=ret$neval, nregions=ret$nregions,
+                   ifail=ret$fail, value=ret$integral,
+                abs.error=ret$error,
+                   neval=ret$neval,
+                   prob=ret$prob, message=mess,
+                call=match.call())
+    attr(retour, "class") = c("cuba")
+    return(retour)
+  } # End divonne
+
+