Mailman 3 May 2001 - NumPy-Discussion

Numpy bug?
by Karl Bellve May 11, 2001

May 11, 2001

Since I didn't get a response, I will post again with a more appropiate subject heading. If you run the following code, you get an error. It appears that Numpy doesn't like to be restarted. If you take out the lines "import_array();" and "PyRun_SimpleString("from Numeric import *\n");" the code loops fine. for (int x = 0; x < 100; x++) { Py_Initialize(); import_array(); PyRun_SimpleString("from Numeric import *\n"); TRACE("%d\n",x); Py_Finalize(); } Here is the output: 0 1 Fatal Python error: UNREF invalid object -- Cheers, Karl Bellve

2 2

Re: [Numpy-discussion] Thoughts on masked arrays
by Michael Haggerty May 11, 2001

May 11, 2001

"Paul F. Dubois" <paul(a)pfdubois.com> writes: > -----Original Message----- > Michael Haggerty wrote > 1. I couldn't find a simple way to tell if all of the cells of a > masked array are unmasked. > ====== > So your test could be if count(x) < product(x.shape): error... > > So your test could be > if make_mask(m.mask(),flag=1) is not None: > error... > > You could also consider if not Numeric.allclose(m.filled(0), m.filled(1)) > or > m.mask() is not None and not Numeric.alltrue(Numeric.ravel(m.mask())): Shouldn't that be m.mask() is not None and Numeric.sometrue(Numeric.ravel(m.mask())) ? ("Proof" that these expressions are nonintuitive.) > Is that enough ways to do it? (TM) (:-> Frankly, it's too many ways to do it, none of them obvious to the writer or the reader. This is a simple and useful concept and it should have one obvious implementation. > I'm not against is_unmasked but I'm not sure how much it would get > used and I don't like the name. I hate query methods with side > effects (if you use them in an assert you change the program). In this case the side effect is to change the internal representation of the object without changing its semantics, so I don't find it too objectionable. But omit this optimization if you prefer; the query method would be just as useful even without the side effect. Because of the relationship with filled(), maybe this query function should be called m.isfull(). There should probably also be an isfull(m) function for the same reason that there is a mask(m) function. > A method that replaces the mask with None if possible might make > sense. m.unmask()? m.demask()? m.debride() ? Of these names, I like m.unmask() the best. I assume that it would set m.__mask=None if possible and throw an exception if not. On the other hand, while it would be desirable to have a function equivalent (i.e., unmask(m)), this would be awkward because a function should usually not change its argument. Therefore, I suggest adding a safe analogue of raw_data() that throws an exception if the array has a nontrivial mask and otherwise returns self.__data. E.g. [untested]: class MaskedArray: [...] def data(self): """If no values are masked, return self.__data(). Otherwise raise an exception. """ d = self.__data m = self.__mask if m is not None and Numeric.sometrue(Numeric.ravel(m)): raise MAError, "MaskedArray cannot be converted to array" elif d.iscontiguous(): return d else: return Numeric.array(d, typecode=d.typecode(), copy=1, savespace = d.spacesaver()) def data(a): if isinstance(a, MaskedArray): return m.data() elif isinstance(a, Numeric.ArrayType) and a.iscontiguous(): return a else: return Numeric.array(a) A more obscure name should be chosen since you seem to encourage "from MA import *". > 3. I found the semantics of MA.compress(condition,a,axis=0) to be > inconvenient and inconsistent with those of Numeric.compress. > ====== > It has been an interesting project in that there are hundreds of these > individual little design questions. > Can you propose the semantics you would like in a precise way? Include the > case where the condition has masked values. > ====== In the simple case where the condition has no masked values, I think compress() should simply pick slices out according to condition, without regard to which cells of x are masked. When condition is masked, I don't think that there is a sensible interpretation for compress() because a "masked" value in condition means you don't know whether that slice of x should be included or not. Since you can't have an output array of indeterminate shape, I would throw an exception if condition is masked. Here is my attempt [untested]: def compress(condition, x, dimension=-1): # data function is defined above (throws exception if condition is masked): c = data(condition) if mask(x) is None: mask = None else: mask=Numeric.compress(condition, mask(x), dimension) return array(Numeric.compress(condition, filled(x), dimension), mask=mask) Yours, Michael -- Michael Haggerty mhagger(a)alum.mit.edu

2 1

20.1.0a3 -- for those who like pickles
by Paul F. Dubois May 10, 2001

May 10, 2001

I put a source tarball up for 20.1.0a3, also in CVS as of now. This contains: DL_EXPORT changes Konrad's fix to import_array A new MA that pickles A new package "Properties", that replaces MA.activeattr. Properties allows you to have properties such as MA/Numeric's shape attribute, that look like attributes but actually invoke functions on get, set, and/or delete. It can also be used to make attributes unwriteable or undeleteable. It doesn't have anything to do with Numeric but I use it in MA and kinds and don't have much else of an idea where to put it.

1 0

Modified header files
by Konrad Hinsen May 10, 2001

May 10, 2001

Recently we had a discussion about how to use NumPy arrays from extension modules with multiple source files, on various platforms. The modified header files that are attached to this message provide a (hopefully!) general solution. In fact, I propose to make them part of the official distribution, unless there are objections. If used like before, these header files give exactly the same result as the ones in NumPy 20.0.0. However, they permit to define the name of the C API pointer array and make it globally visible. Under the condition that the chosen name is unique, this should not create problems under any platform, no matter if static or dynamic linking is used. To use NumPy features in multiple source file extension modules, you have to write #define PY_ARRAY_UNIQUE_SYMBOL PyArrayXXX #include "Numeric/arrayobject.h" in the main source file (the one that contains the init function) and #define PY_ARRAY_UNIQUE_SYMBOL PyArrayXXX #define NO_IMPORT_ARRAY #include "Numeric/arrayobject.h" in all the others. The symbol you choose instead of PyArrayXXX should contain both the name of the imported module (array) and the name of the importing module (whatever your module is called) in order to be unique with a reasonably high probability. The same applies to the Ufunc module, just replace "array" by "ufunc" in the example. I have also applied the "static" correction to the Ufunc header file, there is no reason not to do it. Konrad. -- ------------------------------------------------------------------------------- Konrad Hinsen | E-Mail: hinsen(a)cnrs-orleans.fr Centre de Biophysique Moleculaire (CNRS) | Tel.: +33-2.38.25.56.24 Rue Charles Sadron | Fax: +33-2.38.63.15.17 45071 Orleans Cedex 2 | Deutsch/Esperanto/English/ France | Nederlands/Francais -------------------------------------------------------------------------------

2 1

UNREF
by Karl Bellve May 9, 2001

May 9, 2001

Occasionally, if I reinitialize Python and rerun a script, I get the following error: FATAL Python error: UNREF invalid object from my code I do the following: /* alot of other crap and then */ PyRun_SimpleString("From Numeric import *\n"); eventually in multiarraymodule.c array_set_string_function() Py_INCREF(Py_None) which eventually calls (according to my stack trace) where the error occurs PyArray_SetStringFunction() Py_XDECREF(PyArray_ReprFunction) Is this a Numpy problem, Python problem, or my problem? I can run the same python script several times, and then I get the error. Something isn't being cleaned up properly around Py_Initialize() Py_Finalize() pair. I am using Python 2.1b1 and Numeric 18_1 with Boost -- Cheers, Karl Bellve

1 1

Re: [Numpy-discussion] Can this be done efficiently with array functions?
by Chris Barker May 4, 2001

May 4, 2001

Janko and Paul, Thanks for your suggestions. I had come up with something like Janko's solution, but not as cleanly. It seems that Paul's suggestion is essentially the same thing. Maybe I'll do some tests to see which is faster. By the way, I tried Janko's code, and got an error on "a == 0" (I'm running Python 2.0). I looked at it again, and then looked at the release notes for Python 2.1, and saw that rich comparisons are working..Glory be! I'll be upgrading real soon, just for that! thanks again, -Chris PS: anyone have a response for my Int32 question? Why isn't it a "long" on Linux or windows on Intel? Janko Hauser wrote: > truely_max=100 > > a=array([[10, 0, 5, 0,], > [ 0, 0, 0, 0,], > [ 0, 5,15, 0,], > [ 0, 0, 0, 0,], > [ 0, 3, 1,25,], > [ 0, 0, 0, 0,], > [ 4, 7, 2,12,], > [ 0, 0, 0, 0,]]) > > am=minimum.reduce(where(a==0,truely_max,a),1) > am=where(am==truely_max,0,am) > print am > [5 0 5 0 1 0 2 0] > "Paul F. Dubois" wrote: > but there doesn't seem to be > a minimum.reduce() in MA. > > --> I'll add this to my list. That will be great if you get a chance. > Meanwhile you could do something like this, > where x is your array: > y = MA.masked_equal(x,0,copy=0) > Numeric.minimum.reduce(y.filled(100000), 1) > If x is floating point use masked_value rather than masked_equal. > --- > --> There are two problems. There is a hit, even in the case of no mask. > It is hard to quantify because it depends on what operations you do. > Much of that could be eliminated if MA was in C but it is Python. > It would have driven me mad to write it in C, however. WOW, I'm impressed, I had assumed that some C was involved. Nice job. > More seriously, there are a few semantic issues where it isn't clear > what you mean if a mask is present. I'm sure these could be resolved. It would be nice if we could have in mind to move in that direction. I'm not going ot be able to do any of the coding, so that's the last you'll here from me about it. -Chris -- Christopher Barker, Ph.D. ChrisHBarker(a)home.net --- --- --- http://members.home.net/barkerlohmann ---@@ -----@@ -----@@ ------@@@ ------@@@ ------@@@ Oil Spill Modeling ------ @ ------ @ ------ @ Water Resources Engineering ------- --------- -------- Coastal and Fluvial Hydrodynamics -------------------------------------- ------------------------------------------------------------------------

1 0

[Q]Best way for an array operation?
by Daehyok Shin May 4, 2001

May 4, 2001

What is the best Numpy way for the following work? for i in range(len(x)): if x[i] > 0: y[i] = v[i] z[i] = u[i]+2 Daehyok Shin (Peter)

5 5

curious FFFT timing
by Eric Hagemann May 4, 2001

May 4, 2001

I have an application that needs to perform close 1/2 million FFT's (actually inverse FFTs) I was stymied by the time it was taking so I created a simple benchmark program to measure raw fft speed - hoping to extrapolate and see if the FFT algorithm was the time consuming portion of the algorithm. The result I obtained surprised me Inverse FFT ( 4096) Number of seconds 5.608000 Time per FFT 0.005608 Forward FFT ( 4096) Number of seconds 3.164000 Time per FFT 0.003164 Inverse FFT ( 4095) Number of seconds 8.222000 Time per FFT 0.008222 Forward FFT ( 4095) Number of seconds 5.468000 Time per FFT 0.005468 Inverse FFT ( 8192) Number of seconds 12.578000 Time per FFT 0.012578 Forward FFT ( 8192) Number of seconds 7.551000 Time per FFT 0.007551 Inverse FFTs were taking almost twice as long as forward FFTs ! Bottom line I found that in the inverse code the multiplication of 1/N at the end was the culprit. When I removed the /n from the line return _raw_fft(a, n, axis, fftpack.cffti, fftpack.cfftb, _fft_cache)/n the code timing was more along the line of what was expected. Any one want to speculate on the timing difference for the 4095 vice 4096 long transforms ? Since 4095 is not a power of two I would have expected a greater time difference (DFT vice FFT) N*N / N*log2(N) = N/log2(N) = 4096/12 = 341.33 which is >> than the 1.72 seen above Cheers Eric ------------------------------------ Code Snippet --------------------------------------------- import time from Numeric import * from FFT import * def f_fft(number,length): a = arange(float(length)) start_time = time.time() for i in xrange(number): b = fft(a) durat = time.time() - start_time durat_per = durat / number print "Forward FFT (%10d) Number of seconds %12.6f Time per FFT %12.6f" % (length,durat,durat_per) def i_fft(number,length): a = arange(float(length)) start_time = time.time() for i in xrange(number): b = inverse_fft(a) durat = time.time() - start_time durat_per = durat / number print "Inverse FFT (%10d) Number of seconds %12.6f Time per FFT %12.6f" % (length,durat,durat_per) i_fft(1000,4096) f_fft(1000,4096) i_fft(1000,4095) f_fft(1000,4095) i_fft(1000,8192) f_fft(1000,8192)

2 2

For Python 6.1b2 you need a new Pyfort
by Paul F. Dubois May 3, 2001

May 3, 2001

There is a new Pyfort release 6.3 that lets Pyfort work with Python 6.1b2. Please get it from http://sourceforge.net/projects/pyfortran if you are using 6.1b2. This Pyfort will work just fine with any Python version it used to work with.

2 1

static pointer
by Paul F. Dubois May 3, 2001

May 3, 2001

Mr. Verveer wrote: "I recently installed this version and it unfortunately broke my own extension C code. I found that my code did not link properly anymore because of the addition of the static declaration tot the API pointer (as noted in changes.txt). Since my extensions consist of several source files, this led to problems when importing my extensions. If I understand correctly, the addition of the static keyword means that the import_array() function can now only be used by extensions that are defined in a single C source file. I can easily work around this problem, by defining my own api pointer, and using a modified import_array() in my code, but I thought I should draw your attention to this problem in any case. Maybe another solution can be found that does not need this pointer to be declared static?" -- I tried this change to see whether it had any bad consequences. It did not appear to have bad consequences in that Numpy itself and all its packages built ok, and we were formerly not following the instructions in the Python documentation and now we are. And, OS-X now worked. So far so good. If I understand you correctly, this does not work properly if the extension is built in pieces. I think I see why. (Frankly, linking and dynamic loading are not something I know much about; as a user I just flounder around.) This all started because of the need to have access to the C API of one dynamic extension from another. Is there a way of institutionalizing your work-around so that others can use it? Frankly, I'm out of ideas. The community is going to have to rescue me on this one. Paul P.S. Wait. I have an idea. Konrad was the one who put all this in. He should be punished for his good deed...(:->

3 3