Re: [Numpy-discussion] rcond in polyfit

2006-10-15 Thread Charles R Harris 
On 10/14/06, A. M. Archibald <[EMAIL PROTECTED]> wrote:
On 14/10/06, Charles R Harris <[EMAIL PROTECTED]> wrote: 
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=lnk&kid=120709&bid=263057&dat=121642___
Numpy-discussion mailing list
Numpy-discussion@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/numpy-discussion

Re: [Numpy-discussion] rcond in polyfit

2006-10-14 Thread Charles R Harris 
On 10/14/06, Tim Hochberg <[EMAIL PROTECTED]> wrote:
Charles R Harris wrote:>>> On 10/14/06, *A. M. Archibald* <[EMAIL PROTECTED]> [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=lnk&kid=120709&bid=263057&dat=121642___
Numpy-discussion mailing list
Numpy-discussion@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/numpy-discussion

Re: [Numpy-discussion] rcond in polyfit

2006-10-14 Thread Tim Hochberg 
Charles R Harris wrote:
>
>
> On 10/14/06, *A. M. Archibald* <[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=lnk&kid=120709&bid=263057&dat=121642
___
Numpy-discussion mailing list
Numpy-discussion@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/numpy-discussion

Re: [Numpy-discussion] rcond in polyfit

2006-10-14 Thread Charles R Harris 
On 10/14/06, A. M. Archibald <[EMAIL PROTECTED]> wrote:
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:> 
Numerical Recipes (http://www.nrbook.com/a/bookcpdf/c15-4.pdf
 )recommend setting rcond to the number of data points times machineepsilon (which of course is different for single/double). We shoulddefinitely warn the user if any singular value is below s[0]*rcond (as
that means that there is effectively a useless basis function, up toroundoff).Well, that would work. On the other hand, it seems overly pessimistic from my experiments here. What seems to be a better guide is rank reduction. For instance, I can do a perfectly decent fit to Greg's data with single precision, degree 10,  and ~1000 data points with rcond = ~5e-7 (effectively single precision precision). Degree 11 blows up entirely even though it loses rank. It is also true that rcond=1e-3 fails for degree 11 even though rank is strongly reduced. What looks to be taking place is roundoff error in evaluating the polynomial in single precision in the presence of higher order terms that still belong to the reduced basis functions and rank reduction is a good indicator of this. Bear in mind that I am now normalizing x by dividing it by its largest element, which futzes with the condition number. The condition number of the unscaled fit doesn't bear thinking about.
Hmmm, I wonder if we have a dictionary of precisions indexed by dtype somewhere? 
I don't get the impression that the warnings module is much tested; Ihad similar headaches.I see some folks bitchin' at it when I google. There are remarkably few hits, though. For the moment I am going with the smaller rcond numbers and raising an error on rank reduction. I suspect something similar should be done for pinv.
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=lnk&kid=120709&bid=263057&dat=121642___
Numpy-discussion mailing list
Numpy-discussion@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/numpy-discussion

Re: [Numpy-discussion] rcond in polyfit

2006-10-14 Thread A. M. Archibald 
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:
>
> 
>
> > > 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=lnk&kid=120709&bid=263057&dat=121642
___
Numpy-discussion mailing list
Numpy-discussion@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/numpy-discussion

Re: [Numpy-discussion] rcond in polyfit

2006-10-14 Thread Charles R Harris 
On 10/13/06, A. M. Archibald <[EMAIL PROTECTED]> wrote:
On 13/10/06, Tim Hochberg <[EMAIL PROTECTED]> wrote:> Charles R Harris wrote: 
> 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=lnk&kid=120709&bid=263057&dat=121642___
Numpy-discussion mailing list
Numpy-discussion@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/numpy-discussion

Re: [Numpy-discussion] rcond in polyfit

2006-10-13 Thread A. M. Archibald 
On 13/10/06, Tim Hochberg <[EMAIL PROTECTED]> wrote:
> Charles R Harris wrote:

> > 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 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).

The warnings module is frustratingly non-reentrant, unfortunately,
which makes writing tests very awkward.

> > 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.

This is a very good idea. scipy.integrate.quad returns a pair (result,
error_estimate) and every time I use it I trip over that. (Perhaps if
I were a fine upstanding numerical analyst I would be checking the
error estimate every time, but it is a pain.) Another option would be
a "full_output" optional argument.

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=lnk&kid=120709&bid=263057&dat=121642
___
Numpy-discussion mailing list
Numpy-discussion@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/numpy-discussion

Re: [Numpy-discussion] rcond in polyfit

2006-10-13 Thread Tim Hochberg 
Charles R Harris wrote:
>
>
> On 10/13/06, *A. M. Archibald* <[EMAIL PROTECTED] 
> > wrote:
>
> On 12/10/06, Charles R Harris <[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=lnk&kid=120709&bid=263057&dat=121642
___
Numpy-discussion mailing list
Numpy-discussion@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/numpy-discussion

Re: [Numpy-discussion] rcond in polyfit

2006-10-13 Thread A. M. Archibald 
On 13/10/06, Charles R Harris <[EMAIL PROTECTED]> wrote:
>
>
>
> On 10/13/06, A. M. Archibald <[EMAIL PROTECTED]> wrote:
> > On 12/10/06, Charles R Harris <[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? I would
> also have to modify lstsq so it returns the degree of the fit which would
> mess up the current  interface.

Python's warnings module is a decent solution for providing this
information. Goodness-of-fit worries me less than ill-conditioning -
users are going to expect the curve to deviate from their function
(and an easy reliable way to get goodness of fit is
sqrt(sum(abs(f(xs)-polynomial(xs))**2)); this is certain to take into
account any roundoff errors introduced anywhere). But they may well
have no idea they should be worried about the condition number of some
matrix they've never heard of.

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=lnk&kid=120709&bid=263057&dat=121642
___
Numpy-discussion mailing list
Numpy-discussion@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/numpy-discussion

Re: [Numpy-discussion] rcond in polyfit

2006-10-13 Thread Charles R Harris 
On 10/13/06, A. M. Archibald <[EMAIL PROTECTED]> wrote:
On 12/10/06, Charles R Harris <[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 theysupplied 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? I would also have to modify lstsq so it returns the degree of the fit which would mess up the current  interface. 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=lnk&kid=120709&bid=263057&dat=121642___
Numpy-discussion mailing list
Numpy-discussion@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/numpy-discussion

Re: [Numpy-discussion] rcond in polyfit

2006-10-13 Thread A. M. Archibald 
On 12/10/06, Charles R Harris <[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.

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=lnk&kid=120709&bid=263057&dat=121642
___
Numpy-discussion mailing list
Numpy-discussion@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/numpy-discussion

[Numpy-discussion] rcond in polyfit

2006-10-12 Thread Charles R Harris 
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.
ChuckPS, what is the best way of converting arbitrary arrays to the appropriate c data types float and double?
-
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=lnk&kid=120709&bid=263057&dat=121642___
Numpy-discussion mailing list
Numpy-discussion@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/numpy-discussion