[SciPy-User] Integration over Voronoi cells

[Friedrich Romstedt]

2011/11/24 Eraldo Pomponi <eraldo.pomponi at 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.