search

MagneticResonanceImaging / MRIReco.jl

Julia Package for MRI Reconstruction
https://magneticresonanceimaging.github.io/MRIReco.jl/latest/
Other
85 stars 22 forks source link

Error when using CGNR #110

Closed aTrotier closed 1 year ago

aTrotier commented 1 year ago

This error does not happens with admm / fista solver. I am using ComplexF32 data, not sure it is related.

nested task error: MethodError: no method matching MRIReco.CompositeNormalOp(::SparsityOperators.FFTOp{ComplexF32}, ::SparsityOperators.NormalOp{LinearOperators.LinearOperator{Int64, Int64, LinearOperators.var"#144#146"{Vector{Int64}}, LinearOperators.var"#145#147"{Vector{Int64}}, LinearOperators.var"#145#147"{Vector{Int64}}, Vector{Int64}}, LinearOperators.LinearOperator{ComplexF64, Int64, LinearOperators.var"#134#136"{Vector{ComplexF64}}, LinearOperators.var"#134#136"{Vector{ComplexF64}}, LinearOperators.var"#135#137"{typeof(conj), Vector{ComplexF64}}, Vector{ComplexF64}}, Vector{Int64}}) Closest candidates are: MRIReco.CompositeNormalOp(::S, ::U, ::V) where {S, U, V} at ~/.julia/packages/MRIReco/ehYO5/src/Operators/Composition.jl:89

tknopp commented 1 year ago

Could you provide the script that generates the error?

ADMM and Fista do not use the NormalOp right now.

aTrotier commented 1 year ago

I get that with a custom MP2RAGE sequence. I'll try to generate the error with simulated data

tknopp commented 1 year ago

Ah okay, would you then probably provide a longer stacktrace? I probably can see the issue from that. And is this on all packages dev-ed? If MRIReco is dev-ed so does SparsityOperators at the moment.

tknopp commented 1 year ago

This commit will fix the issue: https://github.com/MagneticResonanceImaging/MRIReco.jl/commit/abf335c871cea52f288b2d7566c0737a9edcb1d1

But I am actually not sure if that is the only issue since I wonder, why our tests did to catch that case. So it would be definitely good to get this case into our test suite.

aTrotier commented 1 year ago

on my side this test is failling :


function testCSRecoMultiCoilCGNR(;N=32,redFac=1.1,type = ComplexF32)
    x = shepp_logan(N)

    # coil sensitivites
    sensMaps = zeros(ComplexF64,N*N,2,1)
    sensMaps[1:floor(Int64, N*N/2),1,1] .= 1.0
    sensMaps[floor(Int64, N*N/2)+1:end,2,1] .= 1.0

    # convert to type
    x = convert.(type,x)
    sensMaps = convert.(type,sensMaps)

    # simulation
    params = Dict{Symbol, Any}()
    params[:simulation] = "fast"
    params[:trajName] = "Cartesian"
    params[:numProfiles] = floor(Int64, N)
    params[:senseMaps] = reshape(sensMaps,N,N,1,2)
    params[:numSamplingPerProfile] = N

    acqData = simulation(x,params)
    Random.seed!(1234)
    acqData = MRIReco.sample_kspace(acqData, redFac, "poisson", calsize=5)

    # reco
    params[:reco] = "multiCoil"
    params[:reconSize] = (N,N)
    params[:senseMaps] = reshape(sensMaps,N,N,1,2)
    params[:regularization] = "L2"       # regularization
    params[:λ] = 1.e-3
    params[:solver] = "cgnr"
    params[:iterations] = 1
    params[:ρ] = 1.0e-1
    params[:absTol] = 1.e-5
    params[:relTol] = 1.e-4

    x_approx = vec(reconstruction(acqData,params))
    @test (norm(vec(x)-x_approx)/norm(vec(x))) < 1e-1
end

testCSRecoMultiCoilCGNR(N=N,redFac = 1.1, type=ComplexF64)
testCSRecoMultiCoilCGNR(N=N,redFac = 1.1, type=ComplexF32)

I'll try with your fix

aTrotier commented 1 year ago

new error :) I guess we can add the test


      nested task error: TaskFailedException

          nested task error: InexactError: Int64(0.6068992520447697 - 0.1552710088505263im)
          Stacktrace:
            [1] Real
              @ ./complex.jl:44 [inlined]
            [2] convert
              @ ./number.jl:7 [inlined]
            [3] setindex!
              @ ./array.jl:966 [inlined]
            [4] _unsafe_copyto!(dest::Vector{Int64}, doffs::Int64, src::Vector{ComplexF64}, soffs::Int64, n::Int64)
              @ Base ./array.jl:253
            [5] unsafe_copyto!
              @ ./array.jl:307 [inlined]
            [6] _copyto_impl!
              @ ./array.jl:331 [inlined]
            [7] copyto!
              @ ./array.jl:317 [inlined]
            [8] copyto!
              @ ./array.jl:343 [inlined]
            [9] copyto!
              @ ./broadcast.jl:954 [inlined]
           [10] copyto!
              @ ./broadcast.jl:913 [inlined]
           [11] materialize!
              @ ./broadcast.jl:871 [inlined]
           [12] materialize!
              @ ./broadcast.jl:868 [inlined]
           [13] mulRestrict!(res::Vector{Int64}, I::Vector{Int64}, v::Vector{ComplexF64}, α::ComplexF64, β::ComplexF64)
              @ LinearOperators ~/.julia/packages/LinearOperators/O1SvF/src/special-operators.jl:160
           [14] #144
              @ ~/.julia/packages/LinearOperators/O1SvF/src/special-operators.jl:182 [inlined]
           [15] mul!(res::Vector{Int64}, op::LinearOperators.LinearOperator{Int64, Int64, LinearOperators.var"#144#146"{Vector{Int64}}, LinearOperators.var"#145#147"{Vector{Int64}}, LinearOperators.var"#145#147"{Vector{Int64}}, Vector{Int64}}, v::Vector{ComplexF64}, α::ComplexF64, β::ComplexF64)
              @ LinearOperators ~/.julia/packages/LinearOperators/O1SvF/src/operations.jl:29
           [16] mul!
              @ ~/.julia/packages/LinearOperators/O1SvF/src/operations.jl:36 [inlined]
           [17] mul!
              @ ~/.julia/packages/SparsityOperators/m5Lmk/src/NormalOp.jl:35 [inlined]
           [18] mul!
              @ ~/.julia/dev/MRIReco/src/Operators/Composition.jl:107 [inlined]
           [19] macro expansion
              @ ~/.julia/dev/MRIReco/src/Operators/Operators.jl:185 [inlined]
           [20] (::MRIReco.var"#157#158"{Vector{MRIReco.CompositeNormalOp{SparsityOperators.FFTOp{ComplexF64}, SparsityOperators.NormalOp{LinearOperators.LinearOperator{Int64, Int64, LinearOperators.var"#144#146"{Vector{Int64}}, LinearOperators.var"#145#147"{Vector{Int64}}, LinearOperators.var"#145#147"{Vector{Int64}}, Vector{Int64}}, LinearOperators.LinearOperator{ComplexF64, Int64, LinearOperators.var"#134#136"{Vector{ComplexF64}}, LinearOperators.var"#134#136"{Vector{ComplexF64}}, LinearOperators.var"#135#137"{typeof(conj), Vector{ComplexF64}}, Vector{ComplexF64}}, Vector{Int64}}, Vector{ComplexF64}}}, Vector{Int64}, Vector{ComplexF64}, Vector{ComplexF64}, Int64})()
              @ MRIReco ./threadingconstructs.jl:258

      ...and 1 more exception.
tknopp commented 1 year ago

Perfect thanks, I will later fix that. With the test that will be possible.

nbwuzhe commented 1 year ago

@alexjaffray and I are also encountering exactly same two issues. We fixed the first one as abf335c but counldn't find which part triggered this conversion from ComplexF32 to Int64. Would appreciate it if that could be investigated, thanks.

tknopp commented 1 year ago

I have a fix locally and push it in some minutes.

By the way, I wonder why you are doing CS with CGNR, which will not give good results.

tknopp commented 1 year ago

And for the records, the issue is

julia> eltype(SamplingOp([1,2,3,4],(4,4)))
Int64

which seems to be technically correct. But it prevents us from generating a proper temporary vector. My fix will circumvent that but I need to think how to tackle that more generally.

alexjaffray commented 1 year ago

By the way, I wonder why you are doing CS with CGNR, which will not give good results.

@nbwuzhe and I aren't doing CS, we get the error when we try to do a routine Cartesian trajectory reconstruction

tknopp commented 1 year ago

Hm, okay, then I am not sure why it generated the SamplingOp. Once I got my fix in, it would be great if you could test your case.

alexjaffray commented 1 year ago

Sounds good @tknopp.

I'll test the fix once it's posted and let you know how it goes.

aTrotier commented 1 year ago

I have a fix locally and push it in some minutes.

By the way, I wonder why you are doing CS with CGNR, which will not give good results.

I used the CGNR solver to perform the zero-filling reconstruction while combining the data with coil sensitivity with iteration = 1. For the real reconstruction, I am using the ADMM solver + Wavelet reg.

tknopp commented 1 year ago

Thats a valid use case of course. The fix is committed: https://github.com/MagneticResonanceImaging/MRIReco.jl/commit/5296327cbb16dbb9ae88921fd2a002e63e96d8a0

CI is hopefully happy soon. Needed to make a release of SparsityOperators.jl

nbwuzhe commented 1 year ago

Hi @tknopp thanks for the quick fix. Just had a trial on our Cartesian data and seems the error is still the same. Would appreciate some extra leads to the source of this error:

` ERROR: TaskFailedException

nested task error: TaskFailedException

    nested task error: InexactError: Int64(-2.2006506f-7 + 1.1055942f-7im)
    Stacktrace:
      [1] Real
        @ .\complex.jl:37 [inlined]
      [2] convert
        @ .\number.jl:7 [inlined]
      [3] setindex!
        @ .\array.jl:843 [inlined]
      [4] _unsafe_copyto!(dest::Vector{Int64}, doffs::Int64, src::Vector{ComplexF32}, soffs::Int64, n::Int64)
        @ Base .\array.jl:235
      [5] unsafe_copyto!
        @ .\array.jl:289 [inlined]
      [6] _copyto_impl!
        @ .\array.jl:313 [inlined]
      [7] copyto!
        @ .\array.jl:299 [inlined]
      [8] copyto!
        @ .\array.jl:325 [inlined]
      [9] copyto!
        @ .\broadcast.jl:977 [inlined]
     [10] copyto!
        @ .\broadcast.jl:936 [inlined]
     [11] materialize!
        @ .\broadcast.jl:894 [inlined]
     [12] materialize!
        @ .\broadcast.jl:891 [inlined]
     [13] mulRestrict!(res::Vector{Int64}, I::Vector{Int64}, v::Vector{ComplexF32}, α::ComplexF32, β::ComplexF32)
        @ LinearOperators C:\Users\TimWu\.julia\packages\LinearOperators\O1SvF\src\special-operators.jl:160
     [14] #144
        @ C:\Users\TimWu\.julia\packages\LinearOperators\O1SvF\src\special-operators.jl:182 [inlined]
     [15] mul!(res::Vector{Int64}, op::LinearOperators.LinearOperator{Int64, Int64, LinearOperators.var"#144#146"{Vector{Int64}}, LinearOperators.var"#145#147"{Vector{Int64}}, LinearOperators.var"#145#147"{Vector{Int64}}, Vector{Int64}}, v::Vector{ComplexF32}, α::ComplexF32, β::ComplexF32)
        @ LinearOperators C:\Users\TimWu\.julia\packages\LinearOperators\O1SvF\src\operations.jl:29
     [16] mul!
        @ C:\Users\TimWu\.julia\packages\LinearOperators\O1SvF\src\operations.jl:36 [inlined]
     [17] mul!
        @ C:\Users\TimWu\.julia\packages\SparsityOperators\m5Lmk\src\NormalOp.jl:35 [inlined]
     [18] mul!
        @ C:\Users\TimWu\.julia\dev\MRIReco\src\Operators\Composition.jl:108 [inlined]
     [19] macro expansion
        @ C:\Users\TimWu\.julia\dev\MRIReco\src\Operators\Operators.jl:188 [inlined]
     [20] (::MRIReco.var"#157#158"{Vector{MRIReco.CompositeNormalOp{SparsityOperators.FFTOp{ComplexF32}, SparsityOperators.NormalOp{LinearOperators.LinearOperator{Int64, Int64, LinearOperators.var"#144#146"{Vector{Int64}}, LinearOperators.var"#145#147"{Vector{Int64}}, LinearOperators.var"#145#147"{Vector{Int64}}, Vector{Int64}}, LinearOperators.LinearOperator{ComplexF32, Int64, LinearOperators.var"#134#136"{Vector{ComplexF32}}, LinearOperators.var"#134#136"{Vector{ComplexF32}}, LinearOperators.var"#135#137"{typeof(conj), Vector{ComplexF32}}, Vector{ComplexF32}}, Vector{Int64}}, Vector{ComplexF32}}}, Vector{Int64}, Vector{ComplexF32}, Vector{ComplexF32}, Int64})()
        @ MRIReco .\threadingconstructs.jl:169

...and 19 more exceptions.

Stacktrace:
  [1] sync_end(c::Channel{Any})
    @ Base .\task.jl:369
  [2] macro expansion
    @ .\task.jl:388 [inlined]
  [3] _produ_diagnormalop(ops::Vector{MRIReco.CompositeNormalOp{SparsityOperators.FFTOp{ComplexF32}, SparsityOperators.NormalOp{LinearOperators.LinearOperator{Int64, Int64, LinearOperators.var"#144#146"{Vector{Int64}}, LinearOperators.var"#145#147"{Vector{Int64}}, LinearOperators.var"#145#147"{Vector{Int64}}, Vector{Int64}}, LinearOperators.LinearOperator{ComplexF32, Int64, LinearOperators.var"#134#136"{Vector{ComplexF32}}, LinearOperators.var"#134#136"{Vector{ComplexF32}}, LinearOperators.var"#135#137"{typeof(conj), Vector{ComplexF32}}, Vector{ComplexF32}}, Vector{Int64}}, Vector{ComplexF32}}}, idx::Vector{Int64}, x::Vector{ComplexF32}, y::Vector{ComplexF32})
    @ MRIReco C:\Users\TimWu\.julia\dev\MRIReco\src\Operators\Operators.jl:186
  [4] mul!
    @ C:\Users\TimWu\.julia\dev\MRIReco\src\Operators\Operators.jl:175 [inlined]
  [5] mul!(y::Vector{ComplexF32}, S::MRIReco.CompositeNormalOp{LinearOperators.LinearOperator{ComplexF32, Int64, MRIReco.var"#61#63"{Matrix{ComplexF32}, Int64, Int64, Int64}, Nothing, MRIReco.var"#62#64"{Int64, Matrix{ComplexF32}, Int64, Int64}, Vector{ComplexF32}}, MRIReco.DiagNormalOp{Vector{MRIReco.CompositeOp{ComplexF32, LinearOperators.LinearOperator{Int64, Int64, LinearOperators.var"#144#146"{Vector{Int64}}, LinearOperators.var"#145#147"{Vector{Int64}}, LinearOperators.var"#145#147"{Vector{Int64}}, Vector{Int64}}, SparsityOperators.FFTOp{ComplexF32}}}, Vector{MRIReco.CompositeNormalOp{SparsityOperators.FFTOp{ComplexF32}, SparsityOperators.NormalOp{LinearOperators.LinearOperator{Int64, Int64, LinearOperators.var"#144#146"{Vector{Int64}}, LinearOperators.var"#145#147"{Vector{Int64}}, LinearOperators.var"#145#147"{Vector{Int64}}, Vector{Int64}}, LinearOperators.LinearOperator{ComplexF32, Int64, LinearOperators.var"#134#136"{Vector{ComplexF32}}, LinearOperators.var"#134#136"{Vector{ComplexF32}}, LinearOperators.var"#135#137"{typeof(conj), Vector{ComplexF32}}, Vector{ComplexF32}}, Vector{Int64}}, Vector{ComplexF32}}}, ComplexF32}, Vector{ComplexF32}}, x::Vector{ComplexF32})
    @ MRIReco C:\Users\TimWu\.julia\dev\MRIReco\src\Operators\Composition.jl:108
  [6] iterate(solver::RegularizedLeastSquares.CGNR{Vector{ComplexF32}, ComplexF32, Nothing}, iteration::Int64)
    @ RegularizedLeastSquares C:\Users\TimWu\.julia\packages\RegularizedLeastSquares\OA9GQ\src\CGNR.jl:172
  [7] iterate
    @ C:\Users\TimWu\.julia\packages\RegularizedLeastSquares\OA9GQ\src\CGNR.jl:167 [inlined]
  [8] iterate
    @ .\iterators.jl:159 [inlined]
  [9] iterate
    @ .\iterators.jl:158 [inlined]
 [10] solve(solver::RegularizedLeastSquares.CGNR{Vector{ComplexF32}, ComplexF32, Nothing}, u::Vector{ComplexF32}; S::MRIReco.CompositeOp{ComplexF32, LinearOperators.LinearOperator{ComplexF32, Int64, LinearOperators.var"#134#136"{Vector{ComplexF32}}, LinearOperators.var"#134#136"{Vector{ComplexF32}}, LinearOperators.var"#135#137"{typeof(conj), Vector{ComplexF32}}, Vector{ComplexF32}}, MRIReco.CompositeOp{ComplexF32, MRIReco.DiagOp{ComplexF32}, LinearOperators.LinearOperator{ComplexF32, Int64, MRIReco.var"#61#63"{Matrix{ComplexF32}, Int64, Int64, Int64}, Nothing, MRIReco.var"#62#64"{Int64, Matrix{ComplexF32}, Int64, Int64}, Vector{ComplexF32}}}}, startVector::Vector{ComplexF32}, weights::Vector{ComplexF32}, solverInfo::SolverInfo{ComplexF32}, kargs::Base.Iterators.Pairs{Symbol, Any, NTuple{10, Symbol}, NamedTuple{(:normalizeReg, :solver, :senseMaps,

:regularization, :sparseTrafoName, :λ, :reco, :ρ, :iterations, :reconSize), Tuple{Bool, String, Array{ComplexF32, 4}, String, String, Float64, String, Float64, Int64, Tuple{Int64, Int64}}}}) @ RegularizedLeastSquares C:\Users\TimWu.julia\packages\RegularizedLeastSquares\OA9GQ\src\CGNR.jl:153 [11] macro expansion @ C:\Users\TimWu.julia\dev\MRIReco\src\Reconstruction\IterativeReconstruction.jl:196 [inlined] [12] (::MRIReco.var"#215#216"{AcquisitionData{Float32}, Tuple{Int64, Int64}, Vector{Regularization}, Vector{Vector{ComplexF32}}, String, Array{ComplexF32, 4}, Nothing, Dict{Symbol, Any}, Dict{Any, Any}, Int64, Int64, Int64, Int64})() @ MRIReco .\threadingconstructs.jl:169

...and 14 more exceptions.

Stacktrace: [1] sync_end(c::Channel{Any}) @ Base .\task.jl:369 [2] macro expansion @ .\task.jl:388 [inlined] [3] reconstruction_multiCoil(acqData::AcquisitionData{Float32}, reconSize::Tuple{Int64, Int64}, reg::Vector{Regularization}, sparseTrafo::Nothing, weights::Vector{Vector{ComplexF32}}, solvername::String, senseMaps::Array{ComplexF32, 4}, normalize::Bool, encodingOps::Nothing, params::Dict{Symbol, Any}) @ MRIReco C:\Users\TimWu.julia\dev\MRIReco\src\Reconstruction\IterativeReconstruction.jl:181 [4] reconstruction(acqData::AcquisitionData{Float32}, recoParams::Dict{Symbol, Any}) @ MRIReco C:\Users\TimWu.julia\dev\MRIReco\src\Reconstruction\Reconstruction.jl:43 [5] top-level scope @ .\timing.jl:210 [inlined] [6] top-level scope @ d:\SpiralDiffusion\GIRFReco\recon\CartesianRecon_Mar2022_Human.jl:0`

tknopp commented 1 year ago

Have you updated SparsityOperators?

nbwuzhe commented 1 year ago

Thanks @tknopp after upgrading all depending packages both Cartesian and spiral data were passed as in the old version.

alexjaffray commented 1 year ago

@nbwuzhe and I have updated SparsityOperators and things seem to be working now. Thanks @tknopp!

Total allocations are a bit higher than in v0.5.0 (albeit not problematically so) and things run and reconstructions look as they did before. Thanks for the quick fixes! We are still trying to figure out how to efficiently put our data in a test to include into MRIReco to catch these things more quickly, have been using v0.5.0 for a while now and decided just recently to try to update to the newest developments.

Is there a plan for formalizing the changes on master into a new release anytime soon?

tknopp commented 1 year ago

releases are triggered.