Mailman 3 February 2000 - NumPy-Discussion

An Experiment in code-cleanup.
by Travis Oliphant Feb. 10, 2000

Feb. 10, 2000

I wanted to let users of the community know (so they can help if they want, or offer criticism or comments) that over the next several months I will be experimenting with a branch of the main Numerical source tree and endeavoring to "clean-up" the code for Numerical Python. I have in mind a few (in my opinion minor) alterations to the current code-base which necessitates a branch. Guido has made some good suggestions for improving the code base, and both David Ascher and Paul Dubois have expressed concerns over the current state of the source code and given suggestions as to how to improve it. That said, I should emphasize that my work is not authorized, or endorsed, by any of the people mentioned above. It is simply my little experiment. My intent is not to re-create Numerical Python --- I like most of the current functionality --- but to merely, clean-up the code, comment it, and change the underlying structure just a bit and add some features I want. One goal I have is to create something that can go into Python 1.7 at some future point, so this incarnation of Numerical Python may not be completely C-source compatible with current Numerical Python (but it will be close). This means C-extensions that access the underlying structure of the current arrayobject may need some alterations to use this experimental branch if it every becomes useful. I don't know how long this will take me. I'm not promising anything. The purpose of this announcement is just to invite interested parties into the discussion. These are the (somewhat negotiable) directions I will be pursuing. 1) Still written in C but heavily (in my opinion) commented. 2) Addition of bit-types and unsigned integer types. 3) Facility for memory-mapped dataspace in arrays. 4) Slices become copies with the addition of methods for current strict referencing behavior. 5) Handling of sliceobjects which consist of sequences of indices (so that setting and getting elements of arrays using their index is possible). 6) Rank-0 arrays will not be autoconverted to Python scalars, but will still behave as Python scalars whenever Python allows general scalar-like objects in it's operations. Methods will allow the user-controlled conversion to the Python scalars. 7) Addition of attributes so that different users can configure aspects of the math behavior, to their hearts content. If their is anyone interested in helping in this "unofficial branch work" let me know and we'll see about setting up someplace to work. Be warned, however, that I like actual code or code-templates more than just great ideas (truly great ideas are never turned away however ;-) ) If something I do benefits the current NumPy source in a non-invasive, backwards compatible way, I will try to place it in the current CVS tree, but that won't be a priority, as my time does have limitations, and I'm scratching my own itch at this point. Best regards, Travis Oliphant

10 18

digest
by Les Schaffer Feb. 9, 2000

Feb. 9, 2000

just switched over from matrix-sig to numpy-discussion. in the process i changed to the digest version and got my first issue. is it possible to distribute the digests properly formatted as multipart/digests as per rfc822 and company? having such a formatted digest makes it very easy when using an emailer like VM in emacs: VM automatically displays the digest as a virtual folder, allowing one to browse all the posts in a given digest very quickly and easily. don't know whether the other (lacklusters) emailers out there will handle it so nicely, but i don't think the extra required markers will interfere with your reading of the digests at all. highly recommended. i'd be glad to work with whoever has control over this to ensure that the proper markers get placed into the digests. les schaffer

2 1

RE: [Matrix-SIG] An Experiment in code-cleanup.
by Beausoleil, Raymond Feb. 9, 2000

Feb. 9, 2000

From: Konrad Hinsen [mailto:hinsen@cnrs-orleans.fr] > > (1) The extent of the support for LAPACK. Do we want to stick > > with LAPACK Lite? > > There has been a full LAPACK interface for a long while, of which > LAPACK Lite is just the subset that is needed for supporting the > high-level routines in the module LinearAlgebra. I seem to have lost > the URL to the full version, but it's on my disk, so I can put it > onto my FTP server if there is a need. Yes, I'd like to get a copy! You can simply e-mail it to me, if you'd prefer. > > (2) The storage format. If we've still got row-ordered matrices > > under the hood, and we want to use native LAPACK libraries that > > were compiled using column-major format, then we'll have to be > > careful to set all of the flags correctly. This isn't going to > > be a big deal, _unless_ NumPy will support more of LAPACK when a > > native library is available. Then, of course, there ... > > The low-level interface routines don't take care of this. It's the > high-level Python code (module LinearAlgebra) that sets the > transposition argument correctly. That looks like a good compromise > to me. I'll have to look at this more carefully. Due to my relative lack of Python experience, I hacked the C code so that Fortran routines could be called instead, producing the expected results. > > (3) Through the judicious use of header files with compiler- > > dependent flags, we could accommodate the various naming > > conventions used when the FORTRAN libraries were compiled (e.g., > > sgetrf_ or SGETRF). > > That's already done! Where? Even in the latest f2c'd source code that I downloaded from SourceForge, I see all names written using the lower-case-trailing-underscore convention (e.g., dgeqrf_). The Intel MKL was compiled from Fortran source using the upper-case-no-underscore convention (e.g., DGEQRF). If I replace dgeqrf_ with DGEQRF in dlapack_lite.c (and a few other tweaks), then the subsequent link with the IMKL succeeds. ============================ Ray Beausoleil Hewlett-Packard Laboratories mailto:beausol@hpl.hp.com Vox: 425-883-6648 Fax: 425-883-2535 HP Telnet: 957-4951 ============================

2 1

Upcasting
by Pitts, Todd A., Ph.D. Feb. 9, 2000

Feb. 9, 2000

Here are my two cents worth on the subject... Most of what has been said in this thread (at least what I have read) I find very valuable. Apparently, many people have been thinking about the subject. I view this problem as inherent in a language without an Lvalue (like C) that allows a very explicit and clear definition from the programmer's point of view as to the size of the container you are going put things in. The language in many cases simply returns an object to you and has made some decision as to what you "needed" or "wanted". Of course, this is one of the things that makes languages such an Numerical Python, Matlab, IDL, etc. very nice for protyping and investigating. In many cases this decision will be adequate or acceptable. In some, quite simply, it will not be. At this point the programmer has to have a good means of managing this decision himself. If memory is not a constraint, I can think of very few situations (none, actually) where I would choose to go with something other than the Numerical Python default of double. In general, that is what you get when creating python arrays unless you make some effort to obtain some other type. However, in some important (read "cases that affect the author") situations memory is a very critical constraint. Typically, in Numerical Python I use 4-byte floats. In fact, one of the reasons I use Numerical Python is because I don't *need* doubles and matlab for example is really only setup to work gracefully with doubles. I do *need* to conserve memory as I deal with very large data sets. It seems the question we are discussing is not really "what *should* be done in terms of casting?" but "what provides good enough decisions much of the time *and* a gracefull way to manage the decisions when "good enough" no longer applies to you?" Currently, this is not a trivial thing to manage. Reading in a 100 MB data set and multiplying by the python scalar 2 produces a 200 MB data set. I manage this by wrapping the 2 in an array. This happens, of course, all the time. Having to do this once is not a big deal, but doing everywhere in code that uses floats makes for cluttered code -- not something which I expect to have to write in an otherwise very elegant and concise language. Also, I often find myself trudging through code looking for the subtlety that converted my floats to doubles, doubled my memory usage and then caused subsequent float only routines to error out. To those who are constrained to use floats this is awkward and time consuming. To those who are not I would say -- use doubles. The flag that causes an array to be a "space saving array" seems to be a temporary fix (that doesn't mean it was a bad idea -- just that it feels messy and effectively adds complexity that shouldn't be there). It also, mearly postpones the problem as I understand it -- what happens when I multiply two space saving arrays? We simply will never get away from situations where we have to manage the interaction ourselves and so we should be careful not to make that management so awkward (and currently I think it is awkward) that the floats, bytes, shorts, etc. become marginalized in their utility. My suggestion is to go with the rule that a simple hirearchy (in which downcasting is the rule) longs integers shorts bytes cardinals booleans doubles complex doubles <--- default floats complex floats for the most part makes good decisions: Principally because people who are not constrained to conserve memory will use the larger, default types all the time and not wince. They don't *need* floats or bytes. If anyone gives them a float a simple astype('d') or astype('D') to make sure it becomes a double lets them go on their way. Types like integers and bytes are effectively treated as being precise. If you are constrained to conserve memory by staying with floats or bytes instead of just reading things in from disk and making them doubles it will not be so awkward to manage the types in large programs. If I use someones code and they have a scalar anywhere in it at some point, even if I (or they) cast the output, memory usage swells at least for intermediate calculations. Effectively, python *has* 4-byte floats but programming with them is awkward. This means, of course, that multiplying a float array by a double array produces a float. Multiplying a double array by anything above it produces a double. etc. For my work, if I have a float anywhere in the calculation I don't believe precision beyond that in the output so getting a float back is reasonable. I know that some operations produce "more precision" and so I would cast the array if I needed to take advantage of that. Perhaps the downcasting is *not* the way to go. However, I definately think the current awkwardness should be eliminated. I hope my comments will not be percieved as being critical of the original language designers. I find python to be very useful or I wouldn't have bothered to make the comments at all. -Todd Pitts

1 0

Re: [Matrix-SIG] Re: Matrix-SIG digest, Vol 1 #364 - 9 msgs
by Greg Kochanski Feb. 9, 2000

Feb. 9, 2000

"Andrew P. Mullhaupt" wrote: > > It's most important that the rules be simple, and (preferably) close > > to common languages. I'd suggest C. > > That is a good example of a language which has a pretty weird history on > this particular matter. True. The only real advantage of C is that so many people are used to it. Don't forget the human element. FORTRAN would also be a reasonable choice. There's a big cost to learning a new language; if it gets too big, people simply won't use Python. > > Since automatic typecasting only buys a small improvement > > in ease of use, I'd want to be extremely sure that it doesn't cause > > many problems. > > Au contraire. It is a huge win. Try writing a "generic" function with six > arguments which can sensibly be integers, or single or double precision > variables. If you have to test variables to see what they are, then you have > to essentially write a table driven typecaster. If, as in Fortran, you have > to write different functions for different argument types then you have the > dangerous programming practice of having several different pieces of code > which do essentially the same computation. While that's nice to say, it doesn't really translate completely to practice. A lot of functions don't make sense with arbitrary objects; and some require subtle changes. For instance, who wants a matrix inversion function that operates on integers, using integer division inside? Lots of functions have DOUBLE_EPS or FLOAT_EPS embedded inside them. One has to change the small number when you change the data type. I'll grant you that running things with both doubles or floats is often useful. I'd be happy with automatic upcasting among them. I'd be moderately happy with upcasting among the integers. I really don't see any crying need to mix integers with floating point numbers. I'd like some examples to make me believe that mixing ints and floats is a 'huge win'.

1 0

Re: Numpy-discussion digest, Vol 1 #10 - 10 msgs
by Greg Kochanski Feb. 9, 2000

Feb. 9, 2000

> From: "Andrew P. Mullhaupt" <amullhau(a)zen-pharaohs.com> > > The upcasting rule thus ensures that > > > > 1) No precision is lost accidentally. > > More or less. > > More precisely, it depends on what you call an accident. What happens when > you add the IEEE single precision floating point value 1.0 to the 32-bit > integer 2^30? A _lot_ of people don't expect to get the IEEE single > precision floating point value 2.0^30, but that is what happens in some > languages. Is that an "upcast"? Would the 32 bit integer 2^30 make more > sense? Now what about the case where the 32 bit integer is signed and adding > one to it will "wrap around" if the value remains an integer? Because these > two examples might make double precision or a wider integer (if available) > seem the correct answer, suppose it's only one element of a gigantic array? > Let's now talk about complex values.... > It's most important that the rules be simple, and (preferably) close to common languages. I'd suggest C. In my book, anyone who carelessly mixes floats and ints deserves whatever punishment the language metes out. I've done numeric work in languages where casting was by request _only_ (e.g., Limbo, for Inferno), and I found, to my surprise, that automatic type casting these type casting is only a mild convenience. Writing code with manual typecasting is surprisingly easy. Since automatic typecasting only buys a small improvement in ease of use, I'd want to be extremely sure that it doesn't cause many problems. It's very easy to write some complicated set of rules that wastes more time (in the form of unexpected, untraceable bugs) than it saves. By the way, automatic downcasting has a hidden problems if python is ever set to trap underflow errors. I had a program that would randomly crash every 10th (or so) time I ran it with a large dataset (1000x1000 linear algebra). After days of hair-pulling, I found that the matrix was being converted from double to float at one step, and about 1 in 10,000,000 of the entries was too small to represent as a single precision number. That very rare event would underflow, be trapped, and crash the program with a floating point exception.

1 0

Come take an informal survey.
by Travis Oliphant Feb. 8, 2000

Feb. 8, 2000

In an effort to try and get data about what users' attitudes are toward Numerical Python, I'm conducting a survey at sourceforge.net If you would like to participate in the survey, please go to http://www.sourceforge.net, log-in with your sourceforge id and go to the numpy page: http://sourceforge.net/project/?group_id=1369 In the Public Survey section there is a short survey you can fill out. Thank you, Travis Oliphant NumPy Developer

1 0

Re: [Matrix-SIG] An Experiment in code-cleanup.
by Paul Barrett Feb. 8, 2000

Feb. 8, 2000

Travis Oliphant writes: > > > > 1) The re-use of temporary arrays -- to conserve memory. > > Please elaborate about this request. When Python evaluates the expression: >>> Y = B*X + A where A, B, X, and Y are all arrays, B*X creates a temporary array, T. A new array, Y, will be created to hold the result of T + A, and T will be deleted. If T and Y have the same shape and typecode, then instead of creating Y, T can be re-used to conserve memory. > > > > 2) A copy-on-write option -- to enhance performance. > > > > I need more explanation of this as well. This would be an advanced feature of arrays that use memory-mapping or access their arrays from disk. It is similar to the secondary cache of a CPU. The data is held in memory until a write request is made. > > > > 3) The initialization of arrays by default -- to help novices. > > What kind of initialization are you taking about (we have zeros and ones > and random already). For mixed-type (or object) arrays containing strings, zeros() and ones() would be confusing. Therefore by default, integer and floating types are initialized to 0 and string types to ' ', and the option would be available to not initialize the array for performance. > > > > 4) The creation of a standard API -- which I guess is assumed, if it > > is to be part of the Python standard distribution. > > Any suggestions as to what needs to be changed in the already somewhat > standard API. No, not exactly. But the last time I looked, I thought some improvements could be made to it. > > > > 5) The inclusion of IEEE support. > > This was supposed to be there from the beginning, but it didn't get > finished. Jim's original idea was to have two math modules, one which > checked and gave error's for 1/0 and another that returned IEEE inf for > 1/0. > > The current umath does both with different types which is annoying. When I last spoke to Jim about this at IPC6, I was under the impression that IEEE support was not fully implemented and much work still needed to be done. Has this situation changed since then? > > > > And > > > > 6) Enhanced support for mixed-types or objects. > > > > This last issue is very import to me and the astronomical community, > > since we routinely store data as (multi-dimensional) arrays of fixed > > length records or C-structures. A current deficiency of NumPy is that > > the object typecode does not work with the fromstring() method, so > > importing arrays of records from a binary file is just not possible. > > I've been developing my own C-extension type to handle this situation > > and have come to realize that my record type is really just a > > generalization of NumPy's types. > > > I would like to see the code for your generalized type which would help me > see if there were some relatively painless way the two could be merged. recordmodule.c is part of my PyFITS module for dealing with FITS files. You can find it here: ftp://ra.stsci.edu/pub/barrett/PyFITS_0.3.tgz I use NumPy to access fixed-type arrays and the record type for accessing mixed-type arrays. A common example is accessing the second element of a mixed-type (ie. an object) from the entire array. This returns a record type with a single element, which is equivalent to a NumPy array of fixed type. Therefore users expect this object to be a NumPy array and it isn't. They have to convert it to one. > > two C-extension types merged. I think this enhancement can be done > > with minimal change to the current NumPy behavior and minor changes to > > the typecode system. > > If you already see how to do it, then great. Note that NumPy already has some support for an Object type. It has been proposed that it be removed, because it is not well supported and hence few people use it. I have the contrary opinion and feel we should enhance the Object type and make it much more usable. If you don't need it, then you don't have to use it. This enhancement really shouldn't get in the way of those who only use fixed-type arrays. So what changes to NumPy are needed? 1) Instead of a typecode (or in addition to the typecode for backward compatibility), I suggest an optional format keyword, which can be used to specify the mixed-type or object format. Namely, format = 'i, f, s10', where 'i' is an integer type, 'f' a floating point type, and s10 is a string of 10 characters. 2) Array access will be the same as it is now. For example # Create a 10x10 mixed-type array. A = array((10, 10), format = 'i, f, 10s') # Create a 10x10 fixed-type array. B = array((10, 10), typecode = 'i') # Print a 5x5 subarray of mixed-type. print A[:5,:5] # Print a 5x5 subarray of fixed-type print B[:5,:5] # Or # (Note that the 3rd index is optional for fixed-type arrays, it # always defaults to 0.) print B[:5,:5,0] # Print the second element of the mixed-type of the entire array. # Note that this is now an array of fixed-type. print A[:,:,1] The major thorn that I see at this point is how to reconcile the behavior of numbers and strings during operations. But I don't see this as an intractable problem. I actually believe this enhancement will encourage us to create a better and more generic multi-dimensional array module by concentrating on the behavioral aspects of this extension type. Note that J, which NumPy is base upon, allows such mixed-types. -- Dr. Paul Barrett Space Telescope Science Institute Phone: 410-516-6714 DESD/DPT FAX: 410-516-8615 Baltimore, MD 21218

2 1

Re: A proposal for dot product
by John Travers Feb. 8, 2000

Feb. 8, 2000

This message was sent from Geocrawler.com by "John Travers" <j.c.travers(a)durham.ac.uk> Be sure to reply to that address. If the above was implemented, I would be very happy indeed. As a maths student, I use NumPy a lot. And get infuriated with the current implementation. John. Geocrawler.com - The Knowledge Archive

1 0

Re: Numpy-discussion digest, Vol 1 #7 - 1 msg
by Travis Oliphant Feb. 4, 2000

Feb. 4, 2000

> I'm having a problem with PyArray_Check. If I just call > PyArray_Check(args) I don't have a problem, but if I try to assign the > result to anything, etc., it crashes (due to acces violation). So, for > example the code at the end of this note doesn't work, yet I know an array > is being passed and I can, for example, calculate its trace correctly if I > type cast it as a PyArrayObject*. > > Also, a more general question: is this the recommended way to input numpy > arrays when using swig, or do most people find it easier to use more > elaborate typemaps, or something else? I have some experience with SWIG but it is not my favorite method to use Numerical Python with C, since you have so little control over how things get allocated. Your problem is probably due to the fact that you do not run import_array() in the module header. There is a typemap in SWIG that let's you put commands to run at module initialization. Try this in your *.i file. %init %{ import_array(); %} This may help. Best, Travis

1 0