Give system matrix as argument to least squares parameters

docs/src/literate/tutorials/lowlevel.jl

@@ -29,7 +29,7 @@ S, grid = getSF(bSF, freqs, sparseTrafo, Kaczmarz)
 typeof.([S, grid])
 
 # Now we can configur a low-level reconstruction:
-cLow = reconstruct("LowLevel", u; S = S, iterations = params[:iterations], reg = params[:reg])
+cLow = reconstruct("LowLevel", S, u; iterations = params[:iterations], reg = params[:reg])
 # Note that the low-level reconstruction returns a matrix without any metadata unlike the other reconstructions.
 # The second dimension of the result matrix are the frames. To compare and plot our data we have to reshape it:
 sliceLow = sliceLow = reshape(cLow[:, 1], Tuple(grid.shape))

src/AlgorithmInterface.jl

@@ -122,6 +122,16 @@ function reconstruct(name::AbstractString, data::Union{MPIFile, AbstractArray},
   plan = loadRecoPlan(name, cache, modules; kwargs...)
   return reconstruct(build(plan), data)
 end
+function reconstruct(name::AbstractString, S, data::Union{MPIFile, AbstractArray}, cache::Bool = false, modules = getRecoPlanModules(); kwargs...)
+  plan = loadRecoPlan(name, cache, modules; kwargs...)
+  algo = build(plan)
+  result = nothing
+  lock(algo) do
+    put!(algo, S, data)
+    result = take!(algo)
+  end
+  return result
+end
 # Load plan with RecoCache consideration
 function loadRecoPlan(name::AbstractString, cache::Bool, modules; kwargs...)
   planfile = planpath(name)

src/Algorithms/HandsFreeAlgorithms/HandsFreeLeastSquares.jl

@@ -11,18 +11,17 @@ Base.@kwdef struct HandsFreeSolverParameters <: AbstractSolverParameters{Kaczmar
   flattenIters::Bool=false
 end
 
-function process(t::Type{<:AbstractMPIRecoAlgorithm}, params::LeastSquaresParameters{Kaczmarz, O, SF, R, SL, W}, u::AbstractArray, snr::AbstractVector) where {O, SF, R, SL <: HandsFreeSolverParameters, W}
+function process(t::Type{<:AbstractMPIRecoAlgorithm}, params::LeastSquaresParameters{Kaczmarz, O, R, SL, W}, S, u::AbstractArray, snr::AbstractVector) where {O, SF, R, SL <: HandsFreeSolverParameters, W}
 
   solverParams = params.solverParams
-  N = size(params.S, 2)
-  M = div(length(params.S), N)
+  N = size(S, 2)
+  M = div(length(S), N)
   L = size(u)[end]
   u = reshape(u, M, L)
   c = zeros(N, L)
 
   reg, _ = prepareRegularization([L2Regularization(real(eltype(u))(solverParams.startλ))], params)
 
-  S = params.S
   if !isnothing(params.weights)
     S = ProdOp(WeightingOp(params.weights), S)
     for l = 1:L

src/Algorithms/HandsFreeAlgorithms/HandsFreeSinglePatch.jl

@@ -25,11 +25,11 @@ function process(algo::SinglePatchReconstructionAlgorithm, params::SinglePatchHa
 
   B = getLinearOperator(algo, params)
 
-  solver = LeastSquaresParameters(op = B, S = algo.S, reg = L2Regularization[], solverParams = params.solverParams, weights = weights)
+  solver = LeastSquaresParameters(op = B, reg = L2Regularization[], solverParams = params.solverParams, weights = weights)
 
   snr = real(eltype(algo.S)).(vec(MPIFiles.getCalibSNR(algo.sf)[algo.freqs, 1]))
 
-  result = process(algo, solver, u, snr)
+  result = process(algo, solver, algo.S, u, snr)
 
   return gridresult(result, algo.grid, algo.sf)
 end

src/Algorithms/LowLevelAlgorithm.jl

@@ -12,6 +12,6 @@ Base.isready(algo::LowLevelReconstructionAlgorithm) = isready(algo.output)
 Base.wait(algo::LowLevelReconstructionAlgorithm) = wait(algo.output)
 AbstractImageReconstruction.take!(algo::LowLevelReconstructionAlgorithm) = Base.take!(algo.output)
 
-function AbstractImageReconstruction.put!(algo::LowLevelReconstructionAlgorithm, u::AbstractArray)
-  put!(algo.output, process(algo, algo.params, u))
+function AbstractImageReconstruction.put!(algo::LowLevelReconstructionAlgorithm, S, u::AbstractArray)
+  put!(algo.output, process(algo, algo.params, S, u))
 end

src/Algorithms/MultiPatchAlgorithms/MultiPatchAlgorithm.jl

@@ -135,9 +135,9 @@ end
 function process(algo::MultiPatchReconstructionAlgorithm, params::MultiPatchReconstructionParameter, u::AbstractArray)
   weights = process(algo, params.weightingParams, u, WeightingType(params.weightingParams))
 
-  solver = LeastSquaresParameters(S = algo.ffOp, reg = params.reg, solverParams = params.solverParams, weights = weights)
+  solver = LeastSquaresParameters(reg = params.reg, solverParams = params.solverParams, weights = weights)
 
-  result = process(algo, solver, u)
+  result = process(algo, solver, algo.ffOp, u)
 
   return gridresult(result, algo.ffOp.grid, algo.sf)
 end

src/Algorithms/MultiPatchAlgorithms/MultiPatchPeriodicMotion.jl

@@ -91,9 +91,9 @@ function process(algoT::Type{<:MultiPatchReconstructionAlgorithm},
 end
 
 function process(algo::MultiPatchReconstructionAlgorithm, params::PeriodicMotionReconstructionParameter, u::Array)
-  solver = LeastSquaresParameters(S = algo.ffOp, reg = params.reg, solverParams = params.solverParams, weights = algo.weights)
+  solver = LeastSquaresParameters(reg = params.reg, solverParams = params.solverParams, weights = algo.weights)
 
-  result = process(algo, solver, u)
+  result = process(algo, solver, algo.ffOp, u)
 
   return gridresult(result, algo.ffOp.grid, algo.sf)
 end

src/Algorithms/SinglePatchAlgorithms/SinglePatchAlgorithm.jl

@@ -64,9 +64,9 @@ function process(algo::SinglePatchReconstructionAlgorithm, params::SinglePatchRe
 
   B = getLinearOperator(algo, params)
 
-  solver = LeastSquaresParameters(op = B, S = algo.S, reg = params.reg, solverParams = params.solverParams, weights = weights)
+  solver = LeastSquaresParameters(op = B, reg = params.reg, solverParams = params.solverParams, weights = weights)
 
-  result = process(algo, solver, u)
+  result = process(algo, solver, algo.S, u)
 
   return gridresult(result, algo.grid, algo.sf)
 end

src/Algorithms/SinglePatchAlgorithms/SinglePatchBGEstimationAlgorithm.jl

@@ -75,9 +75,9 @@ function process(algo::SinglePatchBGEstimationAlgorithm, params::SinglePatchBGEs
   solverParams = fromKwargs(typeof(params.solverParams); toKwargs(params.solverParams, overwrite = Dict{Symbol, Any}(:normalizeReg => NoNormalization()))...)
   solverParams = ConstraintMaskedSolverParameters(;constraintMask = constraintMask, params = params.solverParams)
 
-  solver = LeastSquaresParameters(solver = Kaczmarz, S = G, reg = [reg], solverParams = solverParams)
+  solver = LeastSquaresParameters(solver = Kaczmarz, reg = [reg], solverParams = solverParams)
 
-  temp = process(algo, solver, u)
+  temp = process(algo, solver, G, u)
 
   result = zeros(eltype(temp), N, size(temp, 2))
 

src/Algorithms/SinglePatchAlgorithms/SinglePatchTemporalRegularizationAlgorithm.jl

@@ -85,9 +85,9 @@ function process(algo::SinglePatchTemporalRegularizationAlgorithm, params::Singl
   solverParams = fromKwargs(typeof(params.solverParams); toKwargs(params.solverParams, overwrite = Dict{Symbol, Any}(:normalizeReg => NoNormalization()))...)
   solverParams = ConstraintMaskedSolverParameters(;constraintMask = constraintMask, params = params.solverParams)
 
-  solver = LeastSquaresParameters(solver = Kaczmarz, S = op, reg = [reg], solverParams = solverParams)
+  solver = LeastSquaresParameters(solver = Kaczmarz, reg = [reg], solverParams = solverParams)
 
-  temp = process(algo, solver, u)
+  temp = process(algo, solver, op, u)
   temp = real.( reshape(temp[1:(NSub*J),:],NSub,J) ./ sqrt(λ) ) 
 
   cInterp = similar(temp, size(c,1), op.L)

src/LeastSquares.jl

@@ -3,9 +3,8 @@ export LeastSquaresParameters
 abstract type AbstractSolverParameters{AbstractLinearSolver} <: AbstractMPIRecoParameters end
 
 export LeastSquaresParameters
-Base.@kwdef struct LeastSquaresParameters{L<:AbstractLinearSolver, O, M, R<:AbstractRegularization, P<:AbstractSolverParameters{L}, W} <: AbstractMPIRecoParameters
+Base.@kwdef struct LeastSquaresParameters{L<:AbstractLinearSolver, O, R<:AbstractRegularization, P<:AbstractSolverParameters{L}, W} <: AbstractMPIRecoParameters
   op::O = nothing
-  S::M
   reg::Vector{R} 
   solverParams::P
   weights::W = nothing
@@ -55,18 +54,17 @@ Base.propertynames(params::RecoPlan{ElaborateSolverParameters}) = union([:solver
 
 getSolverKwargs(::Type{SL}) where SL <: AbstractLinearSolver = intersect(union(Base.kwarg_decl.(methods(SL))...), fieldnames(ElaborateSolverParameters))
 
-function process(t::Type{<:AbstractMPIRecoAlgorithm}, params::LeastSquaresParameters{SL}, u::AbstractArray) where SL
+function process(t::Type{<:AbstractMPIRecoAlgorithm}, params::LeastSquaresParameters{SL}, S, u::AbstractArray) where SL
 
-  N = size(params.S, 2)
-  M = div(length(params.S), N)
+  N = size(S, 2)
+  M = div(length(S), N)
   L = size(u)[end]
   u = reshape(u, M, L)
   c = zeros(N, L)
 
   reg, args = prepareRegularization(params.reg, params)
   args[:reg] = reg
 
-  S = params.S
   if !isnothing(params.weights)
     S = ProdOp(WeightingOp(params.weights), S)
     u = params.weights.*u

test/LowLevel.jl

@@ -20,7 +20,7 @@
 
     S = high.S
     u = process(high, high.params.pre, b, high.freqs)
-    cLow = reconstruct("LowLevel", u; S = S, iterations = params[:iterations], reg = params[:reg], solver = params[:solver])
+    cLow = reconstruct("LowLevel", S, u;, iterations = params[:iterations], reg = params[:reg], solver = params[:solver])
 
     cLow = reshape(cLow, size(cHigh))
     @test isapprox(cHigh.data.data, cLow)
@@ -33,7 +33,7 @@
     S = high.S
     u = process(high, high.params.pre, b, high.freqs)
     weights = high.weights
-    cLow = reconstruct("LowLevel", u; S = S, iterations = params[:iterations], reg = params[:reg], solver = params[:solver], weights = weights)
+    cLow = reconstruct("LowLevel", S, u; iterations = params[:iterations], reg = params[:reg], solver = params[:solver], weights = weights)
 
     cLow = reshape(cLow, size(cHigh))
     @test isapprox(cHigh.data.data, cLow)
@@ -63,7 +63,7 @@
     S = high.S
     u = process(high, high.params.pre, b, high.freqs)
     op = MPIReco.getLinearOperator(high, high.params.reco)
-    cLow = reconstruct("LowLevel", u; S = S, op = op, iterations = params[:iterations], reg = params[:reg], solver = params[:solver])
+    cLow = reconstruct("LowLevel", S, u; op = op, iterations = params[:iterations], reg = params[:reg], solver = params[:solver])
 
     cLow = reshape(cLow, size(cHigh))
     @test isapprox(cHigh.data.data, cLow)
@@ -77,7 +77,7 @@
     u = process(high, high.params.pre, b, high.freqs)
     op = MPIReco.getLinearOperator(high, high.params.reco)
     weights = high.weights
-    cLow = reconstruct("LowLevel", u; S = S, op = op, iterations = params[:iterations], reg = params[:reg], solver = params[:solver], weights = weights)
+    cLow = reconstruct("LowLevel", S, u; op = op, iterations = params[:iterations], reg = params[:reg], solver = params[:solver], weights = weights)
 
     cLow = reshape(cLow, size(cHigh))
     @test isapprox(cHigh.data.data, cLow)
@@ -105,7 +105,7 @@
 
     S = copy(high.ffOp)
     u = process(high, high.params.pre, b, high.freqs)
-    cLow = reconstruct("LowLevel", u; S = S, iterations = params[:iterations], reg = params[:reg], solver = params[:solver])
+    cLow = reconstruct("LowLevel", S, u; iterations = params[:iterations], reg = params[:reg], solver = params[:solver])
 
     cLow = reshape(cLow, size(cHigh))
     @test isapprox(cHigh.data.data, cLow)
@@ -118,7 +118,7 @@
     S = copy(high.ffOp)
     u = process(high, high.params.pre, b, high.freqs)
     weights = high.weights
-    cLow = reconstruct("LowLevel", u; S = S, iterations = params[:iterations], reg = params[:reg], solver = params[:solver], weights = weights)
+    cLow = reconstruct("LowLevel", S, u; iterations = params[:iterations], reg = params[:reg], solver = params[:solver], weights = weights)
 
     cLow = reshape(cLow, size(cHigh))
     @test isapprox(cHigh.data.data, cLow)