Hi folks, I'm working on the integration of a function like: K(r,t) = 1/(2piDr)exp[-r^2/Dt] [1] over Voronoi cells (r is the distance from the point at which is associated the cell). I googled a lot and I found this two useful hints: http://stackoverflow.com/questions/5941113/looking-for-python-package-for-nu... http://mathforum.org/kb/message.jspa?messageID=4963570&tstart=0 but I'm still not able to understand how I should do this integration. I have the function that returns the segments ( list of [(x_star,y_start),(x_sop,y,stop)] ) necessary to construct the Voronoi cells associated to a set of points. Could someone suggest how to proceed ? There's also a numerical problem connected with the integration due to the singularity in r==0. Could you suggest which is a reasonable stable integration method available in scipy that could handle the function [1]? Cheers, Eraldo
2011/11/23 Eraldo Pomponi <eraldo.pomponi@gmail.com>:
Hi folks, I'm working on the integration of a function like: K(r,t) = 1/(2piDr)exp[-r^2/Dt] [1] over Voronoi cells (r is the distance from the point at which is associated the cell). I googled a lot and I found this two useful hints: http://stackoverflow.com/questions/5941113/looking-for-python-package-for-nu... http://mathforum.org/kb/message.jspa?messageID=4963570&tstart=0 but I'm still not able to understand how I should do this integration. I have the function that returns the segments ( list of [(x_star,y_start),(x_sop,y,stop)] ) necessary to construct the Voronoi cells associated to a set of points. Could someone suggest how to proceed ? There's also a numerical problem connected with the integration due to the singularity in r==0. Could you suggest which is a reasonable stable integration method available in scipy that could handle the function [1]?
Well, AISI, your function is radially symmetric, so the integration over 1/r can be done analytically within a small disk where the Gaussian is approximately one. Since: 2 \pi r * [1 / (2 \pi D r)] is just 1/D constantly. So you're left with [the integration up to] R_0/D when R_0 is the radius of that disk (where the integration is done analytically). [No more divergence, since you cut out the centre region.] You could also just integrate with the full Gaussian, since the kernel just integrates up the lengthes of the disk perimeters, if you understand. All you need to do, AISI atm, is to calculate a function that gives you the perimeter length of the circle which is inside your Voronoi cell. This is zero at r = 0 and will stay bounded for all finite radii. So integration should be fairly straightforward. You could even just do a sum over a 1D grid for the radius, since the function varies slowly, this should be fast and easy (both in runtime as in coding time). You might not even need scipy for this. Additionally, the Gaussian just introduces a suppression of radii which are farther outside. It makes the complete function bounded even for an infinitely large Voronoi cell :-) The kernel 1/r makes things easy and convenient, instead of making things troublesome: 1) because it exhibits rotational invariance; 2) because it converges nicely in 2D for an integral over a circle line I might overlook something obvious. Is the function you gave really the full function or only the "kernel", the weighting function? Friedrich P.S.: You divide your tringles resulting from the centre point and the boundary lines into two parts, separated by the closest point on the boundary, which is always in the bounds of that boundary line. Then you just have a perimeter length which is linear in r until you touch the line, and then it'll be a littel more complicated. I think you'll figure it out. P.P.S.: You're really lucky that your kernel is 1/r and not 1/r^2, iirc ;-)
Dear Friedrich, First of all let me thank you for your suggestions. It took a little bit to get into them but I still have some BIG doubts.
Well, AISI, your function is radially symmetric, so the integration over 1/r can be done analytically within a small disk where the Gaussian is approximately one. Since:
Not so clear why the Gaussian is approximately one (it is true just when r is small enough)
2 \pi r * [1 / (2 \pi D r)] is just 1/D constantly. So you're left with [the integration up to] R_0/D when R_0 is the radius of that disk (where the integration is done analytically). [No more divergence, since you cut out the centre region.] You could also just integrate with the full Gaussian, since the kernel just integrates up the lengthes of the disk perimeters, if you understand.
Let me rewrite the function (it is just the kernel not the full function): K(r,t) = 1/(2piDr)exp[-r^2/Dt] [1] This function is rotational invariant [a] so a clever way to integrate it is to consider a disk of radius R0. The perimeter of this disk is equal to: l=2*pi*R0 [2] The term 1/(2*pi*D*r) gives a constant contribution along the path eq.2, equal to 1/D so what is left to do is to integrate the term: exp[-r^2/Dt] in the interval [R0/D,R0/D+delta] and multiply by the length of the chose path, i.e. eq.2. and sum on a 1D grid in the interval [0,R0]. ------>>>>> Is it correct? <<<<<----------- The term exp[-r^2/Dt] can be analytically integrated (removing t for convenience) : g(r) = sqrt(pi)*erf(r*sqrt(1/D))/(2*sqrt(1/D)) [3] Its defined integral in [R0/D,R0/D+delta], multiplied by the length of the choose path, give us: g'(r) = 2*pi*r * [ g(r/D + delta) - g(r/D)] [4] summing in the grid k=1..R0, we obtain the integral of eq.1 in [0,R0]: integral([1]) = sum([g'(k) for k in np.arange(0,R0,0.001)] ) What do you think ?
All you need to do, AISI atm, is to calculate a function that gives
you the perimeter length of the circle which is inside your Voronoi cell. This is zero at r = 0 and will stay bounded for all finite radii. So integration should be fairly straightforward. You could even just do a sum over a 1D grid for the radius, since the function varies slowly, this should be fast and easy (both in runtime as in coding time). You might not even need scipy for this.
Additionally, the Gaussian just introduces a suppression of radii
which are farther outside. It makes the complete function bounded even for an infinitely large Voronoi cell :-)
The kernel 1/r makes things easy and convenient, instead of making things troublesome:
1) because it exhibits rotational invariance; 2) because it converges nicely in 2D for an integral over a circle line
I might overlook something obvious. Is the function you gave really the full function or only the "kernel", the weighting function?
Friedrich
P.S.: You divide your tringles resulting from the centre point and the boundary lines into two parts, separated by the closest point on the boundary, which is always in the bounds of that boundary line.
I didn't understood what you mean here. Sorry
you just have a perimeter length which is linear in r until you touch the line, and then it'll be a littel more complicated. I think you'll figure it out.
Yes , I hope so. When the path in eq.2 is no more circular, the integration along it is more complicated. A bit difficult to figure out how to do that programmatically but this is the next problem.
P.P.S.: You're really lucky that your kernel is 1/r and not 1/r^2, iirc ;-)
I guess so but I think that I will fully understand this sentence just when my doubts will be cleared. Thanks a lot for your help. Cheers, Eraldo
2011/11/24 Eraldo Pomponi <eraldo.pomponi@gmail.com>:
Dear Friedrich, First of all let me thank you for your suggestions. It took a little bit to get into them but I still have some BIG doubts.
:-) Apparently I don't see any doubts on your side, rather a way to go.
Well, AISI, your function is radially symmetric, so the integration over 1/r can be done analytically within a small disk where the Gaussian is approximately one. Since:
Not so clear why the Gaussian is approximately one (it is true just when r is small enough)
"r is small enough" = "within a small disk" :-)
2 \pi r * [1 / (2 \pi D r)] is just 1/D constantly. So you're left with [the integration up to] R_0/D when R_0 is the radius of that disk (where the integration is done analytically). [No more divergence, since you cut out the centre region.] You could also just integrate with the full Gaussian, since the kernel just integrates up the lengthes of the disk perimeters, if you understand.
Let me rewrite the function (it is just the kernel not the full function):
That's a pity. Still you can average your operand undergoing the kernel multiplication just on the path length. That's what your integral does.
K(r,t) = 1/(2piDr)exp[-r^2/Dt] [1] This function is rotational invariant [a] so a clever way to integrate it is to consider a disk of radius R0. The perimeter of this disk is equal to: l=2*pi*R0 [2] The term 1/(2*pi*D*r) gives a constant contribution along the path eq.2, equal to 1/D so what is left to do is to integrate the term: exp[-r^2/Dt]
Exactly. You got this. :-)
in the interval [R0/D,R0/D+delta] and multiply by the length of the chose path, i.e. eq.2. and sum on a 1D grid in the interval [0,R0]. ------>>>>> Is it correct? <<<<<-----------
Not quite. But nearly. Let's rewrite the remaining Gaussian like this: G = exp(-r^2/(2 k^2)) where: (2 k^2) = Dt . Just to bring it into the well-known standard Gaussian form. The std deviation of your Gaussian (1 sigma range) is then k. So what you get is: I(R_1) = \int_0^{R_1} G(r) 1/D \mathrm{d}r . The integral over the angle along the circle line is already consumed! It walked into the 1/D term!
The term exp[-r^2/Dt] can be analytically integrated (removing t for convenience) : g(r) = sqrt(pi)*erf(r*sqrt(1/D))/(2*sqrt(1/D)) [3]
Well, afair, there is no closed form for the integral? Prove me wrong. I think erf (error function) is just the "name" for what the integral is, no? When we cannot write it down in a "closed form" we just invent a new primitive function to make it closed by definition. As I said, the closed form for erf just does not pop into my mind. Of course, you have some points of erf(x) defined, when you say erf is the integral starting at -oo. In fact, I(R_1) "is" just erf, you can just read it off the equation. In general, one more comment, it's IMHO not a good idea to just suppress parameters like t, it simplifies too much at the cost of unit inconstency. A physicist speaking :-) Better introduce another quantity and avoid confusion. Then you can safely say "for lim k to 0" and stuff. ... I'm not sure on the erf Eq. (3), but it's a straightforward linear variable substitution, speaking algebra. Just pull out the 1/D, transform r s.t. k = 1, and I think you should be able to apply erf.
Its defined integral in [R0/D,R0/D+delta], multiplied by the length of the
Still, I have no clue why you have a /D in the span.
choose path, give us: g'(r) = 2*pi*r * [ g(r/D + delta) - g(r/D)] [4] summing in the grid k=1..R0, we obtain the integral of eq.1 in [0,R0]: integral([1]) = sum([g'(k) for k in np.arange(0,R0,0.001)] ) What do you think ?
Looks rather wrong. There might be some kernel of truth aside of that you did some integration I don't have the inclination to find the error in. That's up to you. But the 2 pi r should be consumed by the angular integral, you just integrate a constant! (The constant is 1/D) (And you integrate it over a radius quantity.) I would consider the following approach to understand the problem: 1. Figure out how to just integrate the kernel, without the "kerneled" function. And without trinangular bounds, just on the full circle. 2. Then introduce an angular, but constant bound on the angular integral (i.e. not \int_0^{2 \pi} but \int_Q where Q is a subset of [0, 2 pi]). 3. Next, you will be able to vary the Q with r, i.e. Q(r). You will see that the full integral (withour kerneled function) is just the integral over the measure length of Q with r. 4. And then, you can safely introduce the kerneled function (is there a term for "kerneled"?). AISI, it's just the average of the kerneled function over Q, integrated with r. AISI, from (4.), the full integral *is just the average of the kerneled function, with a weighting s.t. all circles have in net the same weight.* It's stunningly simple. Is it true? And forget about the small disk. You don't need it at all, your kernel cancels nicely with the r factor from the integration (path length).
All you need to do, AISI atm, is to calculate a function that gives you the perimeter length of the circle which is inside your Voronoi cell. This is zero at r = 0 and will stay bounded for all finite radii. So integration should be fairly straightforward. You could even just do a sum over a 1D grid for the radius, since the function varies slowly, this should be fast and easy (both in runtime as in coding time). You might not even need scipy for this.
Additionally, the Gaussian just introduces a suppression of radii which are farther outside. It makes the complete function bounded even for an infinitely large Voronoi cell :-)
The kernel 1/r makes things easy and convenient, instead of making things troublesome:
1) because it exhibits rotational invariance; 2) because it converges nicely in 2D for an integral over a circle line
I might overlook something obvious. Is the function you gave really the full function or only the "kernel", the weighting function?
Friedrich
P.S.: You divide your tringles resulting from the centre point and the boundary lines into two parts, separated by the closest point on the boundary, which is always in the bounds of that boundary line.
I didn't understood what you mean here. Sorry
you just have a perimeter length which is linear in r until you touch the line, and then it'll be a littel more complicated. I think you'll figure it out.
Yes , I hope so. When the path in eq.2 is no more circular, the integration along it is more complicated. A bit difficult to figure out how to do that programmatically but this is the next problem.
See above. Skip for now, I'd say. :-)
P.P.S.: You're really lucky that your kernel is 1/r and not 1/r^2, iirc ;-)
I guess so but I think that I will fully understand this sentence just when my doubts will be cleared. Thanks a lot for your help.
:-) Good luck, and please don't believe anything I say, just believe what you figure out yourself! Happy sciencing! Friedrich :-) P.S.: What's your kerneled function? I.e. does it have some symmetries (or not). And why integrating over Voronoi cells? Since the integral, thanks to the Gaussian, converges also on full domain, I believe if the Voronoi cells are introduced artificially you can drop them and do the full integral without the hassle do track the boundaries. Just a shot in the dark.
participants (2)
-
Eraldo Pomponi -
Friedrich Romstedt