Lluis.Hurtado at uv.es
2015-Mar-25 11:55 UTC
[R] Fast evaluation of functions in 3D domains
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)" soperhaps 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 sumof 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 (likea> > quadrature scheme) and then summing up the output array, but I would like toknow 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).
On 25/03/2015 7:55 AM, Lluis.Hurtado at 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.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). > > ______________________________________________ > 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. >
On Mar 25, 2015, at 4:55 AM, <Lluis.Hurtado at uv.es> <Lluis.Hurtado at 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.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 >>> > > -- > Llu?s Hurtado-Gil > Observatori Astron?mic. Universitat de Val?ncia.David Winsemius Alameda, CA, USA
On Mar 25, 2015, at 8:45 AM, David Winsemius wrote:> > On Mar 25, 2015, at 4:55 AM, <Lluis.Hurtado at uv.es> <Lluis.Hurtado at 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)snipped quite a bit.> > 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:I meant to type: dV = 4*pi*r^2*dr:> >> f <- function(x) 4*pi*x^2*4.0/((x/5)*(1+(x/5))^2)But function definition was correct.>> 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.snipped David Winsemius Alameda, CA, USA