Will Julia likely ever allow negative indexing of arrays

[Mark Sherlock]

Hi,

I work in computational physics. The main reason we all use Fortran in my area is because it allows arrays to have negative indices. This is very useful when solving some partial differential equations (in e.g. plasma physics, astrophysics, fluid mechanics).

I and my colleagues frequently consider alternative languages but in the end never change due to the headaches involved regarding this. Since Julia seems to be focused on computational science,
I am wondering how likely it is that this would ever be implemented, and/or how we could encourage the developers to do this?

In all other areas Julia looks fantastic for our needs!

[Matt Bauman]

There has been a lot of discussion about this in the past few weeks

https://groups.google.com/d/msg/julia-users/ScwXMfQIBGs/wD1HTXeZBQAJ

https://groups.google.com/d/msg/julia-users/fNisYpMdZ6o/DvFaQi_ZBAAJ

TL;DR, yes, it is possible, but it takes some care since it's violating a fairly well-entrenched assumption about how arrays behave.

[Tom Breloff]

Do you need the bracket notarion 'x[-5]'? This would be best implemented as a package with explicit get/set, as Matt implied... As otherwise you risk some tricky bugs. Also if you're implementing "array-like" types, I would definitely use 0.4+. 

[Daniel Carrera]

I don't know anything about PDEs for plasma physics, but you made me curious. What do you use negative indices for? Are arrays supposed to wrap around?

You can currently use the `end` keyword to write `foo[end-5]`, but I assume that this is not useful to you. Right?

[John Gibson]

Probably it's for something like indexing by Fourier wavenumber.  You represent a real-valued periodic function as a linear combination of Fourier modes exp(2 pi i k x/L) where k varies from -K to K-1. The primary representation in data is an array of complex coefficients for those modes. It might be convenient to index that array using k over the same range. 

Personally, I'm ok with using 0-based or 1-based indices to index into the array and translating into Fourier wavenumber as necessary. In C++ I just have a couple array index -> wavenumber and wavenumber -> array index member functions in my Fourier-expansion classes. 

John

[Scott Jones]

Could you elaborate on what sort of tricky bugs that would cause?  Thanks!

[Kevin Squire]

An older attempt at this can be found here:

https://gist.github.com/alsam/8283205

[Mark Sherlock]

Thanks for all the replies, I'll check the links provided.
Just to give a summary of the reasons.... As John Gibson pointed out Fourier expansions would be one example.

When solving PDE's we use "ghost cells" to provide boundary conditions. A physical variable, let's say density, could be represented in space by a 3D array,
with each spatial index (ix,iy,iz) taking on values like ix=1..nx, iy=1..ny, iz=1..nz etc. If the boundary is, say, reflective, then you need to set the value of the density in
the part of space you can't "see", e.g. in the regions ix=-2..0 and ix=nx+1..nx+3

Obviously this can be handled by offsetting each array index by some amount. The problem is mainly that what we have written down on paper, in books and papers on
numerical analysis etc are numerical equations which naturally count from 1 to n along some dimension. Couple that to the fact that many codes rely on arrays that
consist of mixtures of real space (3D say) and special function expansions, with Fourier expansion being one example.
I for example use a 7-dimensional array to describe some problems (3 dimensions for real space, and 3 dimensions for velocity space). Real space is not expanded,
so it counts from 1..n naturally, while velocity space has 2 dimensions expanded with one type of series (counting from 0..n) and the other dimension expanded in
a different series (counting from -n..n). The final seventh dimension naturally has only 2 indices 0 and 1 (which accounts nicely for the real and imaginary parts of the function).

Within each array loop, it is nice to be able to use if statements that correspond to what you have written down on paper, e.g.:

if(in>0) then blah
if(im<1) then something else

or something like:

do i=0,1
do im=-mmax,mmax
do in=max(im,0),nmax
do ir=1,nr
do ix=1,nx
do iy=1,ny
do iz=2,nz-1

blah

enddo
enddo
enddo
enddo
enddo
enddo
enddo

(incidentally, imagine how many tabs would be required to write the above in python!).

Within a code, there may be many such loops and structures with different looping conditions.

When you have a code with 20000-40000 lines that you want to update to a new language, having to change all of these to account for offsetting is too much work for anyone to seriously consider.
Not to mention the fact you can no longer easily code up the difference equations you have written on paper (and possible spent months creating). That in a nutshell is why nobody is switching from
Fortran (even the young guys are writing new codes in Fortran). There are other nice things about Fortran arrays, such as the ability to quickly set all elements of an array to zero with:

a=0.0

I hope that sort of makes it clear the type of problems we deal with?!

Cheers

[Luke Stagner]

Nice to see more plasma physicists using Julia.

a=0.0

You could just do

a[:] = 0.0

to set all the elements to zero or you could do

a=zeros(n)

In regards to the negative indexing I think that while negative indexing may not be a part of base you should be able to replicate the effect with a package e.g. https://github.com/simonster/TwoBasedIndexing.jl

[Mark Sherlock]

Yeah thanks, two-based indexing is not particularly useful on its own...I think until arbitrary indexing is either in the standard release or very simple to implement with a library,
most people I mention won't be motivated to switch....

[data.pu...@gmail.com]

What about negative indices for reverse selection for example x[-[1, 3]] is everything but the first and third element?

DP

[Milan Bouchet-Valat]

Le lundi 05 octobre 2015 à 09:12 -0700, data.pu...@gmail.com a
écrit :

> What about negative indices for reverse selection for example x[-[1,
> 3]] is everything but the first and third element?

The plan is to use a special "negated index" type instead, which would
be clearer and more flexible (e.g. with non-numeric indices). See
https://github.com/JuliaLang/julia/issues/1032


Regards

[Gray Calhoun]

Sorry that this reply is a little late.... In addition to the other replies, this seems like the sort of problem that might be handled by macros: the macro can rewrite a more natural indexing strategy into the array's index and can automatically add boilerplate that checks if the array is the right size, etc. Writing a `@fortran_indexing` macro and putting it in front of your first loop might be enough to convert some code directly.

As an example, I wrote a toy macro for a presentation that does something that might be similar: automatically generate loops to generate random variables from an ARMA time series model. I.e., when you write

@loop_ts y[t+1] = 0.8y[t] + 0.02y[t-2] + e[t+1]

it automatically replaces it with the code

for _t in 3:(length(y) - 1)
    y[_t+1] = 0.8y[_t] + 0.02y[_t-2] + e[_t+1]
end

Like I said, it's just a toy, but I've put the code in a gist if the idea useful:

https://gist.github.com/grayclhn/d948a9d11eb05339ebaa

(I wrote this code a year ago, but it still seems to work fine.) It would be pretty straightforward to extend the code to create they `y` vector automatically, generate e[t+1] from a specified distribution, etc.

Macros are described in the metaprogramming section of the user's guide:

http://julia.readthedocs.org/en/latest/manual/metaprogramming/

--Gray