R help - Mar 2015 - Fast evaluation of functions in 3D domains

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)

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

On 23/03/15 23:44, 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.
I gather that you mean that you want to evaluate the *integral* of the 
function over the window.

What do you mean by "non analytical"?  You have, it would appear, some
way of calculating the function at an arbitrary point (x,y,z) in 
3-space, which makes it at least in some some sense "analytical".
> 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)
>
> 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).
There is nothing very sophisticated about the way integral.im() works. 
Just look at the code.  The operative line of code is:

    a <- with(f, sum(v, na.rm = TRUE) * xstep * ystep)
> 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.
It sounds, at least vaguely, as though you are imitating in 3D what 
integral.im() does in 2D.  So, given that you are doing it right, you 
are doing as well as integral.im() does.

The calculations should be pretty fast since sum() is pretty fast.  The 
bulk of the computational burden is likely to be the evaluation of the 
function at all of the voxel centres.

Note that this is a very unsophisticated method of numerical quadrature. 
  The integrand is being approximated by a step function.  The method 
seems to work well enough for the uses that we put it to in spatstat. 
Whether it works well enough for your purposes depends on what your 
purposes are.

cheers,

Rolf Turner

-- 
Rolf Turner
Technical Editor ANZJS
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276
Home phone: +64-9-480-4619

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
> 
> ______________________________________________
> R-help at 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-guide.html
> and provide commented, minimal, self-contained, reproducible code.
David Winsemius
Alameda, CA, USA

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 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
> > 
> > ______________________________________________
> > R-help at 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).