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

[Kacper Kowalik]

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.

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.

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
)

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

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 900 bytes
Desc: OpenPGP digital signature
URL: <http://lists.spacepope.org/pipermail/yt-dev-spacepope.org/attachments/20130215/20c5a2d0/attachment.pgp>