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

[Kacper Kowalik]

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

-------------- 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/20130218/28fe99ae/attachment.sig>