[yt-dev] Avoiding roundoff errors in comparison of floats

[Matthew Turk]

Hi Kacper,

On Mon, Feb 18, 2013 at 4:28 PM, Kacper Kowalik <xarthisius.kk at gmail.com> wrote:
> On 18.02.2013 21:02, Matthew Turk wrote:
>> Hi Kacper,
>>
>> Sorry this email has languished so long -- I meant to get to it, but
>> have not had the chance to give it the thought it needs until now.
>>
>> On Fri, Feb 15, 2013 at 12:36 PM, Kacper Kowalik
>> <xarthisius.kk at gmail.com> wrote:
>>> Hi All!
>>> I've been recently struck by roundoff errors in yt. AFAIR some variants
>>> of such bugs have been already reported from time to time on #yt or
>>> mailing lists and I'd like to propose a general solution that should
>>> solve this.
>>
>> Yes, they have.  And our solutions have not been terribly good.  For
>> one example, we have (in the past) done reconstruction of hierarchies
>> to lock against integer positions.
>>
>>>
>>> The problem can be seen here: http://imgur.com/oIaarPR
>>> It's a dense blob (\rho = 2) moving diagonally across uniform medium
>>> (\rho = 1). As you can see there are blueish stripes on the few block
>>> boundaries. They're caused by yt thinking that blocks on n+1 refinement
>>> level overlap with those regions in blocks on "n" level. It all happens
>>> here:
>>>
>>> (yt/data_objects/object_finding_mixin.py)
>>> ObjectFindingMixin::get_box_grids
>>> ...
>>> grid_i = np.where(
>>>   (np.all(self.grid_right_edge > left_edge, axis=1)
>>>  & np.all(self.grid_left_edge < right_edge, axis=1)) == True
>>> )
>>>
>>> The problem is that comparing floats is *not* accurate. What should be
>>> done: instead of "a > b", we should do e.g. "(a - b) > eps", where eps
>>> is small number grater than machine precision.
>>
>> Well, I think there are two problems here.  The first is that this is
>> asymmetric.  It's possible in this case to have a slice, if exactly
>> set on the boundary, that does not intersect anything.  The second is
>> that the floating point comparison can make this so.
>>
>>>
>>> I've fixed this particular bit as follows:
>>>
>>> ...
>>> dx = self.get_smallest_dx() * 0.1
>>> grid_i = np.where(
>>>   (np.all(self.grid_right_edge > (left_edge + dx), axis=1)
>>>  & np.all(self.grid_left_edge < (right_edge - dx), axis=1)) == True
>>> )
>>
>> I have to think this through, but doesn't this set it up such that if
>> you have a slice right on the edge of grid boundaries (for instance,
>> in a symmetric domain this would be at 0.0) you can have it
>> intersecting both grids that touch it?
>>
>> My inclination is to instead insert a preferred direction, for edge cases.
>
> I think that I got it wrong but following changes should be better and
> give you *preferred* direction. Let's consider two grids [-1, 0] and [0,
> 1] and slice through 0.0. In AMRSliceBase:
>
> def _get_list_of_grids(self):
>   goodI = ((self.pf.h.grid_right_edge[:,self.axis] > self.coord)
>         &  (self.pf.h.grid_left_edge[:,self.axis] <= self.coord ))
>   self._grids = self.pf.h.grids[goodI] # Using sources not hierarchy
>
> As it's now we would have:
> grid(0): (0 > 0) & (-1 <= 0) which can be either true or false
> grid(1): (1 > 0) & (0 <= 0) which can be either true or false
>
> if we instead did:
> def _get_list_of_grids(self):
>   goodI = ((self.pf.h.grid_right_edge[:,self.axis] - self.coord >
> np.finfo(np.float64).eps)
>         &  (self.pf.h.grid_left_edge[:,self.axis] - self.coord <=
> np.finfo(np.float64).eps))
>   self._grids = self.pf.h.grids[goodI] # Using sources not hierarchy
>
> we would have:
> grid(0): (0 - 0 > eps) & (-1 - 0 <= eps) always false
> grid(1): (1 - 0 > eps) & (0 - 0 <= eps) always true
>
> What do you think?

I like it.  I think we should have an epsilon set in the hierarchy,
though, rather than re-calculating it.

One thing I should note is that all of the geometry handling has been
completely changed in 3.0.  It still does similar operations, but they
are defined in several routines:

 * select_grid(left_edge[3], right_edge[3], level ..)
 * select_cell(pos[3], dds[3] ..)

This means all of the epsilon handling code will need to be ported
over, *but*, I think that this will be much easier in the new way.
I'm rather happy with the new selection code as well, because it's
much more centralized.  Is there a chance you could add a unit test in
your PR that checks if a slice that runs on a boundary of a fake pf
selects the correct grids?

Here's the source for the selection routines in 3.0, for anyone interested:

https://bitbucket.org/yt_analysis/yt-3.0/src/162feb61f27ac820077d7bc9bd81e92cb5a414c2/yt/geometry/selection_routines.pxd?at=yt-3.0
https://bitbucket.org/yt_analysis/yt-3.0/src/162feb61f27ac820077d7bc9bd81e92cb5a414c2/yt/geometry/selection_routines.pyx?at=yt-3.0

The slice selector (a good example) is on line 651 of the second one.
Another good example is RegionSelector on line 453.

-Matt

> Kacper
>
>
> _______________________________________________
> yt-dev mailing list
> yt-dev at lists.spacepope.org
> http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org
>