From 4eba737782f7d3f52cfe5c8535fcb2098b21d8b2 Mon Sep 17 00:00:00 2001 From: nHackel Date: Fri, 26 Jul 2024 11:08:42 +0200 Subject: [PATCH 01/10] Add RowNorm and Composite Weighting --- src/Weighting.jl | 24 +++++++++++++++++++++--- 1 file changed, 21 insertions(+), 3 deletions(-) diff --git a/src/Weighting.jl b/src/Weighting.jl index d34f569..124ddf7 100644 --- a/src/Weighting.jl +++ b/src/Weighting.jl @@ -5,25 +5,43 @@ abstract type AbstractWeightingParameters <: AbstractMPIRecoParameters end export NoWeightingParameters struct NoWeightingParameters <: AbstractWeightingParameters end -process(::Type{<:AbstractMPIRecoAlgorithm}, params::NoWeightingParameters, data) = nothing +process(::Type{<:AbstractMPIRecoAlgorithm}, params::NoWeightingParameters, args...) = nothing export ChannelWeightingParameters Base.@kwdef struct ChannelWeightingParameters <: AbstractWeightingParameters channelWeights::Vector{Float64} = [1.0, 1.0, 1.0] end -process(::Type{<:AbstractMPIRecoAlgorithm}, params::ChannelWeightingParameters, freqs::Vector{CartesianIndex{2}}) = map(x-> params.channelWeights[x[2]], freqs) +process(::Type{<:AbstractMPIRecoAlgorithm}, params::ChannelWeightingParameters, freqs::Vector{CartesianIndex{2}}, args...) = map(x-> params.channelWeights[x[2]], freqs) export WhiteningWeightingParameters Base.@kwdef struct WhiteningWeightingParameters <: AbstractWeightingParameters whiteningMeas::MPIFile tfCorrection::Bool = false end -function process(::Type{<:AbstractMPIRecoAlgorithm}, params::WhiteningWeightingParameters, freqs::Vector{CartesianIndex{2}}) +function process(::Type{<:AbstractMPIRecoAlgorithm}, params::WhiteningWeightingParameters, freqs::Vector{CartesianIndex{2}}, args...) u_bg = getMeasurementsFD(params.whiteningMeas, false, frequencies=freqs, frames=measBGFrameIdx(params.whiteningMeas), bgCorrection = false, tfCorrection=false) bg_std = std(u_bg, dims=3) weights = minimum(abs.(vec(bg_std))) ./ abs.(vec(bg_std)) return weights end + +export RowNormWeightingParameters +Base.@kwdef struct RowNormWeightingParameters <: AbstractWeightingParameters + # NOP +end +function process(::Type{<:AbstractMPIRecoAlgorithm}, params::RowNormWeightingParameters, freqs, op, args...) + weights = map(r -> 1/sqrt(rownorm²(op, r)), 1:size(op, 1)) + return weights +end + +export CompositeWeightingParameters +Base.@kwdef struct CompositeWeightingParameters{WS} <: AbstractWeightingParameters where WS <: AbstractWeightingParameters + weightingParameters::Vector{WS} +end +function process(algoT::Type{<:AbstractMPIRecoAlgorithm}, params::CompositeWeightingParameters, args...) + weights = map(p -> process(algoT, p, args...), params.weightingParameters) + return reduce(.*, weights) +end #= baremodule WeightingType None = 0 From ef3dca09df076ba9378fdeb1cb53489652af139e Mon Sep 17 00:00:00 2001 From: nHackel Date: Fri, 26 Jul 2024 11:09:01 +0200 Subject: [PATCH 02/10] Adapt SinglePatch Reco to weighting interface changes --- .../SinglePatchAlgorithms/SinglePatchAlgorithm.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/Algorithms/SinglePatchAlgorithms/SinglePatchAlgorithm.jl b/src/Algorithms/SinglePatchAlgorithms/SinglePatchAlgorithm.jl index 2caccd2..54a903a 100644 --- a/src/Algorithms/SinglePatchAlgorithms/SinglePatchAlgorithm.jl +++ b/src/Algorithms/SinglePatchAlgorithms/SinglePatchAlgorithm.jl @@ -62,7 +62,13 @@ function process(algo::SinglePatchReconstructionAlgorithm, params::SinglePatchRe return gridresult(result, algo.grid, algo.sf) end -process(algo::SinglePatchReconstructionAlgorithm, params::Union{W, ProcessResultCache{W}}, u) where {W <: Union{ChannelWeightingParameters, WhiteningWeightingParameters}} = map(real(eltype(algo.S)), process(typeof(algo), params, algo.freqs)) +function process(algo::SinglePatchReconstructionAlgorithm, params::Union{W, ProcessResultCache{W}}, u) where {W <: AbstractWeightingParameters} + result = process(typeof(algo), params, algo.freqs, algo.sf) + if !isnothing(result) + result = map(real(eltype(algo.S)), result) + end + return nothing +end function getLinearOperator(algo::SinglePatchReconstructionAlgorithm, params::SinglePatchReconstructionParameter{<:DenseSystemMatixLoadingParameter, S}) where {S} return nothing From cf4f6040f80d072fd9ddf673bb047eab83674dc7 Mon Sep 17 00:00:00 2001 From: nHackel Date: Fri, 26 Jul 2024 11:09:43 +0200 Subject: [PATCH 03/10] Add SinglePatch rownorm and composite weighting tests --- test/Reconstruction.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/test/Reconstruction.jl b/test/Reconstruction.jl index b0b1a84..0ec7295 100644 --- a/test/Reconstruction.jl +++ b/test/Reconstruction.jl @@ -128,4 +128,12 @@ c8 = reconstruct("SinglePatch", b; params...) exportImage(joinpath(imgdir, "Reconstruction8.png"), Array(c8[1,:,:,1,1])) @test compareImg("Reconstruction8.png") + + c9a = reconstruct("SinglePatch", b; params..., reg = [L2Regularization(0.0)], weightingParams = NoWeightingParameters()) + c9b = reconstruct("SinglePatch", b; params..., reg = [L2Regularization(0.0)], weightingParams = RowNormWeightingParameters()) + @test isapprox(arraydata(c9a), arraydata(c9b)) + + c10a = reconstruct("SinglePatch", b; params..., reg = [L2Regularization(0.0)], weightingParams = WhiteningWeightingParameters(whiteningMeas = bSF)) + c10b = reconstruct("SinglePatch", b; params..., reg = [L2Regularization(0.0)], weightingParams = CompositeWeightingParameters([WhiteningWeightingParameters(whiteningMeas = bSF), RowNormWeightingParameters()])) + @test isapprox(arraydata(c9a), arraydata(c9b)) end From add9276586ac51d8ac6663b2862add8180ca8d58 Mon Sep 17 00:00:00 2001 From: nHackel Date: Fri, 26 Jul 2024 14:27:31 +0200 Subject: [PATCH 04/10] Add Weighting to MultiPatchRecos and fix bug in SinglePatch weighting-process --- .../MultiPatchAlgorithms/MultiPatchAlgorithm.jl | 16 +++++++++++++--- .../SinglePatchAlgorithm.jl | 2 +- 2 files changed, 14 insertions(+), 4 deletions(-) diff --git a/src/Algorithms/MultiPatchAlgorithms/MultiPatchAlgorithm.jl b/src/Algorithms/MultiPatchAlgorithms/MultiPatchAlgorithm.jl index f5d5c48..aaaaa8a 100644 --- a/src/Algorithms/MultiPatchAlgorithms/MultiPatchAlgorithm.jl +++ b/src/Algorithms/MultiPatchAlgorithms/MultiPatchAlgorithm.jl @@ -1,5 +1,5 @@ export MultiPatchReconstructionAlgorithm, MultiPatchReconstructionParameter -Base.@kwdef struct MultiPatchReconstructionParameter{matT <: AbstractArray,F<:AbstractFrequencyFilterParameter,O<:AbstractMultiPatchOperatorParameter, S<:AbstractSolverParameters, FF<:AbstractFocusFieldPositions, FFSF<:AbstractFocusFieldPositions, R <: AbstractRegularization} <: AbstractMultiPatchReconstructionParameters +Base.@kwdef struct MultiPatchReconstructionParameter{matT <: AbstractArray,F<:AbstractFrequencyFilterParameter,O<:AbstractMultiPatchOperatorParameter, S<:AbstractSolverParameters, FF<:AbstractFocusFieldPositions, FFSF<:AbstractFocusFieldPositions, R <: AbstractRegularization, W<:AbstractWeightingParameters} <: AbstractMultiPatchReconstructionParameters arrayType::Type{matT} = Array # File sf::MultiMPIFile @@ -9,7 +9,7 @@ Base.@kwdef struct MultiPatchReconstructionParameter{matT <: AbstractArray,F<:Ab ffPosSF::FFSF = DefaultFocusFieldPositions() solverParams::S reg::Vector{R} = AbstractRegularization[] - # weightingType::WeightingType = WeightingType.None + weightingParams::Union{W, ProcessResultCache{W}} = NoWeightingParameters() end Base.@kwdef mutable struct MultiPatchReconstructionAlgorithm{P, matT <: AbstractArray} <: AbstractMultiPatchReconstructionAlgorithm where {P<:AbstractMultiPatchAlgorithmParameters} @@ -121,9 +121,19 @@ function process(::Type{<:MultiPatchReconstructionAlgorithm}, params::ExternalPr end function process(algo::MultiPatchReconstructionAlgorithm, params::MultiPatchReconstructionParameter, u::AbstractArray) - solver = LeastSquaresParameters(S = algo.ffOp, reg = params.reg, solverParams = params.solverParams) + weights = adapt(algo.arrayType, process(algo, params.weightingParams, u)) + + solver = LeastSquaresParameters(S = algo.ffOp, reg = params.reg, solverParams = params.solverParams, weights = weights) result = process(algo, solver, u) return gridresult(result, algo.ffOp.grid, algo.sf) +end + +function process(algo::MultiPatchReconstructionAlgorithm, params::Union{W, ProcessResultCache{W}}, u) where {W <: AbstractWeightingParameters} + result = process(typeof(algo), params, algo.freqs, algo.ffOp) + if !isnothing(result) + result = map(real(eltype(algo.ffOp)), result) + end + return result end \ No newline at end of file diff --git a/src/Algorithms/SinglePatchAlgorithms/SinglePatchAlgorithm.jl b/src/Algorithms/SinglePatchAlgorithms/SinglePatchAlgorithm.jl index 54a903a..181e202 100644 --- a/src/Algorithms/SinglePatchAlgorithms/SinglePatchAlgorithm.jl +++ b/src/Algorithms/SinglePatchAlgorithms/SinglePatchAlgorithm.jl @@ -67,7 +67,7 @@ function process(algo::SinglePatchReconstructionAlgorithm, params::Union{W, Proc if !isnothing(result) result = map(real(eltype(algo.S)), result) end - return nothing + return result end function getLinearOperator(algo::SinglePatchReconstructionAlgorithm, params::SinglePatchReconstructionParameter{<:DenseSystemMatixLoadingParameter, S}) where {S} From f55f540201ca5ac69f8121af5c8c467f95f19078 Mon Sep 17 00:00:00 2001 From: nHackel Date: Fri, 26 Jul 2024 14:28:28 +0200 Subject: [PATCH 05/10] Add hash for MutliPatchOperators --- ext/MPIRecoKernelAbstractionsExt/MultiPatch.jl | 15 +++++++++++++++ src/MultiPatch.jl | 7 +++++++ 2 files changed, 22 insertions(+) diff --git a/ext/MPIRecoKernelAbstractionsExt/MultiPatch.jl b/ext/MPIRecoKernelAbstractionsExt/MultiPatch.jl index f99af94..a1d8af4 100644 --- a/ext/MPIRecoKernelAbstractionsExt/MultiPatch.jl +++ b/ext/MPIRecoKernelAbstractionsExt/MultiPatch.jl @@ -140,4 +140,19 @@ function RegularizedLeastSquares.kaczmarz_update!(op::DenseMultiPatchOperator{T, kernel(x, op.S, row, beta, op.xcc, op.xss, op.sign, Int32(div(op.M, op.nPatches)), op.RowToPatch, op.patchToSMIdx; ndrange = size(op.xss, 1)) synchronize(backend) return x +end + +function Base.hash(op::DenseMultiPatchOperator{T, V}, h::UInt64) where {T, V <: AbstractGPUArray{T}} + @warn "Hashing of GPU DenseMultiPatchOperator is inefficient" + h = hash(typeof(op), h) + h = @allowscalar hash(op.S, h) + h = hash(op.grid, h) + h = hash(op.N, h) + h = hash(op.M, h) + h = @allowscalar hash(op.RowToPatch, h) + h = @allowscalar hash(op.xcc, h) + h = @allowscalar hash(op.xss, h) + h = @allowscalar hash(op.sign, h) + h = hash(op.nPatches, h) + h = @allowscalar hash(op.patchToSMIdx, h) end \ No newline at end of file diff --git a/src/MultiPatch.jl b/src/MultiPatch.jl index 6f40c2b..5a85256 100644 --- a/src/MultiPatch.jl +++ b/src/MultiPatch.jl @@ -111,6 +111,13 @@ Adapt.adapt_structure(::Type{Array}, op::AbstractMultiPatchOperator) = op LinearOperators.storage_type(op::MultiPatchOperator) = LinearOperators.storage_type(first(op.S)) LinearOperators.storage_type(op::DenseMultiPatchOperator) = typeof(similar(op.S, 0)) +function Base.hash(op::AbstractMultiPatchOperator, h::UInt64) + h = hash(typeof(op), h) + for field in fieldnames(typeof(op)) + h = hash(getfield(op, field), h) + end + return h +end function Base.convert(::Type{DenseMultiPatchOperator}, op::MultiPatchOperator) S = stack(op.S) From e299ed833d8f62bfde780ad06d99c6c7971a4910 Mon Sep 17 00:00:00 2001 From: nHackel Date: Mon, 29 Jul 2024 12:03:39 +0200 Subject: [PATCH 06/10] Differentiate between systemmatrix and measurement based weighting This allows for easier caching. In particular we can cache/compute weighting before the SM is moved to the GPU. Since GPU arrays don't implement a hash function this allows for performant hash computations that dont have to rely on @allowscalar --- .../MultiPatchAlgorithm.jl | 35 +++++++++++-------- .../MultiPatchPeriodicMotion.jl | 25 +++++++------ .../SinglePatchAlgorithm.jl | 34 ++++++++++-------- .../SinglePatchAlgorithms.jl | 2 +- .../SinglePatchBGEstimationAlgorithm.jl | 11 +++--- ...glePatchTemporalRegularizationAlgorithm.jl | 8 ++--- src/Weighting.jl | 14 ++++++++ 7 files changed, 79 insertions(+), 50 deletions(-) diff --git a/src/Algorithms/MultiPatchAlgorithms/MultiPatchAlgorithm.jl b/src/Algorithms/MultiPatchAlgorithms/MultiPatchAlgorithm.jl index aaaaa8a..203c442 100644 --- a/src/Algorithms/MultiPatchAlgorithms/MultiPatchAlgorithm.jl +++ b/src/Algorithms/MultiPatchAlgorithms/MultiPatchAlgorithm.jl @@ -1,6 +1,6 @@ export MultiPatchReconstructionAlgorithm, MultiPatchReconstructionParameter -Base.@kwdef struct MultiPatchReconstructionParameter{matT <: AbstractArray,F<:AbstractFrequencyFilterParameter,O<:AbstractMultiPatchOperatorParameter, S<:AbstractSolverParameters, FF<:AbstractFocusFieldPositions, FFSF<:AbstractFocusFieldPositions, R <: AbstractRegularization, W<:AbstractWeightingParameters} <: AbstractMultiPatchReconstructionParameters - arrayType::Type{matT} = Array +Base.@kwdef struct MultiPatchReconstructionParameter{arrT <: AbstractArray,F<:AbstractFrequencyFilterParameter,O<:AbstractMultiPatchOperatorParameter, S<:AbstractSolverParameters, FF<:AbstractFocusFieldPositions, FFSF<:AbstractFocusFieldPositions, R <: AbstractRegularization, W<:AbstractWeightingParameters} <: AbstractMultiPatchReconstructionParameters + arrayType::Type{arrT} = Array # File sf::MultiMPIFile freqFilter::F @@ -12,12 +12,13 @@ Base.@kwdef struct MultiPatchReconstructionParameter{matT <: AbstractArray,F<:Ab weightingParams::Union{W, ProcessResultCache{W}} = NoWeightingParameters() end -Base.@kwdef mutable struct MultiPatchReconstructionAlgorithm{P, matT <: AbstractArray} <: AbstractMultiPatchReconstructionAlgorithm where {P<:AbstractMultiPatchAlgorithmParameters} +Base.@kwdef mutable struct MultiPatchReconstructionAlgorithm{P, arrT <: AbstractArray, vecT <: AbstractArray} <: AbstractMultiPatchReconstructionAlgorithm where {P<:AbstractMultiPatchAlgorithmParameters} params::P # Could also do reconstruction progress meter here opParams::Union{AbstractMultiPatchOperatorParameter, ProcessResultCache{<:AbstractMultiPatchOperatorParameter},Nothing} = nothing sf::MultiMPIFile - arrayType::Type{matT} + weights::Union{Nothing, vecT} = nothing + arrayType::Type{arrT} ffOp::Union{Nothing, AbstractMultiPatchOperator} ffPos::Union{Nothing,AbstractArray} ffPosSF::Union{Nothing,AbstractArray} @@ -40,7 +41,7 @@ function MultiPatchReconstructionAlgorithm(params::MultiPatchParameters{<:Abstra ffPosSF = [vec(ffPos(SF))[l] for l=1:L, SF in reco.sf] end - return MultiPatchReconstructionAlgorithm(params, reco.opParams, reco.sf, reco.arrayType, nothing, ffPos_, ffPosSF, freqs, Channel{Any}(Inf)) + return MultiPatchReconstructionAlgorithm{typeof(params), reco.arrayType, typeof(reco.arrayType{Float32}(undef, 0))}(params, reco.opParams, reco.sf, nothing, reco.arrayType, nothing, ffPos_, ffPosSF, freqs, Channel{Any}(Inf)) end recoAlgorithmTypes(::Type{MultiPatchReconstruction}) = SystemMatrixBasedAlgorithm() AbstractImageReconstruction.parameter(algo::MultiPatchReconstructionAlgorithm) = algo.origParam @@ -50,7 +51,7 @@ AbstractImageReconstruction.take!(algo::MultiPatchReconstructionAlgorithm) = Bas function AbstractImageReconstruction.put!(algo::MultiPatchReconstructionAlgorithm, data::MPIFile) #consistenceCheck(algo.sf, data) - algo.ffOp = process(algo, algo.opParams, data, algo.freqs) + algo.ffOp, algo.weights = process(algo, algo.opParams, data, algo.freqs, algo.params.reco.weightingParams) result = process(algo, algo.params, data, algo.freqs) @@ -65,7 +66,7 @@ function AbstractImageReconstruction.put!(algo::MultiPatchReconstructionAlgorith Base.put!(algo.output, result) end -function process(algo::MultiPatchReconstructionAlgorithm, params::Union{OP, ProcessResultCache{OP}}, f::MPIFile, frequencies::Vector{CartesianIndex{2}}) where OP <: AbstractMultiPatchOperatorParameter +function process(algo::MultiPatchReconstructionAlgorithm, params::Union{OP, ProcessResultCache{OP}}, f::MPIFile, frequencies::Vector{CartesianIndex{2}}, weightingParams) where OP <: AbstractMultiPatchOperatorParameter ffPos_ = ffPos(f) periodsSortedbyFFPos = unflattenOffsetFieldShift(ffPos_) idxFirstPeriod = getindex.(periodsSortedbyFFPos,1) @@ -77,7 +78,11 @@ function process(algo::MultiPatchReconstructionAlgorithm, params::Union{OP, Proc end result = process(typeof(algo), params, algo.sf, frequencies, gradient, ffPos_, algo.ffPosSF) - return adapt(algo.arrayType, result) + # Kinda of hacky. MultiPatch parameters don't map nicely to the SinglePatch inspired pre, reco, post structure + # Have to create weights before ffop is (potentially) moved to GPU, as GPU arrays don't have efficient hash implementations + # Which makes this process expensive to cache + weights = process(typeof(algo), weightingParams, frequencies, result, nothing, algo.arrayType) + return adapt(algo.arrayType, result), weights end function process(algo::MultiPatchReconstructionAlgorithm, params::Union{A, ProcessResultCache{<:A}}, f::MPIFile, args...) where A <: AbstractMPIPreProcessingParameters @@ -121,7 +126,7 @@ function process(::Type{<:MultiPatchReconstructionAlgorithm}, params::ExternalPr end function process(algo::MultiPatchReconstructionAlgorithm, params::MultiPatchReconstructionParameter, u::AbstractArray) - weights = adapt(algo.arrayType, process(algo, params.weightingParams, u)) + weights = process(algo, params.weightingParams, u, WeightingType(params.weightingParams)) solver = LeastSquaresParameters(S = algo.ffOp, reg = params.reg, solverParams = params.solverParams, weights = weights) @@ -130,10 +135,10 @@ function process(algo::MultiPatchReconstructionAlgorithm, params::MultiPatchReco return gridresult(result, algo.ffOp.grid, algo.sf) end -function process(algo::MultiPatchReconstructionAlgorithm, params::Union{W, ProcessResultCache{W}}, u) where {W <: AbstractWeightingParameters} - result = process(typeof(algo), params, algo.freqs, algo.ffOp) - if !isnothing(result) - result = map(real(eltype(algo.ffOp)), result) - end - return result +function process(algo::MultiPatchReconstructionAlgorithm, params::Union{W, ProcessResultCache{W}}, u, ::MeasurementBasedWeighting) where W<:AbstractWeightingParameters + return process(typeof(algo), params, algo.freqs, algo.ffOp, u, algo.arrayType) +end + +function process(algo::MultiPatchReconstructionAlgorithm, params::Union{W, ProcessResultCache{W}}, u, ::SystemMatrixBasedWeighting) where W<:AbstractWeightingParameters + return algo.weights end \ No newline at end of file diff --git a/src/Algorithms/MultiPatchAlgorithms/MultiPatchPeriodicMotion.jl b/src/Algorithms/MultiPatchAlgorithms/MultiPatchPeriodicMotion.jl index 2415f7a..af55a9d 100644 --- a/src/Algorithms/MultiPatchAlgorithms/MultiPatchPeriodicMotion.jl +++ b/src/Algorithms/MultiPatchAlgorithms/MultiPatchPeriodicMotion.jl @@ -1,5 +1,5 @@ export PeriodicMotionPreProcessing, PeriodicMotionReconstructionParameter -Base.@kwdef struct PeriodicMotionPreProcessing{BG<:AbstractMPIBackgroundCorrectionParameters} <: AbstractMPIPreProcessingParameters{BG} +Base.@kwdef struct PeriodicMotionPreProcessing{BG<:AbstractMPIBackgroundCorrectionParameters, W <: AbstractWeightingParameters} <: AbstractMPIPreProcessingParameters{BG} # Periodic Motion frames::Union{Nothing, UnitRange{Int64}, Vector{Int64}} = nothing alpha::Float64 = 3.0 @@ -10,21 +10,21 @@ Base.@kwdef struct PeriodicMotionPreProcessing{BG<:AbstractMPIBackgroundCorrecti sf::MultiMPIFile tfCorrection::Bool = false bgParams::BG = NoBackgroundCorrectionParameters() - # weightingType::WeightingType = WeightingType.None + weightingParams::Union{W, ProcessResultCache{W}} = NoWeightingParameters() end -Base.@kwdef struct PeriodicMotionReconstructionParameter{F<:AbstractFrequencyFilterParameter, S<:AbstractSolverParameters} <: AbstractMultiPatchReconstructionParameters +Base.@kwdef struct PeriodicMotionReconstructionParameter{F<:AbstractFrequencyFilterParameter, S<:AbstractSolverParameters, R <: AbstractRegularization, arrT <: AbstractArray} <: AbstractMultiPatchReconstructionParameters sf::MultiMPIFile freqFilter::F solverParams::S - λ::Float32 - # weightingType::WeightingType = WeightingType.None + reg::Vector{R} = AbstractRegularization[] + arrayType::Type{arrT} = Array end function MultiPatchReconstructionAlgorithm(params::MultiPatchParameters{<:PeriodicMotionPreProcessing,<:PeriodicMotionReconstructionParameter,<:AbstractMPIPostProcessingParameters}) reco = params.reco freqs = process(MultiPatchReconstructionAlgorithm, reco.freqFilter, reco.sf) - return MultiPatchReconstructionAlgorithm(params, nothing, reco.sf, Array, nothing, nothing, nothing, freqs, Channel{Any}(Inf)) + return MultiPatchReconstructionAlgorithm{typeof(params), reco.arrayType, typeof(reco.arrayType{Float32}(undef, 0))}(params, nothing, reco.sf, nothing, reco.arrayType, nothing, nothing, nothing, freqs, Channel{Any}(Inf)) end function AbstractImageReconstruction.put!(algo::MultiPatchReconstructionAlgorithm{MultiPatchParameters{PT, R, T}}, data::MPIFile) where {R, T, PT <: PeriodicMotionPreProcessing} @@ -43,8 +43,9 @@ end function process(algo::MultiPatchReconstructionAlgorithm, params::Union{OP, ProcessResultCache{OP}}, f::MPIFile, frequencies::Union{Vector{CartesianIndex{2}}, Nothing} = nothing) where OP <: PeriodicMotionPreProcessing - uReco, ffOp = process(typeof(algo), params, f, algo.sf, frequencies) + uReco, ffOp, weights = process(typeof(algo), params, f, algo.sf, frequencies) algo.ffOp = adapt(algo.arrayType, ffOp) + algo.weights = adapt(algo.arrayType, weights) return adapt(algo.arrayType, uReco) end @@ -67,26 +68,28 @@ function process(algoT::Type{<:MultiPatchReconstructionAlgorithm}, params::Perio resortedInd[i,:] = unflattenOffsetFieldShift(ffPos_)[i][1:p] end + # Can't adapt data here because it might be used in background correction ffOp = MultiPatchOperator(sf, frequencies, #indFFPos=resortedInd[:,1], unused keyword FFPos=ffPos_[:,resortedInd[:,1]], mapping=mapping, FFPosSF=ffPos_[:,resortedInd[:,1]], bgCorrection = false, tfCorrection = params.tfCorrection) - return uReco, ffOp + weights = process(algoT, params.weightingParams, frequencies, ffOp, nothing, Array) + return uReco, ffOp, weights end function process(algoT::Type{<:MultiPatchReconstructionAlgorithm}, params::PeriodicMotionPreProcessing{SimpleExternalBackgroundCorrectionParameters}, f::MPIFile, sf::MPIFile, frequencies::Union{Vector{CartesianIndex{2}}, Nothing} = nothing) # Foreground fgParams = fromKwargs(PeriodicMotionPreProcessing; toKwargs(params)..., bgParams = NoBackgroundCorrectionParameters()) - result, ffOp = process(algoT, fgParams, f, sf, frequencies) + result, ffOp, weights = process(algoT, fgParams, f, sf, frequencies) # Background bgParams = fromKwargs(ExternalPreProcessedBackgroundCorrectionParameters; toKwargs(params)..., bgParams = params.bgParams, spectralLeakageCorrection=true) - return process(algoT, bgParams, result, frequencies), ffOp + return process(algoT, bgParams, result, frequencies), ffOp, weights end function process(algo::MultiPatchReconstructionAlgorithm, params::PeriodicMotionReconstructionParameter, u::Array) - solver = LeastSquaresParameters(S = algo.ffOp, reg = [L2Regularization(params.λ)], solverParams = params.solverParams) + solver = LeastSquaresParameters(S = algo.ffOp, reg = params.reg, solverParams = params.solverParams, weights = algo.weights) result = process(algo, solver, u) diff --git a/src/Algorithms/SinglePatchAlgorithms/SinglePatchAlgorithm.jl b/src/Algorithms/SinglePatchAlgorithms/SinglePatchAlgorithm.jl index 181e202..039484b 100644 --- a/src/Algorithms/SinglePatchAlgorithms/SinglePatchAlgorithm.jl +++ b/src/Algorithms/SinglePatchAlgorithms/SinglePatchAlgorithm.jl @@ -1,21 +1,22 @@ Base.@kwdef struct SinglePatchReconstructionParameter{L<:AbstractSystemMatrixLoadingParameter, SL<:AbstractLinearSolver, - matT <: AbstractArray, SP<:AbstractSolverParameters{SL}, R<:AbstractRegularization, W<:AbstractWeightingParameters} <: AbstractSinglePatchReconstructionParameters + arrT <: AbstractArray, SP<:AbstractSolverParameters{SL}, R<:AbstractRegularization, W<:AbstractWeightingParameters} <: AbstractSinglePatchReconstructionParameters # File sf::MPIFile sfLoad::Union{L, ProcessResultCache{L}} - arrayType::Type{matT} = Array + arrayType::Type{arrT} = Array # Solver solverParams::SP reg::Vector{R} = AbstractRegularization[] weightingParams::Union{W, ProcessResultCache{W}} = NoWeightingParameters() end -Base.@kwdef mutable struct SinglePatchReconstructionAlgorithm{P, matT <: AbstractArray} <: AbstractSinglePatchReconstructionAlgorithm where {P<:AbstractSinglePatchAlgorithmParameters} +Base.@kwdef mutable struct SinglePatchReconstructionAlgorithm{P, SM, arrT <: AbstractArray, vecT <: arrT} <: AbstractSinglePatchReconstructionAlgorithm where {P<:AbstractSinglePatchAlgorithmParameters} params::P # Could also do reconstruction progress meter here sf::Union{MPIFile, Vector{MPIFile}} - S::AbstractArray - arrayType::Type{matT} + S::SM + weights::Union{Nothing, vecT} = nothing + arrayType::Type{arrT} grid::RegularGridPositions freqs::Vector{CartesianIndex{2}} output::Channel{Any} @@ -26,16 +27,20 @@ function SinglePatchReconstruction(params::SinglePatchParameters{<:AbstractMPIPr end function SinglePatchReconstructionAlgorithm(params::SinglePatchParameters{<:AbstractMPIPreProcessingParameters, R, PT}) where {R<:AbstractSinglePatchReconstructionParameters, PT <:AbstractMPIPostProcessingParameters} freqs, S, grid, arrayType = prepareSystemMatrix(params.reco) - return SinglePatchReconstructionAlgorithm(params, params.reco.sf, S, arrayType, grid, freqs, Channel{Any}(Inf)) + weights = prepareWeights(params.reco, freqs, S) + return SinglePatchReconstructionAlgorithm{typeof(params), typeof(S), arrayType, typeof(arrayType{real(eltype(S))}(undef, 0))}(params, params.reco.sf, S, weights, arrayType, grid, freqs, Channel{Any}(Inf)) end recoAlgorithmTypes(::Type{SinglePatchReconstruction}) = SystemMatrixBasedAlgorithm() AbstractImageReconstruction.parameter(algo::SinglePatchReconstructionAlgorithm) = algo.params function prepareSystemMatrix(reco::SinglePatchReconstructionParameter{L,S}) where {L<:AbstractSystemMatrixLoadingParameter, S<:AbstractLinearSolver} - freqs, sf, grid = process(AbstractMPIRecoAlgorithm, reco.sfLoad, reco.sf, S, reco.arrayType) + freqs, sf, grid = process(AbstractMPIRecoAlgorithm, reco.sfLoad, reco.sf, nothing, reco.arrayType) return freqs, sf, grid, reco.arrayType end +function prepareWeights(reco::SinglePatchReconstructionParameter{L,S,arrT,SP,R,W}, freqs, sf) where {L, S, arrT, SP, R, W<:AbstractWeightingParameters} + return process(AbstractMPIRecoAlgorithm, reco.weightingParams, freqs, sf, nothing, reco.arrayType) +end AbstractImageReconstruction.take!(algo::SinglePatchReconstructionAlgorithm) = Base.take!(algo.output) @@ -51,7 +56,7 @@ end function process(algo::SinglePatchReconstructionAlgorithm, params::SinglePatchReconstructionParameter, u) - weights = adapt(algo.arrayType, process(algo, params.weightingParams, u)) + weights = process(algo, params.weightingParams, u, WeightingType(params.weightingParams)) B = getLinearOperator(algo, params) @@ -62,12 +67,13 @@ function process(algo::SinglePatchReconstructionAlgorithm, params::SinglePatchRe return gridresult(result, algo.grid, algo.sf) end -function process(algo::SinglePatchReconstructionAlgorithm, params::Union{W, ProcessResultCache{W}}, u) where {W <: AbstractWeightingParameters} - result = process(typeof(algo), params, algo.freqs, algo.sf) - if !isnothing(result) - result = map(real(eltype(algo.S)), result) - end - return result +function process(algo::SinglePatchReconstructionAlgorithm, params::Union{W, ProcessResultCache{W}}, u, ::MeasurementBasedWeighting) where W<:AbstractWeightingParameters + return process(typeof(algo), params, algo.freqs, algo.S, u, algo.arrayType) +end + + +function process(algo::SinglePatchReconstructionAlgorithm, params::Union{W, ProcessResultCache{W}}, u, ::SystemMatrixBasedWeighting) where W<:AbstractWeightingParameters + return algo.weights end function getLinearOperator(algo::SinglePatchReconstructionAlgorithm, params::SinglePatchReconstructionParameter{<:DenseSystemMatixLoadingParameter, S}) where {S} diff --git a/src/Algorithms/SinglePatchAlgorithms/SinglePatchAlgorithms.jl b/src/Algorithms/SinglePatchAlgorithms/SinglePatchAlgorithms.jl index f1f0c3a..df86405 100644 --- a/src/Algorithms/SinglePatchAlgorithms/SinglePatchAlgorithms.jl +++ b/src/Algorithms/SinglePatchAlgorithms/SinglePatchAlgorithms.jl @@ -19,7 +19,7 @@ function process(algo::T, params::SinglePatchParameters, data::MPIFile, frequenc end function AbstractImageReconstruction.put!(algo::AbstractSinglePatchReconstructionAlgorithm, data) - consistenceCheck(algo.sf, data) + #consistenceCheck(algo.sf, data) result = process(algo, algo.params, data, algo.freqs) diff --git a/src/Algorithms/SinglePatchAlgorithms/SinglePatchBGEstimationAlgorithm.jl b/src/Algorithms/SinglePatchAlgorithms/SinglePatchBGEstimationAlgorithm.jl index a63f09a..2305beb 100644 --- a/src/Algorithms/SinglePatchAlgorithms/SinglePatchBGEstimationAlgorithm.jl +++ b/src/Algorithms/SinglePatchAlgorithms/SinglePatchBGEstimationAlgorithm.jl @@ -1,10 +1,10 @@ export SinglePatchBGEstimationAlgorithm, SinglePatchBGEstimationReconstructionParameter Base.@kwdef struct SinglePatchBGEstimationReconstructionParameter{L<:DenseSystemMatixLoadingParameter, - SP<:AbstractSolverParameters, matT <: AbstractArray} <: AbstractSinglePatchReconstructionParameters + SP<:AbstractSolverParameters, arrT <: AbstractArray} <: AbstractSinglePatchReconstructionParameters # File sf::MPIFile sfLoad::Union{L, ProcessResultCache{L}} - arrayType::Type{matT} = Array + arrayType::Type{arrT} = Array # Solver solverParams::SP λ::Float32 @@ -13,12 +13,13 @@ Base.@kwdef struct SinglePatchBGEstimationReconstructionParameter{L<:DenseSystem bgDict::BGDictParameter end -Base.@kwdef mutable struct SinglePatchBGEstimationAlgorithm{P, matT <: AbstractArray} <: AbstractSinglePatchReconstructionAlgorithm where {P<:AbstractSinglePatchAlgorithmParameters} +Base.@kwdef mutable struct SinglePatchBGEstimationAlgorithm{P, arrT <: AbstractArray} <: AbstractSinglePatchReconstructionAlgorithm where {P<:AbstractSinglePatchAlgorithmParameters} params::P # Could also do reconstruction progress meter here sf::Union{MPIFile,Vector{MPIFile}} S::AbstractArray - arrayType::Type{matT} + #weights::Union{Nothing, vecT} = nothing + arrayType::Type{arrT} bgDict::AbstractArray grid::RegularGridPositions freqs::Vector{CartesianIndex{2}} @@ -30,7 +31,7 @@ function SinglePatchReconstruction(params::SinglePatchParameters{<:AbstractMPIPr end function SinglePatchBGEstimationAlgorithm(params::SinglePatchParameters{<:AbstractMPIPreProcessingParameters,R,PT}) where {R<:SinglePatchBGEstimationReconstructionParameter,PT<:AbstractMPIPostProcessingParameters} freqs, S, grid, arrayType = prepareSystemMatrix(params.reco) - return SinglePatchBGEstimationAlgorithm(params, params.reco.sf, S, arrayType, process(SinglePatchBGEstimationAlgorithm, freqs, params.reco.bgDict), grid, freqs, Channel{Any}(Inf)) + return SinglePatchBGEstimationAlgorithm(params, params.reco.sf, S, nothing, arrayType, process(SinglePatchBGEstimationAlgorithm, freqs, params.reco.bgDict), grid, freqs, Channel{Any}(Inf)) end recoAlgorithmTypes(::Type{SinglePatchBGEstimationAlgorithm}) = SystemMatrixBasedAlgorithm() AbstractImageReconstruction.parameter(algo::SinglePatchBGEstimationAlgorithm) = algo.origParam diff --git a/src/Algorithms/SinglePatchAlgorithms/SinglePatchTemporalRegularizationAlgorithm.jl b/src/Algorithms/SinglePatchAlgorithms/SinglePatchTemporalRegularizationAlgorithm.jl index fad1b81..b7e0d07 100644 --- a/src/Algorithms/SinglePatchAlgorithms/SinglePatchTemporalRegularizationAlgorithm.jl +++ b/src/Algorithms/SinglePatchAlgorithms/SinglePatchTemporalRegularizationAlgorithm.jl @@ -1,10 +1,10 @@ export SinglePatchTemporalRegularizationAlgorithm, SinglePatchTemporalRegularizationReconstructionParameter Base.@kwdef struct SinglePatchTemporalRegularizationReconstructionParameter{L<:DenseSystemMatixLoadingParameter, - SP<:AbstractSolverParameters, matT <: AbstractArray} <: AbstractSinglePatchReconstructionParameters + SP<:AbstractSolverParameters, arrT <: AbstractArray} <: AbstractSinglePatchReconstructionParameters # File sf::MPIFile sfLoad::Union{L, ProcessResultCache{L}} - arrayType::Type{matT} = Array + arrayType::Type{arrT} = Array # Solver solverParams::SP λ::Float32 @@ -15,11 +15,11 @@ Base.@kwdef struct SinglePatchTemporalRegularizationReconstructionParameter{L<:D bgDict::BGDictParameter end -Base.@kwdef mutable struct SinglePatchTemporalRegularizationAlgorithm{P, matT <: AbstractArray} <: AbstractSinglePatchReconstructionAlgorithm where {P<:AbstractSinglePatchAlgorithmParameters} +Base.@kwdef mutable struct SinglePatchTemporalRegularizationAlgorithm{P, arrT <: AbstractArray} <: AbstractSinglePatchReconstructionAlgorithm where {P<:AbstractSinglePatchAlgorithmParameters} params::P sf::Union{MPIFile,Vector{MPIFile}} S::AbstractArray - arrayType::Type{matT} + arrayType::Type{arrT} bgDict::AbstractArray idxFG::Union{Nothing, UnitRange{Int64}, Vector{Int64}} = nothing idxFG::Union{Nothing, UnitRange{Int64}, Vector{Int64}} = nothing diff --git a/src/Weighting.jl b/src/Weighting.jl index 124ddf7..a977488 100644 --- a/src/Weighting.jl +++ b/src/Weighting.jl @@ -2,6 +2,20 @@ export getWeights, WeightingType, setFreqToZero export AbstractWeightingParameters abstract type AbstractWeightingParameters <: AbstractMPIRecoParameters end +function process(type::Type{<:AbstractMPIRecoAlgorithm}, params::AbstractWeightingParameters, freqs, op = nothing, u = nothing, arrayType = Array) + result = process(type, params, freqs, op, u) + if !isnothing(result) + result = map(real(eltype(algo.S)), result) + end + return adapt(arrayType, result) +end + +abstract type WeightingType end +struct MeasurementBasedWeighting <: WeightingType end +struct SystemMatrixBasedWeighting <: WeightingType end + +WeightingType(::AbstractWeightingParameters) = SystemMatrixBasedWeighting() +WeightingType(cache::ProcessResultCache) = WeightingType(cache.param) export NoWeightingParameters struct NoWeightingParameters <: AbstractWeightingParameters end From 42122bffe1416621cbb63e67ef55c23ed62de2b2 Mon Sep 17 00:00:00 2001 From: nHackel Date: Mon, 29 Jul 2024 16:51:28 +0200 Subject: [PATCH 07/10] Implement more efficient weighted normalization for MultiPatch op --- .../MPIRecoKernelAbstractionsExt.jl | 2 +- .../MultiPatch.jl | 59 +++++++++++++++++++ 2 files changed, 60 insertions(+), 1 deletion(-) diff --git a/ext/MPIRecoKernelAbstractionsExt/MPIRecoKernelAbstractionsExt.jl b/ext/MPIRecoKernelAbstractionsExt/MPIRecoKernelAbstractionsExt.jl index e3d43e2..4eed007 100644 --- a/ext/MPIRecoKernelAbstractionsExt/MPIRecoKernelAbstractionsExt.jl +++ b/ext/MPIRecoKernelAbstractionsExt/MPIRecoKernelAbstractionsExt.jl @@ -1,6 +1,6 @@ module MPIRecoKernelAbstractionsExt -using MPIReco, MPIReco.Adapt, MPIReco.LinearAlgebra, MPIReco.RegularizedLeastSquares +using MPIReco, MPIReco.Adapt, MPIReco.LinearAlgebra, MPIReco.RegularizedLeastSquares, MPIReco.LinearOperatorCollection using KernelAbstractions, GPUArrays using KernelAbstractions.Extras: @unroll using Atomix diff --git a/ext/MPIRecoKernelAbstractionsExt/MultiPatch.jl b/ext/MPIRecoKernelAbstractionsExt/MultiPatch.jl index a1d8af4..05157aa 100644 --- a/ext/MPIRecoKernelAbstractionsExt/MultiPatch.jl +++ b/ext/MPIRecoKernelAbstractionsExt/MultiPatch.jl @@ -142,6 +142,65 @@ function RegularizedLeastSquares.kaczmarz_update!(op::DenseMultiPatchOperator{T, return x end +function RegularizedLeastSquares.normalize(::SystemMatrixBasedNormalization, op::OP, x) where {T, V <: AbstractGPUArray{T}, OP <: DenseMultiPatchOperator{T, V}} + weights = one(real(eltype(op))) + energy = normalize_dense_op(op, weights) + return norm(energy)^2/size(op, 2) +end + +function RegularizedLeastSquares.normalize(::SystemMatrixBasedNormalization, prod::ProdOp{T, <:WeightingOp, OP}, x) where {T, V <: AbstractGPUArray{T}, OP <: DenseMultiPatchOperator{T, V}} + op = prod.B + weights = prod.A.weights + energy = normalize_dense_op(op, weights) + return norm(energy)^2/size(prod, 2) +end + +function normalize_dense_op(op::DenseMultiPatchOperator{T, V}, weights) where {T, V <: AbstractGPUArray{T}} + backend = get_backend(op.S) + kernel = normalize_kernel!(backend, 256) + energy = KernelAbstractions.zeros(backend, real(eltype(op)), size(op, 1)) + kernel(energy, weights, op.S, op.xss, op.sign, Int32(div(op.M, op.nPatches)), op.RowToPatch, op.patchToSMIdx; ndrange = (256, size(op, 1))) + synchronize(backend) + return energy +end + +# The normalization kernels are structured the same as the mul!-kernel. The multiplication with x is replaced by abs2 for the rownorm² +@kernel cpu = false inbounds = true function normalize_kernel!(energy, weights, @Const(S), @Const(xss), @Const(signs), @Const(M), @Const(RowToPatch), @Const(patchToSMIdx)) + # Each group/block handles a single row of the operator + operator_row = @index(Group, Linear) # k + patch = RowToPatch[operator_row] # p + patch_row = mod1(operator_row, M) # j + smIdx = patchToSMIdx[patch] + sign = eltype(energy)(signs[patch_row, smIdx]) + grid_stride = prod(@groupsize()) + N = Int32(size(xss, 1)) + + localIdx = @index(Local, Linear) + shared = @localmem eltype(energy) grid_stride + shared[localIdx] = zero(eltype(energy)) + + @unroll for i = localIdx:grid_stride:N + shared[localIdx] = shared[localIdx] + abs2(sign * S[patch_row, xss[i, patch], smIdx]) + end + @synchronize + + @private s = div(min(grid_stride, N), Int32(2)) + while s > Int32(0) + if localIdx <= s + shared[localIdx] = shared[localIdx] + shared[localIdx + s] + end + s >>= 1 + @synchronize + end + + if localIdx == 1 + energy[operator_row] = sqrt(get_kernel_weights(weights, operator_row)^2 * shared[localIdx]) + end +end + +@inline get_kernel_weights(weights::AbstractArray, operator_row) = weights[operator_row] +@inline get_kernel_weights(weights::Number, operator_row) = weights + function Base.hash(op::DenseMultiPatchOperator{T, V}, h::UInt64) where {T, V <: AbstractGPUArray{T}} @warn "Hashing of GPU DenseMultiPatchOperator is inefficient" h = hash(typeof(op), h) From 07163d32486319c43ea87e7236802d608c4f6ef3 Mon Sep 17 00:00:00 2001 From: nHackel Date: Tue, 30 Jul 2024 15:24:11 +0200 Subject: [PATCH 08/10] Improve adapt-GPU array performance for MultiPatch Operator --- ext/MPIRecoKernelAbstractionsExt/MultiPatch.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/ext/MPIRecoKernelAbstractionsExt/MultiPatch.jl b/ext/MPIRecoKernelAbstractionsExt/MultiPatch.jl index 05157aa..793d2a0 100644 --- a/ext/MPIRecoKernelAbstractionsExt/MultiPatch.jl +++ b/ext/MPIRecoKernelAbstractionsExt/MultiPatch.jl @@ -5,10 +5,10 @@ function Adapt.adapt_structure(::Type{arrT}, op::MultiPatchOperator) where {arrT # Ideally we create a DenseMultiPatchOperator on the GPU if validSMs && validXCC && validXSS - S = adapt(arrT, stack(op.S)) + S = stack(adapt.(arrT, op.S)) # We want to use Int32 for better GPU performance - xcc = Int32.(adapt(arrT, stack(op.xcc))) - xss = Int32.(adapt(arrT, stack(op.xss))) + xcc = Int32.(stack(adapt.(arrT, op.xcc))) + xss = Int32.(stack(adapt.(arrT, op.xss))) sign = Int32.(adapt(arrT, op.sign)) RowToPatch = Int32.(adapt(arrT, op.RowToPatch)) patchToSMIdx = Int32.(adapt(arrT, op.patchToSMIdx)) From 1aea871f8c7c9d5fa6673c560df2a997a617a24c Mon Sep 17 00:00:00 2001 From: nHackel Date: Tue, 30 Jul 2024 15:44:26 +0200 Subject: [PATCH 09/10] Add rownorm weighted reconstruction tests for MultiPatch operator --- test/ReconstructionGPU.jl | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/test/ReconstructionGPU.jl b/test/ReconstructionGPU.jl index 8044605..98f8f62 100644 --- a/test/ReconstructionGPU.jl +++ b/test/ReconstructionGPU.jl @@ -103,6 +103,10 @@ for arrayType in arrayTypes c2 = reconstruct("MultiPatch", b; params..., arrayType = arrayType) @test isapprox(arraydata(c1), arraydata(c2), rtol = 0.02) + c3 = reconstruct("MultiPatch", b; params..., weightingParams = RowNormWeightingParameters()) + c4 = reconstruct("MultiPatch", b; params..., weightingParams = RowNormWeightingParameters(), arrayType = arrayType) + @test isapprox(arraydata(c3), arraydata(c4), rtol = 0.02) + # Test Kaczmarz since it uses specific functions and not just mul! bSF = MultiMPIFile([joinpath(datadir, "calibrations", "12.mdf")]) dirs = ["1.mdf", "2.mdf", "3.mdf", "4.mdf"] @@ -118,9 +122,9 @@ for arrayType in arrayTypes params[:reg] = [L2Regularization(0.0f0)] params[:tfCorrection] = false params[:solver] = Kaczmarz - c3 = reconstruct("MultiPatch", b; params...) - c4 = reconstruct("MultiPatch", b; params..., arrayType = arrayType) - @test isapprox(arraydata(c3), arraydata(c4)) + c5 = reconstruct("MultiPatch", b; params...) + c6 = reconstruct("MultiPatch", b; params..., arrayType = arrayType) + @test isapprox(arraydata(c5), arraydata(c6)) end end \ No newline at end of file From 5b813e9857ac209b9d4a93b6d6b65411bdb4bf11 Mon Sep 17 00:00:00 2001 From: nHackel Date: Tue, 30 Jul 2024 17:45:11 +0200 Subject: [PATCH 10/10] Update prepareSystemMatrix function to include linear solver parameter for single patch algo --- src/Algorithms/SinglePatchAlgorithms/SinglePatchAlgorithm.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Algorithms/SinglePatchAlgorithms/SinglePatchAlgorithm.jl b/src/Algorithms/SinglePatchAlgorithms/SinglePatchAlgorithm.jl index 039484b..2dc74c1 100644 --- a/src/Algorithms/SinglePatchAlgorithms/SinglePatchAlgorithm.jl +++ b/src/Algorithms/SinglePatchAlgorithms/SinglePatchAlgorithm.jl @@ -34,7 +34,7 @@ recoAlgorithmTypes(::Type{SinglePatchReconstruction}) = SystemMatrixBasedAlgorit AbstractImageReconstruction.parameter(algo::SinglePatchReconstructionAlgorithm) = algo.params function prepareSystemMatrix(reco::SinglePatchReconstructionParameter{L,S}) where {L<:AbstractSystemMatrixLoadingParameter, S<:AbstractLinearSolver} - freqs, sf, grid = process(AbstractMPIRecoAlgorithm, reco.sfLoad, reco.sf, nothing, reco.arrayType) + freqs, sf, grid = process(AbstractMPIRecoAlgorithm, reco.sfLoad, reco.sf, S, reco.arrayType) return freqs, sf, grid, reco.arrayType end