Type Stability of Normal Operator

[JakobAsslaender]

I am wondering if we should also re-write the way we construct the normal operator depending on whether we want to use a Toeplitz Operator or not. Unfortunately, the function

MRIReco.jl/src/Operators/NFFTOp.jl

Lines 105 to 111 in 355a5c8

	function SparsityOperators.normalOperator(S::NFFTOp, W::UniformScaling=I)
	if S.toeplitz
	return Toeplitz_NormalOp(S,W)
	else
	return NormalOp(S,W)
	end
	end

is not type stable:

julia> @code_warntype SparsityOperators.normalOperator(F_nfft)
Variables
  #self#::Core.Const(SparsityOperators.normalOperator)
  S::NFFTOp{ComplexF64}

Body::Union{NormalOp{NFFTOp{ComplexF64}, UniformScaling{Bool}}, MRIReco.Toeplitz_NormalOp{_A, _B, UniformScaling{Bool}} where {_A, _B}}
1 ─ %1 = (#self#)(S, MRIReco.I)::Union{NormalOp{NFFTOp{ComplexF64}, UniformScaling{Bool}}, MRIReco.Toeplitz_NormalOp{_A, _B, UniformScaling{Bool}} where {_A, _B}}
└──      return %1

Any ideas on how we could do this?

[tknopp]

Type instabilities are not always evil and actually it is a feature of Julia that you can mix dynamic and static code. In this particular case I don't think that it will be an issue. The question is what happens with the return type of normalOperator. If you pass it to a solver the code within the solver will be fast because there was a function barrier.

Basically a reconstruction chain on a high-level will always be dynamic, i.e. the operator types being constructed will depend on the input parameters (e.g. non equidistant sampling points -> NFFT, otherwise -> FFT). One the operators a constructed all low-level code should be fast.

The danger is if one has type-unstable structs and passes them from high to low level.