[Numpy-discussion] Catching and dealing with floating point errors

[Skipper Seabold]

On Mon, Nov 8, 2010 at 4:04 PM, Bruce Southey <bsouthey at gmail.com> wrote:
> On 11/08/2010 02:52 PM, Skipper Seabold wrote:
>> On Mon, Nov 8, 2010 at 3:42 PM, Bruce Southey<bsouthey at gmail.com>  wrote:
>>> On 11/08/2010 02:17 PM, Skipper Seabold wrote:
>>>> On Mon, Nov 8, 2010 at 3:14 PM, Skipper Seabold<jsseabold at gmail.com>    wrote:
>>>>> I am doing some optimizations on random samples.  In a small number of
>>>>> cases, the objective is not well-defined for a given sample (it's not
>>>>> possible to tell beforehand and hopefully won't happen much in
>>>>> practice).  What is the most numpythonic way to handle this?  It
>>>>> doesn't look like I can use np.seterrcall in this case (without
>>>>> ignoring its actual intent).  Here's a toy example of the method I
>>>>> have come up with.
>>>>>
>>>>> import numpy as np
>>>>>
>>>>> def reset_seterr(d):
>>>>>      """
>>>>>      Helper function to reset FP error-handling to user's original settings
>>>>>      """
>>>>>      for action in [i+'='+"'"+d[i]+"'" for i in d]:
>>>>>          exec(action)
>>>>>      np.seterr(over=over, divide=divide, invalid=invalid, under=under)
>>>>>
>>>> It just occurred to me that this is unsafe.  Better options for
>>>> resetting seterr?
>>>>
>>>>> def log_random_sample(X):
>>>>>      """
>>>>>      Toy example to catch a FP error, re-sample, and return objective
>>>>>      """
>>>>>      d = np.seterr() # get original values to reset
>>>>>      np.seterr('raise') # set to raise on fp error in order to catch
>>>>>      try:
>>>>>          ret = np.log(X)
>>>>>          reset_seterr(d)
>>>>>          return ret
>>>>>      except:
>>>>>          lb,ub = -1,1  # includes bad domain to test recursion
>>>>>          X = np.random.uniform(lb,ub)
>>>>>          reset_seterr(d)
>>>>>          return log_random_sample(X)
>>>>>
>>>>> lb,ub = 0,0
>>>>> orig_setting = np.seterr()
>>>>> X = np.random.uniform(lb,ub)
>>>>> log_random_sample(X)
>>>>> assert(orig_setting == np.seterr())
>>>>>
>>>>> This seems to work, but I'm not sure it's as transparent as it could
>>>>> be.  If it is, then maybe it will be useful to others.
>>>>>
>>>>> Skipper
>>>>>
>>>> _______________________________________________
>>>> NumPy-Discussion mailing list
>>>> NumPy-Discussion at scipy.org
>>>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>> What do you mean by 'floating point error'?
>>> For example, log of zero is not what I would consider a 'floating point
>>> error'.
>>>
>>> In this case, if you are after a log distribution, then you should be
>>> ensuring that the lower bound to the np.random.uniform() is always
>>> greater than zero. That is, if lb<= zero then you *know* you have a
>>> problem at the very start.
>>>
>>>
>> Just a toy example to get a similar error.  I call x<= 0 on purpose here.
>>
>>
> I was aware of that.
>

Ah, ok.  I don't mean to be short, just busy.

> Messing about warnings is not what I consider Pythonic because you
> should be fixing the source of the problem. In this case, your sampling
> must be greater than zero. If you are sampling from a distribution, then
> that should be built into the call otherwise your samples will not be
> from the requested distribution.
>

Basically, it looks like a small sample issue with an estimator.  I'm
not sure about the theory yet (or the underlying numerical issue), but
I've confirmed that the solution also breaks down using several
different solvers with a constrained version of the primal in GAMS to
ensure that it's not just a domain error or numerical
underflow/overflow.  So at this point I just want to catch the warning
and resample.  Am going to explore the "bad" cases further at a later
time.

Skipper