Mailman 3 January 2006 - NumPy-Discussion

PyArray_ContiguousFromObject & PyDECREF
by Michiel Jan Laurens de Hoon 16 May '11

16 May '11

Dear pythoneers, In Chapter 12 ("Writing a C extension") in the Numerical Python manual, in the first example PyArray_ContiguousFromObject is used together with PyDECREF, while in the last example PyArray_ContiguousFromObject is used but the arrays are not PyDECREF'ed before the function returns. From Chapter 13 ("C API Reference") it seems that a PyDECREF is needed in both cases. Or am I missing something here? Is PyDECREF not needed for some reason in the last example? Thanks in advance, Michiel de Hoon U Tokyo First example: ============== static PyObject * trace(PyObject *self, PyObject *args) { PyObject *input; PyArrayObject *array; double sum; int i, n; if (!PyArg_ParseTuple(args, "O", &input)) return NULL; array = (PyArrayObject *) PyArray_ContiguousFromObject(input, PyArray_DOUBLE, 2, 2); if (array == NULL) return NULL; n = array->dimensions[0]; if (n > array->dimensions[1]) n = array->dimensions[1]; sum = 0.; for (i = 0; i < n; i++) sum += *(double *)(array->data + i*array->strides[0] + i*array->strides[1]); Py_DECREF(array); return PyFloat_FromDouble(sum); } Last example: ============= static PyObject * matrix_vector(PyObject *self, PyObject *args) { PyObject *input1, *input2; PyArrayObject *matrix, *vector, *result; int dimensions[1]; double factor[1]; double real_zero[1] = {0.}; long int_one[1] = {1}; long dim0[1], dim1[1]; extern dgemv_(char *trans, long *m, long *n, double *alpha, double *a, long *lda, double *x, long *incx, double *beta, double *Y, long *incy); if (!PyArg_ParseTuple(args, "dOO", factor, &input1, &input2)) return NULL; matrix = (PyArrayObject *) PyArray_ContiguousFromObject(input1, PyArray_DOUBLE, 2, 2); if (matrix == NULL) return NULL; vector = (PyArrayObject *) PyArray_ContiguousFromObject(input2, PyArray_DOUBLE, 1, 1); if (vector == NULL) return NULL; if (matrix->dimensions[1] != vector->dimensions[0]) { PyErr_SetString(PyExc_ValueError, "array dimensions are not compatible"); return NULL; } dimensions[0] = matrix->dimensions[0]; result = (PyArrayObject *)PyArray_FromDims(1, dimensions, PyArray_DOUBLE); if (result == NULL) return NULL; dim0[0] = (long)matrix->dimensions[0]; dim1[0] = (long)matrix->dimensions[1]; dgemv_("T", dim1, dim0, factor, (double *)matrix->data, dim1, (double *)vector->data, int_one, real_zero, (double *)result->data, int_one); return PyArray_Return(result); } -- Michiel de Hoon, Assistant Professor University of Tokyo, Institute of Medical Science Human Genome Center 4-6-1 Shirokane-dai, Minato-ku Tokyo 108-8639 Japan http://bonsai.ims.u-tokyo.ac.jp/~mdehoon

3 2

Numpy x Matlab: some synthetic benchmarks
by Paulo Jose da Silva e Silva 11 Feb '07

11 Feb '07

Hello, Travis asked me to benchmark numpy versus matlab in some basic linear algebra operations. Here are the resuts for matrices/vectors of dimensions 5, 50 and 500: Operation x'*y x*y' A*x A*B A'*x Half 2in2 Dimension 5 Array 0.94 0.7 0.22 0.28 1.12 0.98 1.1 Matrix 7.06 1.57 0.66 0.79 1.6 3.11 4.56 Matlab 1.88 0.44 0.41 0.35 0.37 1.2 0.98 Dimension 50 Array 9.74 3.09 0.56 18.12 13.93 4.2 4.33 Matrix 81.99 3.81 1.04 19.13 14.58 6.3 7.88 Matlab 16.98 1.94 1.07 17.86 0.73 1.57 1.77 Dimension 500 Array 1.2 8.97 2.03 166.59 20.34 3.99 4.31 Matrix 17.95 9.09 2.07 166.62 20.67 4.11 4.45 Matlab 2.09 6.07 2.17 169.45 2.1 2.56 3.06 Obs: The operation Half is actually A*x using only the lower half of the matrix and vector. The operation 2in2 is A*x using only the even indexes. Of course there are many repetitions of the same operation: 100000 for dim 5 and 50 and 1000 for dim 500. The inner product is number of repetitions is multiplied by dimension (it is very fast). The software is numpy svn version 1926 Matlab 6.5.0.180913a Release 13 (Jun 18 2002) Both softwares are using the *same* BLAS and LAPACK (ATLAS for sse). As you can see, numpy array looks very competitive. The matrix class in numpy has too much overhead for small dimension though. This overhead is very small for medium size arrays. Looking at the results above (specially the small dimensions ones, for higher dimensions the main computations are being performed by the same BLAS) I believe we can say: 1) Numpy array is faster on usual operations but outerproduct (I believe the reason is that the dot function uses the regular matrix multiplication to compute outer-products, instead of using a special function. This can "easily" changes). In particular numpy was faster in matrix times vector operations, which is the most usual in numerical linear algebra. 2) Any operation that involves transpose suffers a very big penalty in numpy. Compare A'*x and A*x, it is 10 times slower. In contrast Matlab deals with transpose quite well. Travis is already aware of this and it can be probably solved. 3) When using subarrays, numpy is a slower. The difference seems acceptable. Travis, can this be improved? Best, Paulo Obs: Latter on (in a couple of days) I may present less synthetic benchmarks (a QR factorization and a Modified Cholesky).

13 27

RE: [Numpy-user] possible error with isarrtype
by O'Keefe, Michael 07 Feb '06

07 Feb '06

> Now, the data-type object itself is the preferred way. I need to figure > out what to do with these functions now and see if they even have any > use anymore. Thanks Travis. With regard to what to do with the functions, I only ended up trying them because they were there. From my point of view, I'd love to learn how to code the "right" way. To that end, any means of informing the user of "deprecated" functions (or just taking them out) would be fine with me. However, I realize some may have legacy code and so I would defer to their opinions. > What use were you making of them? I was just migrating some working scipy/numeric code I'd been using to the newer numpy (0.9.4), python (2.4.2), and scipy (0.4.4). When I ran my unit-tests, one line that gave me problems was: import scipy as sp ... assert( type(deltaTime_sec)==float ) This line of code was there to inform me if I was accidentally passing in an array type instead of a scalar. This no longer worked in the new environment. The type of deltaTime_sec came out as float64_arrtype which did not equal type float causing the assertion to trip which I believe to be caused by the following: >>> a=sp.array([1.],float) >>> type(a) <type 'numpy.ndarray'> >>> type(a[0]) # somewhere I assign deltaTime_sec = someArray[0] <type 'float64_arrtype'> I ended up fixing it with the following code which seems to do the trick (though I'm always keen to learn better ways of doing things): assert( sp.size(deltaTime_sec)==1 ) Thanks! Michael

6 8

numpy can't use ACML
by George Nurser 01 Feb '06

01 Feb '06

I'm sure this is well known, but I just realized that numpy cannot use the _dotblas.so that it makes when it is compiled with ACML. This is because ACML only has the fortran blas libraries, not cblas. numpy will find the acml libraries and use them to make a _dotblas.so without complaining, if you have [blas] blas_libs = acml language = f77 in your site.cfg. But attempting to import this _dotblas.so gives errors .... so numpy never actually uses it. >>> import _dotblas.so Traceback (most recent call last): File "<stdin>", line 1, in ? ImportError: ./_dotblas.so: undefined symbol: cblas_zaxpy nm _dotblas.so gives a whole stream of undefined cblas_xxxx symbols. So what are the options? Forget about ACML? Find an optimized cblas for numpy _dotblas but use the ACML flapack for scipy? Persuade somebody to write a scipy-style ?f2py interface to generate _dotblas.so using the fblas? George Nurser.

2 2

Numerical Mathematics Consortium vs. scipy and numpy
by Zdeněk Hurák 01 Feb '06

01 Feb '06

Hello all, the other day I ran across the following link ----> http://www.nmconsortium.org/. Are you, Python numerics guys, aware of this? I admit I have absolutely no idea about the actual usefullness and relevancy of this consortium and their drafts. But it might be useful just to be aware... Important players in scientific computation are involved. Best regards, Zdenek Hurak

3 5

Английский на уровне носителя языка!
by Выучить английский - просто 31 Jan '06

31 Jan '06

Не учится или не совершенствуется английский? Курсы английского не дают желаемого результата? Нет способностей к языкам? Используйте возможность всего за 3 месяца научиться грамотному английскому языку! В России существует множество различных языковых школ, курсов, методик, но: к сожалению, язык либо не учится, либо не совершенствуется до должного уровня. Даже после изучения языка в престижных отечественных ВУЗах, люди совершают грубейшие ошибки, зачастую ведущие к абсолютному непониманию их речи носителями языка. Предлагаем Вам дойти до такого уровня владения языком, который бы позволил чувствовать себя, как рыба в воде, - свободно говорить на грамотном английском и избавиться от страха общения. Кроме того, Вы легко сможете довести свой словарный запас до уровня носителя языка. Уникальность методики и ее жизнеспособность признана носителями языка. Методика включает в себя постановку абсолютно правильного произношения, позволяет обучить самого необучаемого взрослого и уже через 3 месяца свободно и грамотно говорить и писать. В процессе обучения Вы научитесь: воспринимать речь на слух освоите устойчивые навыки чтения, в том числе вслух и понимания текста, а также качественной имитации носителей языка восстановите навыки эффективного самообучения улучшите свою память освоите эффективные способы расширения словарного запаса и совершенствования навыков устной и письменной речи освоите соответствие знаков транскрипции реальным звукам увеличите скорость чтения увеличите скорость пропускания символьной и звуковой информации через мозг ускорите скорость мышления Автор методики - наш соотечественник, живущий в США, широко известный в СНГ и на Западе журналист - педагогический обозреватель ряда международных средств массовой информации, лингвист-полиглот, cоветник Дня Земли и ряда других организаций ООН, российского и американского правительств по решению проблем гарантированного переобучения и повышения квалификации взрослых. Информация об интеллектуальных достижениях автора этого курса и его системе обучения занесена во всемирное издание Книги рекордов Гиннесса, международную энциклопедию Who's Who of Professionals, российскую Книгу феноменальных достижений человека Диво и др. подобные издания. С 9 по 23 декабря 2005 г. автор провел очные блиц-курсы английского языка и интеллектуальной машинописи для высшего руководства (Совета директоров) крупнейшего предприятия Европы - ОАО Камаз. Заочные курсы - Вы можете приобрести заочный курс изучения английского языка на CD и DVD. Вам нужно будет заниматься всего по 1-2 часа каждый день в течение 2-3 месяцев. Очно-заочные курсы - 16 или 40 часов обучения английскому непосредственно автором методики и 4 часа интеллектуальной машинописи + каждый студент получает заочный курс для самостоятельных занятий. Харизма автора и большой опыт обучения языку сделают обучение невероятно легким и интересным. Корпоративное обучение - Успех развития любой компании в первую очередь зависит от квалификации персонала. В наше время международные деловые отношения выходят на новый уровень. Общение с зарубежными партнерами через переводчика говорит о низком уровне специалистов. Более того, во всем мире отношение к партнерам, говорящим грамотно на иностранном языке, на порядок выше. Если Вы заинтересованы в обучении Вашего персонала, а, особенно, топ-менеджмента, обращайтесь к нам Для желающих, возможна организация выездных семинаров (обучение + отдых за рубежом). Позвоните нам: Все виды курсов: Тел. +7 (926) 197-99-05 Заочные курсы: Тел. +7(495) 291-72-64 (после 18.00), +7 (916) 855-17-34 Напишите: englishbest(a)netscape.net

1 0

gfortran, numpy build
by Matthew Brett 31 Jan '06

31 Jan '06

Hi, I have the latest version of Mandrake, which has no g77, only gfortran - and needed the following patch to compile (with some added g77, gfortran outputs): Index: numpy/distutils/fcompiler/gnu.py =================================================================== --- numpy/distutils/fcompiler/gnu.py (revision 2033) +++ numpy/distutils/fcompiler/gnu.py (working copy) @@ -20,8 +20,13 @@ # GNU Fortran (GCC) 3.3.3 (Debian 20040401) # GNU Fortran 0.5.25 20010319 (prerelease) # Redhat: GNU Fortran (GCC 3.2.2 20030222 (Red Hat Linux 3.2.2-5)) 3.2.2 20030222 (Red Hat Linux 3.2.2-5) + # Mandrake: GNU Fortran (GCC) 3.3.6 (Mandrakelinux release 3.3.6-2mdk) + # + # 'gfortran --version' results + # Redhat: GNU Fortran 95 (GCC 4.0.2 20051125 (Red Hat 4.0.2-8)) + # Mandrake: GNU Fortran 95 (GCC 4.0.1 (4.0.1-5mdk for Mandriva Linux release 2006.0)) - for fc_exe in map(find_executable,['g77','f77']): + for fc_exe in map(find_executable,['g77','f77','gfortran']): if os.path.isfile(fc_exe): break executables = { @@ -41,7 +46,10 @@ if os.name != 'nt' and sys.platform!='cygwin': pic_flags = ['-fPIC'] - g2c = 'g2c' + if os.path.basename(fc_exe) == 'gfortran': + g2c = 'gfortran' + else: + g2c = 'g2c' #def get_linker_so(self): # # win32 linking should be handled by standard linker Can I suggest this a patch to svn numpy? Thanks, Matthew

3 3

bdist_egg trips up over non-extension extensions
by Robert Kern 31 Jan '06

31 Jan '06

In particular, numpy can't build into an egg currently because numpy/__config__.py is created using Configuration.add_extension() and a generator function. In order to make extension modules usable from zipped eggs, setuptools generates a pure Python stub loader that will extract the extension module to a .hidden directory and load it from there. In this case, since the Extension isn't actually an extension, the stub loader is broken, because it expects a numpy/__config__.so. Can we find another way to generate Python code instead of abusing Extension this way? -- Robert Kern robert.kern(a)gmail.com "In the fields of hell where the grass grows high Are the graves of dreams allowed to die." -- Richard Harter

4 6

Re: [Numpy-user] possible error with isarrtype
by Travis Oliphant 31 Jan '06

31 Jan '06

O'Keefe, Michael wrote: >I was just updating my code to work with python 2.4.2 and numpy 0.9.4 and observed the following behavior with the function isarrtype. Just thought I would report it (see below). > > > Thanks for this report. I've fixed the problem, but you should probably not be using these functions anymore. These were put in place before the data-type objects materialized. In the past, what is now the typeobject of the array scalars was thought to be the basic way to determine the data-type of the array. Now, the data-type object itself is the preferred way. I need to figure out what to do with these functions now and see if they even have any use anymore. What use were you making of them? Best, -Travis

1 0

A couple of errors in numpy/records.py
by Francesc Altet 31 Jan '06

31 Jan '06

Hi, I've detected a couple of issues in records.py. For the first one, the next should fix it: Index: numpy/core/records.py =================================================================== --- numpy/core/records.py (revision 2033) +++ numpy/core/records.py (working copy) @@ -410,7 +410,8 @@ offset=offset) elif isinstance(obj, sb.ndarray): res = obj.view(recarray) - if issubclass(res.dtype, nt.void): + if issubclass(res.dtype.type, nt.void): res.dtype = sb.dtype((record, res.dtype)) + return res else: raise ValueError("Unknown input type") Now, another curious thing (sorry, no patch this time): In [6]: for i in numpy.rec.array([(1,2)], formats='u4,u4'): print i ...: (1L, 2L) [No problem with this] In [7]: for i in numpy.rec.array([(1,2)], formats='u4,u4')[0]: print i ...: --------------------------------------------------------------------------- exceptions.KeyError Traceback (most recent call last) /home/faltet/computacio.nobackup/hdf5-1.7.51/test/<console> /usr/lib/python2.3/site-packages/numpy/core/records.py in __getitem__(self, obj) 126 127 def __getitem__(self, obj): --> 128 return self.getfield(*(self.dtype.fields[obj][:2])) 129 130 def __setitem__(self, obj, val): KeyError: 0 I'd expect something like 1L 2L Cheers, -- >0,0< Francesc Altet http://www.carabos.com/ V V Cárabos Coop. V. Enjoy Data "-"

2 1