Index broken in ManifoldsBase.jl 0.13

[kellertuer]

See https://github.com/JuliaManifolds/Manopt.jl/issues/129

the indexing in this line https://github.com/JuliaManifolds/Manopt.jl/blob/e7725773ffc97dcf9276649223ee9d17dbc936bb/src/functions/Jacobi_fields.jl#L222-L224

is broken for multidimensional indices (i.e. tuples i). That did work with ManifoldsBase.jl 0.12

Here is an smaller example to reproduce the error

using Manifolds
M = PowerManifold(Sphere(2),NestedPowerRepresentation(), 2,2)
p = reshape([
        [0., 0., 1.], [1., 0., 0.], 1/sqrt(2) * [1.,1.,0.], 1/sqrt(2) * [0.,1.,1.]
    ], 2, 2)
A = collect(Manifolds.get_iterator(M))
p[M,A[1]] # does not work but worked earlier
p[M,A[1]...] # does work

To me it seems the corresponding get index/setindex are missing (for tuples)? But even if I add that we seem to collect once too much and the _read at

https://github.com/JuliaManifolds/ManifoldsBase.jl/blob/bd55280f43fb8482726b61350f38443029df69a2/src/PowerManifold.jl#L1013-L1020

(which seems to obtain a Tuple of a Tuple if I just add tuples to getindex/setindex!. fails. So I am a little lost and call @mateuszbaran the Power manifold expert to the rescue.

[mateuszbaran]

I'm not sure why it worked before without splatting i, it isn't supposed to work. jacobi_field is incorrectly implemented (it should use the _read/_write scheme).

[mateuszbaran]

This seems correctly implemented: https://github.com/JuliaManifolds/ManifoldDiff.jl/blob/main/src/Jacobi_fields.jl#L159-L188 .

[kellertuer]

So read(M, p, rep, i) works but a simple p[M,i] for the same index from a for loop we provide would not? I am asking before rewriting all places where I use that style would be a lot of work and this difference also basically renders the p[M,i] stuff useless if it can not be used with the iterators we provide?

[kellertuer]

This seems correctly implemented: https://github.com/JuliaManifolds/ManifoldDiff.jl/blob/main/src/Jacobi_fields.jl#L159-L188 .

Oh I misread this as it is _in_correctly implemented. Strange that this is so different.

Still. Why do we have p[M,i] if it does not work with our iterators (which worked before) and why should we use the more complicated _read which is even interal? I had hoped that our indexing would work with our iterators.

[kellertuer]

Ah and one main difference is, that the adjoint one, that does work, actually also uses i... and not just I, I feel that is a little counterintuitive to use.

But as soon as I have understood, I will try to unify it in Manopt then. My favourite way, however, would be something easy to use to the end user as well.

[mateuszbaran]

p[M,i...] should be OK as a substitute for _read (at least I don't immediately see why it wouldn't be). p[M,i...] for setting an element should also be more or less OK. Splatting is there for consistency with array indexing:

julia> A = randn(3, 3)
3×3 Matrix{Float64}:
 -0.01988   -0.764316   0.981408
  0.218246  -0.544568  -0.513191
  1.52856    1.02999   -0.670607

julia> A[(2,2)]
ERROR: ArgumentError: invalid index: (2, 2) of type Tuple{Int64, Int64}
Stacktrace:
 [1] to_index(i::Tuple{Int64, Int64})
   @ Base ./indices.jl:300
 [2] to_index(A::Matrix{Float64}, i::Tuple{Int64, Int64})
   @ Base ./indices.jl:277
 [3] to_indices
   @ ./indices.jl:333 [inlined]
 [4] to_indices
   @ ./indices.jl:325 [inlined]
 [5] getindex(A::Matrix{Float64}, I::Tuple{Int64, Int64})
   @ Base ./abstractarray.jl:1218
 [6] top-level scope
   @ REPL[55]:1

julia> A[2,2]
-0.5445681210073285

julia> A[CartesianIndex(2,2)]
-0.5445681210073285

I prefer _read/_write because it's specifically designed to be as fast as possible without making it completely unmaintainable but I think indexing should generally work.

[kellertuer]

Ok I just checked the splatting with normal arrays as well, and you are right we should not support tuple indexing then.

I will think about it for a moment, since I also think the p[M,i...] I might use then is more readable, and also getindex/setindex internally use _read/_write.

Thanks for the checks and explanation.

[mateuszbaran]

Yes, p[M,i...] should be fine. Also note that I've linked the ManifoldDiff.jl implementation of Jacobi fields. I might have made some changes since copying it from Manopt.jl, at some point we should remove this duplication.

[kellertuer]

I will check. Yes once Diff is registered we should remove that duplication, though I would still have to check the dependencies then (currently Manopt loads a little bit of stuff is manifolds is loaded but does not depend on it).

[kellertuer]

It seems most dependencies are stuff that should move anyways?

[kellertuer]

(this is getting off topic since the issue is basically resolved).

[mateuszbaran]

It seems most dependencies are stuff that should move anyways?

Which dependencies?

(this is getting off topic since the issue is basically resolved).

Yes, the duplication should be a separate issue :slightly_smiling_face: .

[kellertuer]

Mainly Manifolds, but it seems that might be resolved when Diff is registered (and loaded by both Manifolds and Manopt).

[kellertuer]

I'll close here, I think once ManifoldDiff is released and there is no dependency to Manifolds anymore (because the backends have been moved over) then it is also no problem to remove the duplication with Manopt.

[mateuszbaran]

Removing the Manifolds.jl dependency in ManifoldDiff.jl is actually complicated and would require making some additional interface packages. I'd prefer to avoid having to do that until ManifoldDiff.jl is at least somewhat mature.