mul! is failing with mul!(::Vector, ::Matrix, ::Matrix)

[phajy]

This follow on from #77. If you have loaded the data as described there and try to fit with, e.g, an XSPEC power law you get the following error.

I'm not suer what the problem is but it seems to be related to the fact that we want to get a vector by multiplying a sparse matrix by another matrix. This should be fine and has worked in the past. Could be some issue with the NuSTAR responses being sparse matrices? Although when I tested mul! by hand it works fine with sparse matrices. Could be related to automatic differentiation (?)

I will investigate but am posting this issue here as a starting point.

julia> model = XS_PowerLaw()
┌ XS_PowerLaw
│    K  => (1 ± 0.1)
│    a  => (1 ± 0.1)
└ 

julia> prob = FittingProblem(model => data)
┌ FittingProblem:
│   Models:
│     . XS_PowerLaw{FitParam{Float64}}((1 ± 0.1), (1 ± 0.1))
│   Data:
│     . NuStarData[obs_id=60601017002]
└ 

julia> result = fit(prob, LevenbergMarquadt())
ERROR: MethodError: no method matching mul!(::Vector{Float64}, ::SparseArrays.SparseMatrixCSC{Float64, Int64}, ::Matrix{Float64}, ::Bool, ::Bool)

Closest candidates are:
  mul!(::ChainRulesCore.AbstractThunk, ::Any, ::Any, ::Any, ::Any)
   @ ChainRulesCore ~/.julia/packages/ChainRulesCore/DOJhJ/src/tangent_types/thunks.jl:99
  mul!(::Any, ::ChainRulesCore.AbstractThunk, ::Any, ::Any, ::Any)
   @ ChainRulesCore ~/.julia/packages/ChainRulesCore/DOJhJ/src/tangent_types/thunks.jl:109
  mul!(::Any, ::Any, ::ChainRulesCore.AbstractThunk, ::Any, ::Any)
   @ ChainRulesCore ~/.julia/packages/ChainRulesCore/DOJhJ/src/tangent_types/thunks.jl:110
  ...

Stacktrace:
  [1] mul!
    @ ~/.julia/juliaup/julia-1.10.1+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/matmul.jl:237 [inlined]
  [2] (::SpectralFitting.var"#_transformer!!#58"{…})(output::Vector{…}, energy::Vector{…}, flux::SubArray{…})
    @ SpectralFitting ~/.julia/packages/SpectralFitting/lBW1r/src/datasets/spectraldata.jl:140
  [3] _invoke_and_transform!(cache::SpectralFitting.SpectralCache{…}, domain::Vector{…}, params::Vector{…})
    @ SpectralFitting ~/.julia/packages/SpectralFitting/lBW1r/src/fitting/cache.jl:61
  [4] (::SpectralFitting.var"#f!!#89"{FittingConfig{…}})(domain::Vector{Float64}, parameters::Vector{Float64})
    @ SpectralFitting ~/.julia/packages/SpectralFitting/lBW1r/src/fitting/cache.jl:104
  [5] (::LsqFit.var"#32#34"{SpectralFitting.var"#f!!#89"{…}, Vector{…}, Vector{…}, Vector{…}})(p::Vector{Float64})
    @ LsqFit ~/.julia/packages/LsqFit/OglWj/src/curve_fit.jl:186
  [6] lmfit(f::LsqFit.var"#32#34"{…}, p0::Vector{…}, wt::Vector{…}; autodiff::Symbol, kwargs::@Kwargs{…})
    @ LsqFit ~/.julia/packages/LsqFit/OglWj/src/curve_fit.jl:69
  [7] curve_fit(model::SpectralFitting.var"#f!!#89"{…}, xdata::Vector{…}, ydata::Vector{…}, wt::Vector{…}, p0::Vector{…}; inplace::Bool, kwargs::@Kwargs{…})
    @ LsqFit ~/.julia/packages/LsqFit/OglWj/src/curve_fit.jl:187
  [8] _lsq_fit(f::Function, x::Vector{…}, y::Vector{…}, cov::Vector{…}, parameters::Vector{…}, alg::LevenbergMarquadt{…}; verbose::Bool, max_iter::Int64, kwargs::@Kwargs{…})
    @ SpectralFitting ~/.julia/packages/SpectralFitting/lBW1r/src/fitting/methods.jl:22
  [9] _lsq_fit
    @ ~/.julia/packages/SpectralFitting/lBW1r/src/fitting/methods.jl:11 [inlined]
 [10] fit(prob::FittingProblem{…}, alg::LevenbergMarquadt{…}; verbose::Bool, max_iter::Int64, kwargs::@Kwargs{})
    @ SpectralFitting ~/.julia/packages/SpectralFitting/lBW1r/src/fitting/methods.jl:60
 [11] fit(prob::FittingProblem{FittableMultiModel{…}, FittableMultiDataset{…}, Vector{…}}, alg::LevenbergMarquadt{Float64})
    @ SpectralFitting ~/.julia/packages/SpectralFitting/lBW1r/src/fitting/methods.jl:52
 [12] top-level scope
    @ REPL[6]:1
Some type information was truncated. Use `show(err)` to see complete types.

[fjebaker]

Yes I encountered this bug too but I assumed it was due to uncommited changes. Hmmm... I couldn't work out what was causing it, but it's something since Julia 1.10, as the CI running 1.9 doesn't encounter it.

[fjebaker]

I think JuliaSparse/SparseArrays.jl#469 could be related

[fjebaker]

The problematic line mul!(output, R, f) can be made to work if it's replaced with mul!(output, R, vec(f)), which in theory shouldn't be causing allocations since f is a contiguous single column matrix...

[phajy]

Yes, everything seems to be working if we use vec(f) as you suggest. f is a Matrix{Float64}.

[fjebaker]

Yeah I get all tests passing with vec but it allocates 80 bytes when it should allocate 0. Also it's not unique to SparseMatrix; it seems even with Matrix it fails to find a method

[fjebaker]

MWE:

mat = diagm([1.0, 2.0, 3.0, 4.0, 5.0])
out = zeros(Float64, size(mat, 1))
f = (2 * ones(Float64, size(mat, 1)))
f2 = ones(Float64, (length(f), 1))
f2 .= f
mul!(out, mat, f2)

[fjebaker]

Looking at the signature

mul!(C::AbstractMatrix, A::AbstractVecOrMat, B::AbstractVecOrMat, α::Number, β::Number)

If you change the output array to also be a Matrix, then this works. I don't know why C is not AbstractVecOrMat?

[fjebaker]

This is odd.

From my understanding this is okay: the memory layout should be exactly the same, since Julia is column major, but the wrapping type may differ. Since f2 is a 1 x 5 matrix it is a 2d array, but in memory a vector of 5 elements, and vec(f2) is now a 5 element vector, also a 5 element vector in memory.

[fjebaker]

Examining all signatures with methods(mul!):

 [134] mul!(y::AbstractVector, A::AbstractVecOrMat, x::AbstractVector, alpha::Number, beta::Number)
     @ ~/.julia/juliaup/julia-1.10.1+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/matmul.jl:66
 [135] mul!(C::AbstractMatrix, A::AbstractVecOrMat, B::AbstractVecOrMat, α::Number, β::Number)
     @ ~/.julia/juliaup/julia-1.10.1+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/matmul.jl:263

This feels like a bug?