I would like to be able to use a “diagonal coordinate system” to access array items the same way as with a “column/row system”. In other words I would like to access the elements using indices and slices which specify the diagonal number (the default left-right or the right-left one for anti-diagonals) and the element offset in this diagonal counted from a “root” which I expect to be at the edge of the array. Is there a way to accomplish this in Julia?
Probably. You either need a view that is essentially something that wraps around an existing array or a new index type.
Could you provide an example to provide more specificity?
2 Likes
The issue with the “diagonal coordinate system” is that the range for possible values for the second coordinate depends on the value provided as the first one. So there is no simple formula used in shaping an array as the “width” is in this coordinate system variable. Only the height is constant and equal to numberOfRows+numberOfColumns-1. The width varies between 1 and min(numberOfRows, numberOfColumns). I suppose that it needs support for providing a function to an array as parameter which calculates the flat index out of given indices. This way any possible coordinate system can be established. Python numpy for example does not allow this supporting only strides with constant values.
In other words the question asks which kinds of access to array elements are supported by the native array data structure? Is there a possibility to define a custom function for interpretation of to the array passed indices?
Are you interested in arrays in general, or only in matrices (two-dimensional arrays)?
Yes, as @mkitti explained, this is straightforward.
1 Like
julia> using LinearAlgebra
julia> A = reshape(collect(1:25), 5, 5)
5×5 Matrix{Int64}:
1 6 11 16 21
2 7 12 17 22
3 8 13 18 23
4 9 14 19 24
5 10 15 20 25
julia> struct DiagView{A}
parent::A
end
julia> Base.getindex(v::DiagView, d, o) = diag(v.parent, d)[o]
julia> v = DiagView(A)
DiagView{Matrix{Int64}}([1 6 … 16 21; 2 7 … 17 22; … ; 4 9 … 19 24; 5 10 … 20 25])
julia> v[0,:]
5-element Vector{Int64}:
1
7
13
19
25
julia> v[0,3]
13
julia> v[-2,3]
15
Is this what you had in mind?
6 Likes
May you please explain in detail the provided lines of code line by line?
It seems to do what I am after, but I want to understand how it works and how the Julia arrays differ from Python numpy arrays making this possible?
Looks like he’s defined a custom data structure/type that simply wraps around a matrix. Then he defines a method to tell Julia how to index into this custom type using the diag function built into the LinearAlgebra library. So, you construct regular matrix, use it to define an instance of this new custom type, and then indexing should work as defined.
1 Like
OK … this works for diagonals … may you be so kind to provide the code which demonstrates the same behavior, but addressing right-left(anti-)diagonals? And how to make the indexing the diagonals start with one, and not with zero and positive/negative values? Just simply using diag(v.parent, d - size(v.parent)[1], d)[o] ?
Would the new defined View work as usual in all array operations except getting values using the indices? Or need the underlying array be used for array operations and the View only for getting values? How to assign new values given the new indices? It seems that it needs a further definition …
julia> Base.setindex!(v::DiagView, value, d, o) =
v.parent[diagind(v.parent, d)[o]] = value
julia> A = reshape(collect(1:25), 5, 5)
5×5 Matrix{Int64}:
1 6 11 16 21
2 7 12 17 22
3 8 13 18 23
4 9 14 19 24
5 10 15 20 25
julia> v = DiagView(A)
DiagView{Matrix{Int64}}([1 6 … 16 21; 2 7 … 17 22; … ; 4 9 … 19 24; 5 10 … 20 25])
julia> v[0,3] = -1
-1
julia> A
5×5 Matrix{Int64}:
1 6 11 16 21
2 7 12 17 22
3 8 -1 18 23
4 9 14 19 24
5 10 15 20 25
julia> v[3,2] = -2
-2
julia> A
5×5 Matrix{Int64}:
1 6 11 16 21
2 7 12 17 -2
3 8 -1 18 23
4 9 14 19 24
5 10 15 20 25
I’ll leave the rest as an exercise to the reader.
5 Likes
Even better than a custom array here would be to hook into to_indices, I would think.
1 Like
May you please provide some details about what you mean?
1 Like
Do I understand it right that Julia Vector, Array and Matrix data structures consist of flat values data, shape information and pointers to setindex! and getindex functions which calculate the flat index to flat data values to get or set element values? And the difference to Python numpy arrays is, that Numpy arrays do not provide function pointer, but only strides values which can be manipulated from the default values set by shaping the array in order to achieve effects like transposition, flipping left/right and upside/down? So Numpy function for calculating the indices is always the same, but taking strides as parameter, where Julia function can be any function, so it is possible to define a custom one for defining a “diagonal coordinate system” or any other one?
I’m not sure what you ask for here? Many “array operations” are defined by indexing the array, many are not. If you want to index your array both with cartesian indices and with diagonal indices, there must be some way to distinguish them. That is, mkitti’s code above redefines the indices on the DiagView to be your diagonal indices. To use cartesian indices you must use the underlying array.
As an alternative, since you only want diagonal indexing when you fetch content yourself, it is possible to encapsulate the index instead, e.g. like:
struct DI
d::Int
end
Base.getindex(A::Matrix, idx::DI, o) = @views A[diagind(A, idx.d)][o]
# or
# Base.getindex(A::Matrix, idx::DI, o) = diag(A, idx.d)[o]
# or specialize for scalar index
Base.getindex(A::Matrix, idx::DI, o::Integer) = diag(A, idx.d)[o]
A[DI(2), :]
A[DI(1), [1,3]]
A[DI(1), 2]
Alternatively you can build an entire new array type, subtype of AbstractArray. That is some work to make it perform smoothly w.r.t. broadcasting, slicing etc.
As new to Julia I am myself not sure what I am asking for … just in the process of building understanding what is what and why, checking if it is helpful in coming up with custom code, grasping existing one and express new ideas in terms of Julia syntax. In other words I am trying to grasp what you are speaking about deducting from provided code what the underlying mechanism could be to arrive at some understanding. Currently I don’t really know what the code means, even if it is possible for me to customize it without understanding the details what it does and why. In other words the more self-explaining code examples and explanations are there, the higher the chance I will get the right understanding helpful to grasp the “spirit of Julia”.
Ok.
So, one of the main ideas of julia is that of multiple dispatch. Another feature is that many basic operations are defined in julia itself. For instance, array indexing. Whenever you write something like A[i, j] it is translated by the compiler to a function call getindex(A, i, j). This function, which lives in the module Base, which is always available, can be overloaded with your own types.
So in @mkitti’s code above, a new type (struct) DiagView is constructed. Then he overloads Base.getindex so that his code is called whenever a DiagView is indexed.
My code above does something similar. I create a new type DI which just holds an Int. Then I overload Base.getindex for every Matrix indexed by my new type DI as the first index. So whenever a Matrix is indexed with its first index being a DI, my code is called. It uses a standard function diagind to create an ordinary index into A to create a view into A, which is then indexed by o. Then I also overload the special case that o is an Integer (not a Colon or Vector), so that it return the element of the Matrix, not a view into the matrix.
There is a more general way to do this, that’s the reference to to_indices.
Julia has a help-system, you can write e.g. ?to_indices to access it. There is also a nice macro called @less (and @edit, and @which) which you can use to inspect the source code for a function call. You can e.g. use @less A[2, 3] to see the source of the getindex method being called for this indexing operation.
2 Likes
There is no standard function for finding anti-diagonals, as far as I know. There is a way to create a new matrix with the columns in reverse order: reverse(A, dims=2). However, if you look at the diagind function, you’ll see that it just computes a range for indexing as a linear index (in column major order).
julia> A = reshape(collect(1:25), 5, 5)
julia> di = diagind(A, 2)
11:6:23
julia> A[di]
3-element Vector{Int64}:
11
17
23
With some effort one can create a similar antidiagind function. A method of the diagind function computes cartesian indices:
julia> collect(diagind(A, 2, IndexCartesian()))
3-element Vector{CartesianIndex{2}}:
CartesianIndex(1, 3)
CartesianIndex(2, 4)
CartesianIndex(3, 5)
This one should be fairly easy to use to compute anti diagonal indices.
1 Like