add CPMG-CPMG and SR-CPMG

docs/src/types_structs.md

@@ -20,6 +20,8 @@ SR
 CPMG
 PFG
 IRCPMG
+SRCPMG
+CPMGCPMG
 ```
 
 ## Inversion solvers

example_data/spinsolve_txt/CPMGCPMG/acqu.par

@@ -0,0 +1,53 @@
+90Amplitude1H = -6
+180Amplitude1H = 0
+accumulate = "yes"
+acqTime = 0.032
+autoPhase = "yes"
+b1Freq1H = 44.081126099999999
+bandwidth = 1000
+bigDelta = 200
+dataDirectory = "C:\Users\Lab\Desktop\Yiheng Xiao\2026\02\09"
+dummyEchoes = 0
+dummyScans = 1
+duration = 140.254
+duration = 5248.1844068998471
+durSpoil = 2
+dwellTime = 1
+echoShift = 0
+echoTime = 135
+experiment = "T2T2"
+expName = "260209-045705 T2T2 (80%-water-200)"
+flatFilter = "no"
+gSpoil = 100
+incExpNr = "yes"
+logspace = "yes"
+maxEcho = 200
+minEcho = 0.1
+minEchoSteps = 8
+nrEchoes = 512
+nrEchoSteps = 40
+nrPnts = 32
+nrScans = 32
+nucleus = "1H-1H"
+percentageCompleted = 100
+percentageCompleted = 100
+phaseMethod = "minsd"
+pulseLength1H = 18
+repTime = 4000
+rxChannel = "1H"
+rxGain = 28
+rxPhase = 0
+saveData = "true"
+smoothFactor = 1
+softwareVersion = "2.01.18"
+specID = "SPA1722"
+specType = "P43Grad"
+useEndDelay = 1
+usePhaseCycle = "yes"
+useStartDelay = 0
+x_maximum = 1000
+x_minimum = 1
+x_steps = 50
+y_maximum = 1000
+y_minimum = 1
+y_steps = 50

ext/gui_2D.jl

@@ -83,20 +83,24 @@ function Makie.plot!(fig::Union{Makie.Figure,Makie.GridPosition}, res::NMRInvers
     x = res.X_indirect
     y = res.X_direct
 
-    xlbl = if seq == IRCPMG
+    xlbl = if seq in (IRCPMG, SRCPMG)
         L"T_1 \, \textrm{(s)}"
     elseif seq == PFGCPMG
         L"D \, \textrm{(m^2/s)}"
+    elseif seq == CPMGCPMG
+        L"T_{2A} \,\textrm{(s)}"
     end
 
-    ylbl = if seq in [IRCPMG, PFGCPMG]
+    ylbl = if seq in (IRCPMG, SRCPMG, PFGCPMG)
         L"T_2 \,\textrm{(s)}"
+    elseif seq == CPMGCPMG
+        L"T_{2B} \,\textrm{(s)}"
     end
 
-    x_low = seq == IRCPMG ? exp10(floor(log10(min(x[1],y[1])))) : exp10(floor(log10(x[1])))
-    x_high = seq == IRCPMG ? exp10(ceil(log10(max(x[end],y[end])))) : exp10(ceil(log10(x[end])))
-    y_low = seq == IRCPMG ? exp10(floor(log10(min(x[1],y[1])))) : exp10(floor(log10(y[1])))
-    y_high = seq == IRCPMG ? exp10(ceil(log10(max(x[end],y[end])))) : exp10(ceil(log10(y[end])))
+    x_low = seq in (IRCPMG, SRCPMG, CPMGCPMG) ? exp10(floor(log10(min(x[1],y[1])))) : exp10(floor(log10(x[1])))
+    x_high = seq in (IRCPMG, SRCPMG, CPMGCPMG) ? exp10(ceil(log10(max(x[end],y[end])))) : exp10(ceil(log10(x[end])))
+    y_low = seq in (IRCPMG, SRCPMG, CPMGCPMG) ? exp10(floor(log10(min(x[1],y[1])))) : exp10(floor(log10(y[1])))
+    y_high = seq in (IRCPMG, SRCPMG, CPMGCPMG) ? exp10(ceil(log10(max(x[end],y[end])))) : exp10(ceil(log10(y[end])))
 
     gr = fig[1:10, 1:10] = GridLayout()
 
@@ -203,7 +207,7 @@ function draw_on_axes(
     markers = collect('a':'z')
     colors = cgrad(:tab10)
 
-    if res.seq == IRCPMG
+    if res.seq in (IRCPMG, SRCPMG, CPMGCPMG)
         plot_diagonal(axmain,x,y)
     end
 
@@ -220,11 +224,14 @@ function draw_on_axes(
         xc = xdist ⋅ x / sum(spo)
         yc = ydist ⋅ y / sum(spo)
 
-        lbl = if res.seq == IRCPMG
+        lbl = if res.seq in (IRCPMG, SRCPMG)
             "T₁/T₂ = $(round(wa_indir[i]/wa_dir[i] , digits=1)) \nVolume = $(round(volumes[i] *100 ,digits = 1))%"
         elseif res.seq == PFGCPMG
             "D = $(round(wa_indir[i], sigdigits=3)), T₂ = $(round(wa_dir[i], sigdigits=3))\nVolume = $(round(volumes[i] *100 ,digits = 1))%"
 
+        elseif res.seq == CPMGCPMG
+            "T₂ₐ = $(round(wa_indir[i], sigdigits=3)), T₂ᵦ = $(round(wa_dir[i], sigdigits=3))\nVolume = $(round(volumes[i] *100 ,digits = 1))%"
+
         end
 
         alphas = 0.83
@@ -342,14 +349,18 @@ function Makie.plot(res::NMRInversions.inv_out_2D)
     coord_label = Observable("")
     coord_label[] = "Hover mouse over plot to view coordinates."
 
-    xlbl = if res.seq == IRCPMG
+    xlbl = if res.seq in (IRCPMG, SRCPMG)
         "T₁"
     elseif res.seq == PFGCPMG
         "D"
+    elseif res.seq == CPMGCPMG
+        "T₂ₐ"
     end
 
-    ylbl = if res.seq in [IRCPMG, PFGCPMG]
+    ylbl = if res.seq in (IRCPMG, SRCPMG, PFGCPMG)
         "T₂"
+    elseif res.seq == CPMGCPMG
+        "T₂ᵦ"
     end
 
     on(events(gui).mouseposition) do _

ext/gui_data.jl

@@ -91,7 +91,7 @@ function Makie.plot(data::NMRInversions.input2D)
     ax_dir.xlabel = "time (s)"
     ax_full.xlabel = "time (s)"
 
-    if data.seq == IRCPMG
+    if data.seq in (IRCPMG, SRCPMG, CPMGCPMG)
         ax_indir.xlabel = "time (s)"
         ax_full.ylabel = "time (s)"
     elseif data.seq == PFGCPMG

src/import_data.jl

@@ -135,8 +135,8 @@ function read_tnt_de0_times(filename, len::Int)
         data_length = read(io,Int32) # how much data (in bytes)
         skip(io, data_length)
 
-        readuntil(io, "de0:2")
-        readuntil(io, "de0:2")
+        readuntil(io, "de0:2") # 1st instance is not relevant
+        readuntil(io, "de0:2") # 2nd instance contains data
         read(io,Int32) # discard this (empty bytes)
         x = Vector{Float64}(undef,0)
         while !eof(io)
@@ -147,7 +147,15 @@ function read_tnt_de0_times(filename, len::Int)
             end
         end
         if length(x) > 1
-            return x[1:len]
+            if length(x) >= len
+                return x[1:len]
+            else
+                extra_points = len - length(x)
+                extrapolated = x[end] .* ( (x[end] / x[end-1]) .^ (1:extra_points) )
+                @warn "de0:2 points not enough, extrapolating logarithmically \
+                to match size of the data."
+                return [x ; extrapolated]
+            end
         else
             throw("Could not read de0 times.")
         end
@@ -195,7 +203,7 @@ function import_tnt(seq::Type{<:Union{pulse_sequence1D, pulse_sequence2D}},
         x = read_tnt_echo_times(filename, length(y), echotime)
         return autophase(input1D(seq, x, y))
 
-    elseif seq == IRCPMG
+    elseif seq in (IRCPMG, CPMGCPMG, SRCPMG)
 
         data_mat = reshape(y , :, header["actual_points_2d"]) 
         x_dir = read_tnt_echo_times(filename, size(data_mat,1), echotime)
@@ -242,6 +250,8 @@ function import_spinsolve(files=pick_multi_file(pwd()))
         seq = PFG
     elseif exp == "DT2"
         seq = PFGCPMG
+    elseif exp == "T2T2"
+        seq = CPMGCPMG
     else
         error("Unrecognised pulse sequence")
     end
@@ -257,6 +267,10 @@ function import_spinsolve(files=pick_multi_file(pwd()))
     elseif seq in [PFGCPMG]
 
         return spinsolve_read_PFGCPMG(acqufile, datafile)
+
+    elseif seq in [CPMGCPMG]
+
+        return spinsolve_read_CPMGCPMG(acqufile, datafile)
     end
 
 end
@@ -327,7 +341,40 @@ function spinsolve_read_PFGCPMG(acqufile, datafile)
 
 end
 
+function spinsolve_read_CPMGCPMG(acqufile, datafile)
+
+    n_echoes = parse(Int32 ,read_acqu(acqufile, "nrEchoes"))
+    t_echo = parse(Float64, read_acqu(acqufile, "echoTime")) * 1e-6
+    τ_steps = parse(Int32, read_acqu(acqufile, "nrEchoSteps"))
+    τ_min = parse(Float64, read_acqu(acqufile, "minEcho")) * 1e-3
+    τ_max = parse(Float64, read_acqu(acqufile, "maxEcho")) * 1e-3
+
+    Raw = readdlm(datafile, ' ')
 
+    if size(Raw, 2) == 1
+        Raw = readdlm(datafile, ',')
+    end
+
+    Data = collect(transpose(complex.(Raw[:, 1:2:end], Raw[:, 2:2:end])))
+
+    ## Make time arrays
+    # Time array in direct dimension
+    t_direct = collect(1:n_echoes) * t_echo
+
+    # Time array in direct dimension
+    if  read_acqu(acqufile, "logspace") == "yes" # if log spacing is selected, do log array
+
+        t_indirect = exp10.(range(log10(τ_min), log10(τ_max), τ_steps))
+    else                   # otherwise, do a linear array
+
+        t_indirect = collect(range(τ_min, τ_max, τ_steps))
+    end
+
+    # make it so that it's multiples of echotime
+    #=map!(x -> round(x / t_echo) * t_echo, t_indirect)=#
+
+    return input2D(CPMGCPMG, t_direct, t_indirect, Data)
+end
 export import_geospec
 """
     import_geospec(dir)

src/inversion_functions.jl

@@ -150,10 +150,10 @@ function invert(
     end
 
     if isa(lims2, Type{<:pulse_sequence2D})
-        if lims2 == IRCPMG
-            X_indirect = collect(NMRInversions.logrange(minimum(x_indirect)/7, maximum(x_indirect)*7, 64)) 
-        elseif lims2 == PFGCPMG
+        if lims2 == PFGCPMG
             X_indirect = collect(NMRInversions.logrange(minimum(x_indirect)*7, maximum(x_indirect)*10, 64)) 
+        else
+            X_indirect = collect(NMRInversions.logrange(minimum(x_indirect)/7, maximum(x_indirect)*7, 64)) 
         end
     elseif isa(lims2, Tuple)
         X_indirect = exp10.(range(lims2...))

src/kernels.jl

@@ -139,6 +139,12 @@ function create_kernel(seq::Type{<:pulse_sequence2D},
     elseif seq == PFGCPMG
         K_dir = create_kernel(CPMG, x_direct, X_direct, y = vec(Data[:,1]))
         K_indir = create_kernel(PFG, x_indirect, X_indirect, y = vec(Data[1,:]))
+    elseif seq == SRCPMG
+        K_dir = create_kernel(CPMG, x_direct, X_direct, y = vec(Data[:,end]))
+        K_indir = create_kernel(SR, x_indirect, X_indirect, y = vec(Data[1,:]))
+    elseif seq == CPMGCPMG
+        K_dir = create_kernel(CPMG, x_direct, X_direct, y = vec(Data[:,1]))
+        K_indir = create_kernel(CPMG, x_indirect, X_indirect, y = vec(Data[1,:]))
     else
         error("2D inversion not yet implemented for $(seq)")
     end
@@ -158,7 +164,18 @@ function create_kernel(seq::Type{<:pulse_sequence2D},
 
     g̃ = diag(usv_indir.U[:, si]' * G' * usv_dir.U[:, sj])
     Ũ₀ = Array{Float64}(undef, 0, 0) # no such thing as U in this case, it is absorbed by g 
-    Ṽ₀ = repeat(usv_dir.V[:, sj], size(usv_indir.V, 1), 1) .* reshape(repeat(usv_indir.V[:, si]', size(usv_dir.V, 1), 1), ñ, size(usv_indir.V,1)*size(usv_dir.V,1) )'
+
+    Ṽ₀ = zeros(size(usv_indir.V, 1) * size(usv_dir.V, 1), ñ)
+    for k in 1:ñ
+        # Get the specific vectors for this significant component
+        v_indir = usv_indir.V[:, si[k]]
+        v_dir   = usv_dir.V[:, sj[k]]
+
+        # The Kronecker product combines them into a single basis vector
+        # This represents the joint "space" of both dimensions
+        Ṽ₀[:, k] = kron(v_indir, v_dir)
+    end
+
     K̃₀ = Diagonal(s̃) * Ṽ₀'
 
     return svd_kernel_struct(K̃₀, g̃, Ũ₀, s̃, Ṽ₀)

src/misc.jl

@@ -102,9 +102,9 @@ function autophase(data::input2D; rotation::Real=0)
 
         d = real.(y_phased[1,1] - y_phased[1,end])
 
-        if data.seq in (PFGCPMG,) && d < 0
+        if data.seq in (PFGCPMG, CPMGCPMG) && d < 0
             return autophase(data, rotation = pi)
-        elseif data.seq in (IRCPMG,) && d > 0
+        elseif data.seq in (IRCPMG, SRCPMG) && d > 0
             return autophase(data, rotation = pi)
         end
     end
@@ -225,7 +225,7 @@ function weighted_averages(r::inv_out_2D; silent::Bool = false)
         volumes[i] = sum(spo)/sum(z)
 
         dir_lbl = ["<T₂>","s"]
-        ind_lbl = if r.seq == IRCPMG
+        ind_lbl = if r.seq in (IRCPMG, SRCPMG)
             ["<T₁>","s"]
         elseif r.seq == PFGCPMG
             ["<D>", "m²/s"]
@@ -235,7 +235,7 @@ function weighted_averages(r::inv_out_2D; silent::Bool = false)
             display("Selection $(collect('a':'z')[i]) :")
             display(ind_lbl[1] * " = $(round(wa_indir[i], sigdigits=4)) "*ind_lbl[2])
             display(dir_lbl[1] * " = $(round(wa_dir[i], sigdigits=4)) "*dir_lbl[2])
-            if r.seq == IRCPMG
+            if r.seq in (IRCPMG, SRCPMG)
                 display("T₁/T₂ = $(round(wa_indir[i]/wa_dir[i], sigdigits=2)) ")
             end
             display("Volume = $(round(volumes[i], sigdigits=4) * 100) %")

src/types.jl

@@ -1,7 +1,7 @@
 ## The following are custom types, defined mainly for multiple dispatch purposes
 
 export pulse_sequence1D, pulse_sequence2D
-export IR, SR, CPMG, PFG, IRCPMG, PFGCPMG
+export IR, SR, CPMG, PFG, IRCPMG, SRCPMG, PFGCPMG, CPMGCPMG
 export alpha_optimizer, regularization_solver
 
 abstract type pulse_sequence1D end
@@ -44,12 +44,24 @@ It can be used wherever the `seq` argument is required.
 """
 struct IRCPMG <: pulse_sequence2D end
 
+"""
+Saturation recovery - CPMG pulse sequence for 2D relaxation experiments (T1-T2).
+The direct dimension is the T2, and the indirect dimension is the T1 acquisition times.
+It can be used wherever the `seq` argument is required.
+"""
+struct SRCPMG <: pulse_sequence2D end
+
 """
 PFG - CPMG pulse sequence for 2D D-T2 experiments.
 It can be used wherever the `seq` argument is required.
 """
 struct PFGCPMG <: pulse_sequence2D end
 
+"""
+CPMG - CPMG pulse sequence for 2D T2-T2 experiments.
+It can be used wherever the `seq` argument is required.
+"""
+struct CPMGCPMG <: pulse_sequence2D end
 
 export regularization_solver
 

test/runtests.jl

@@ -31,6 +31,34 @@ function test_artificial_data(seq::Type{<:pulse_sequence1D}, SNR = 100 ; kwargs.
     #=return f=#
 end
 
+function T2T2_data()
+
+    x_direct = range(1:5000,500) .* (250 * 1e-6)
+    x_indirect = logrange(1,5000,32) .* (250 * 1e-6)
+    X = exp10.(range(-5, 1, 128)) # T range
+    f = [0.5exp.(-(x+0.7)^2 / 0.2) + exp.(-(x - 1.3)^2 / 1.9) for x in range(-5, 5, length(X))]
+
+    # T2a = 5e-1
+    # T2b = 5e-2
+    # A = [exp.(-x_direct ./ T2a)  exp.(-x_direct ./ T2b)]
+    # B = [exp.(-x_indirect ./ T2a)  exp.(-x_indirect ./ T2b)]
+    # K = [0.5  0.0 ; 
+    #     0.1  0.4]
+
+    A = create_kernel(CPMG, x_direct, X)
+    B = create_kernel(CPMG, x_indirect, X)
+    K = Diagonal(f) # off diagonal peaks are for exchange
+
+    G = A * K * B'
+
+    noise_real = (maximum(G) * 0.01) .* randn(size(G)...)
+    noise_imag = (maximum(G) * 0.01) .* randn(size(G)...)
+
+    data = input2D(CPMGCPMG, x_direct, x_indirect, complex.(G .+ noise_real, noise_imag))
+
+    return data
+end
+
 function test_artificial_data_2D() #throws errors, needs fixing
 
     x_direct = exp10.(range(log10(1e-4), log10(5), 1024)) # acquisition range