Matrix contains Infs or NaNs when fitting

[phajy]

We've run into a problem fitting the SWIFT BAT data where the minimisation algorithm runs into an invalid matrix at some point. Here's an example using the SWIFT BAT dataset. The datasets you need to reproduce these are as follows (in a different repository)

Data: https://github.com/phajy/SWIFT_J0909/blob/main/data/swift_bat/swift_bat_157month.pha Response: https://github.com/phajy/SWIFT_J0909/blob/main/data/swift_bat/swiftbat_survey_full_157m.rsp

swift_bat_path = "./data/swift_bat/swift_bat_157month.pha"
SpectralFitting.read_fits_header(swift_bat_path; hdu = 3)
SpectralFitting.OGIP.read_paths_from_spectrum(swift_bat_path)
swift_bat_data = SpectralFitting.OGIPDataset(swift_bat_path; rmf_matrix_index = 2, rmf_energy_index = 3)
drop_bad_channels!(swift_bat_data)
regroup!(swift_bat_data)
normalize!(swift_bat_data)

plot(swift_bat_data, xaxis = :log, xrange=[1,200], yaxis = :log, yrange=[1e-8,1e-5])
model = XS_PowerLaw()
prob = FittingProblem(model => swift_bat_data)
result = fit(prob, LevenbergMarquadt())

gives the following error

ERROR: ArgumentError: matrix contains Infs or NaNs
Stacktrace:
  [1] chkfinite
    @ ~/.julia/juliaup/julia-1.10.1+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/lapack.jl:86 [inlined]
  [2] getrf!(A::Matrix{Float64}; check::Bool)
    @ LinearAlgebra.LAPACK ~/.julia/juliaup/julia-1.10.1+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/lapack.jl:559
  [3] getrf!
    @ ~/.julia/juliaup/julia-1.10.1+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/lapack.jl:557 [inlined]
  [4] #lu!#158
    @ ~/.julia/juliaup/julia-1.10.1+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/lu.jl:82 [inlined]
  [5] lu!
    @ ~/.julia/juliaup/julia-1.10.1+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/lu.jl:81 [inlined]
  [6] #lu#164
    @ ~/.julia/juliaup/julia-1.10.1+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/lu.jl:300 [inlined]
  [7] lu (repeats 2 times)
    @ ~/.julia/juliaup/julia-1.10.1+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/lu.jl:299 [inlined]
  [8] \(A::Matrix{Float64}, B::Vector{Float64})
    @ LinearAlgebra ~/.julia/juliaup/julia-1.10.1+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/generic.jl:1124
  [9] levenberg_marquardt(df::NLSolversBase.OnceDifferentiable{…}, initial_x::Vector{…}; x_tol::Float64, g_tol::Float64, maxIter::Int64, maxTime::Float64, lambda::Float64, tau::Float64, lambda_increase::Float64, lambda_decrease::Float64, min_step_quality::Float64, good_step_quality::Float64, show_trace::Bool, store_trace::Bool, lower::Vector{…}, upper::Vector{…}, avv!::Nothing)
    @ LsqFit ~/.julia/packages/LsqFit/OglWj/src/levenberg_marquardt.jl:200
 [10] levenberg_marquardt
    @ ~/.julia/packages/LsqFit/OglWj/src/levenberg_marquardt.jl:79 [inlined]
 [11] #lmfit#15
    @ ~/.julia/packages/LsqFit/OglWj/src/curve_fit.jl:82 [inlined]
 [12] lmfit
    @ ~/.julia/packages/LsqFit/OglWj/src/curve_fit.jl:75 [inlined]
 [13] lmfit(f::LsqFit.var"#32#34"{…}, p0::Vector{…}, wt::Vector{…}; autodiff::Symbol, kwargs::@Kwargs{…})
    @ LsqFit ~/.julia/packages/LsqFit/OglWj/src/curve_fit.jl:72
 [14] lmfit
    @ ~/.julia/packages/LsqFit/OglWj/src/curve_fit.jl:54 [inlined]
 [15] curve_fit(model::SpectralFitting.var"#f!!#90"{…}, xdata::Vector{…}, ydata::Vector{…}, wt::Vector{…}, p0::Vector{…}; inplace::Bool, kwargs::@Kwargs{…})
    @ LsqFit ~/.julia/packages/LsqFit/OglWj/src/curve_fit.jl:187
 [16] _lsq_fit(f::Function, x::Vector{…}, y::Vector{…}, cov::Vector{…}, parameters::Vector{…}, alg::LevenbergMarquadt{…}; verbose::Bool, max_iter::Int64, kwargs::@Kwargs{…})
    @ SpectralFitting ~/Documents/GitHub/SWIFT_J0909/dev/SpectralFitting/src/fitting/methods.jl:22
 [17] _lsq_fit
    @ ~/Documents/GitHub/SWIFT_J0909/dev/SpectralFitting/src/fitting/methods.jl:11 [inlined]
 [18] #fit#120
    @ ~/Documents/GitHub/SWIFT_J0909/dev/SpectralFitting/src/fitting/methods.jl:60 [inlined]
 [19] fit(prob::FittingProblem{FittableMultiModel{…}, FittableMultiDataset{…}, Vector{…}}, alg::LevenbergMarquadt{Float64})
    @ SpectralFitting ~/Documents/GitHub/SWIFT_J0909/dev/SpectralFitting/src/fitting/methods.jl:52
 [20] top-level scope
    @ REPL[9]:1
Some type information was truncated. Use `show(err)` to see complete types.

Perhaps the model explores an invalid region of parameter space of becomes zero? I need to do some more investigation but thought I'd start a thread here.

This does work for other datasets, e.g., NuSTAR.

[fjebaker]

I see the problem:

┌ OGIPDataset:
│   . Object              : [no object]
│   . Observation ID      : [no observation id]
│   . Exposure ID         : [no exposure id]
│ SpectralData with 7 active channels:
│   . Chn. E (min/max)    : (14.0, 195.0)
│   . Masked channels     : 0 / 8
│   . Model domain size   : 205
│   . Domain (min/max)    : (10.0, 9000.0)
│ Primary Spectrum:
│  Spectrum: SWIFT[BAT]
│   Units                 : counts s^-1 keV^-1
│   . Exposure time       : 2.704453e8
│   . Channels            : 8
│   . Data (min/max)      : (0.0, 8.4029e-07)
│   . Grouped             : yes
│   . Bad channels        : no
│ Response:
│  Response Matrix (8 x 204) channels:
│   . Chn. E (min/max)    : (14.0, 195.0)
│   . Domain E (min/max)  : (10.0, 9000.0)
│ Background: missing
│ Ancillary: missing

Specifically

│   . Data (min/max)      : (0.0, 8.4029e-07)

There are some bins with zero counts and zero uncertainty, so when it tries to do a matrix inverse it seems to end up with a singular matrix and the thing blows up. I think you'd need to exclude those with a mask

[phajy]

I bumped into this same problem again and it wasn't obvious what was wrong. Then I remembered this thread! Made me think that perhaps we should report an error if there are data bins with zero counts and uncertainty.

UPDATE: turned out to be a dodgy dataset, but might still be worth letting the user know. (This issue is low priority!)