Avoiding roundoff errors in comparison of floats
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
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@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@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org
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@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? Kacper
Hi Kacper, On Mon, Feb 18, 2013 at 4:28 PM, Kacper Kowalik <xarthisius.kk@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@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/162feb61f27ac820077d7bc9bd81e92... https://bitbucket.org/yt_analysis/yt-3.0/src/162feb61f27ac820077d7bc9bd81e92... 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@lists.spacepope.org http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org
participants (2)
-
Kacper Kowalik -
Matthew Turk