Re: [Numpy-discussion] rcond in polyfit
On 10/14/06, A. M. Archibald [EMAIL PROTECTED] wrote: On 14/10/06, Charles R Harris [EMAIL PROTECTED] wrote:snip I don't get the impression that the warnings module is much tested; Ihad similar headaches.Turns out to be a rather simple bug (feature?) in warnings, where it shortcircuits the filters by doing a quick look in a registry , which I assume means that the warning has been previously disabled. Commenting out one line makes warnings behave as expected. I also note that clearing the filter doesn't clear the warn once registry, which is probably some sort of bug also waiting to happen. Looking at the warning code leaves me feeling somewhat disappointed: few comments and little discussion of intended behaviour. Chuck - Using Tomcat but need to do more? Need to support web services, security? Get stuff done quickly with pre-integrated technology to make your job easier Download IBM WebSphere Application Server v.1.0.1 based on Apache Geronimo http://sel.as-us.falkag.net/sel?cmd=lnkkid=120709bid=263057dat=121642___ Numpy-discussion mailing list Numpy-discussion@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/numpy-discussion
Re: [Numpy-discussion] rcond in polyfit
On 10/13/06, A. M. Archibald [EMAIL PROTECTED] wrote:
On 13/10/06, Tim Hochberg [EMAIL PROTECTED] wrote: Charles R Harris wrote:snip
On the other hand if error handling is set to 'raise', then a FloatingPointError is issued. Is a FloatingPointWarning in order to mirror the FloatingPointError? And if so, would it be appropriate to use
for condition number?I submitted a patchto use warnings for several functions in scipy awhile ago, and the approach I took was to create a ScipyWarning, fromwhich more specific warnings were derived (IntegrationWarning, for
example). That was perhaps a bit short-sighted.I'd suggest a FloatingPointWarning as a base class, withIllConditionedMatrix as a subclass (it should include the conditionnumber, but probably not the matrix itself unless it's small, as
debugging information).Let's pin this down a bit. Numpy seems to need its own warning classes so we can control the printing of the warnings. For the polyfit function I think it should *always* warn by default because it is likely to be used interactively and one can fool around with the degree of the fit. For instance, I am now scaling the the input x vector so norm_inf(x) == 1 and this seems to work pretty well for lots of stuff with rcond=-1 (about 1e-16) and a warning on rank reduction. As long as the rank stays the same things seem to work ok, up to fits of degree 21 on the test data that started this conversation.
BTW, how does one turn warnings back on? If I do a warnings.simplefilter('always', mywarn)things work fine. Following this by warnings.simplefilter('once', mywarn)
does what is supposed to do. Once again issuing warnings.simplefilter('always', mywarn)fails to have any effect. Resetwarnings doesn't help. Hmmm...Chuck
-
Using Tomcat but need to do more? Need to support web services, security?
Get stuff done quickly with pre-integrated technology to make your job easier
Download IBM WebSphere Application Server v.1.0.1 based on Apache Geronimo
http://sel.as-us.falkag.net/sel?cmd=lnkkid=120709bid=263057dat=121642___
Numpy-discussion mailing list
Numpy-discussion@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/numpy-discussion
Re: [Numpy-discussion] rcond in polyfit
On 14/10/06, Charles R Harris [EMAIL PROTECTED] wrote:
On 10/13/06, A. M. Archibald [EMAIL PROTECTED] wrote:
On 13/10/06, Tim Hochberg [EMAIL PROTECTED] wrote:
Charles R Harris wrote:
snip
On the other hand if error handling is set to 'raise', then a
FloatingPointError is issued. Is a FloatingPointWarning in order to
mirror the FloatingPointError? And if so, would it be appropriate to use
for condition number?
I submitted a patchto use warnings for several functions in scipy a
while ago, and the approach I took was to create a ScipyWarning, from
which more specific warnings were derived (IntegrationWarning, for
example). That was perhaps a bit short-sighted.
I'd suggest a FloatingPointWarning as a base class, with
IllConditionedMatrix as a subclass (it should include the condition
number, but probably not the matrix itself unless it's small, as
debugging information).
Let's pin this down a bit. Numpy seems to need its own warning classes so we
can control the printing of the warnings. For the polyfit function I think
it should *always* warn by default because it is likely to be used
interactively and one can fool around with the degree of the fit. For
instance, I am now scaling the the input x vector so norm_inf(x) == 1 and
this seems to work pretty well for lots of stuff with rcond=-1 (about 1e-16)
and a warning on rank reduction. As long as the rank stays the same things
seem to work ok, up to fits of degree 21 on the test data that started this
conversation.
Numerical Recipes (http://www.nrbook.com/a/bookcpdf/c15-4.pdf )
recommend setting rcond to the number of data points times machine
epsilon (which of course is different for single/double). We should
definitely warn the user if any singular value is below s[0]*rcond (as
that means that there is effectively a useless basis function, up to
roundoff).
I'm not sure how to control the default warnings setting (once vs.
always); it's definitely not possible using the standard API to save
the warnings state and restore it later. One might be able to push
such a change into the warnings module by including it in a
ContextManager.
ipython should probably reset all the show once warnings every time
it shows an interactive prompt. I suppose more accurately, it should
do that only for warnings the user hasn't given instructions about.
That way you'll get a warning about bad polynomial fits every time you
run a command that contains one, but if your function runs thousands
of fits you don't drown in warnings.
BTW, how does one turn warnings back on? If I do a
warnings.simplefilter('always', mywarn)
things work fine. Following this by
warnings.simplefilter('once', mywarn)
does what is supposed to do. Once again issuing
warnings.simplefilter('always', mywarn)
fails to have any effect. Resetwarnings doesn't help. Hmmm...
I don't get the impression that the warnings module is much tested; I
had similar headaches.
A. M. Archibald
-
Using Tomcat but need to do more? Need to support web services, security?
Get stuff done quickly with pre-integrated technology to make your job easier
Download IBM WebSphere Application Server v.1.0.1 based on Apache Geronimo
http://sel.as-us.falkag.net/sel?cmd=lnkkid=120709bid=263057dat=121642
___
Numpy-discussion mailing list
Numpy-discussion@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/numpy-discussion
Re: [Numpy-discussion] rcond in polyfit
Charles R Harris wrote: On 10/14/06, *A. M. Archibald* [EMAIL PROTECTED] mailto:[EMAIL PROTECTED] wrote: [SNIP] Hmmm, I wonder if we have a dictionary of precisions indexed by dtype somewhere? Here's some code I stole from somewhere for computing EPS. It would easy enough to generate the dictionary you are looking for at startup using this. I can't recall the pedigree of this code though, so caveat emptor: def bits_of_precision(dtype): one = np.array([1], dtype) i = 0 while not np.alltrue(one + (one / 2.**i) == one): i += 1 return i - 1 EPSS = 1.0 / 2**bits_of_precision(float) * 10 # XXX safety factor It's sorta old and translated from numpy too, so it could probably be rewritten in better style. -tim [SNIP] - Using Tomcat but need to do more? Need to support web services, security? Get stuff done quickly with pre-integrated technology to make your job easier Download IBM WebSphere Application Server v.1.0.1 based on Apache Geronimo http://sel.as-us.falkag.net/sel?cmd=lnkkid=120709bid=263057dat=121642 ___ Numpy-discussion mailing list Numpy-discussion@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/numpy-discussion
Re: [Numpy-discussion] rcond in polyfit
On 10/14/06, Tim Hochberg [EMAIL PROTECTED] wrote: Charles R Harris wrote: On 10/14/06, *A. M. Archibald* [EMAIL PROTECTED] mailto:[EMAIL PROTECTED] wrote:[SNIP] Hmmm, I wonder if we have a dictionary of precisions indexed by dtype somewhere?Here's some code I stole from somewhere for computing EPS. It would easy enough to generate the dictionary you are looking for at startup usingthis. I can't recall the pedigree of this code though, so caveat emptor:def bits_of_precision(dtype): dtype) i = 0while not np.alltrue(one + (one / 2.**i) == one):i += 1return i - 1Yep, that works. There is a version in my zeros module in scipy also, it's just that ISTR seeing similar code in numpy somewhere and I was hoping Travis would tell me where ;) Grep is my friend, I guess, ... ah, here it is In [140]: np.MachAr(np.single).epsOut[140]: 1.1920928955078125e-07In [141]: np.MachAr(np.double).epsOut[141]: 2.2204460492503131e-16Chuck - Using Tomcat but need to do more? Need to support web services, security? Get stuff done quickly with pre-integrated technology to make your job easier Download IBM WebSphere Application Server v.1.0.1 based on Apache Geronimo http://sel.as-us.falkag.net/sel?cmd=lnkkid=120709bid=263057dat=121642___ Numpy-discussion mailing list Numpy-discussion@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/numpy-discussion
Re: [Numpy-discussion] rcond in polyfit
Charles R Harris wrote:
On 10/13/06, *A. M. Archibald* [EMAIL PROTECTED]
mailto:[EMAIL PROTECTED] wrote:
On 12/10/06, Charles R Harris [EMAIL PROTECTED]
mailto:[EMAIL PROTECTED] wrote:
Hi all,
I note that polyfit looks like it should work for single and
double, real
and complex, polynomials. On the otherhand, the default rcond
needs to
depend on the underlying precision. On the other, other hand,
all the svd
computations are done with dgelsd or zgelsd, i.e., double
precision. Even so
problems can arise from inherent errors of the input data if it
is single
precision to start with. I also think the final degree of the
fit should be
available somewhere if wanted, as it is an indication of what is
going on.
Sooo, any suggestions as to what to do? My initial impulse would
be to set
rcond=1e-6 for single, 1e-14 for double, make rcond a keyword,
and kick the
can down the road on returning the actual degree of the fit.
I'd also be inclined to output a warning (which the user can ignore,
read or trap as necessary) if the condition number is too bad or they
supplied an rcond that is too small for the precision of their data.
That sounds good, but how to do it? Should I raise an exception?
Use the warnings framework:
import warnings
warnings.warn(condition number is BAD)
__main__:1: UserWarning: condition number is BAD
The user can turn warnings on or off or turned in exceptions based on a
variety of criteria. Look for the warnings filter in the docs.
Which brings up a question: do we want to have a FloatingPointWarning or
some such? Currently, if you use set the error handling to warn using
seterr a runtime warning is issued:
np.seterr(all='warn')
{'over': 'ignore', 'divide': 'ignore', 'invalid': 'ignore', 'under':
'ignore'}
np.arange(1) / 0
__main__:1: RuntimeWarning: divide by zero encountered in divide
On the other hand if error handling is set to 'raise', then a
FloatingPointError is issued. Is a FloatingPointWarning in order to
mirror the FloatingPointError? And if so, would it be appropriate to use
for condition number?
I would also have to modify lstsq so it returns the degree of the fit
which would mess up the current interface.
One approach would be to write lstsqcond (or a better name) that returns
both the fit and the condition number. listsq could then be just a
wrapper over that which dumped the condition number. IIRC, the
condition number is available, but we're not returning it.
-tim
-
Using Tomcat but need to do more? Need to support web services, security?
Get stuff done quickly with pre-integrated technology to make your job easier
Download IBM WebSphere Application Server v.1.0.1 based on Apache Geronimo
http://sel.as-us.falkag.net/sel?cmd=lnkkid=120709bid=263057dat=121642
___
Numpy-discussion mailing list
Numpy-discussion@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/numpy-discussion
