85 views

Skip to first unread message

Apr 20, 2017, 5:10:29 PM4/20/17

to sage-devel

polynomials with numpy coefficients strike again, this time if you build Sage with clang:

sage: import numpy as np

sage: x=polygen(RR) # this is problematic

sage: y=polygen(RDF) # this works fine

sage: np.float('1.5')*x

1.50000000000000*x

sage: np.float('1.5')*y

1.5*x

sage: np.float32('1.5')*x

/usr/home/dima/Sage/sage/src/bin/sage-ipython:1: RuntimeWarning: invalid value encountered in multiply

#!/usr/bin/env python

1.50000000000000*x

sage: np.float64('1.5')*x

/usr/home/dima/Sage/sage/src/bin/sage-ipython:1: RuntimeWarning: invalid value encountered in multiply

#!/usr/bin/env python

1.50000000000000*x

sage: np.float64('1.5')*y

1.5*x

sage: np.float32('1.5')*y

1.5*x

Finding out where exactly in numpy this warning happens is not easy (one can set up a numpy interrupt handler, but because it happens in compiled code, error message is almost meaningless)

but it looks as if it might be a coercion problem.

Any ideas where to look?

Apr 21, 2017, 4:11:37 AM4/21/17

to sage-...@googlegroups.com

Dima Pasechnik wrote:

> but it looks as if it might be a coercion problem.

> Any ideas where to look?

Not really, but it does look like the common parent
> but it looks as if it might be a coercion problem.

> Any ideas where to look?

discovered by the coercion system is incorrect:

sage: import numpy as np

sage: b = np.float32('1.5')

sage: get_coercion_model().common_parent(b, polygen(RR))

Univariate Polynomial Ring in x over Real Field with 53 bits of

precision

sage: RR.coerce(a)

1.50000000000000

sage: RR.coerce(b)

...

TypeError: no canonical coercion from <type 'numpy.float32'> to Real

Field with 53 bits of precision

sage: get_coercion_model().common_parent(b, RR)

<type 'numpy.float32'>

--

Marc

Apr 21, 2017, 5:10:12 AM4/21/17

to sage-devel, ma...@mezzarobba.net

This does not look like this is the root cause of the problem.

Namely, note that also

sage: a128 = np.float128('1.5')

sage: RR.coerce(a128)

# throws the same as above TypeError, but

sage: a128*x

# does not print the numpy warning.

what is also somewhat puzzling is

that while

sage: np.float32('1.5')*x

# does print the numpy warning

sage: x*np.float32('1.5')

# just works (i.e. does not print the numpy warning)

So this is a subtle combination of an apparent bug in numpy (or in clang) with

some Sage coercion weirdness.

Apparently numpy's floats are coerced into RDF for the purpose of

dealing with polynomials, see py_scalar_to_element in structure/coerce.pyx, but apparently this only happens

for

x*np.float32('1.5'), and not for np.float32('1.5')*x ?

Dima

--

Marc

Apr 21, 2017, 10:46:03 AM4/21/17

to sage-devel, ma...@mezzarobba.net

Quick guess (don't have time to look at this soon): the issue could be in the action discovery code, which may try creating elements of numpy.float32 using weird inputs. This could cause a difference between right and left multiplication. One way to check if this is the source is to try converting each of the polynomials from RR['x'].some_elements() into numpy.float32 and see if you get the same error.

--

You received this message because you are subscribed to the Google Groups "sage-devel" group.

To unsubscribe from this group and stop receiving emails from it, send an email to sage-devel+unsubscribe@googlegroups.com.

To post to this group, send email to sage-...@googlegroups.com.

Visit this group at https://groups.google.com/group/sage-devel.

For more options, visit https://groups.google.com/d/optout.

Apr 21, 2017, 11:29:26 AM4/21/17

to sage-devel, ma...@mezzarobba.net

On Friday, April 21, 2017 at 3:46:03 PM UTC+1, David Roe wrote:

Quick guess (don't have time to look at this soon): the issue could be in the action discovery code, which may try creating elements of numpy.float32 using weird inputs. This could cause a difference between right and left multiplication. One way to check if this is the source is to try converting each of the polynomials from RR['x'].some_elements() into numpy.float32 and see if you get the same error.

no, this does not give anything that looks like the error in question (no errors on degree 0 terms,

"ValueError: setting an array element with a sequence." on degree >0 terms)

To unsubscribe from this group and stop receiving emails from it, send an email to sage-devel+...@googlegroups.com.

Apr 21, 2017, 9:10:48 PM4/21/17

to sage-devel, ma...@mezzarobba.net

On Fri, Apr 21, 2017 at 11:29 AM, Dima Pasechnik <dim...@gmail.com> wrote:

On Friday, April 21, 2017 at 3:46:03 PM UTC+1, David Roe wrote:Quick guess (don't have time to look at this soon): the issue could be in the action discovery code, which may try creating elements of numpy.float32 using weird inputs. This could cause a difference between right and left multiplication. One way to check if this is the source is to try converting each of the polynomials from RR['x'].some_elements() into numpy.float32 and see if you get the same error.no, this does not give anything that looks like the error in question (no errors on degree 0 terms,"ValueError: setting an array element with a sequence." on degree >0 terms)

Sorry, I don't have any other ideas, and can't play around with it easily since the standard Sage build isn't producing these errors.

David

To unsubscribe from this group and stop receiving emails from it, send an email to sage-devel+unsubscribe@googlegroups.com.

Apr 22, 2017, 4:22:21 AM4/22/17

to sage-devel

it is not hard to build with clang on OSX, see #12426.

and this issue is causing the majority of doctest failures in this setting, as well as on clang/FreeBSD.(#22679).

(where the building is not as smooth yet).

I suppose we will have to merge these clang-enabling changes one way or another soon.

Apr 22, 2017, 5:01:52 AM4/22/17

to sage-devel, ma...@mezzarobba.net, roed...@gmail.com

On Saturday, April 22, 2017 at 2:10:48 AM UTC+1, David Roe wrote:

On Fri, Apr 21, 2017 at 11:29 AM, Dima Pasechnik <dim...@gmail.com> wrote:

On Friday, April 21, 2017 at 3:46:03 PM UTC+1, David Roe wrote:no, this does not give anything that looks like the error in question (no errors on degree 0 terms,"ValueError: setting an array element with a sequence." on degree >0 terms)Sorry, I don't have any other ideas, and can't play around with it easily since the standard Sage build isn't producing these errors.David

Is there a "normal" way to trace all the coersion/conversion steps for a given Sage expression?

If not, there should be one, IMHO (otherwise it's always a huge guessing game for this sort of issues)

Apr 22, 2017, 5:47:48 AM4/22/17

to sage-devel, ma...@mezzarobba.net, roed...@gmail.com

There is this:

sage: import numpy as np

sage: cm = sage.structure.element.coercion_model

sage: cm.explain(polygen(RR), np.float32('1.5'), operator.mul)

Action discovered.

Right scalar multiplication by Real Double Field on Univariate Polynomial Ring in x over Real Field with 53 bits of precision

with precomposition on right by Native morphism:

From: Set of Python objects of type 'numpy.float32'

To: Real Double Field

Result lives in Univariate Polynomial Ring in x over Real Field with 53 bits of precision

Univariate Polynomial Ring in x over Real Field with 53 bits of precision

Apr 22, 2017, 5:50:31 AM4/22/17

to sage-...@googlegroups.com

On 2017-04-22 11:01, Dima Pasechnik wrote:

> Is there a "normal" way to trace all the coersion/conversion steps for a

> given Sage expression?

> If not, there should be one, IMHO (otherwise it's always a huge guessing

> game for this sort of issues)

What helps is .record_exceptions and .coercion_traceback:
> Is there a "normal" way to trace all the coersion/conversion steps for a

> given Sage expression?

> If not, there should be one, IMHO (otherwise it's always a huge guessing

> game for this sort of issues)

http://doc.sagemath.org/html/en/reference/structure/sage/structure/element.html#sage.structure.element.coercion_traceback

Daniel

Apr 22, 2017, 6:22:59 PM4/22/17

to sage-devel, ma...@mezzarobba.net, roed...@gmail.com

On Saturday, April 22, 2017 at 10:47:48 AM UTC+1, Volker Braun wrote:

There is this:sage: import numpy as npsage: cm = sage.structure.element.coercion_modelsage: cm.explain(polygen(RR), np.float32('1.5'), operator.mul)Action discovered.Right scalar multiplication by Real Double Field on Univariate Polynomial Ring in x over Real Field with 53 bits of precisionwith precomposition on right by Native morphism:From: Set of Python objects of type 'numpy.float32'To: Real Double FieldResult lives in Univariate Polynomial Ring in x over Real Field with 53 bits of precisionUnivariate Polynomial Ring in x over Real Field with 53 bits of precision

Thanks! It seems that cm.explain() lies in this case, i.e., the action is actually computed in a different way!!!

Indeed:

sage: import numpy as np

sage: x=polygen(RealField())

sage: np.float32('1.5')*x

/usr/home/dima/Sage/sage/src/bin/sage-ipython:1: RuntimeWarning: invalid value encountered in multiply

#!/usr/bin/env python

1.50000000000000*x

sage: cm = sage.structure.element.get_coercion_model()

sage: cm.explain(np.float32('1.5'), x, operator.mul)

Action discovered.

Left scalar multiplication by Real Double Field on Univariate Polynomial Ring in x over Real Field with 53 bits of precision

with precomposition on left by Native morphism:

From: Set of Python objects of type 'numpy.float32'

To: Real Double Field

Result lives in Univariate Polynomial Ring in x over Real Field with 53 bits of precision

Univariate Polynomial Ring in x over Real Field with 53 bits of precision

sage: RDF(np.float32('1.5'))*x

1.50000000000000*x

that is, cm.explain tells me that np.float32('1.5')*x is computed by doing

RDF(np.float32('1.5'))*x

But this is not the case, as RDF(np.float32('1.5'))*x does not produce any numpy warnings, whereas

np.float32('1.5')*x does.

Then I went digging for the place which is supposed to do apply RDF(), and tried to break it:

--- a/src/sage/structure/coerce.pyx

+++ b/src/sage/structure/coerce.pyx

@@ -231,7 +231,8 @@ cpdef py_scalar_to_element(x):

return Integer(x)

elif isinstance(x, numpy.floating):

from sage.rings.real_double import RDF

- return RDF(x)

+ # return RDF(x)

+ return RDF(42)

elif isinstance(x, numpy.complexfloating):

from sage.rings.complex_double import CDF

return CDF(x)

but to no avail, 1.5 does not get replaced by 42.0 in the answer.

It seems that this part of coersion is broken in the sense that it does something rogue,

not prescribed by the framework---and it has nothing to do with clang (I also tried this "patch" in

a usual configuration).

Reply all

Reply to author

Forward

0 new messages

Search

Clear search

Close search

Google apps

Main menu