On 25/03/2015 7:55 AM, 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.
I don't know how well the cubature package would deal with the
discontinuity, but a simple way to integrate over a sphere would be to
set the function value to zero outside of it. Just change
r <- 4.0/((d/5)*(1+(d/5))^2)
to
r <- ifelse(d < 40, 4.0/((d/5)*(1+(d/5))^2), 0)
Duncan Murdoch
>
> 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
>>>
>>> ______________________________________________
>>> 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-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).
>
> ______________________________________________
> 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.
>
______________________________________________
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.
