A 100x speedup with unsafe Python

From everything I'm seeing, it follows that Matlab & Fortran have been decimated by Python and C/C++ around 1990s and 2010s, respectively. Of course I could be wrong; any evidence of their still strong position will be greatly appreciated, doubly so evidence that this position is due to stride issues.

(Of course Python provides wrappers to Fortran code, eg FITPACK, but this code specifically was mostly written in the 80s, with small updates in this century, and is probably used more thru the numpy wrappers, stride conversion issues and all, than thru its native Fortran interface)

But format conversion is inevitable because C and FORTRAN has a different axis ordering anyway, isn't it?

What kind of support would you have hoped for?

I don't use Python because of type safety, I use it because it's easy to write. I couldn't give a single fuck about some missed pointers that weren't cleaned up in the 2.5 seconds my program ran. I'm not deploying public code.

I tend to prototype in Python and then just rewrite the whole thing in C with classes if I need the 10000x speedup.

Seems to me if you could do most of the work in Python and then just make the critical loop unsafe and 100x faster then that would certainly have some appeal.

The "unsafe" in the title appears to be used in the sense of "exposing memory layout details", but not in the sense of "direct unbounded access to memory". It's probably not the memory safety issue you are thinking of.

Thing about academia (like computer science and the like) is, you need a programming language you can get across easily, quickly and can get people to produce satisfying results in no time.

Back then it might have been java, then python, till the next language comes around.

The java thing was so apparent, that when you looked at c++ code, you were able to spot the former java user pretty easily. :)

About python, sometime ago, I was looking into arudino programming in c and found someone offering a library (I think) in c (or c++, can't remember). What i remember is, that this person stopped working on it, because they got many requests and questions regarding to be able to use their library in python.

while python is eminently criticable, yosef kreinin has designed his own cpu instruction set and gotten it fabbed in a shipping hardware product, and is the author of the c++fqa, which is by far the best-informed criticism of c++; criticizing him with 'this generation is clueless about anything that isn't Python' seems off the mark

"Why would anyone use" -> "when" usage generally means lots of use cases are being ignored / swept under the rug.

1 billion people exist. Each has a unique opinion / viewpoint.

Why would anyone use unsafe Python when Python doesn't have other features such as type safety and all the debugging tools that have been built for C and C++ over time is beyond me.

C and C++ 'type safety' is barely there (compared to more sane languages like OCaml or Rust or even Java etc). As to why would anyone do that? The question is 'why not?' It's fun.

Is this supposed to be a joke?

Have a look at the linked article to see what meaning of 'unsafe' the author has in mind.

ctypes, c extensions, and CFFI are unsafe.

Python doesn't have an unsafe keyword like Rustlang.

In Python, you can dereference a null pointer and cause a Segmentation fault given a ctypes import.

lancedb/lance and/or pandas' dtype_backend="pyarrow" might work with pygame

As others commented "safety" is a heavily overloaded term, for instance there's "type safety" [1] and "memory safety" [2], and then there's "static vs dynamic typing" [3], and "weak vs strong typing" [4], so talking about types and safety offered by a language can be very nuanced.

I suspect when you say "unsafe by design" you may be referring to the dynamic type checking aspect of Python, although Python has supported for a while type annotations, that can be statically checked by linter-like tools like PyRight.

--

1: https://en.wikipedia.org/wiki/Type_safety

2: https://en.wikipedia.org/wiki/Memory_safety

3: https://en.wikipedia.org/wiki/Type_system#Type_checking

4: https://en.wikipedia.org/wiki/Strong_and_weak_typing

It is possible that my brain had been damaged by convention, but it looks much tidier to have your vector multiplicand on the right of the matrix when doing a matrix vector multiplication y=Ax, where A is the matrix and x is the vector, y is the result.

So, x and y must be columns. So, it is nice if our language has column-order as the default, that way a vector is also just a good old 1d array.

If we did y=xA, x and y would be rows, the actual math would be the same but… I dunno, isn’t it just hideous? Ax is like a sentence, it should start with an upper case!

It also fits well with how we tend to learn math. First we learn about functions, f(x). Then we learn that a function can be an operator. And a matrix can represent a linear operator. So, we tend to have a long history of stuff sitting to the left of x getting ready to transform it.

Numpy was originally designed for mathematicians and scientists, and in those domains convention often lines up better with column major order. For example, vectors are a single column, and the column index comes before the row index in many notations. So using column major order meant familar formulas and algorithms (including translating from fortran code which is column major, or matlab for that matter) are easy to translate to numpy code.

Also, numpy is built on some fortran libraries, like BLAS and LAPACK that assume a column major order. If it took row major input, it would need to transpose matrices in some cases change the order, or rewrite those libraries to use row major form.

Probably because Fortran stores matrices and other multidimensional arrays in column order. Traditionally most numerical computation software was written in Fortran and numpy calls those under the hood. Storing in row order would have meant copying the data to column major order and back for any call to Fortran.

Hang on, doesn’t numpy use C (row major) array ordering by default? checks docs Seems like it does. However, numpy array indexing also follows the maths convention where 2D arrays are indexed by row, column (as does C, by the way), so to access pixel x, y you need to say im[y, x]. And the image libraries where I have toyed with with numpy integration (only Pillow, to be honest) seem to work just fine like this — a row of pixels is stored contiguously in memory. So I don’t quite see why the author claims that numpy stores a column of pixels contiguously, but I have only glanced at the article, so quite probably I have missed something.

I think in math columns as vectors is a more common representation. Especially if we talk about matrix multiplication and linear systems.

Closed source and not relevant.

That article fails to touch on the fundamental issue with its title "The Secrets of Colour Interpolation": RGB values are a nonlinear function of the light emitted (because we are better at distinguishing dark colours, so it's better to allow representing more of those), so to interpolate properly you need to invert that function, interpolate, then reapply. The difference that makes to colour gradients is really striking.

While this is a valid point, cv2.resize doesn't actually implement color conversion in a way addressing this issue; at least in my testing I get identical results whether I interpret the data as RGB or BGR. So if you want to use cv2.resize, AFAIK you can count on it not caring which channel is which. And if you need fast resizing, you're quite likely to settle for the straightforward interpolation implemented by cv2.resize.

Even if the RGB components correspond to sRGB, to linearize you apply the same non-linear function to each component value, independently. So even if you do the interpolation in a linear colorspace, the order of the sRGB components does not matter.

You can exploit memory unsafety in cpython using only built-in methods, no imports needed at all: https://github.com/DavidBuchanan314/unsafe-python/

I think you'll have a problem with BGR data absent ctypes, since you have a numpy array with the base pointing at the first red channel value, and you want an array with the base pointing 2 bytes before the first red channel value. This is almost definitely "unsafe"?.. or does numpy have functions knowing that due to the negative z stride, the red value before the blue value is within the bounds of the original array? I somehow doubt it though it would be very impressive. And I think the last A value is hopelessly out of reach absent ctypes since a safe API has no information that this last value is within the array bounds; more strictly speaking all the alpha values are out of bounds according to the shape and strides of the BGRA array.

I think what you want is a combination of swapaxes and flip (reshape does not turn rows into columns, it turns m input rows into n output rows), but yeah.

Actually let me include a little figure

  original  reshape  swapaxes

  abc       ab       ad
  def       cd       be
            ef       cf

25 line pull request doing what? Supporting this format efficiently is probably more LOC and unlikely to get merged as few need this and it complicates the code. Doing the thing I did in Python inside OpenCV C++ code instead (reinterpreting the data) is not really possible since you have less knowledge about the input data at this point (like, you're going to assume it's an RGBA image and access past the area defined as the allocated array data range for the last A value? The user gives a 3D array and you blithely work on the 4th dimension?)

And what about the other examples I measured, including pure numpy code? More patches to submit?

You need to modify base pointer in this example, specifically to point 2 bytes before the current base pointer (moving back from the first red pixel value to the first blue value.) I don't think you can do it without ctypes, maybe I'm wrong.

What does Cython do better than numbs except static compilation? Honest q, I know little about both

I dunno if it's "more than that" since they're not directly comparable? What you describe sounds more like numba or Cython than what TFA describes which is a different use case?..

Image byte order, axis directions, coordinate system handedness when in 3D... after enough trying, you eventually figure out the order of looping at any given stage of your program, and then it Just Works, and then you never touch it again.

I just measured this with the np.ascontiguousarray call (inside the loop of course, if you do it outside the loop when i2 is assigned to, then ofc you don't measure the overhead of np.ascontiguousarray) and the runtime is about the same as without this call, which I would expect given everything TFA explains. So you still have the 100x slowdown.

TFA links to a GitHub repo with code resizing non-zero images, which you could use both to easily benchmark your suggested change, and to check whether the output is still correct.