update readme/docstrings/docs;

restrict GCV to m <= n and update tests;
test for `SVDValsWorkspace`

README.md

@@ -33,7 +33,7 @@ $ julia --project=@decaes -e 'import Pkg; Pkg.add("DECAES"); Pkg.build("DECAES")
 This will do two things:
 
 1. Add DECAES.jl to a named Julia project environment separate from your global environment
-2. Build the `decaes` executable at `~/.julia/bin` for running DECAES from the command line
+2. Build the `decaes` launcher script at `~/.julia/bin` for running DECAES from the command line
 
 DECAES can then be run from the command line via `decaes <COMMAND LINE ARGS>`, provided `~/.julia/bin` is added to your `PATH`.
 Run `decaes --help` for available arguments.

api/decaes.m

@@ -49,11 +49,11 @@
 %     * Lastly, we indicate that the regularization parameters should be
 %       saved using the --SaveRegParam flag
 %
-%       decaes image.nii.gz --output results --mask image_mask.mat --T2map --T2part --TE 7e-3 --nT2 60 --T2Range 7e-3 2.0 --SPWin 7e-3 25e-3 --MPWin 25e-3 200.0e-3 --Reg lcurve --SaveRegParam
+%       decaes image.nii.gz --output results --mask image_mask.mat --T2map --T2part --TE 7e-3 --nT2 40 --T2Range 7e-3 2.0 --SPWin 7e-3 25e-3 --MPWin 25e-3 200.0e-3 --Reg lcurve --SaveRegParam
 %
 %   Run the same command using function syntax:
 %
-%       decaes('image.nii.gz', '--output', 'results', '--mask', 'image_mask.mat', '--T2map', '--T2part', '--TE', 7e-3, '--nT2', 60, '--T2Range', [7e-3, 2.0], '--SPWin', [7e-3, 25e-3], '--MPWin', [25e-3, 200.0e-3], '--Reg', 'lcurve', '--SaveRegParam')
+%       decaes('image.nii.gz', '--output', 'results', '--mask', 'image_mask.mat', '--T2map', '--T2part', '--TE', 7e-3, '--nT2', 40, '--T2Range', [7e-3, 2.0], '--SPWin', [7e-3, 25e-3], '--MPWin', [25e-3, 200.0e-3], '--Reg', 'lcurve', '--SaveRegParam')
 %
 %   Create a settings file called 'settings.txt' containing the settings
 %   from the above example (note: only one value or flag per line):
@@ -68,7 +68,7 @@
 %       --TE
 %       7e-3
 %       --nT2
-%       60
+%       40
 %       --T2Range
 %       7e-3
 %       2.0

docs/src/cli.md

@@ -6,7 +6,7 @@ DECAES provides a command line interface (CLI) for calling the main analysis fun
 
 There are two equivalent ways to use the [command line interface (CLI)](https://jondeuce.github.io/DECAES.jl/dev/cli), assuming DECAES is already [installed](@ref installation):
 
-**1. (Recommended) `decaes` launcher:** Use the executable `~/.julia/bin/decaes` which comes with DECAES:
+**1. (Recommended) `decaes` launcher:** Use the script `~/.julia/bin/decaes` which comes with DECAES:
 
 ```bash
 $ decaes <COMMAND LINE ARGS>

docs/src/index.md

@@ -40,7 +40,7 @@ Using Julia v1.9 or later you can install DECAES as follows:
 $ julia --project=@decaes -e 'import Pkg; Pkg.add("DECAES"); Pkg.build("DECAES")'
 ```
 
-This will add DECAES.jl to a named Julia project environment separate from your global environment, and build the `decaes` executable at `~/.julia/bin` for running DECAES from the command line.
+This will add DECAES.jl to a named Julia project environment separate from your global environment, and build the `decaes` launcher script at `~/.julia/bin` for running DECAES from the command line.
 
 ## [Updating DECAES](@id updating)
 

docs/src/ref.md

@@ -29,10 +29,10 @@ main
 ```@docs
 lsqnonneg
 lsqnonneg_tikh
-lsqnonneg_mdp
-lsqnonneg_chi2
-lsqnonneg_gcv
 lsqnonneg_lcurve
+lsqnonneg_gcv
+lsqnonneg_chi2
+lsqnonneg_mdp
 lcurve_corner
 ```
 

src/T2mapSEcorr.jl

@@ -138,9 +138,11 @@ See also:
 
   - [`T2partSEcorr`](@ref)
   - [`lsqnonneg`](@ref)
-  - [`lsqnonneg_chi2`](@ref)
-  - [`lsqnonneg_gcv`](@ref)
+  - [`lsqnonneg_tikh`](@ref)
   - [`lsqnonneg_lcurve`](@ref)
+  - [`lsqnonneg_gcv`](@ref)
+  - [`lsqnonneg_chi2`](@ref)
+  - [`lsqnonneg_mdp`](@ref)
   - [`EPGdecaycurve`](@ref)
 """
 T2mapSEcorr(image::Array{T, 4}; kwargs...) where {T} = T2mapSEcorr(image, T2mapOptions(image; kwargs...))
@@ -424,6 +426,7 @@ end
 # T2-distribution fitting
 # =========================================================
 abstract type RegularizationMethod end
+struct NoRegularization <: RegularizationMethod end
 struct LCurve <: RegularizationMethod end
 struct GCV <: RegularizationMethod end
 struct ChiSquared{T} <: RegularizationMethod
@@ -433,24 +436,23 @@ end
 struct MDP{T} <: RegularizationMethod
     NoiseLevel::T
 end
-struct NoRegularization <: RegularizationMethod end
 
 function regularization_method(o::T2mapOptions)
     reg =
+        o.Reg == "none"   ? NoRegularization() : # Fit T2 distribution using unregularized NNLS
         o.Reg == "lcurve" ? LCurve() : # Fit T2 distribution using L-curve-based regularized NNLS
         o.Reg == "gcv"    ? GCV() : # Fit T2 distribution using GCV-based regularized NNLS
         o.Reg == "chi2"   ? ChiSquared(o.Chi2Factor, o.legacy) : # Fit T2 distribution using chi2-based regularized NNLS
         o.Reg == "mdp"    ? MDP(o.NoiseLevel) : # Fit T2 distribution using Morizov discrepancy principle-based regularized NNLS
-        o.Reg == "none"   ? NoRegularization() : # Fit T2 distribution using unregularized NNLS
         error("Unrecognized regularization method: $(o.Reg)")
     return reg
 end
 
+nnls_workspace(::NoRegularization, decay_basis::AbstractMatrix{T}, decay_data::AbstractVector{T}) where {T} = lsqnonneg_work(decay_basis, decay_data)
 nnls_workspace(::LCurve, decay_basis::AbstractMatrix{T}, decay_data::AbstractVector{T}) where {T} = lsqnonneg_lcurve_work(decay_basis, decay_data)
 nnls_workspace(::GCV, decay_basis::AbstractMatrix{T}, decay_data::AbstractVector{T}) where {T} = lsqnonneg_gcv_work(decay_basis, decay_data)
 nnls_workspace(::ChiSquared, decay_basis::AbstractMatrix{T}, decay_data::AbstractVector{T}) where {T} = lsqnonneg_chi2_work(decay_basis, decay_data)
 nnls_workspace(::MDP, decay_basis::AbstractMatrix{T}, decay_data::AbstractVector{T}) where {T} = lsqnonneg_mdp_work(decay_basis, decay_data)
-nnls_workspace(::NoRegularization, decay_basis::AbstractMatrix{T}, decay_data::AbstractVector{T}) where {T} = lsqnonneg_work(decay_basis, decay_data)
 
 struct T2DistWorkspace{Reg, T, W}
     reg::Reg
@@ -470,6 +472,12 @@ function T2DistWorkspace(reg::RegularizationMethod, decay_basis::Matrix{T}, deca
     return T2DistWorkspace(reg, nnls_work, decay_basis, decay_data, decay_scale, μ, χ²fact)
 end
 
+function T2_distribution!(t2work::T2DistWorkspace{NoRegularization, T}) where {T}
+    (; nnls_work, μ, χ²fact) = t2work
+    μ[], χ²fact[] = zero(T), one(T)
+    return lsqnonneg!(nnls_work)
+end
+
 function T2_distribution!(t2work::T2DistWorkspace{LCurve, T}) where {T}
     (; nnls_work, μ, χ²fact) = t2work
     x, μ[], χ²fact[] = lsqnonneg_lcurve!(nnls_work)
@@ -496,12 +504,6 @@ function T2_distribution!(t2work::T2DistWorkspace{MDP{T}, T}) where {T}
     return x
 end
 
-function T2_distribution!(t2work::T2DistWorkspace{NoRegularization, T}) where {T}
-    (; nnls_work, μ, χ²fact) = t2work
-    μ[], χ²fact[] = zero(T), one(T)
-    return lsqnonneg!(nnls_work)
-end
-
 T2_distribution(t2work::T2DistWorkspace) = solution(t2work.nnls_work)
 
 # =========================================================

src/lsqnonneg.jl

@@ -1017,6 +1017,7 @@ struct NNLSGCVRegProblem{T, TA <: AbstractMatrix{T}, Tb <: AbstractVector{T}, W0
 end
 function NNLSGCVRegProblem(A::AbstractMatrix{T}, b::AbstractVector{T}) where {T}
     m, n = size(A)
+    @assert m >= n "Generalized Cross-Validation (GCV) regularization requires at least as many rows as columns, but got m = $m < n = $n"
     svd_work = SVDValsWorkspace(A) # workspace for computing singular values
     nnls_prob = NNLSProblem(A, b)
     nnls_prob_smooth_cache = NNLSTikhonovRegProblemCache(A, b)

src/types.jl

@@ -63,9 +63,9 @@ See also:
     nRefAnglesMin::Int = !legacy ? min(5, nRefAngles) : nRefAngles
     @assert 2 <= nRefAnglesMin <= nRefAngles "Minimum number of angles to check during flip angle optimization must be in the range [2, nRefAngles], but nRefAngles = $nRefAngles and nRefAnglesMin = $nRefAnglesMin."
 
-    "Regularization routine to use. One of \"none\", \"chi2\", \"gcv\", or \"lcurve\", representing no regularization, `Chi2Factor`-based Tikhonov regularization, Generalized Cross-Validation method, or L-Curve based regularization, respectively."
+    "Regularization routine to use. One of \"none\", \"lcurve\", \"gcv\", \"chi2\", or \"mdp\", representing no regularization, the L-Curve method, the Generalized Cross-Validation method, `Chi2Factor`-based Tikhonov regularization, or the Morozov discrepancy principle, respectively."
     Reg::String
-    @assert Reg ∈ ("none", "lcurve", "chi2", "gcv", "mdp")
+    @assert Reg ∈ ("none", "lcurve", "gcv", "chi2", "mdp")
 
     "Constraint on ``\\chi^2`` used for regularization when `Reg == \"chi2\"`."
     Chi2Factor::Union{T, Nothing} = nothing

src/utils.jl

@@ -53,12 +53,12 @@ end
 ####
 
 function with_singlethreaded_blas(f)
-    nblasthreads = LinearAlgebra.BLAS.get_num_threads()
+    nblasthreads = BLAS.get_num_threads()
     try
-        LinearAlgebra.BLAS.set_num_threads(1) # Prevent BLAS from stealing julia threads
+        BLAS.set_num_threads(1) # Prevent BLAS from stealing julia threads
         f()
     finally
-        LinearAlgebra.BLAS.set_num_threads(nblasthreads) # Reset BLAS threads
+        BLAS.set_num_threads(nblasthreads) # Reset BLAS threads
     end
 end
 
@@ -95,7 +95,7 @@ end
 
 function LinearAlgebra.svdvals!(work::SVDValsWorkspace, A::AbstractMatrix)
     copyto!(work.A, A)
-    return LinearAlgebra.svdvals!(work)
+    return svdvals!(work)
 end
 
 # See: https://github.com/JuliaLang/julia/blob/64de065a183ac70bb049f7f9e30d790f8845dd2b/stdlib/LinearAlgebra/src/lapack.jl#L1590
@@ -545,7 +545,7 @@ function mock_t2map_opts(::Type{T} = Float64; kwargs...) where {T}
         nT2 = 40,
         Reg = "lcurve",
         Chi2Factor = 1.02,
-        NoiseLevel = 1e-3, # SNR = 60
+        NoiseLevel = exp10(-T(get(kwargs, :SNR, 60.0)) / 20),
         t2map_kwargs...,
     )
 end