expfit fixed, commit version 0.9.0

ext/gui_ext.jl

@@ -55,7 +55,7 @@ function Makie.plot!(fig::Union{Makie.Figure,Makie.GridPosition}, res::expfit_st
         ax = Axis(fig[1, 1], xlabel="time (s)", ylabel="Signal (a.u.)")
 
     elseif res[1].seq in [NMRInversions.PFG]
-        ax = Axis(fig[1, 1], xlabel="b factor", ylabel="Signal (a.u.)")
+        ax = Axis(fig[1, 1], xlabel="b factor (s/m² e-9)", ylabel="Signal (a.u.)")
     end
 
     for (i, r) in enumerate(res)
@@ -107,15 +107,15 @@ function Makie.plot!(fig::Union{Makie.Figure,Makie.GridPosition}, res::NMRInvers
 
     # Make axes
     if res[1].seq in [NMRInversions.IR]
-        ax1 = Axis(fig[:, 1], xlabel="time", ylabel="Signal")
-        ax2 = Axis(fig[:, 2], xlabel="T₁", xscale=log10)
+        ax1 = Axis(fig[:, 1], xlabel="time (s)", ylabel="Signal (a.u.)")
+        ax2 = Axis(fig[:, 2], xlabel="T₁ (s)", xscale=log10)
 
     elseif res[1].seq in [NMRInversions.CPMG]
-        ax1 = Axis(fig[:, 1], xlabel="time", ylabel="Signal")
-        ax2 = Axis(fig[:, 2], xlabel="T₂", xscale=log10)
+        ax1 = Axis(fig[:, 1], xlabel="time (s)", ylabel="Signal (a.u.)")
+        ax2 = Axis(fig[:, 2], xlabel="T₂ (s)", xscale=log10)
 
     elseif res[1].seq in [NMRInversions.PFG]
-        ax1 = Axis(fig[:, 1], xlabel="b factor", ylabel="Signal")
+        ax1 = Axis(fig[:, 1], xlabel="b factor (s/m² e-9)", ylabel="Signal (a.u.)")
         ax2 = Axis(fig[:, 2], xlabel="D (m²/s)", xscale=log10)
     end
 

src/exp_fits.jl

@@ -61,10 +61,22 @@ Arguments:
 - `x` : acquisition x parameter (time for relaxation or b-factor for diffusion).
 - `y` : acquisition y parameter (magnetization).
 - `solver` : optimization solver (default is BFGS).
+
+
+Note that initial conditions are automatically determined. 
+If you get NaN for any of the resulting parameters, 
+try changing the initial conditions by calling the funtion using
+the `u0` argument instead of `n`.
 """
 function expfit(n::Int, seq::Type{<:NMRInversions.pulse_sequence1D}, x::Vector, y::Vector; kwargs...)
-    u0 = (ones(2 * n))
-    u0[2:2:end] .= [10.0^(1-x)  for x in 1:n]
+
+    # Make an educated guess of u0
+    if seq == PFG
+        u0 = [v for l in 1:n for v in ( abs(y[1])/(3*l -2) ,maximum(x)/100 * 10.0^(-(l/2 -0.5)) ) ]
+    else 
+        u0 = [v for l in 1:n for v in ( abs(y[1])/(3*l -2) ,maximum(x)/5 * 10.0^(-(l/2)) ) ]
+    end
+
     return expfit(u0, seq, x, y; kwargs...)
 end
 

src/finding_alpha.jl

@@ -45,10 +45,12 @@ function solve_gcv(svds::svd_kernel_struct, solver::Union{regularization_solver,
     α = []
     φ = []
     f_star = []
+    count = 0
 
     done = false
-    while ~done
+    while ~done && count <= 20
 
+        count += 1
         display("Testing α = $(round(αₙ,sigdigits=3))")
         f, r = solve_regularization(svds.K, svds.g, αₙ, solver)
 
@@ -63,6 +65,9 @@ function solve_gcv(svds::svd_kernel_struct, solver::Union{regularization_solver,
             done = true
 
             display("The optimal α is $(round(α[end-1],digits=3))")
+            if count == 20
+                @warn("Reached maximum iterations looking for alpha, results might be inaccurate")
+            end
         end
 
     end

src/inversion_functions.jl

@@ -32,12 +32,18 @@ function invert(seq::Type{<:pulse_sequence1D}, x::AbstractArray, y::Vector;
         X = lims
     elseif isa(lims, Type{<:pulse_sequence1D})
         if lims == PFG
-            X = exp10.(range(-10, -7, 128))
+            X = exp10.(range(-11, -8, 128))
         else
             X = exp10.(range(-5, 1, 128))
         end
     end
 
+
+    # Change scale to match bfactor, which is s/m²e-9, undo below to go back to SI
+    if seq == PFG
+        X .= X .* 1e9
+    end
+
     α = 1 #placeholder, will be replaced below 
 
     if isa(alpha, Real)
@@ -62,6 +68,11 @@ function invert(seq::Type{<:pulse_sequence1D}, x::AbstractArray, y::Vector;
 
     isreal(y) ? SNR = NaN : SNR = calc_snr(y)
 
+
+    if seq == PFG
+        X .= X ./ 1e9
+    end
+
     weighted_average =  f'X / sum(f)
 
     return inv_out_1D(seq, x, y, x_fit, y_fit, X, f, r, SNR, α, weighted_average)

test/runtests.jl

@@ -187,12 +187,23 @@ function test_phase_correction(plots=false)
     return abs(2π - (ϕd + ϕc)) < 0.05
 end
 
+function test_expfit()
+
+    x = [range(0.001, 3, 32)...]
+    u = [3.0, 0.2 ,4.0 ,0.05]
+    y = mexp(CPMG, u, x) + 0.01 .* randn(length(x))
+    data = input1D(CPMG, x, y)
+    results = expfit(2,data)
+
+    return sum((sort(u) .- sort(results.u)) .^ 2) < 0.5
+end
 
 
 @testset "NMRInversions.jl" begin
     # Write your tests here.
     @test test1D(IR)
     @test testT1T2()
+    @test test_expfit()
     @test test_phase_correction()
 
 end