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

[Matthew Turk]

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.

>
> 10% of the smallest cell size seems to be good approximation of eps at
> this point. As a result I've got: http://imgur.com/vCGQUKv
> Still one artifact remained. It was caused by exact the same issue in
> other place in the code
> (utilities/lib/geometry_utils.pyx:get_box_grids_level)
> Using similar trick there yielded: http://imgur.com/wpOCGaf \o/
>
> Sorry for the long mail, but I wanted to share this so we could inspect
> other place that could suffer from similar issues too.
> Cheers,
> Kacper
>
> P.s. Nice animations are available here:
> http://tinyurl.com/aw2tbof
> http://tinyurl.com/amp7noc
> http://tinyurl.com/arnqmfw
>
>
> _______________________________________________
> yt-dev mailing list
> yt-dev at lists.spacepope.org
> http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org
>