Indexing incorrectness with ℝⁿ2ℝⁿ_function

[BambOoxX]

Hi @baggepinnen . I am not entirely sure this is a bug or a misunderstanding on my part, so I am asking this here instead of the discourse.
I tried to perform QR decompositions of matrices of Particles but I ran into some issues, in particular, but not only with complex particles and with the Xⁿ2Xⁿ_functions.

Let's start with a tallmatrix

H = (10± 0.1) .*rand(ComplexF64,10,5)

The R factor of the QR decomposition of the mean of H gives a 5x5 upper triangular matrix

qr(mean_object(H)).R

And applying the same to H gives agreeing results in terms of mean at least.

qr(H).R

However, in my actual application I run into the "Comparison of uncertain values using comparison mode safe failed" error. So I figured I should use [ℂⁿ2ℂⁿ_function to get a proper decomposition. Let me know if it's not intended for that purpose. But I guessed it was a proper use given the docs.
When applying this function with

MonteCarloMeasurements.ℂⁿ2ℂⁿ_function(x->Matrix(qr(x).R),H)

The results inside the matrix seem to agree with the direct qr call, but the indexing is way off.

So this seems to be a silent incorrectness. Moreover, trying ℝⁿ2ℝⁿ_function on the real part returned an error.

MonteCarloMeasurements.ℝⁿ2ℝⁿ_function(x->Matrix(qr(x).R),real(H))

So it seems both methods have some kind of indexing issue, but only the real version catches it.

Here is an example with some numeric results.

julia> H = (10± 0.1) .*rand(ComplexF64,10,5)
10×5 Matrix{Complex{Particles{Float64, 2000}}}:
     (8.0 ± 0.08)+(5.53 ± 0.055)im    (2.41 ± 0.024)+(3.98 ± 0.04)im   (5.72 ± 0.057)+(9.83 ± 0.098)im    (7.11 ± 0.071)+(0.893 ± 0.0089)im    (5.84 ± 0.058)+(1.84 ± 0.018)im
   (3.66 ± 0.037)+(7.17 ± 0.072)im    (4.47 ± 0.045)+(8.37 ± 0.084)im  (1.28 ± 0.013)+(5.86 ± 0.059)im    (3.47 ± 0.035)+(9.36 ± 0.094)im      (4.46 ± 0.045)+(5.49 ± 0.055)im
    (3.1 ± 0.031)+(8.47 ± 0.085)im     (2.04 ± 0.02)+(7.99 ± 0.08)im   (2.83 ± 0.028)+(4.67 ± 0.047)im    (5.57 ± 0.056)+(2.03 ± 0.02)im     (0.275 ± 0.0027)+(9.41 ± 0.094)im
   (7.73 ± 0.077)+(9.06 ± 0.091)im    (9.63 ± 0.096)+(1.87 ± 0.019)im   (3.05 ± 0.03)+(7.57 ± 0.076)im    (6.21 ± 0.062)+(9.97 ± 0.1)im        (5.72 ± 0.057)+(0.887 ± 0.0089)im
   (2.41 ± 0.024)+(5.61 ± 0.056)im    (1.92 ± 0.019)+(6.72 ± 0.067)im   (9.9 ± 0.099)+(6.09 ± 0.061)im    (8.89 ± 0.089)+(6.93 ± 0.069)im      (7.36 ± 0.074)+(6.25 ± 0.062)im
   (9.29 ± 0.093)+(7.06 ± 0.071)im    (9.94 ± 0.099)+(2.62 ± 0.026)im  (2.14 ± 0.021)+(9.56 ± 0.096)im     (7.99 ± 0.08)+(8.06 ± 0.081)im     (0.71 ± 0.0071)+(0.773 ± 0.0077)im
    (7.97 ± 0.08)+(4.14 ± 0.041)im    (8.32 ± 0.083)+(4.0 ± 0.04)im    (8.79 ± 0.088)+(7.1 ± 0.071)im     (2.28 ± 0.023)+(8.86 ± 0.089)im      (7.27 ± 0.073)+(6.86 ± 0.069)im
   (3.84 ± 0.038)+(6.17 ± 0.062)im    (3.45 ± 0.035)+(3.64 ± 0.036)im  (1.58 ± 0.016)+(6.61 ± 0.066)im    (8.41 ± 0.084)+(5.79 ± 0.058)im       (9.7 ± 0.097)+(3.9 ± 0.039)im
    (8.8 ± 0.088)+(3.49 ± 0.035)im    (3.15 ± 0.031)+(3.38 ± 0.034)im  (4.32 ± 0.043)+(0.518 ± 0.0052)im  (9.32 ± 0.093)+(5.47 ± 0.055)im       (6.7 ± 0.067)+(5.8 ± 0.058)im
 (0.332 ± 0.0033)+(7.59 ± 0.076)im  (0.855 ± 0.0085)+(1.16 ± 0.012)im  (6.27 ± 0.063)+(8.32 ± 0.083)im    (7.95 ± 0.079)+(8.08 ± 0.081)im      (1.67 ± 0.017)+(7.5 ± 0.075)im

julia> qr(mean_object(H)).R
5×5 Matrix{ComplexF64}:
 -28.926+0.0im  -20.8581+2.96627im   -23.652-3.24835im   -27.4408+2.3586im   -20.5518+2.5752im
     0.0+0.0im  -10.9714+0.0im       1.98025+0.453082im   2.37197-4.43032im   -3.1989+4.50696im
     0.0+0.0im       0.0+0.0im      -14.8365+0.0im       -6.97618-1.25375im  -4.25821-3.5999im
     0.0+0.0im       0.0+0.0im           0.0+0.0im        13.3929+0.0im       5.11278-0.901183im
     0.0+0.0im       0.0+0.0im           0.0+0.0im            0.0+0.0im      -11.5312+0.0im

julia> qr(H).R
5×5 Matrix{Complex{Particles{Float64, 2000}}}:
 (-28.9 ± 0.29)+(-0.0)im     (-20.9 ± 0.21)+(2.97 ± 0.03)im    (-23.7 ± 0.24)-(3.25 ± 0.032)im     (-27.4 ± 0.27)+(2.36 ± 0.024)im    (-20.6 ± 0.21)+(2.58 ± 0.026)im
            (0.0)+(0.0)im  (-11.0 ± 0.11)+(-0.0)im              (1.98 ± 0.02)+(0.453 ± 0.0045)im   (2.37 ± 0.024)-(4.43 ± 0.044)im    (-3.2 ± 0.032)+(4.51 ± 0.045)im
            (0.0)+(0.0)im             (0.0)+(0.0)im          (-14.8 ± 0.15)+(-0.0)im               (-6.98 ± 0.07)-(1.25 ± 0.013)im   (-4.26 ± 0.043)-(3.6 ± 0.036)im
            (0.0)+(0.0)im             (0.0)+(0.0)im                     (0.0)+(0.0)im             (13.4 ± 0.13)+(-0.0)im              (5.11 ± 0.051)-(0.901 ± 0.009)im
            (0.0)+(0.0)im             (0.0)+(0.0)im                     (0.0)+(0.0)im                       (0.0)+(0.0)im           (-11.5 ± 0.12)+(-0.0)im

julia> MonteCarloMeasurements.ℂⁿ2ℂⁿ_function(x->Matrix(qr(x).R),H)
10×5 Matrix{Complex{Particles{Float64, 2000}}}:
 (-28.9 ± 0.29)+(0.0)im          (-23.7 ± 0.24)-(3.25 ± 0.032)im     (-20.6 ± 0.21)+(2.58 ± 0.026)im   (0.0)+(0.0)im  (0.0)+(0.0)im
          (0.0)+(0.0)im           (1.98 ± 0.02)+(0.453 ± 0.0045)im   (-3.2 ± 0.032)+(4.51 ± 0.045)im   (0.0)+(0.0)im  (0.0)+(0.0)im
          (0.0)+(0.0)im          (-14.8 ± 0.15)+(0.0)im             (-4.26 ± 0.043)-(3.6 ± 0.036)im    (0.0)+(0.0)im  (0.0)+(0.0)im
          (0.0)+(0.0)im                   (0.0)+(0.0)im              (5.11 ± 0.051)-(0.901 ± 0.009)im  (0.0)+(0.0)im  (0.0)+(0.0)im
          (0.0)+(0.0)im                   (0.0)+(0.0)im              (-11.5 ± 0.12)+(0.0)im            (0.0)+(0.0)im  (0.0)+(0.0)im
 (-20.9 ± 0.21)+(2.97 ± 0.03)im  (-27.4 ± 0.27)+(2.36 ± 0.024)im              (0.0)+(0.0)im            (0.0)+(0.0)im  (0.0)+(0.0)im
 (-11.0 ± 0.11)+(0.0)im          (2.37 ± 0.024)-(4.43 ± 0.044)im              (0.0)+(0.0)im            (0.0)+(0.0)im  (0.0)+(0.0)im
          (0.0)+(0.0)im          (-6.98 ± 0.07)-(1.25 ± 0.013)im              (0.0)+(0.0)im            (0.0)+(0.0)im  (0.0)+(0.0)im
          (0.0)+(0.0)im           (13.4 ± 0.13)+(0.0)im                       (0.0)+(0.0)im            (0.0)+(0.0)im  (0.0)+(0.0)im
          (0.0)+(0.0)im                   (0.0)+(0.0)im                       (0.0)+(0.0)im            (0.0)+(0.0)im  (0.0)+(0.0)im

julia> MonteCarloMeasurements.ℝⁿ2ℝⁿ_function(x->Matrix(qr(x).R),real(H))
ERROR: BoundsError: attempt to access 5×5 Matrix{Float64} at index [26]
Stacktrace:
  [1] getindex
    @ .\essentials.jl:13 [inlined]
  [2] _broadcast_getindex_evalf
    @ .\broadcast.jl:709 [inlined]
  [3] _broadcast_getindex
    @ .\broadcast.jl:682 [inlined]
  [4] getindex
    @ .\broadcast.jl:636 [inlined]
  [5] macro expansion
    @ .\broadcast.jl:1004 [inlined]
  [6] macro expansion
    @ .\simdloop.jl:77 [inlined]
  [7] copyto!
    @ .\broadcast.jl:1003 [inlined]
  [8] copyto!
    @ .\broadcast.jl:956 [inlined]
  [9] copy
    @ .\broadcast.jl:928 [inlined]
 [10] materialize
    @ .\broadcast.jl:903 [inlined]
 [11] _finish_individuals(::Type{Particles}, ::Val{2000}, individuals::Vector{Matrix{…}}, p::Matrix{Particles{…}})
    @ MonteCarloMeasurements D:\bamboo\.julia\packages\MonteCarloMeasurements\hHq6h\src\particles.jl:226
 [12] ℝⁿ2ℝⁿ_function(f::var"#75#76", p::Matrix{Particles{Float64, 2000}})
    @ MonteCarloMeasurements D:\bamboo\.julia\packages\MonteCarloMeasurements\hHq6h\src\particles.jl:245
 [13] top-level scope
    @ REPL[151]:1
Some type information was truncated. Use `show(err)` to see complete types.

[baggepinnen]

The issue is that ℂⁿ2ℂⁿ_function assumes that the dimensions of the input and output are the same (ℂⁿ -> ℂⁿ). It should probably have a check and error if this is not the case

[BambOoxX]

Ah, ok. I tried that syntax

R = MonteCarloMeasurements.ℂⁿ2ℂⁿ_function(x->Matrix(qr(x)).R,H)

because

MonteCarloMeasurements.ℂⁿ2ℂⁿ_function(x->qr(x),H).R

gave an error about keys, I am assuming because a Factorization object is returned

julia> MonteCarloMeasurements.ℂⁿ2ℂⁿ_function(x->qr(x),H)
ERROR: MethodError: no method matching keys(::LinearAlgebra.QRCompactWY{ComplexF64, Matrix{ComplexF64}, Matrix{ComplexF64}})

Closest candidates are:
  keys(::Pkg.Types.Manifest)
   @ Pkg D:\bamboo\AppData\Local\Programs\Julia-1.10.4\share\julia\stdlib\v1.10\Pkg\src\Types.jl:318
  keys(::Core.SimpleVector)
   @ Base essentials.jl:781
  keys(::RegexMatch)
   @ Base regex.jl:245
  ...

Stacktrace:
 [1] eachindex(itrs::LinearAlgebra.QRCompactWY{ComplexF64, Matrix{ComplexF64}, Matrix{ComplexF64}})
   @ Base .\abstractarray.jl:318
 [2] ℂⁿ2ℂⁿ_function(f::var"#7#8", R::Matrix{Complex{Particles{Float64, 2000}}})
   @ MonteCarloMeasurements D:\bamboo\.julia\dev\MonteCarloMeasurements\src\complex.jl:118
 [3] top-level scope
   @ REPL[8]:1

[baggepinnen]

The key to understand the capabilities of this function is in the name, ℂⁿ2ℂⁿ, it takes complex-valued arrays of dimension n in expects the same type of object in return, matrices (10,5) -> (5,5) do not satisfy this, and neither does (10,5) -> factoization.

In #168 I made it slightly more flexible in order so support the (10,5) -> (5,5) though, so hopefully that works in the latest release (out in a few mintes)