#20438: Wrapper_cdf.call_c() drops imaginary part
-------------------------------------------------+-------------------------
Reporter: cswiercz | Owner:
Type: defect | Status: new
Priority: major | Milestone: sage-7.2
Component: cython | Resolution:
Keywords: fast_callable, interpreters, | Merged in:
wrapper_cdf | Reviewers:
Authors: cswiercz | Work issues:
Report Upstream: N/A | Commit:
Branch: | Stopgaps:
Dependencies: |
-------------------------------------------------+-------------------------
Description changed by jdemeyer:
Old description:
> A `Wrapper_cdf` object is created by `fast_callable(f, vars=[],
> domain=None)` when the domain is set to `CDF`. The standard method of
> evaluating the function represented by the resulting object by simple
> calling it behaves normally.
>
> However, the Cython-accesible method `Wrapper_cdf.call_c(double_complex*
> out, double_complex* in)` drops the imaginary component of the input and
> output. A demonstration of this behavior:
> {{{
> cython("""
> from sage.ext.interpreters.wrapper_cdf cimport Wrapper_cdf
>
> # cannot do below: it redefines double_complex from the
> # auto-generated headers created by ext.interpreters.py. generates
> # cython error:
> #
> # "Cannot assign type 'double complex [2]' to 'double_complex *'"
> #
> #ctypedef complex double_complex "double complex"
>
> # Instead, use the (incorrect) definition as given in the auto-gen
> # wrappers. This ctypedefs "double_complex" as a "double"! That is,
> # all complex parts of calculations are lost when invoking
> # Wrapper_cdf.call_c()!
> #
> # The following is the self-proclaimed "hack" in the auto-gened file
> # ext/interpreters/wrapper_cdf.pyx
> cdef extern from "complex.h":
> ctypedef double double_complex "double complex" # cannot just change
> "double" to "complex"
> cdef double creal(double_complex)
> cdef double cimag(double_complex)
> cdef double_complex _Complex_I
>
> cdef inline complex call_c(Wrapper_cdf f, complex x, complex y):
> # Quickly evaluates a fast_callable function `f` on domain `CDF`
> # (a Wrapper_cdf tupe) with complex arguments `x`, and `y`.
>
> cdef double_complex[2] input
> cdef double_complex* output = [0]
> cdef complex value
> input[0] = <double_complex>x
> input[1] = <double_complex>y
> f.call_c(input, output)
> value = <complex>(output[0])
> return value
>
> from sage.all import QQ, CDF, fast_callable
>
> R = QQ['x,y']
> x,y = R.gens()
> f = y**3 - 2*x**3*y + x**7
> f_fast = fast_callable(f, vars=[x,y], domain=CDF)
>
> cdef complex a = 1. + 1.j
> cdef complex b = 1. - 2.j
>
> cdef complex value_sage = f_fast(a,b)
> cdef complex value_c = call_c(f_fast,a,b)
>
> print 'sage value: ', value_sage
> print 'call_c value:', value_c
> """)
>
> sage value: (-7-18j)
> call_c value: (-7+0j)
> }}}
>
> The source code for `Wrapper_cdf` is generated by
> `sage/src/sage_setup/autogen/interpreters.py`. Within, there are some
> peculiar comments when it comes to defining the custom data type
> `double_complex`. The following is from
> [https://github.com/sagemath/sage/blob/master/src/sage_setup/autogen/interpreters.py#L2533
> Line 2533] where the Cython header for `Wrapper_cdf` is defined.
>
> {{{
> # We need the type double_complex to work around
> # http://trac.cython.org/ticket/869
> # so this is a bit hackish.
> cdef extern from "complex.h":
> ctypedef double double_complex "double complex"
> }}}
>
> But this seems to define `double_complex` as a `double`, which explains
> why the imaginary part is dropped. Since complex externs are not allowed
> in Cython(?) I tried changing this line to
>
> {{{
> ctypedef double complex double_complex "double complex"
> }}}
>
> but the following compilation error is raised
>
> {{{
> [1/3] Cythonizing sage/ext/interpreters/wrapper_cdf.pyx
>
> Error compiling Cython file:
> ------------------------------------------------------------
> ...
> cdef extern from "interp_cdf.h":
> double_complex interp_cdf(double_complex* args,
> double_complex* constants,
> PyObject** py_constants,
> double_complex* stack,
> int* code) except? -1094648119105371
> ^
> ------------------------------------------------------------
>
> sage/ext/interpreters/wrapper_cdf.pyx:64:28: Not allowed in a constant
> expression
> }}}
>
> This part of the wrapper code is generated by
> [https://github.com/sagemath/sage/blob/master/src/sage_setup/autogen/interpreters.py#L3290
> interpreters.py:3290]. Commenting out this `except?` part resolves the
> compile issue as well as the missing imaginary part problems. Attached is
> a patch which makes the appropriate `double_complex` fixes as well as,
> probably erroneously, comments out the except so that the code compiles.
>
> (''Motivation: I would like to directly call the underlying C function of
> a `Wrapper_cdf` for performance purposes.'')
New description:
A `Wrapper_cdf` object is created by `fast_callable(f, vars=[],
domain=None)` when the domain is set to `CDF`. The standard method of
evaluating the function represented by the resulting object by simple
calling it behaves normally.
However, the Cython-accesible method `Wrapper_cdf.call_c(double_complex*
out, double_complex* in)` drops the imaginary component of the input and
output. A demonstration of this behavior:
{{{
cython("""
from sage.ext.interpreters.wrapper_cdf cimport Wrapper_cdf
# The following is the self-proclaimed "hack" in the auto-gened file
# ext/interpreters/wrapper_cdf.pyx
cdef extern from "complex.h":
ctypedef double double_complex "double complex" # cannot just change
"double" to "complex"
cdef double creal(double_complex)
cdef double cimag(double_complex)
cdef double_complex _Complex_I
cdef inline complex call_c(Wrapper_cdf f, complex x, complex y):
# Quickly evaluates a fast_callable function `f` on domain `CDF`
# (a Wrapper_cdf tupe) with complex arguments `x`, and `y`.
cdef double_complex[2] input
cdef double_complex* output = [0]
cdef complex value
input[0] = <double_complex>x
input[1] = <double_complex>y
f.call_c(input, output)
value = <complex>(output[0])
return value
from sage.all import QQ, CDF, fast_callable
R = QQ['x,y']
x,y = R.gens()
f = y**3 - 2*x**3*y + x**7
f_fast = fast_callable(f, vars=[x,y], domain=CDF)
cdef complex a = 1. + 1.j
cdef complex b = 1. - 2.j
cdef complex value_sage = f_fast(a,b)
cdef complex value_c = call_c(f_fast,a,b)
print 'sage value: ', value_sage
print 'call_c value:', value_c
""")
sage value: (-7-18j)
call_c value: (-7+0j)
}}}
The source code for `Wrapper_cdf` is generated by
`sage/src/sage_setup/autogen/interpreters.py`. Within, there are some
peculiar comments when it comes to defining the custom data type
`double_complex`. The following is from
[https://github.com/sagemath/sage/blob/master/src/sage_setup/autogen/interpreters.py#L2533
Line 2533] where the Cython header for `Wrapper_cdf` is defined.
{{{
# We need the type double_complex to work around
# http://trac.cython.org/ticket/869
# so this is a bit hackish.
cdef extern from "complex.h":
ctypedef double double_complex "double complex"
}}}
But this seems to define `double_complex` as a `double`, which explains
why the imaginary part is dropped. Since complex externs are not allowed
in Cython(?) I tried changing this line to
{{{
ctypedef double complex double_complex "double complex"
}}}
but the following compilation error is raised
{{{
[1/3] Cythonizing sage/ext/interpreters/wrapper_cdf.pyx
Error compiling Cython file:
------------------------------------------------------------
...
cdef extern from "interp_cdf.h":
double_complex interp_cdf(double_complex* args,
double_complex* constants,
PyObject** py_constants,
double_complex* stack,
int* code) except? -1094648119105371
^
------------------------------------------------------------
sage/ext/interpreters/wrapper_cdf.pyx:64:28: Not allowed in a constant
expression
}}}
This part of the wrapper code is generated by
[https://github.com/sagemath/sage/blob/master/src/sage_setup/autogen/interpreters.py#L3290
interpreters.py:3290]. Commenting out this `except?` part resolves the
compile issue as well as the missing imaginary part problems. Attached is
a patch which makes the appropriate `double_complex` fixes as well as,
probably erroneously, comments out the except so that the code compiles.
(''Motivation: I would like to directly call the underlying C function of
a `Wrapper_cdf` for performance purposes.'')
--
--
Ticket URL: <http://trac.sagemath.org/ticket/20438#comment:4>
Sage <http://www.sagemath.org>
Sage: Creating a Viable Open Source Alternative to Magma, Maple, Mathematica,
and MATLAB
--
You received this message because you are subscribed to the Google Groups
"sage-trac" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To post to this group, send email to [email protected].
Visit this group at https://groups.google.com/group/sage-trac.
For more options, visit https://groups.google.com/d/optout.
