On Mar 25, 2015, at 4:55 AM, <lluis.hurt...@uv.es> <lluis.hurt...@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.hurt...@uv.es> <lluis.hurt...@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
______________________________________________
R-help@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.
