...

3) Facility for memory-mapped dataspace in arrays.

I'd really like to have that...

This is pretty easy to add but it does require some changes to the underlying structure, So you can expect it.

...

4) Slices become copies with the addition of methods for current strict referencing behavior.

This will break a lot of code, and in a way that will be difficult to debug. In fact, this is the only point you mention which would be reason enough for me not to use your modified version; going through all of my code to check what effect this might have sounds like a nightmare.

I know this will be a sticky point. I'm not sure what to do exactly, but the current behavior and implementation makes the semantics for slicing an array using a sequence problematic since I don't see a way to represent a reference to a sequence of indices in the underlying structure of an array. So such slices would have to be copies and not references, which makes for an inconsistent code.

I see the point of having a copying version as well, but why not implement the copying behaviour as methods and leave indexing as it is?

I want to agree with you, but I think we may need to change the behavior eventually so when is it going to happen?

...

5) Handling of sliceobjects which consist of sequences of indices (so that setting and getting elements of arrays using their index is possible).

Sounds good as well...

This facility is already embedded in the underlying structure. My plan is to go with the original idea that Jim Hugunin and Chris Chase had for slice objects. The sliceobject in python is already general enough for this to work.

...

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.

I suspect that full behaviour-compatibility with scalars is impossible, but I am willing to be proven wrong. For example, Python scalars are immutable, arrays aren't. This also means that rank-0 arrays can't be used as keys in dictionaries.

How do you plan to implement mixed arithmetic with scalars? If the return value is a rank-0 array, then a single library returning a rank-0 array somewhere could mess up a program well enough that debugging becomes a nightmare.

Mixed arithmetic in general is another sticky point. I went back and read the discussion of this point which occured 1995-1996. It was very interesting reading and a lot of points were made. Now we have several years of experience and we should apply what we've learned (of course we've all learned different things :-) ). Konrad, you had a lot to say on this point 4 years ago. I've had a long discussion with a colleague who is starting to "get in" to Numerical Python and he has really been annoyed with the current mixed arithmetic rules. The seem to try to outguess the user. The spacesaving concept helps, but it still seem's like a hack to me. I know there are several opinions, so I'll offer mine. We need simple rules that are easy to teach a newcomer. Right now the rule is farily simple in that coercion always proceeds up. But, mixed arithmetic with a float and a double does not produce something with double precision -- yet that's our rule. I think any automatic conversion should go the other way. Konrad, 4 years ago, you talked about unexpected losses of precision if this were allowed to happen, but I couldn't understand how. To me, it is unexpected to have double precision arrays which are really only carrying single-precision results. My idea of the coercion hierchy is shown below with conversion always happening down when called for. The Python scalars get mapped to the "largest precision" in their category and then normal coercions rules take place. The casual user will never use single precision arrays and so will not even notice they are there unless they request them. If they do request them, they don't want them suddenly changing precision. That is my take anyway. Boolean Character Unsigned long int short Signed long int short Real /* long double */ double float Complex /* __complex__ long double */ __complex__ double __complex__ float Object

...

7) Addition of attributes so that different users can configure aspects of the math behavior, to their hearts content.

You mean global attributes? That could be the end of universally usable library modules, supposing that people actually use them.

I thought I did, but I've changed my mind after reading the discussion in 1995. I don't like global attributes either, so I'm not going there.

...

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

I don't have much time at the moment, but I could still help out with testing etc.

Konrad you were very instrumental in getting NumPy off the ground in the first place and I will always appreciate your input.

On Tue, 8 Feb 2000, Travis Oliphant wrote:

I know there are several opinions, so I'll offer mine. We need simple rules that are easy to teach a newcomer. Right now the rule is farily simple in that coercion always proceeds up. But, mixed arithmetic with a float and a double does not produce something with double precision -- yet that's our rule. I think any automatic conversion should go the other way.

Remark: If you are consistent then you say here that mixed arithmetic with an int and a float/double produces int?! Right? (I hope that I am wrong.)

Boolean Character Unsigned long int short Signed long int short

How about `/* long long */'? Is this left out intentionally?

Real /* long double */ double float

Travis, while you are doing revision on NumPy, could you also estimate the degree of difficulty of introducing column major order arrays? Pearu

Konrad wrote:

In addition to backwards compatibility, there is another argument for keeping indexing behaviour as it is: compatibility with other Python sequence types.

I claim the current Numeric is INconsistent with other Python sequence types:

...

...

x = [1, 2, 3, 4, 5] y = x[2:5] x [1, 2, 3, 4, 5] y [3, 4, 5] y[1] = 7 y [3, 7, 5] x [1, 2, 3, 4, 5]

So, y is a copy of x[2:5], not a reference.

I'd be perfectly happy to have it implemented as a method. Indexing would remain as it is (by reference), and a new method would provide copying behaviour for element extraction and also permit more generalized sequence indices.

I think I can live with that, as long as it _syntactically_ looks like indexing. This is one case where the syntax is more important than functionality. There are things you want to index with indices, etc., and the composition with parenthesis-like (Dyck language) syntax has proved to be one of the few readable ways to do it.

In addition to backwards compatibility, there is another argument for keeping indexing behaviour as it is: compatibility with other Python sequence types. If you have a list of lists, which in many ways behaves like a 2D array, and extract the third element (which is thus a list), then this data is shared with the full nested list.

_Avoiding_ data sharing will eventually be more important that supporting data sharing since memory continues to get cheaper but memory bandwidth and latency do not improve at the same rate. Locality of reference is hard to control when there is a lot of default data sharing, and performance suffers, yet it becomes important on more and more scales as memory systems become more and more hierarchical. Ultimately, the _semantics_ we like will be implemented efficiently by emulating references and copies in code which copies and references as it sees fit and keeps track of which copies are "really" references and which references are really "copies". I've thought this through for the "everything gets copied" languages and it isn't too mentally distressing - you simply reference count fake copies. The "everything is a reference" languages are less clean, but the database people have confronted that problem.

Which reminds me: rank-0 arrays are also incompatible with the nested-list view of arrays.

There are ways out of that trap. Most post-ISO APLs provide examples of how to cope.

...

I know there are several opinions, so I'll offer mine. We need simple rules that are easy to teach a newcomer. Right now the rule is farily simple in that coercion always proceeds up. But, mixed arithmetic

Like in Python (for scalars), C, Fortran, and all other languages that I can think of.

And that is not a bad thing. But which way is "up"? (See example below.)

...

Konrad, 4 years ago, you talked about unexpected losses of precision if this were allowed to happen, but I couldn't understand how. To me, it is unexpected to have double precision arrays which are really only carrying single-precision results.

Most people always hate, and only sometimes detect, when that happens. It specifically contravenes the Geneva conventions on programming mental hygiene.

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.... There's plenty of rough edges like this when you mix numerical types. It's guaranteed that everybody's ox will get gored somewhere.

2) No overflow occurs unless it is unavoidable (the range problem).

...

The casual user will never use single precision arrays and so will not even notice they are there unless they request them. If they do request

There are many ways in which single-precision arrays can creep into a program without a user's attention.

Absolutely.

Suppose you send me some data in a pickled array, which happens to be single-precision. Or I call a library routine that does some internal calculation on huge data arrays, which it keeps at single precision, and (intentionally or by error) returns a single-precision result.

And the worst one is when the accuracy of the result is single precision, but the _type_ of the result is double precision. There is a function in S-plus which does this (without documenting it, of course) and man was that a pain in the neck to sort out. Today, I just found another bug in one of the S-plus functions - turns out that that if you hand a complex triangular matrix and a real right hand side to the triangular solver (backsubstitution) it doesn't cast the right hand side to complex and uses whatever values are subsequent in memoty to the right hand side as if they were part of the vector. Obviously, when testing the function, they didn't try this mixed type case. But interpreters are really convenient for writing code so that you _don't_ have to think about types all the time and do your own casting. Which is why stubbing your head on an unexpected cast is so unlooked for.

I think your "active flag" solution is a rather good solution to the casting problem, because it gives access to a different behaviour in a very explicit way. So unless future experience points out problems, I'd propose to keep it.

Is there a simple way to ensure that no active arrays are ever activated at any time when I use Numerical Python? Later, Andrew Mullhaupt

silently). Consider a function

k0 = 100 k = 99 while k < k0: .. k0 = k k = ...

which refines a calculation until the measure k stops decreasing. This algorithm may terminate when k is a float, but _fail_ when k is a double -- the extra precision may cause the algorithm

I'd call this a buggy implementation. Convergence criteria should be explicit and not rely on the internal representation of data types. Neither Python nor C guarantees you any absolute bounds for precision and value range, and even languages that do (such as Fortran 9x) only promise to give you a data type that is *at least* as big as your specification.

programming is all about. Numerical programmers need to know how big numbers are, and how much significance they have, and optimise calculations accordingly -- sometimes by _using_ the precision of the working types to advantage.

If you care at all about portability, you shouldn't even think about this. Konrad. -- ------------------------------------------------------------------------------- Konrad Hinsen | E-Mail: hinsen@cnrs-orleans.fr Centre de Biophysique Moleculaire (CNRS) | Tel.: +33-2.38.25.55.69 Rue Charles Sadron | Fax: +33-2.38.63.15.17 45071 Orleans Cedex 2 | Deutsch/Esperanto/English/ France | Nederlands/Francais -------------------------------------------------------------------------------

Konrad Hinsen wrote:

...

But back to precision, which is also a popular subject:

but one which even numerical programmers don't seem to understand ...

Some do, some don't.

It is NOT safe to convert floating point from a lower to a higher number of bits.

It is usually safe. Extremely safe. Safe enough that code in which it is _not_ safe is badly designed.

ALL such conversions should be removed for this reason: any conversions should have to be explicit.

I really hope not. A generic function with six different arguments becomes an interesting object in a language without automatic conversions. Usually, a little table driven piece of code has to cast the arguments into conformance, and then multiple versions of similar code are applied.

which refines a calculation until the measure k stops decreasing. This algorithm may terminate when k is a float, but _fail_ when k is a double -- the extra precision may cause the algorithm to perform many useless iterations, in which the precision of the result is in fact _lost_ due to rounding error.

This is a classic bad programming practice and _it_ is what should be eliminated. It is a good, (and required, if you work for me), practice that: 1. All iterations should have termination conditions which are correct; that is, prevent extra iterations. This is typically precision sensitive. But that is simply something that has to be taken into account when writing the termination condition. 2. All iterations should be protected against an unexpectedly large number of iterates taking place. There are examples of iterations which are intrinsically stable in lower precision and not in higher precision (Brun's algorithm) but those are quite rare in practice. (Note that the Fergueson-Forcade algorithm, as implemented by Lenstra, Odlyzko, and others, has completely supplanted any need to use Brun's algorithm as well.) When an algorithm converges because of lack of precision, it is because the rounding error regularizes the problem. This is normally referred to in the trade as "idiot regularization". It is in my experience, invariably better to actually choose a regularization that is specific to the computation than to rely on rounding effects which might be different from machine to machine. In particular, your type of example is in for serious programmer enjoyment hours on Intel or AMD machines, which have 80 bit wide registers for all the floating point arithmetic. Supporting needless machine dependency is not something to argue for, either, since the Cray style floating point arithmetic has a bad error model. Even Cray has been beaten into submission on this, finally releasing IEEE compliant processors, but only just recently.

to put this another way, it is generally bad to keep more digits (bits) or precision than you actually have

I couldn't agree less. The exponential function and inner product accumulation are famous examples of why extra bits are important in intermediate computations. It's almost impossible to have an accurate exponential function without using extra precision - which is one reason why so many machines have extra bits in their FPUs and there is an IEEE "extended" precision type. The storage history effects which result from temporarily increased precision are well understood, mild in that they violate no common error models used in numerical analysis. And for those few cases where testing for equality is needed for debugging purposes, many systems permit you to impose truncation and eliminate storage history issues. Later, Andrew Mullhaupt

it as an example for people who want to learn stuff about mmap. As it stands, there was some similar code I was able to produce at some point. I forget who here has a copy, maybe Konrad, maybe David Ascher.

Later, Andrew Mullhaupt

I did have some of that code, but it was almost 3 years ago and five computers ago. In other words, it's *somewhere*. I'll start a grep, but don't hold your breath... --da