Add example of fast non-copying conversion of numeric arrays to docs

[jlapeyre]

It would be useful to have an example of "Fast non-copying conversion of numeric arrays" somewhere in the documentation. For example, I've never been able to figure out how to do this. I imagine some others have not either.

More specifically, I'd like to pass a numpy array from Python to a Julia function where the Julia function receives the data as an Array with no copying.

[jondeuce]

The closest I've been able to come is PyArray(x; copy = false), but it's still somewhat slow on my machine (~100us, regardless of array size), so I'm not sure if that's the best method. And of course you get a PyArray, not an Array, which may or may not matter for your use case.

[jlapeyre]

Regarding PyArray(x; copy = false); where does x come from? You have to pass it from python to Julia and it already is received on the Julia side as a PyArray, I think with copying. In general, it's very important to get a Julia Array efficiently from the data in a numpy array. First, accessing elements of PyArray is slow. Second, the Julia ecosystem only recognizes PyArray as an AbstractArray

[jondeuce]

Ahh, I've only tried from the Julia side. So for example, after calling a Python function from Julia that returns a numpy.ndarray wrapped in a Julia Py struct, I then convert that Py to a PyArray as above. But I agree that it would be great to wrap (something like a pointer to the numpy memory) directly with Array, even if it were unsafe in some sense (e.g. I assume that python would still have control over freeing the memory?).

[cjdoris]

I'd appreciate it if instead of saying "X is slow" you were to specify exactly what X is, how slow it is (with some benchmark code) and why you expect it should be faster.

I'll address the various points separately.

Speed of constructing a PyArray

np = pyimport("numpy")
a = np.random.randn(1000, 1000)
@benchmark b = PyArray(a)

On my computer, the benchmark is about 145us. I have just made a couple of optimisations to make it 16x faster (9us).

You can shave off 1-2us by fully specifying the array type, in this case

b = PyArray{Float64, 2, true, false, Float64}(a)

The best you can expect with the current method is 2us because this is how long it takes for numpy to construct x.__array_interface__, the object containing all the information to wrap the array. So I'm not too far off optimal for this method.

There is also x.__array_struct__, which provides the same information in a different form and only takes 0.5us to access. So at the absolute most PyArray(x) could be 20x faster.

Speed of access

A PyArray is doing exactly what @jondeuce proposed, namely it wraps a pointer to the data plus some shape and type information. If I do this:

@benchmark c = collect(b::PyArray)  # 3.6ms
@benchmark d = collect(c::Array)    # 2.5ms

so collecting a PyArray is a little slower than collecting an Array. BUT this is entirely because the original numpy matrix x is row-major, changing to b = PyArray(x.T) to instead wrap the column-major transpose, then collecting the PyArray is just as fast as collecting the Array. You often want the transpose anyway, because the dimensions in numpy arrays are usually the reverse of Julia arrays, precisely because of the row-/column-major difference.

Just to check that collect isn't somehow cheating, you can also benchmark this function with pretty much the same results.

function mycopy!(desc, src)
    for i in eachindex(dest, src)
        @inbounds dest[i] = src[i]
    end
end
t = zeros(1000, 1000)
@benchmark mycopy!(t, b)
@benchmark mycopy!(t, c)

The copy argument

As the docs say, the copy argument allows PyArray(x) to wrap a copy of x (by calling x.__array__()) if x itself cannot be wrapped. If x is a numpy array, the copy argument should have no effect because x can be wrapped directly.

To put it another way, passing copy=false will sometimes throw an error where copy=true would not have. But if it doesn't throw, the result is no different: the PyArray is a no-copy wrapper around the data in x.

Back to the OP: No-copy Array

There is no such thing, at least not safely. One could use unsafe_wrap but there is no guarantee the memory will be kept alive, and anyway numpy can use a non-contiguous memory layout that is incompatible with what Array expects.

This is why PyArray exists, and as I've hopefully shown, accessing a PyArray is just as fast as accessing an Array.

Controlling the conversion

If you want more control over conversion into Julia functions, wrap a Julia function with PythonCall.pyfunc. When it is called from Python, the function receives its arguments as unconverted Py objects, so you can convert them as you wish.

Alternatively you can call juliacall.convert to convert the arguments before passing them to Julia.

[jondeuce]

I'd appreciate it if instead of saying "X is slow" you were to specify exactly what X is, how slow it is (with some benchmark code) and why you expect it should be faster.

Apologies! I should have been more clear. Unfortunately I don't have the Python knowledge to ask more concrete questions, rather only the (admittedly naive) intuition that Python function calls take ~1us, and so ~100us was more than the ~5us I naively expected.

I'll address the various points separately.

Thanks very much for the detailed explanations.

I have just made a couple of optimisations to make it 16x faster (9us).

Amazing.

By the way, this is a seriously fabulous package. I've just switched over from PyCall and am really enjoying PythonCall, particularly the clearly specified Python dependencies via CondaPkg.

[cjdoris]

Apologies! I should have been more clear. Unfortunately I don't have the Python knowledge to ask more concrete questions, rather only the (admittedly naive) intuition that Python function calls take ~1us, and so ~100us was more than the ~5us I naively expected.

No worries. Your intuition is about correct - in fact I've improved PyArray(x) some more and it's now more like 2.5us for me, or 0.8us if you fully specify the type (the difference is basically some unavoidable type instability).

By the way, this is a seriously fabulous package.

Glad you're enjoying it!

[jlapeyre]

This is a lot of useful information. I want to work through the comments above with examples on my machine. But I'm away from a computer for 10 days. If it still makes sense to do so after working through the comments I will paste some concrete benchmarks. It's great to see the big performance improvements . I'm looking forward to trying it out.

[jondeuce]

I have a quick question regarding your example above fully specifying the PyArray type:

b = PyArray{Float64, 2, true, false, Float64}(a)

If I'm understanding your code correctly, the false in the PyArray type indicates that numpy arrays do not support fast linear indexing; why is this the case? The only thing I can think of is that row-major storage changes the underlying array strides, which is possibly relevant if Julia requires the unit strides for the IndexLinear() trait.

[cjdoris]

Indeed, you require column-major storage for fast linear indexing, because then the ith element in the array is at offset i * strides(array, 1) in memory.

[cjdoris]

You usually want to pass array.T to Julia for this reason, and because the meaning of the dimensions of an array is usually reversed between Julia and Python (because of the difference in contiguity).

[jondeuce]

Thanks for the quick reply!

Indeed, you require column-major storage for fast linear indexing, because then the ith element in the array is at offset i * strides(array, 1) in memory.

Gotcha, makes sense.

You usually want to pass array.T to Julia for this reason, and because the meaning of the dimensions of an array is usually reversed between Julia and Python (because of the difference in contiguity).

I've been converting array -> PyArray -> PermutedDimsArray, which I think should be equivalent to array.T -> PyArray.

[cjdoris]

I've been converting array -> PyArray -> PermutedDimsArray, which I think should be equivalent to array.T -> PyArray.

Sure, doing it Julia-side is fine too - and works for any array, not just Numpy arrays.

[cjdoris]

I think everything has been resolved here? Feel free to reopen if not.

[orenbenkiki]

I'd still like a concise tutorial about how to use zero-copying between Julia and Python in practice:

• Some sample code showing how a Python ndarray can be passed to a Julia function, and a Julia Array can be passed to a Python function, where these functions do "standard things" with these foreign arrays.

• What happens if the same array is passed back and forth multiple times between languages? Do we get a wrapper added each time? Does this harm performance?

• Can one (with extra effort) create an ndarray in Julia and pass it to Python, or an Array in Python and pass it to Julia, so no wrappers are added at all?

I'm looking at this from the point of view of providing a Julia package to be used from Python - that is, the top-level code is in Python (because of the end users requirements), but the implementation is in Julia (because of performance). I know Julia is supposed to solve the "two language problem", but my end users find it difficult enough to migrate from R to Python, so directly using Julia isn't in the cards for years to come. However the point still stands if one works the other way (have the top level code in Julia invoking Python libraries).

My basic concern is that passing arrays between Python and Julia does not end up using the "natural" type for the data, even though "in theory" this should be possible in most cases. That is, zero copy is well and good, but doesn't help me if I'm forced to then copy the data in order to be able to perform any operations on it.

As a trivial example, m += 1 in Python will fail, if m is an juliacall.ArrayValue (that is, zero-copy-wraps a Julia Array). In general it isn't reasonable to expect everything that works with ndarray to also work with juliacall.ArrayValue. On the other side of things, I expect "everything" in Julia to work with PyArray, as it is an AbstractArray after all, but I have no confidence that I will get the same performance as when using regular Array (e.g., if the memory layout is "well behaved", can I apply Strided.jl to PyArray? Is a PyArray a DenseArray? ...)

I spent the last two days digging into why it isn't possible to do a zero-copy conversion between julia Array and Python's ndarray.

It seems this is mostly due to bad design choices Julia made back in the day, when defining the original Array type. Specifically, conflating the array and the memory buffer, and not having an explicit interface for strided arrays. There are numerous threads going back several years bemoaning these issues and suggesting possible solutions. The fact none of these have been adopted so far indicates they won't be in any relevant time frame, so that's out :-(

So we have a bad case of "accidental complexity" on our hands. Still, is it actually impossible, or just "very messy" to implement this?

• Julia arrays like to be column-major and numpy likes to be row-major. But Julia has PermuteDimsArray, and numpy can deal with column-major data just fine. I personally have a lot of Python code that uses column-major layouts. So this isn't a real hurdle.

• Some Python ndarray representations don't map well to Julia's Array, and an arbitrary Julia AbstractArray can have an utterly bizzare internal representation. So I get it that sometimes one has no option but to use a wrapper like PyArray in Julia and juliacall.ArrayValue in Python to compensate for this. But both Julia and Python are nice dynamically-typed languages. It is perfectly OK to provide Julia code with an Array when possible, and a PyArray otherwise (and likewise provide Python an ndarray when possible, and a juliacall.ArrayValue otherwise).

  I know for a fact one can create a well-behaved, strided, contiguous ndarray wrapping any memory - I've done it myself wrapping memory mapped files, for example - so there's no reason it couldn't directly wrap the memory region of a well-behaved, contiguous Julia Array.

  There are also Julia functions like jl_ptr_to_array which can wrap arbitrary memory as a well-behaved, contiguous Julia Array. I haven't used them myself, though, but on its face, it "should" work for wrapping a memory region we get from an ndarray just as well... right...?

• I understand there's a GC issue here. I'm blissfully ignorant on the details of how the GCs of Julia and Python interact.

  On the Python side it is easier to attach an attribute __thing_to_keep_julia_gc_happy__ to an ndarray object, so presumably things could be made 100% transparent.

  On the Julia side... I'd accept having an option to convert an ndarray to a Julia tuple of (an Array sharing the storage, and also an-opaque-thing-one-must-keep-a-handle-to-to-keep-gc-happy), if that's the price of being able to avoid copying my 100MB matrix when I cross the Julia-Python language boundary. But perhaps this could be a new WrapperForNDArray type that would be sufficiently well-behaved (be a DenseArray, work with Strided.jl, etc.), so it would be transparent after all.

I might be able to put some effort into making this work, if all this makes sense. Of course, I might be missing some inconvenient truth about why this is impossible, even in theory...?