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

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 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? Kacper