Mailman 3 February 2013 - cython-devel

Re: [Cython] Bringing Cython and PyPy closer together
by Stefan Behnel March 29, 2015

March 29, 2015

[copied here from PyPy mailing list] Stefan Behnel, 15.02.2012 12:32: > The current state of the discussion seems to be that PyPy provides ways to > talk to C code, but nothing as complete as CPython's C-API in the sense > that it allows efficient two-way communication between C code and Python > objects. Thus, we need to either improve this or look for alternatives. > > In order to get us more focussed on what can be done and what the > implications are, so that we may eventually be able to decide what should > be done, I started a Wiki page for a PyPy backend CEP (Cython Enhancement > Proposal). > > http://wiki.cython.org/enhancements/pypy The discussion so far makes me rather certain that the most promising short-term solution is to make Cython generate C code that PyPy's cpyext can handle. This should get us a rather broad set of running code somewhat quickly, while requiring the least design-from-scratch type of work in a direction that does not yet allow us to see if it will really make existing code work or not. On top of the basic cpyext interface, it should then be easy to implement obvious optimisations like native C level calls to Cython wrapped functions from PyPy (and potentially also the other direction) and otherwise avoid boxing/unboxing where unnecessary, e.g. for builtins. After all, it all boils down to native code at some point and I'm sure there are various ways to exploit that. Also, going this route will help both projects to get to know each other better. I think that's a required basis if we really aim for designing a more high-level interface at some point. The first steps I see are: - get Cython's test suite to run on PyPy - analyse the failing tests and decide how to fix them - adapt the Cython generated C code accordingly, special casing for PyPy where required Here is a "getting started" guide that tells you how testing works in Cython: http://wiki.cython.org/HackerGuide Once we have the test suite runnable, we can set up a PyPy instance on our CI server to get feed-back on any advances. https://sage.math.washington.edu:8091/hudson/ So, any volunteers or otherwise interested parties to help in getting this to work? Anyone in for financial support? Stefan

2 9

PR on refcounting memoryview buffers
by Sturla Molden March 4, 2013

March 4, 2013

As Stefan suggested, I have posted a PR for a better fix for the issue when MinGW for some reason emits the symbol "__synch_fetch_and_add_4" instead of generating atomic opcode for the __synch_fetch_and_add builtin. The PR is here: https://github.com/cython/cython/pull/185 The discussion probably belongs on this list instead og Cython user: The problem this addresses is when GCC does not use atomic builtins and emits __synch_fetch_and_add_4 and __synch_fetch_and_sub_4 when Cython are internally refcounting memoryview buffers. For some reason it can even happen on x86 and amd64. My PR undos Marks quick fix that always uses PyThread_acquire_lock on MinGW. PyThread_acquire_lock uses a kernel object (semaphore) on Windows and is not very efficient. I want slicing memoryviews to be fast, and that means PyThread_acquire_lock must go. My PR uses Windows API atomic function InterlockedAdd to implement the semantics of __synch_fetch_and_add_4 and __synch_fetch_and_sub_4 instead of using a Python lock. Usually MinGW is configured to compile GNU atomic builtins correctly. I have yet to see a case where it is not. But obviously one user (JF Gallant) has encountered it. I don't think it is a MinGW specific problem, but currently it has only been seen on MinGW and the fix is MinGW specific (well, it should work on Cygwin too). But whenever MinGW does use atomic builtins it just uses them. So it incurs no speed penalty on well-behaved MinGW builds. I took the liberty to use GNU extensions __inline__ and __attribute(always_inline)__. They will make sure the functions always behave like macros. The rationale being that it is GCC specific code so we can assume GNU extensions are available. If we take them away the code should still work, but we have no guarantee the functions will be inlined. I did not use macros because __synch_fetch_and_add is emitted by the preprocessor, and thus GCC will presumably emit __synch_fetch_and_sub_4 after the preprocessing step, which could require __synch_fetch_and_sub_4 to be a function instead of another macro. (I have no way of finding it out since I cannot test for it.) Regarding Linux and OSX: Failure of GCC to use atomic builtins could also happen on other GCC builds though. I don't think it is a MinGW-only issue. It's probably due to how the GCC build was configured. So we should as a safeguard have this for other OSes too. http://developer.apple.com/library/ios/#DOCUMENTATION/System/Conceptual/Man… We probably just need similar code to what I wrote for MinGW. I can write the code, but I don't have a Mac on which to test it. Also we should use OSAtomic* on clang/LLVM, which is now the platform C compiler on OSX. This will avoid PyThread_acquire_lock being the common synch mechanism for refcounting memoryview buffers on OSX. On Linux I am not sure what to suggest if GCC fails to use atomic builtins. I can handcode inline assembly for x86/amd64. I could also use pthreads and pth threads locks. But we could also assume that it never happen and just let the linker fail on __synch_fetch_and_add_4. Sturla

1 2

To Add datetime.pxd to cython.cpython
by ZS March 4, 2013

March 4, 2013

Extended part is in datetime_ex.h: #include "datetime.h" #define PyDateTime_FromDateAndTimeEx(year, month, day, hour, min, sec, usec, tzinfo) \ PyDateTimeAPI->DateTime_FromDateAndTime(year, month, day, hour, \ min, sec, usec, tzinfo, PyDateTimeAPI->DateTimeType) #define PyTime_FromTimeEx(hour, minute, second, usecond, tzinfo) \ PyDateTimeAPI->Time_FromTime(hour, minute, second, usecond, \ tzinfo, PyDateTimeAPI->TimeType) These macros allow to create dattime/time objects with tzinfo. Of course we could do: t = PyTime_FromTime(........) t = t.replace(tzinfo) in absence of that. Zaur Shibzukhov

4 7

About IndexNode and unicode[index]
by ZS March 3, 2013

March 3, 2013

Looking into IndexNode class in ExprNode.py I have seen a possibility for addition of more fast code path for unicode[index] as it done in method `generate_setitem_code` in case of lists. This is files for evaluation of performance difference: #### unicode_index.h /* This is striped version of __Pyx_GetItemInt_Unicode_Fast */ #include "unicodeobject.h" static inline Py_UCS4 unicode_char(PyObject* ustring, Py_ssize_t i); static inline Py_UCS4 unicode_char(PyObject* ustring, Py_ssize_t i) { #if CYTHON_PEP393_ENABLED if (PyUnicode_READY(ustring) < 0) return (Py_UCS4)-1; #endif return __Pyx_PyUnicode_READ_CHAR(ustring, i); } ##### unicode_index.pyx # coding: utf-8 cdef extern from 'unicode_index.h': inline Py_UCS4 unicode_char(unicode ustring, int i) cdef unicode text = u"abcdefghigklmnopqrstuvwxyzabcdefghigklmnopqrstuvwxyz" def f_1(unicode text): cdef int i, j cdef int n = len(text) cdef Py_UCS4 ch for j from 0<=j<=1000000: for i from 0<=i<=n-1: ch = text[i] def f_2(unicode text): cdef int i, j cdef int n = len(text) cdef Py_UCS4 ch for j from 0<=j<=1000000: for i from 0<=i<=n-1: ch = unicode_char(text, i) def test_1(): f_1(text) def test_2(): f_2(text) Timing results: (py33) zbook:mytests $ python3.3 -m timeit -n 100 -r 10 -s "from mytests.unicode_index import test_1" "test_1()" 100 loops, best of 10: 89 msec per loop (py33) zbook:mytests $ python3.3 -m timeit -n 100 -r 10 -s "from mytests.unicode_index import test_2" "test_2()" 100 loops, best of 10: 46.1 msec per loop in setup.py globally: "boundscheck": False "wraparound": False "nonecheck": False Zaur Shibzukhov

4 21

Be more forgiving about memoryview strides
by Sebastian Berg March 1, 2013

March 1, 2013

Hey, Maybe someone here already saw it (I don't have a track account, or I would just create a ticket), but it would be nice if Cython was more forgiving about contiguous requirements on strides. In the future this would make it easier for numpy to go forward with changing the contiguous flags to be more reasonable for its purpose, and second also to allow old (and maybe for the moment remaining) corner cases in numpy to slip past (as well as possibly the same for other programs...). An example is (see also https://github.com/numpy/numpy/issues/2956 and the PR linked there for more details): def add_one(array): cdef double[::1] a = array a[0] += 1. return array giving: >>> add_one(np.ascontiguousarray(np.arange(10.)[::100])) ValueError: Buffer and memoryview are not contiguous in the same dimension. This could easily be changed if MemoryViews check the strides as "can be interpreted as contiguous". That means that if shape[i] == 1, then strides[i] are arbitrary (you can just change them if you like). This is also the case for 0-sized arrays, which are arguably always contiguous, no matter their strides are! Regards, Sebastian PS: A similar thing exists with np.ndarray[...] interface if the user accesses array.strides. They get the arrays strides not the buffers. This is not quite related, but if it would be easy to use the buffer's strides in that case, it may make it easier if we want to change the flags in numpy in the long term, since one could clean up strides for forced contiguous buffer requests.

3 6

Class methods returning C++ class references are not dealt with correctly?
by Yury V. Zaytsev Feb. 28, 2013

Feb. 28, 2013

Hi, I'm sorry if my question would appear to be trivial, but what am I supposed to do, if I want to wrap class methods, that return a reference to another class? >From reading the list, I've gathered that apparently the best strategy of dealing with references is just to not to use them (convert to pointers immediately), because of some scoping rules issues. It works for me for a simple case of POD types, like cdef extern from "test.h": int& foo() cdef int* x = &foo() but in a more complex case, Cython generates incorrect C++ code (first it declares a reference, then assigns to it, which, of course, doesn't even compile): cdef extern from "token.h": cppclass Token: Token(const Datum&) except + cdef extern from "tokenstack.h": cppclass TokenStack: Token& top() except + cdef Token* tok = &self.pEngine.OStack.top() <-> Token *__pyx_v_tok; Token &__pyx_t_5; __pyx_t_5 = __pyx_v_self->pEngine->OStack.top(); __pyx_v_tok = (&__pyx_t_5); I would expect to see this instead: Token *__pyx_v_tok = &__pyx_v_self->pEngine->OStack.top(); Am I doing something wrong? Is there any other way to achieve what I want, other than writing custom C macros? Thanks, -- Sincerely yours, Yury V. Zaytsev

3 6

MemoryViews require writeable arrays?
by Dave Hirschfeld Feb. 28, 2013

Feb. 28, 2013

%%cython cimport cython import numpy as np cimport numpy as np @cython.boundscheck(False) @cython.wraparound(False) @cython.cdivision(True) cpdef double[:] return_one(double[:] x): return np.array([1.0]) In [43]: x = randn(3) ...: return_one(x) Out[43]: <MemoryView of 'ndarray' at 0x8ae14e0> In [44]: x.flags['WRITEABLE'] = False ...: return_one(x) Traceback (most recent call last): File "<ipython-input-44-4fbbd1035d56>", line 2, in <module> return_one(x) File "_cython_magic_7761e77f78c4e321261152684b47c674.pyx", line 11, in _cython_magic_7761e77f78c4e321261152684b47c674.return_one (C:\Users\dhirschfeld\.ipython\cython\_cython_magic_7761e77f78c4e321261152684b47 c674.c:1727) File "stringsource", line 619, in View.MemoryView.memoryview_cwrapper (C:\Users\dhirschfeld\.ipython\cython\_cython_magic_7761e77f78c4e321261152684b47 c674.c:8819) File "stringsource", line 327, in View.MemoryView.memoryview.__cinit__ (C:\Users\dhirschfeld\.ipython\cython\_cython_magic_7761e77f78c4e321261152684b47 c674.c:5594) ValueError: buffer source array is read-only Is this a required restriction? Is there any workaround? The context is calling cython routines using IPython.parallel. IIUC any input arrays sent over zmq are necessarily read-only. As can be seen with the example, even if we don't modify (or use) the input array at all we still get the error. Any help, esp. in regards to a workaround would be greatly appreciated! Thanks, Dave

2 4

MemoryView Casting slow compared to ndarray buffer syntax
by Dave Hirschfeld Feb. 28, 2013

Feb. 28, 2013

%%cython cimport cython import numpy as np cimport numpy as np ctypedef np.float64_t float64_t @cython.boundscheck(False) @cython.wraparound(False) @cython.cdivision(True) def echo_numpy(np.ndarray[float64_t, ndim=1] x): return x @cython.boundscheck(False) @cython.wraparound(False) @cython.cdivision(True) def echo_memview(double[:] x): return np.asarray(x) @cython.boundscheck(False) @cython.wraparound(False) @cython.cdivision(True) def echo_memview_nocast(double[:] x): return x In [19]: %timeit echo_memview(x) ...: %timeit echo_memview_nocast(x) ...: %timeit echo_numpy(x) 10000 loops, best of 3: 38.1 µs per loop 100000 loops, best of 3: 5.58 µs per loop 1000000 loops, best of 3: 749 ns per loop In [20]: 38.1e-6/749e-9 Out[20]: 50.86782376502002 In [21]: 5.58e-6/749e-9 Out[21]: 7.449933244325767 So it seems that the MemoryView is 50x slower than using the ndarray buffer syntax and even 7.5x slower without casting to an array. Is there anything that can be done about this or is it jsut something to be aware of and use each of them in the situations where they perform best? Thanks, Dave

1 0

Non-deterministic behavoiur?
by Dave Hirschfeld Feb. 28, 2013

Feb. 28, 2013

Using the following test code: import numpy as np from lapack import dgelsy from numpy.linalg import lstsq A = np.array([[ 0.12, -8.19, 7.69, -2.26, -4.71], [-6.91, 2.22, -5.12, -9.08, 9.96], [-3.33, -8.94, -6.72, -4.4 , -9.98], [ 3.97, 3.33, -2.74, -7.92, -3.2 ]]) # b = np.array([[ 7.3 , 0.47, -6.28], [ 1.33, 6.58, -3.42], [ 2.68, -1.71, 3.46], [-9.62, -0.79, 0.41]]) # print '#################' print '# ACTUAL RESULT #' print '#################' print lstsq(A, b)[0] print print '#################' print '# DGELSY RESULT #' print '#################' print dgelsy(A, b) I get: ################# # ACTUAL RESULT # ################# [[-0.6858 -0.2431 0.0642] [-0.795 -0.0836 0.2118] [ 0.3767 0.1208 -0.6541] [ 0.2885 -0.2415 0.4176] [ 0.2916 0.3525 -0.3015]] ################# # DGELSY RESULT # ################# [[-0.6858 -0.2431 0.0642] [-0.795 -0.0836 0.2118] [ 0.3767 0.1208 -0.6541] [ 0.2885 -0.2415 0.4176] [ 0.2916 0.3525 -0.3015]] All well and good, however if I type the `tmp` variable as a memview in the following code cdef double[:] res cdef double[:,:] tmp tmp = np.zeros([ldb, nrhs], order='F', dtype=np.float64) tmp[:b.shape[0]] = b res = np.ravel(tmp, order='F') the result changes?!? ################# # DGELSY RESULT # ################# [[-0.7137 -0.2429 0.0876] [-0.2882 -0.0884 -0.2117] [-0.4282 0.1284 0.0185] [ 0.9564 -0.2478 -0.1404] [ 0.3625 0.3519 -0.3607]] Remove the `cdef double[:,:] tmp` and I'm back to the correct result. Does this make any sense? To try and figure out what was going on I put in a couple of debugging print statements: print 'res = ', repr(np.asarray(res)) print 'res.flags = {{{flags}}}'.format(flags=np.asarray(res).flags) Only changing these lines resulted in the same incorrect result ################# # DGELSY RESULT # ################# res = array([ 7.3 , 1.33, 2.68, -9.62, 0., 0.47, 6.58, -1.71, -0.79, 0., -6.28, -3.42, 3.46, 0.41, 0.]) res.flags = { C_CONTIGUOUS : True F_CONTIGUOUS : True OWNDATA : False WRITEABLE : True ALIGNED : True UPDATEIFCOPY : False} [[-0.7137 -0.2429 0.0876] [-0.2882 -0.0884 -0.2117] [-0.4282 0.1284 0.0185] [ 0.9564 -0.2478 -0.1404] [ 0.3625 0.3519 -0.3607]] Removing (only) the print statements again gave me the correct results. So, it seems either typing the array as a memview or printing res will screw up the calculation. The cython code is given below. Any ideas if this is a cython bug or something I'm doing wrong? Thanks, Dave cdef extern from "mkl_lapack.h" nogil: void DGELSY(const MKL_INT* m, const MKL_INT* n, const MKL_INT* nrhs, double* a, const MKL_INT* lda, double* b, const MKL_INT* ldb, MKL_INT* jpvt, const double* rcond, MKL_INT* rank, double* work, const MKL_INT* lwork, MKL_INT* info) @cython.embedsignature(True) @cython.boundscheck(False) @cython.wraparound(False) @cython.cdivision(True) def dgelsy(double[:,:] A, double[:,:] b, double rcond=1e-15, overwrite_A=False, overwrite_b=False): cdef MKL_INT rank, info cdef MKL_INT m = A.shape[0] cdef MKL_INT n = A.shape[1] cdef MKL_INT nrhs = b.shape[1] cdef MKL_INT lda = m cdef MKL_INT ldb = max(m, n) cdef MKL_INT lwork = -1 cdef double worksize = 0 #cdef double[:,:] tmp cdef double[:] res, work cdef MKL_INT[:] jpvt if b.shape[0] != m: message = "b.shape[0] must be equal to A.shape[0].\n" message += "b.shape[0] = {b.shape[0]}\n" message += "A.shape[0] = {A.shape[0]}\n" raise MKL_LAPACK_ERROR(message.format(A=A, b=b)) flags = np.asarray(A).flags if not flags['F_CONTIGUOUS'] or not overwrite_A: A = A.copy_fortran() flags = np.asarray(b).flags if not flags['F_CONTIGUOUS'] or not overwrite_b or b.shape[0] < n: tmp = np.zeros([ldb, nrhs], order='F', dtype=np.float64) tmp[:b.shape[0]] = b res = np.ravel(tmp, order='F') else: res = np.ravel(b, order='F') #print 'res = ', repr(np.asarray(res)) #print 'res.flags = {{{flags}}}'.format(flags=np.asarray(res).flags) jpvt = np.empty(n, dtype=np.int32) DGELSY(&m, &n, &nrhs, &A[0,0], &lda, &res[0], &ldb, &jpvt[0], &rcond, &rank, &worksize, &lwork, &info) if info != 0: message = "Parameter {i} had an illegal value when calling gelsy." raise MKL_LAPACK_ERROR(message.format(i=info)) lwork = int(worksize) work = np.empty(lwork, dtype=np.float64) DGELSY(&m, &n, &nrhs, &A[0,0], &lda, &res[0], &ldb, &jpvt[0], &rcond, &rank, &worksize, &lwork, &info) if info != 0: message = "Parameter {i} had an illegal value when calling gelsy." raise MKL_LAPACK_ERROR(message.format(i=info)) return np.asarray(res).reshape(nrhs, -1).T[:n] #

1 2

Re: [Cython] [cython-users] freelist benchmarks
by Stefan Behnel Feb. 27, 2013

Feb. 27, 2013

mark florisson, 24.02.2013 15:50: > On 24 February 2013 13:52, Stefan Behnel wrote: >> for those who haven't notice my other e-mail, I implemented a new extension >> type decorator "@cython.freelist(N)" that replaces the normal object >> creation and deallocation with a freelist of N recently freed objects. >> Currently, it's only supported for types that do not have a base class (and >> lifting that restriction is not all that easy). > > Very cool! I've been wanting that for a while now :) So did I. > What's the hurdle with base classes? The problem is that the current way types are being instantiated is recursive. The top-most base class calls tp_alloc() and then each step in the hierarchy adds its bit of initialisation. If you want to introduce a freelist into this scheme, then it's still the top-most class that does the allocation, so it would need to manage all freelists of all of its children in order to return the right object struct for a given instantiation request. This cannot really be done at compile time. Imagine a subtype in a different module, for example, for which the code requests a freelist. The compilation of the base type wouldn't even know that it's supposed to manage a freelist at all, only the subtype knows it. There are a couple of ways to deal with this. One is to replicate the freelist in the base type for all subtypes that it finds at runtime. That might actually be the easiest way to do it, but it requires a bit of memory management in order to add a new freelist when a new subtype is found at runtime. It also means that we'd have to find the right freelist before we can get an object from it (or not, if it's empty), which would likely be an operation that's linear with the number of subtypes. And the freelist set would better be bounded in size to prevent users from flooding it with lots and lots of subtypes. Another option would be to split the initialisation up into two functions, one that allocates *and* initialises the instance and one that *only* initialises it. That would allow each hierarchy level to manage its own freelists and to take its own decision about where to get the object from. This approach comes with a couple of tricky details, as CPython doesn't provide support for this. So we'd need to find a way to handle type hierarchies that are implemented across modules. Maybe the best approach would be to let the base type manage everything and just statically limit the maximum number of subtypes for which it provides separate freelists, at a first come, first serve basis. And the freelist selection could be based on the object struct size (tp_basicsize) instead of the specific type. As long as we don't support inheriting from variable size objects (like tuple/bytes), that would cut down the problem quite nicely. I think I should just give it a try at some point. Stefan

4 10