[R] Fast evaluation of functions in 3D domains

[Lluis.Hurtado at uv.es]

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.

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 en uv.es> <Lluis.Hurtado en 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
> > 
> > ______________________________________________
> > R-help en r-project.org mailing list -- To UNSUBSCRIBE and more, see
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide http://www.R-project.org/posting-guidehtml
> > and provide commented, minimal, self-contained, reproducible code.
> 
> David Winsemius
> Alameda, CA, USA
> 
> 
> 


--
Lluís Hurtado-Gil
Observatori Astronòmic. Universitat de València.
Edifici Instituts d'Investigació. Parc Científic.
C/ Catedrático Agustín Escardino n°7
Campus Burjassot-Paterna
46980 Paterna València (Spain).