# Discussion on "Why I no longer recommend Julia" by Yuri Vishnevsky

First of all, please don’t do this here. Trolling and flamebaiting are not welcome

Secondly, most of issues in the blog-post were identified and fixed long before the blogpost was ever posted. The community takes fixing bugs pretty seriously, and the forum is a big part of how they coordinate and approach that work.

27 Likes

I think this brings up some other good points.

Yuri, and others who report similar issues, ran into these problems at least partially because they were developers of packages that try to innovate within the ecosystem. Such package development will naturally find the corner cases.

There are reasons I don’t recommend Julia to some people but bugs are not one. There’s a lot of software I use that is difficult to even get working on windows let alone design novel packaging on top of without creating many bugs

2 Likes

I tend to disagree with the OA on this – I don’t want to blame the entire language/ecosystem because of silent bugs involving OffsetArrays and @inbounds.

I still think OffsetArrays are a bad idea, though, and I don’t see the point for them.

In my own code, I usually need to iterate over a matrix, and I normally do things like:

``````for i ∈ 1:size(A,1), j ∈ 1:size(A,2)
A[i,j] = sin(x[i]) * cos(y[j])
end
``````

which is very succinct and easy to understand… and if I want to traverse the diagonal of `A`, this is simply:

``````for i ∈ 1:size(A,1)
A[i,i] = A[i,i] + 1
end
``````

To support OffsetArrays, I should do:

``````for I ∈ CartesianIndices(A)
A[I[1],I[2]] = sin(x[I[1]]) * cos(y[I[2]])
end
``````

Which is ok. But right now, I don’t know how to traverse the diagonal of a matrix in a way that supports OffsetArrays.

There seems to be an explanation in the OffsetArrays.jl home page about how to use `Diagonal` from the LinearAlgebra package (which looks like it’s converting the OffsetArray to a standard array), but what really caught my eyes was:

Certain libraries, such as LinearAlgebra, require arrays to be indexed from 1

``````ERROR: ArgumentError: offset arrays are not supported but got an array with index other than 1
``````

I believe that would be the quickest/wisest solution, just make more base libraries incompatible with OffsetArrays.

Unless, of course, somebody kindly illustrates what is the real need for OffsetArrays, I think that should be the best solution for most packages.

8 Likes

it’s just a package, there’s no “bad idea”, if you don’t like handle it in your package, don’t. You don’t even have to make a warning, because Julia never made the promise that because the multi dispatch would allow you to run generic code on previously unknown types, it would work correctly logically speaking. One just have to know what they’re using / doing – like in any language

9 Likes

I’m sympathetic to Yuri’s larger point: Julia as a language and ecosystem is permissive by default. You can easily combine packages in ways that are unlikely to work and get mired in the weeds. For example, I can easily imagine attempting to differentiate a distributed SVD of a BlockedArray filled with Unitful Quaternions. Sure, that’s an absurd example and it’s unlikely to work, but I can imagine how it might work and may even be able to construct the problem statement in a few lines of code. And heck, if I’m intrepid, I think I could make it work.

That’s also the promise of Julia: you can (sometimes) take two unrelated packages that know nothing about each other and the combination suddenly provides powerful new functionality. Often with no code changes. And it’s amazing.

But sometimes a package doesn’t quite do what the other expects. It’ll typically just error, but OffsetArrays are a very obvious and salient example where it can result in silent corruption due to `@inbounds`.

52 Likes

That’s fair. I meant that incorporating OffsetArrays as some sort of “standard feature” in the official documentation is probably a bad idea. What the LinearAlgebra package is doing is my preferred solution to ensure correctness and turn all these silent bugs into errors, rather than modifying every loop which starts at 1.

4 Likes

Not to get too off-track, but regarding your question about generic matrix iteration, you can use `axes` to get valid indices for a given array. To index into the diagonal, you could do something like this:

``````julia> v = OffsetArray(zeros(3, 4), OffsetArrays.Origin(0, 0));

julia> for (i, j) in zip(axes(v)...)
v[i, j] += 1
end

julia> v
3×4 OffsetArray(::Matrix{Float64}, 0:2, 0:3) with eltype Float64 with indices 0:2×0:3:
1.0  0.0  0.0  0.0
0.0  1.0  0.0  0.0
0.0  0.0  1.0  0.0
``````
5 Likes

To get the diagonal indices, you can do:

``````julia> using LinearAlgebra

julia> diagind(randn(5,5))
1:6:25
``````
1 Like

That sorta works, but is not the best way to do it. What I’d recommend instead is to use the `axes` function.

``````for i ∈ axes(A,1), j ∈ axes(A,2)
A[i,j] = sin(x[i]) * cos(y[j])
end
``````

and for the diagonal, it’s just

``````@assert axes(A, 1) == axes(A,2)
for i ∈ axes(A,1)
A[i,i] = A[i,i] + 1
end
``````
16 Likes

API · ArrayInterface may be of interest

3 Likes

Coincidentally, today’s interview of Leslie Lamport in Quanta has some interesting points regarding correctness, and computer tools to achieve it.

3 Likes

Thanks, `axes` is a nice function.

It’s funny, though, that `axes(A,1)` seems to be itself a 1-based array.

So, if I want to access all but the boundary points in a matrix built with OffsetArrays, I need to use something like `A[axes(A,1)[2:end-1],axes(A,2)[2:end-1]]`, which would be simply `A[2:end-1,2:end-1]` if we were just sticking with 1-based indexing.

Clearly, supporting OffsetArrays comes at a cost. On the other hand, I see no benefits to it. But maybe we are getting a little off-topic and should perhaps continue this discussion of OffsetArrays on a separate thread.

Returning to the original subject, maybe the general issue here is not so much a “design” problem, but simply supporting too much abstraction in Julia base, like if we were trying to please everyone (i.e. “Julia is 1-based indexed by default, but you can use any indexing you like with AbstractArrays”). I would rather favor a more “opinionated” standard that just sticks to 1-based Arrays.

2 Likes

Except for a handful of mentions, the documentation actually goes out of its way to avoid saying OffsetArrays, even in a page clearly about it. Part of the reason is that the AbstractArray interface isn’t supposed to support OffsetArrays specifically, it’s just supposed to provide generic interface methods for indexing (for example, axes isn’t necessarily 1-based). The page also gives you an out with a couple methods that enforce 1-based indexing. However, this does mean there’s a discrepancy between the documented interface and the existing AbstractArray code in Base and packages like StatsBase that enforce 1-based indexing, and you are right that one way to resolve that is to redefine the interface as 1-based. But it seems that Julia is going in the other direction to facilitate some unnamed algorithms, just not as quickly as OffsetArrays users hoped.

5 Likes

No you don’t need to do that. You can simply write

``````A[begin+1:end-1, begin+1:end-1]
``````

Just like how `end` lowers to `lastindex(A, n)`, `begin` will lower to `firstindex(A, n)` inside an indexing expression.

24 Likes

On the composability of packages that have not met.

It is easier to accomplish for focused interoperability than once was the case.
It takes less time to do well than it would working in other languages I know.

It is hard to accomplish seamlessly
– it takes more code than one expects
– that requires more tests than one would prefer
– it requires that others care enough to raise issues

To do so is gratifying.

9 Likes

For example, relevant to the `axes` discussion:

`indices` is useful. It’s like `axes`, except allows multiple arguments (like `eachindex`) so that it can check that the axes match.

`indices((A,B),(2,1))` is like

``````ax = map(axes, (A,B), (2,1))
fax = first(ax)
all(==(fax), Base.tail(ax)) || throw()
return fax
``````

In behavior, except it’ll also promote static size information.

14 Likes

Also worth bearing in mind that there are also situations when it’s perfectly fine to not use the generic interface methods, which can have a far more intimidating vocabulary (`firstindex`, `axes`, `indices`, `eachindex` vs `LinearIndices` vs `CartesianIndices`) than a more straightforward `1:size(A,1)`. It just so happens that the blog addresses a situation where code should be more generic to compose with something the documented interface nominally supports.

If you’re writing a function that takes in a file and writes to another file, you’re completely free to make any internally allocated and freed arrays entirely 1-based. If you are writing a script that doesn’t need to be very composable, you can simply document its 1-based indexing or any other limitations. A benefit of using more numbers and fewer methods is that it’s easier to share with colleagues who may not have the time to learn Julia enough to recognize the methods. Do keep the self-explanatory syntax and methods like `size`/`begin`/`end`.

7 Likes

I have a quite complicated but completely real use case of OffsetArrays, at least `0`-indexed arrays. The problem I face is a research project for something related to Clifford Algebras.

Basically the researcher that I work with though, along with a PhD student, about a representation of Spinors based on a binary representation with a base over {\mathbb{C}^2}^{\lfloor n/2 \rfloor} that works remarkably well when working with operators of Clifford Algebras.

Given the binary representation, the basis of spinors of an n dimensional space is represented by the integer numbers in the range [0, 2^{\lfloor {n/2}\rfloor}-1], and the coefficients are just complex. So a spinor ends up being just an array with indexes beginning at zero, filled with complex numbers. The beauty of the base they choose is that in the most interesting cases for the work of the researcher, most of the coefficients ends up being zero. So you really only need to use sparse arrays.

You could of course use 1-based index and do all the calculations of the basis putting manually the offsets every time you are going to index, but it is a lot easier (and readable) to just have the indexes start at 0. The researcher is not by any means a programmer, but I promise him that he would be able to write code that would make sense, be fast and easy to understand for him. So I need to use and work with 0 indexed arrays.

23 Likes

in some fields: the Zero-based numbering is more optimal than the 1-based
and there are a lot of examples in the Wikipedia article:

the other problem is the fast and optimal Interoperability with the other 0-based languages.

so no easy answer.

3 Likes