[R] Fast evaluation of functions in 3D domains

[David Winsemius]

On Mar 25, 2015, at 4:55 AM, <Lluis.Hurtado at uv.es> <Lluis.Hurtado at uv.es> wrote:

> Dear all,
> 
> Finally I have tried three different options to integrate a function in a 3D volume. Here 
> I show a test example. The volume is the box [0,100] x [0,100] x [0,100] and the 
> function is
> 
> nfw(d) = 4/((d/5)*(1+(d/5))^2)
> 
> where d is the distance between each point in the box to the point (50,50,40).
> 
> 1-Grid of thick 1 in R (10^6 points)
> 
>> model <- function(x) 
> {
>   d <- sqrt((x[,1]-50)^2 + (x[,2]-50)^2 + (x[,3]-40)^2)
>   r <- 4.0/((d/5)*(1+(d/5))^2)
>   r
> } 
>> sum(model(x))
> [1] 10287.52
>> system.time(sum(model(x)))
>   user  system elapsed 
>  0.052   0.003   0.053 
> 
> 2-Grid with thick 1 in C++ calling from R. Function model_cpp is a function written in 
> C++ reproducing model function as above. (10^6 points)
> 
>> model <- function(x)
>  {
>  	param <- c(50,50,40,5)
>  	.Call('model_cpp',x[,1],x[,2],x[,3],param)
>  }
>> sum(model(x))
> [1] 10287.52
>> system.time(sum(model(x)))
>   user  system elapsed 
>  0.028   0.000   0.028 
> 
> 3-cubature. Mr Tom Jagger kindly proposed to use cubature package:
> http://cran.r-project.org/web/packages/cubature/cubature.pdf
> 
>> model <- function(x) 
> +  {
> +     d <- sqrt((x[1]-50)^2 + (x[2]-50)^2 + (x[3]-40)^2)
> +     r <- 4.0/((d/5)*(1+(d/5))^2)
> +     r
> +   }
>> adaptIntegrate(model, rep(0,3), rep(100,3), tol=1e-4)
> $integral
> [1] 10303.16
> 
> $error
> [1] 1.029888
> 
> $functionEvaluations
> [1] 48609
> 
> $returnCode
> [1] 0
> 
>> system.time(adaptIntegrate(model, rep(0,3), rep(100,3), tol=1e-4))
>   user  system elapsed 
>  0.232   0.002   0.246
> 
> As you can see the second option is the fastest, but the third one is probably more 
> accurate. 
> 
> The function nfw(d) has an analytical primitive when integrated in a sphere If now we 
> reproduce the calculations for cases 1 and 2 in a sphere (named R) of radius 40 
> centered in (50,50,40), the first two methods give me the following result:
> 
> sum(model(R))
> [1] 8204.711
> 
> while the exact solution is
> 
>> 16*pi*(5^3)*(log(5+40) - log(5) - 40/(40+5))
> [1] 8220.516
> 
> However, I can not try the same with cubature since the code is prepared only to be 
> used in hypercubes.

This is the same as one gets with considering this to be a radial symmetric function and transforming the problem to one dimension using the infinitesimal transformation dV = 4*pi*dr:

> f <- function(x) 4*pi*x^2*4.0/((x/5)*(1+(x/5))^2)
> integrate(f, 0, 40 , subdivisions=10000)
8220.516 with absolute error < 0.0012

I wasn't able to achieve convergence to that value using adaptIntegrate using a Boolean variation on Duncan's suggestion:

 model <- function(x) 
   {
      d <- sqrt((x[1]-50)^2 + (x[2]-50)^2 + (x[3]-40)^2)
      r <- 4.0/((d/5)*(1+(d/5))^2)
      dinside <- (d <= 40)*r
    }
> adaptIntegrate(model, rep(0,3), rep(100,3), tol=1e-3, maxEval=1000000)
$integral
[1] 8199.022

$error
[1] 33.33091

$functionEvaluations
[1] 1000065

$returnCode
[1] 0

Without the maxEval the method exceeded my patience. If the real problem is radial symmetric (and it might be noting your domain of investigation) this may offer speed and accuracy advantages.

--
David.

> 
> As I am using non integrable functions I could try to increase the density of the grid 
> and see if I can obtain accurate results before reaching high time costs or study how 
> important is to reach that accuracy, it may be not that important for my algorithm. 
> 
> Anyway, thank you all for you time and ideas.
> 
> Lluís Hurtado
> IFCA
> 
>> 
>> On Mar 23, 2015, at 3:44 AM, <Lluis.Hurtado at uv.es> <Lluis.Hurtado at uv.es> wrote:
>> 
>>> Dear all,
>>> 
>>> I am currently working with the spatstat package with 3D samples. I am trying to 
>>> evaluate a non analytical function over the window that encloses the sample and I 
>>> need to know which is the fastest way of doing it.
>>> 
>>> The function input is a 3 coordinate position in the window (x,y,z) and a list of 
>>> parameters (a,b,c). The output is a numerical value. 
>>> 
>>> n <- function(x,y,z,a,b,c)
>> 
>> Perhaps:
>> 
>> dfrm <- as.data.frame.table(your_volume_matrix) 
>> n.out <- with(dfrm,  mapply( n, x=x, y=y, z=z, MoreArgs=list(a=a,b=b,c=c) ) _
>> dim(n.out) <- dim(your_volume_matrix)
>> 
>> You don't describe the form of this "3 coordinate position in the window (x,y,z)" so 
> perhaps the arguments will need to be extracted. I took a WAG at one approach. If 
> it's not in long-form, you need configure the array indices for either a volume or 
> surface into a dataframe, perhaps with `expand.grid` or `as.data.frame.table`.
>> 
>> You also don't describe the sort of integration you imagine. Why not a simple sum 
> of that result divided by the volume? I cannot imagine any faster procedure .
>> 
>> 
>>> But I need to do it over the whole volume.
>>> 
>>> For 2 dimensions it can be done with
>>> 
>>> A <- as.im(function,window,parameters)
>>> norm <- integral.im(A)
>>> 
>>> For 3 dimensions I have tried to pass an array of a grid covering the window (like 
> a 
>>> quadrature scheme) and then summing up the output array, but I would like to 
> know if 
>>> there is any faster way of integrating the function.
>>> 
>>> Thank you very much,
>>> 
>>> Lluís Hurtado
>>> IFCA
>>> www.ifca.unican.es
>>> 
> 
> --
> Lluís Hurtado-Gil
> Observatori Astronòmic. Universitat de València.

David Winsemius
Alameda, CA, USA