R Under development (unstable) (2024-01-20 r85814 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> #### Tests for "Functionals":  Root finding, Optimization, Integration, etc
> 
> stopifnot(require("Rmpfr"))
Loading required package: Rmpfr
Loading required package: gmp

Attaching package: 'gmp'

The following objects are masked from 'package:base':

    %*%, apply, crossprod, matrix, tcrossprod

C code of R package 'Rmpfr': GMP using 64 bits per limb


Attaching package: 'Rmpfr'

The following object is masked from 'package:gmp':

    outer

The following objects are masked from 'package:stats':

    dbinom, dgamma, dnbinom, dnorm, dpois, dt, pnorm

The following objects are masked from 'package:base':

    cbind, pmax, pmin, rbind

> 
> (f.chk <- system.file("check-tools.R", package="Rmpfr", mustWork=TRUE))
[1] "D:/RCompile/CRANincoming/R-devel/lib/Rmpfr/check-tools.R"
> source(f.chk, keep.source=FALSE)
> ## -> assert.EQ.() showSys.time()
> 
> options(warn = 1)# warnings *immediately*
> (doExtras <- Rmpfr:::doExtras())
[1] FALSE
> 
> ### 1. Integration -----------------------------------------------
> 
> ## Example from Lauren K, June 2014 (~/R/MM/Pkg-ex/Rmpfr/integrateR-LaurenK.R):
> beta0  <- 0.05
> beta1  <- 0.05
> (tau <- sqrt(0.01*0.05*0.95/0.99))# = 0.0219..
[1] 0.02190429
> ##
> Z00  <- 9
> Z01  <- 1
> Z10  <- 18
> Z11  <- 2
> N <- Z00+Z01+Z10+Z11
> 
> integrand <- function(u) {
+     ee.u <- exp(-u^2/2)/(sqrt(2*pi)*tau)
+     b0u <- beta0 + tau*u
+     b1u <- beta1 + b0u # == beta0+beta1+ tau*u
+ 
+     ee.u ^ (Z00+Z01 + Z10+Z11) *
+         (1-b0u)^Z00 * b0u ^ Z01 *
+         (1-b1u)^Z10 * b1u ^ Z11
+ }
> 
> ## MM: note how the integrand() function looks:
> op <- par(mfcol=c(3,1), mgp=c(1.5, .6, 0), mar=.1+c(3,3,0,0))
>   ax1 <- function(a,b) axis(1, at=c(a,b), col=2, col.axis=2, tcl= +3/4, mgp=c(3,0,0))
>   curve(integrand, -5,5, n=1000)
>   cc <- adjustcolor(2, 1/4) ## look closer:
>   ep <- .01; rect(-3, -ep, 0, +ep, col=cc, border=cc); ax1(-3,0)
>   curve(integrand, -3,0, n=1000, ylim = c(-ep,ep))
>   ## but now look really closely :
>   ep <- .001; rect(-3, -ep, -2, +ep, col=cc); ax1(-3,-2)
>   curve(integrand, -3,-2, n=1000, ylim = c(-ep, ep))
> par(op)
> 
> (I1  <- integrate(integrand,lower = -100, upper = 100))
1.396355e+33 with absolute error < 1.6e+29
> (I1. <- integrate(integrand,lower = -100, upper = 100, rel.tol = 1e-14))
1.396355e+33 with absolute error < 8e+18
> 
> showSys.time(I2 <- integrateR(integrand, lower = -100, upper = 100))
Warning in integrateR(integrand, lower = -100, upper = 100) :
  no convergence up to order 13; last relative change = 1.285403e-05
Consider setting 'ord = <n>' (e.g. = 14).
Time    user  system elapsed 
Time    0.01    0.00    0.02 
> I2 ## ... Warning ‘no convergence up to order  13’
Non-convergence message 'no convergence up to order 13; last relative change = 1.285403e-05
Consider setting 'ord = <n>' (e.g. = 14).'
1.3964e+33 with absolute error < 1.7949e+28
> ## Solaris Sparc (2014-06, CRAN checks); thanks Brian: print(I2[1:2], digits=15)
> I2.Solaris <- list(value = 1.3963550396006e+33,  abs.error = 1.79487857486724e+28)
> I.db <-       list(value = 1.39635503960059e+33, abs.error = 1.79487857478077e+28)
> stopifnot(
+     all.equal(I2[1:2], I.db, tol = 1e-10)# Solaris SPARC needs at least 4.8e-11
+ )
> 
> ## Now using high accuracy
> showSys.time(I3 <- integrateR(integrand, lower = mpfr(-100, precBits=256), upper = 100))
Warning in integrateR(integrand, lower = mpfr(-100, precBits = 256), upper = 100) :
  no convergence up to order 13; last relative change = 1.285403e-5
Consider setting 'ord = <n>' (e.g. = 14).
Time    user  system elapsed 
Time    1.60    0.01    1.61 
> ## much slower but not better (and not worse)
> I3
Non-convergence message 'no convergence up to order 13; last relative change = 1.285403e-5
Consider setting 'ord = <n>' (e.g. = 14).'
1.3964e+33 with absolute error < 1.7949e+28
> assert.EQ.(sapply(I3[1:2], asNumeric), unlist(I.db))
Mean relative difference: 3.569039e-15
> ## Really get better when decreasing the integration interval
> ## from  [-100, 100] to  [-10, 10] ... which should give "the same"
> showSys.time(I4 <- integrateR(integrand, lower = mpfr(-10, precBits=256), upper = 10,
+                               ord = 15, verbose=TRUE))
 ord = 15; ==> evaluating integrand at  2^(ord+1)-2 = 65534 locations
n= 1, 2^n=        2 | I = 4.085188566e+34, abs.err =  4.085189e+34
n= 2, 2^n=        4 | I = 8.170377132e+33, abs.err =  3.268151e+34
n= 3, 2^n=        8 | I = 4.712016441e+33, abs.err =  3.458361e+33
n= 4, 2^n=       16 | I = 2.335919214e+33, abs.err =  2.376097e+33
n= 5, 2^n=       32 | I = 1.181962279e+33, abs.err =  1.153957e+33
n= 6, 2^n=       64 | I = 1.220705480e+33, abs.err =  3.874320e+31
n= 7, 2^n=      128 | I = 1.410495454e+33, abs.err =  1.897900e+32
n= 8, 2^n=      256 | I = 1.396201858e+33, abs.err =  1.429360e+31
n= 9, 2^n=      512 | I = 1.396354229e+33, abs.err =  1.523703e+29
n=10, 2^n=     1024 | I = 1.396354970e+33, abs.err =  7.409143e+26
n=11, 2^n=     2048 | I = 1.396354964e+33, abs.err =  5.998797e+24
n=12, 2^n=     4096 | I = 1.396354964e+33, abs.err =  6.619581e+21
n=13, 2^n=     8192 | I = 1.396354964e+33, abs.err =  1.662622e+18
n=14, 2^n=    16384 | I = 1.396354964e+33, abs.err =  1.021861e+14
n=15, 2^n=    32768 | I = 1.396354964e+33, abs.err =   1.561794e+9
Time    user  system elapsed 
Time    6.57    0.11    6.69 
> ## ~ 6.6 sec [lynne 2013]
> I4
1.3964e+33 with absolute error < 1.5618e+9
> 
> ## on the left side, there is "nothing" (and negative, as we know!):
> showSys.time(I0.1 <- integrateR(integrand, lower = mpfr(-1000, precBits=256),
+                                 upper = -10,  ord= 11, verbose=TRUE))
 ord = 11; ==> evaluating integrand at  2^(ord+1)-2 = 4094 locations
n= 1, 2^n=        2 | I = -2.859606864e-613, abs.err = 5.719214e-613
n= 2, 2^n=        4 | I = -1.334483203e-613, abs.err = 1.525124e-613
n= 3, 2^n=        8 | I = -6.566504650e-614, abs.err = 6.778327e-614
n= 4, 2^n=       16 | I = -3.270376825e-614, abs.err = 3.296128e-614
n= 5, 2^n=       32 | I = -1.633589988e-614, abs.err = 1.636787e-614
n= 6, 2^n=       64 | I = -8.165955325e-615, abs.err = 8.169945e-615
n= 7, 2^n=      128 | I = -4.082728442e-615, abs.err = 4.083227e-615
n= 8, 2^n=      256 | I = -2.041333072e-615, abs.err = 2.041395e-615
n= 9, 2^n=      512 | I = -1.020662642e-615, abs.err = 1.020670e-615
n=10, 2^n=     1024 | I = -5.103308345e-616, abs.err = 5.103318e-616
n=11, 2^n=     2048 | I = -2.551653564e-616, abs.err = 2.551655e-616
Time    user  system elapsed 
Time    0.44    0.00    0.44 
> showSys.time(I0.2 <- integrateR(integrand, lower = mpfr(10, precBits=256),
+                                 upper = 1000, ord= 11, verbose=TRUE))
 ord = 11; ==> evaluating integrand at  2^(ord+1)-2 = 4094 locations
n= 1, 2^n=        2 | I = 6.250128073e-618, abs.err = 1.250026e-617
n= 2, 2^n=        4 | I = 2.916726434e-618, abs.err = 3.333402e-618
n= 3, 2^n=        8 | I = 1.435214595e-618, abs.err = 1.481512e-618
n= 4, 2^n=       16 | I = 7.147931510e-619, abs.err = 7.204214e-619
n= 5, 2^n=       32 | I = 3.570472142e-619, abs.err = 3.577459e-619
n= 6, 2^n=       64 | I = 1.784800116e-619, abs.err = 1.785672e-619
n= 7, 2^n=      128 | I = 8.923455870e-620, abs.err = 8.924545e-620
n= 8, 2^n=      256 | I = 4.461659853e-620, abs.err = 4.461796e-620
n= 9, 2^n=      512 | I = 2.230821417e-620, abs.err = 2.230838e-620
n=10, 2^n=     1024 | I = 1.115409645e-620, abs.err = 1.115412e-620
n=11, 2^n=     2048 | I = 5.577046893e-621, abs.err = 5.577050e-621
Time    user  system elapsed 
Time    0.39    0.00    0.39 
> I0.1
-2.5517e-616 with absolute error < 2.5517e-616
> I0.2
5.5770e-621 with absolute error < 5.5770e-621
> I4
1.3964e+33 with absolute error < 1.5618e+9
> 
> ## Integral [-1000, +1000 ] = Int[-1000, -10] + Int[-10, +10] + Int[+10, +1000]:
> 
> I4 $value + I0.1 $value + I0.2 $value
1 'mpfr' number of precision  256   bits 
[1] 1396354963674017387015147490193159.224442394550727848837192551722726424516192161
> ## but this is really the same as just the middle:
> stopifnot(I4 $value + I0.1 $value + I0.2 $value
+           == I4 $value)
> 
> value <- I4$value; delta <- I4$abs.err
> nDig <- -asNumeric(log10(delta/value))
> cat("Correct number of digits: ", round(nDig, 2),"\n",
+     "Integral I =              ", format(I4$value, digits = ceiling(nDig)),
+     " (last change change= ", format(delta, digits = 7),")\n",
+     "integrate(.) =            ", format(I1 $value, digits = 22),"\n",
+     "integrate(., rtol=1e-15)= ", format(I1.$value, digits = 22),"\n", sep="")
Correct number of digits: 23.95
Integral I =              1.39635496367401738701515e+33 (last change change= 1.561794e+9)
integrate(.) =            1.396354963673671429648e+33
integrate(., rtol=1e-15)= 1.396354963674017017894e+33
> 
> 
> 
> 
> ### 2. Root Finding ----------------------------------------------
> 
> ### 3. Optimization / Minimization, .. ---------------------------
> 
> 
> 
> cat('Time elapsed: ', proc.time(),'\n') # "stats"
Time elapsed:  10.1 0.21 10.31 NA NA 
> 
> proc.time()
   user  system elapsed 
  10.10    0.21   10.31