[1] https://numpy.org/doc/stable/user/basics.rec.html

EDIT: I'm always very surprised by the wide acceptance of python/numpy for numerical algorithms, as it seems a notion quite foreign and awkward for this language. Python has excellent native support for string processing and advanced data structures like dictionaries, but those are mostly useless for numerical computation. If you want to do math in Python, you need to write libraries in other languages and provide a Phython interface. This is what scipy does, and it is a very good job, but it does not feel "natural" from the point of view of the language. If Python was a good scientific programming language, the most natural way to implement scipy/numpy and all that would be using Python itself, but this is not at all the case.

To be fair, there's even a standard library `array` module that is a direct wrapper around C arrays, and NumPy shares some functionality with it (namely in regards to buffers, memoryviews, data type encoding etc).

Edit: there's also a built-in buffer type that supports multi-dimensional strides, suboffsets and all that (which, again both `array` and `numpy` make use of in establishing a standardized buffer api).

Paraphrasing Alexander Stepanov, "complexity is part of of the specification" [0]. If your data is not contiguous (as the third figure in the article shows), it is not the same data structure. Thus you cannot implement this data structure in pure Python, period. The "array" module is just a wrapper to C data just like numpy is, it is not a pure Python construction.

[0] Here's the full quote, from http://www.stlport.org/resources/StepanovUSA.html

In 1976, still back in the USSR, I got a very serious case of food poisoning from eating raw fish. While in the hospital, in the state of delirium, I suddenly realized that the ability to add numbers in parallel depends on the fact that addition is associative. (So, putting it simply, STL is the result of a bacterial infection.) In other words, I realized that a parallel reduction algorithm is associated with a semigroup structure type. That is the fundamental point: algorithms are defined on algebraic structures. It took me another couple of years to realize that you have to extend the notion of structure by adding complexity requirements to regular axioms.

[0] https://docs.python.org/3/c-api/buffer.html#c.Py_buffer.stri...

Not at all! The PyBuffer library is all about interfacing to C data. As far as I know, creating a (mutable) contiguous array of floats is impossible in the Python language. This seems to be a deliberate language design decision.

// Pretty much the entirety of CPython is written in C, including lists, dicts and everything else, not much different than the array.

Like I mean sure, maybe it'd be nice to have some syntactic sugar at the language level for tensors (though IMO the slice syntax for numpy.ndarray is already plenty), and maybe it's a microannoyance to have to write `import numpy as np` over and over, but what else is lacking?

To be fair, with modern hardware capabilities, I still wouldn’t be writing out linear algebra algorithms in the language itself even if I were programming in a relatively low-level language like C or C++. I’d be using an implementation of BLAS/LAPACK or some other suitable library, which in turn might well have been written in carefully tuned assembly language on the target platform in order to use any parallel or otherwise specialised operations provided by a CPU or GPU, partition data for optimal cache use, etc.

The situation with numpy and other numerical libraries in Python seems analogous. We are still importing highly optimised code to do the numerical heavy lifting and then using a higher-level language to write the less performance-sensitive logic and glue everything together conveniently.

I disagree strongly with this. Numerical facilities should not be in the stdlib, they should be in the language itself, without need to "import" anything. I can create strings and dictionaries without any imports. I should be able to create a multidimensional array of floats, and perform natural operations with it (like matrix multiplication) without importing anything.

As per the dictionaries in Fortran, there are libraries for that. But I think that not having complex data structures in Fortran is a feature, not a bug. If you actually need these data structures, many Fortran programmers will tell you to just use a different language. You'll never hear such a plainly honest answer from Python programmers; they'll tell you instead to use bizarre libraries with unnatural interfaces, and you'll be forced to multiply matrices using a function "numpy.dot" or an ugly operator "@".

Overloading the language with features that have limited use for most of their users would tie their release cycles and end up a disservice to both users of the language in general and those who use it for heavy number crunching.

As for @, my guess is that Python is running out of ASCII symbols for operators (and unwilling to go full APL with operators). I would imagine that, if NumPy arrays don't support @ (don't they?!), it'd be a desire to be compatible with pre 3.5, which is when the @ infix operator was introduced. It shouldn't be much trouble to implement __matmul__ and __rmatmul__ for them, assuming NumPy's own interfaces are sane.

In fact, one of the common design questions when adding the @-operator to Python was "why do this, when it's not used anywhere in the core language or standard library?" - and the answer was "because it's important for Numpy and the numerical user community".

(and finally, adding __matmul__ and __rmatmul__ is entirely backwards compatible; you can use objects with those methods on Python 3.4 - it's only `A @ B` which is a syntax error)

In Fortran you often need such data structures to keep track of results or to do IO. It’s a pain and I’ve so far not seen a usable library for it.

From dot broadcasting, to multiple dispatch, to parametric types (if you want ), multidimensional array comprehensions. No worrying about 4 different container types (list, np array, pytorch tensor, symbolic etc etc)

Everything works with one array abstraction, including gpu and multi threaded

The way numpy uses python's built in operators makes it feel very native. I really don't find it awkward at all.

I used numeric and numarray in the pre-numpy days, and those did feel more "bolted on". While numpy is really similar to numeric, a lot of little things were fixed during the transition to make numpy very much a native part of python.

Everyone keeps insisting that numpy/scipy/etc are all really "just written in other languages with a python interface", but that's _far_ from being true. Core parts are (e.g. ndarray itself + built-in ufuncs), and a lot more is wrappers around widely-used libraries (think BLAS), but you might be surprised at how much of numpy really is python. Ditto for scipy (though for purely historical reasons, scipy is a dumping ground for X-implemented-in-C stuff). Ditto for things like skimage/sklearn/etc. The bulk of it really is implemented in python. Sure, operations like + aren't, but it might be worth looking through the numpy codebase before stating that it's not written in python.

Next, yes, if you naievely loop over a python array it will be slow, and this is a real limitation of the language + numpy model. Yes, you need to use a different mental model when implementing things. Some problems (e.g. finite-difference-like problems) are best implemented in a way that has poor performance with python+numpy. You do need to drop down to fortran/C/etc (cython is _great_ for those cases) when you have problems that can't be modelled in a particular way. However, that's not the limitation most folks seem to think it is. It's a well-known limitation that there's a ton of support for switching approaches when you hit it (e.g. cython, numba, f2py, etc). Yes, fortran and julia are both nicer in those cases (though my experience with julia is mostly "toy"/learning projects). However, the problems where python+numpy breaks down performance-wise are not anywhere near as common as folks seem to think they are. In my experience, the other types of non-scientific problems are _much_ more common, so it's nice to have a very general purpose language to draw on.

As for matrix multiplication, is it really that bad to do things like x.T.dot(y) instead of x' * y? I really don't understand the argument that the `dot` method of an ndarray is ugly. Also, I actually hate that matlab/et al treat any 2d array as though it's a matrix and fundamentally don't have the concept of a 1d array (row vectors and column vectors are 2d, not 1d). I really do want a 1d array and not a row vector or a column vector most of the time. That's simply impossible in matlab. Also, if you want that behavior in python, it's absolutely possible! Use `np.matrix` instead of `np.array`. Arrays don't behave that way because element-wise operations are more common in practice.

I get the impression you're looking at things from a Julia perspective. Julia is great, but similar to matlab and fortran it falls apart when you start to move outside its core domain. I wouldn't want to write a web service in Julia (though I'm sure it's possible). I do very frequently need to implement web services that do relatively heavy number crunching or deal with things that are fundamentally math-on-big-arrays-of-numbers. Python is _great_ for that type of application. It's also actually pretty good for many "hardcore" number crunching tasks. I've done things where python falls short, sure. However, python + numpy/et al hits a sweet spot that other languages currently don't in what it combines. Domain specific languages are better in some ways, but surprisingly little of a scientific codebase winds up being pure numerical operations. There's an awful lot of other stuff in there, and that's where the scientific python stack is incredibly effective.

See: https://github.com/GenieFramework/Genie.jl

Multiple dispatch + parametric types is a very general and powerful programming model.

It doesn't all apart at all, and anything python is a pycall away.

The "regular product" in math is never intended to be commutative (e.g., the product in a group or in a ring, and that includes the matrix product). If you want to indicate that an operator is commutative, in math, you use always the "+" operator. Thus, it is python that gets it wrong by using a blatantly commutative operator for non-commutative operations such as string concatenation. That python was designed by a mathematician adds sadness to this notational tragedy.

Unfortunately most undergraduates get locked in early with free student versions and continue to use it in industry. Currently I’m trying to expose my students to both Matlab and Python/numpy in a matter of fact manner.

You can do algorithm development in MATLAB or R or Julia or whatever and figure out ways to either translate the logic to Python/C++ later or setup some kind of IPC for production software engineers to call (platform/orchestration stuff is usually Python or Java)...

...OR you can just develop your algorithms in Python and keep code comprehension and platform integration much simpler for everyone

What I've noticed is that the ones who have gained some proficiency in programming can switch to Python in a jiffy. So I don't think it's a life or death choice.

The use of Python in my work group is largely as a problem solving and prototyping tool. The product development teams that we serve have their own software departments, with their own choices of languages and so forth.

Not to mention the IDE is actually very good, and you almost never have to deal with package management etc.

Why no octave? Or julia?

Goes far beyond that

>So, there’s a total of three types of vectors in NumPy: 1D arrays, 2D row vectors, and 2D column vectors.

I think you can have 'vectors' of arbitrarily high dimension, by adding empty axes, or doing, say, .reshape(-1, 1, 1, 1).

I would also tie this point with the rule for broadcasting that dimensions which do not match need to have a value of one. I'm endlessly confused about broadcasting, and find the explicit rule useful.

This is fine, because with high numbers of dimensions you really can't afford to store a dense-matrix representation in RAM and 32 is plenty for any low-dimensional problem.

> Yes, broadcasting is a cool thing, but the docs mostly explains 'how' but not 'why' or 'what for'.

I'd say that, even for the 'how', thinking "OK, I have a [1, 1, 2, 2, 1] and a [3, 2, 1, 2, 5] arrays, which means that element-wise multiplication would result in a [3, 2, 2, 2, 5] array" is very useful in reasoning about broadcasting.

Just relying on simple examples and trying to build an intuition for the general case (which is the route favored in the numpy docs, I'd say) doesn't quite work, at least for me.

> The value of the axis argument is, as a matter of fact, the number of the index in question: The first index is axis=0, the second one is axis=1, and so on

If I'm not wrong - this statement is wrong - the axis does not represent index (or otherwise it would be axis 0 = row, axis 1 = column). The actual idea behind is a little unintuitive, but explained well here: https://aerinykim.medium.com/numpy-sum-axis-intuition-6eb949...

The axis argument gives the index of the axis along which your summate, and which therefore disappears after summation.

There might be a confusion about whether "axis 0" is rows or columns? Its length is the number of rows, but it "points" along columns.

I prefer just calling them axis 0, 1, 2, ... and avoid thinking about rows and columns in numpy, I've found that sometimes avoids confusion.

On the other hand einsum is Crystal clear and not prone to confusion.

> np.einsum("ij -> j", B) # sum along rows to create one column-like array

> np.einsum("ij -> i", B) # sum along columns to create one row-like array

Edit: More Einstein sum fun at https://stackoverflow.com/questions/26089893/understanding-n...

Both of them are of interest and after you were confused once about which one to supply to axis=..., there is no way back to clear the confusion. With einsum there is no confusion.

But why not use whatever suits your problem the best? I use both Python and Julia in my research, and sometimes R when I need some niche statistical packages. "My language is superior to your language" is such a self-limiting mindset.

People are excited about Julia, because frankly it's amazing and a breath of fresh air coming from python. Also, it's really now starting to pick up steam.

Curious though, how are limited by using Julia? Any python or r package is a pycall/rcall away.

Otherwise it's faster, more ergonomic and makes coding fun again

Personally I don't quite get why something like

a[a > 5] = 0

is possible (a is NumPy array).

I think it is confusing to use same name for element wise comparison and the array itself. Are comparison operators overloaded to generate this kind of predicate function or view in other contexts as well or is this some sort of special case that is handled by indexing (__getitem__ call)?

I also found it incredibly confusing the first time I saw this, but like all the numpy magic based on overloading, it’s confusing the first time you see it, but incredibly handy afterwards.

a[lambda i: i > 5] = 0

Probably would expect something like this

np.set(a, lambda i: i > 5, 0)

But that might be too cumbersome when combining values from multiple arrays.

    x[4 5 6] is index access (“at”); also x@4 5 6
    
    x[&x>3] is “x at where x is greater than 3” ; can also be written x@&x>3

[0] https://numpy.org/doc/stable/reference/random/index.html

Anyone have any idea why? I imagine if they could support it they would support it. I wonder if Python doesn't allow this kind of operator overloading?

https://www.python.org/dev/peps/pep-0335/

and deferred PEP for rich comparison chaining specifically

https://www.python.org/dev/peps/pep-0535/

The central issue is that `3<=a<=5` expands to `3<=a and a<=5`, and `and` coerces the variables to booleans in order to work.

>https://www.python.org/dev/peps/pep-0535/

Oh nice! Thank you for sharing. I'm sure the functionality will eventually come.

My response is that, in those cases, go ahead and declare the array you want to catch the result (and it’s shape), then use [:] etc. appropriately to ensure that if the size doesn’t match between computed and expected at the point of assignment, it will fail.

Edit: though now I see it is brought up later on in the matrix section.

For example, if I have a 2-dimensional array and I insert a scalar value by using axis=0, the scalar value needs to be broadcasted to match the column dimension of the array:

  >>> a = np.array([[1,2],[3,4],[5,6]])
  >>> np.insert(a, 1, 11, axis=0)
  array([[ 1,  2],
         [11, 11],
         [ 3,  4],
         [ 5,  6]])

It’s so much more than that. This is not what HN community is about.

3. A thought-provoking text in between, but apparently didn't hit the target in your case.