RR, coercion of numpy.float32() and clang - coercion experts needed

[Dima Pasechnik]

see https://trac.sagemath.org/ticket/22799 :

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?

[Marc Mezzarobba]

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
discovered by the coercion system is incorrect:

sage: import numpy as np

sage: a = np.float('1.5')
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

[Dima Pasechnik]

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

[David Roe]

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.

[Dima Pasechnik]

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.

[David Roe]

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.

[Dima Pasechnik]

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.

[Dima Pasechnik]

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:

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

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)

[Volker Braun]

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

[Daniel Krenn]

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:

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

Daniel

[Dima Pasechnik]

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

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

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