[R] Fast evaluation of functions in 3D domains
Lluis.Hurtado at uv.es
Lluis.Hurtado at uv.es
Wed Mar 25 12:55:00 CET 2015
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).
More information about the R-help
mailing list