Inversion recovery-CPMG experiment (2D ILT)

[1]:

# This file is part of ILTpy examples.
# Author : Dr. Davis Thomas Daniel
# Last updated : 25.08.2025

This example shows steps to invert a \(T_1-T_2\) correlation dataset with data compression. The same workflow may be used for other 2D inversions.

Importing iltpy and loading data

[2]:

# import iltpy
import iltpy as ilt
# other libraries for handling data, plotting
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('../../../examples/latex_style.mplstyle')
from iltpy.output.plotting import iltplot

print(f"ILTpy version: {ilt.__version__}")

ILTpy version: 1.0.0

  • Load already processed data (i.e. after Fourier transformation and routine NMR processing steps)

[3]:

t1t2_corr_data = np.loadtxt('../../../examples/NMR/t1_t2_correlation_2d/t1_t2_correlation_data.txt') ## 2D data
t_recovery = np.loadtxt('../../../examples/NMR/t1_t2_correlation_2d/t_recovery.txt') ## T1 recovery delays
t_cpmg = np.loadtxt('../../../examples/NMR/t1_t2_correlation_2d/t_cpmg.txt') ## cpmg delays

[4]:

t1t2_corr_data.shape,t_cpmg.shape,t_recovery.shape  # CPMG delays along axis 0 and T1 recovery delays along axis 1

[4]:

((256, 22), (256,), (22,))

[5]:

# Plot the data
fig,ax = plt.subplots(figsize=(3,2.5))
a = ax.contourf(t_recovery,t_cpmg,t1t2_corr_data,cmap='jet',levels=500)
ax.set_yscale('log')
ax.set_xscale('log')
ax.set_ylabel(r'$t_2$')
_=ax.set_xlabel(r'$t_1$')
../_images/Gallery_plot_2D_NMR_8_0.png

Data preparation

[6]:

# noise estimation by baseline correction
fig,ax = plt.subplots(ncols=3,figsize=(7,2))
## select a noise region
ax[0].plot(t1t2_corr_data[100:,-1])

## fit a straight line
y = t1t2_corr_data[100:,-1]
x = np.arange(0,len(y))
coeff = np.polyfit(x,y,deg=1)
y_fit = np.polyval(coeff,x=x)

# fit of the noise region
ax[1].plot(x,y)
ax[1].plot(x,y_fit)

# fit of residuals
residuals = y-y_fit
_=ax[2].plot(x,residuals,'k')
../_images/Gallery_plot_2D_NMR_10_0.png 
[7]:

print(f'Standard deviation : {np.std(residuals)}')

Standard deviation : 8.705347107839648

[8]:

# Set noise variance to unity
t1t2_corr_data = t1t2_corr_data/8.705347107839648

Iltpy Workflow

Load data

[9]:

# load data into an IltData object.
t1t2ILT = ilt.iltload(data=t1t2_corr_data,t=[t_cpmg,t_recovery])

[10]:

t1t2ILT.data

[10]:

array([[ 5.29039071e+01,  5.54526280e+01,  5.45551911e+01, ...,
         1.50889063e+01,  1.46252305e+01,  1.64081386e+01],
       [ 5.05885199e+01,  4.96372367e+01,  4.76688584e+01, ...,
         1.29081345e+01,  1.15859108e+01,  1.28213823e+01],
       [ 4.54551806e+01,  4.54821037e+01,  4.64902245e+01, ...,
         1.19179625e+01,  1.22410398e+01,  1.19718087e+01],
       ...,
       [ 1.06196704e+00,  9.48291692e-01,  2.00726727e+00, ...,
         1.05299267e+00,  2.14786572e+00, -5.23504877e-01],
       [ 2.57564400e+00,  1.13974490e+00,  1.19658258e-02, ...,
        -5.32479246e-01,  1.67521561e-01,  1.43889055e+00],
       [ 3.20684130e+00,  2.94359314e+00,  2.77906303e+00, ...,
         1.24743734e+00, -5.86325462e-01,  7.17949546e-02]],
      shape=(256, 22))

Initialization and inversion

  • To reduce computational needs, data can compressed using singular value decomposition along both dimensions.

  • During initialization, set the compress parameter to True for both dimensions.

  • s specifies the number of the singular values.

  • For 2D inversions, another variation of the algorithm which performs faster, is used and is specified by alt_g = 8.

[11]:

# define the parameters as a dict
p_dict = {'compress':[True,True],'s':[12,22],'alt_g':8}
# define tau for both dimensions
tau = [np.logspace(-5,-1,130),np.logspace(-5,5,130)]
# initialize with exponential kernel in both dimensions
t1t2ILT.init(tau,kernel=[ilt.Exponential(),ilt.Exponential()],parameters=p_dict,verbose=2)

Initializing...
Reading fit parameters...
Setting fit parameters...
Initializing fit parameters...
Initializing kernel matrices...
    Initializing kernel matrix in dimension 0...
    SVD compression in dimension 0...
    Initializing kernel matrix in dimension 1...
    SVD compression in dimension 1...
    Kernel initialization done.
Initializing regularization matrices...
   Calculating constant, regularization matrices...
   Calculating diagonal matrices...
   Calculating amplitude regularization...
   Regularization matrices initialization done.
Initialization took 0.092 seconds.

[13]:

t1t2ILT.compress

[13]:

[True, True]

[11]:

t1t2ILT.invert()

Starting iterations ...

100%|██████████| 100/100 [03:30<00:00,  2.11s/it]

Done.

Plot results

[12]:

# ILTpy's plotting function iltplot is used.
# transform_g=True, which is False by default, applies the transformaton :  g = np.sign(g) * np.sqrt(np.abs(g)) and makes it easier to spot weak features.
fig,ax = iltplot(t1t2ILT,transform_g=True,negative_g=True) # for inversion recovery, the negative of g is plotted
ax[1].set_xlabel(r'$log(T_1)$')
_=ax[1].set_ylabel(r'$log(T_2)$')
../_images/Gallery_plot_2D_NMR_23_0.png

Note

  • Inversion recovery curves may exhibit non-zero baselines, which may be interpreted by the algorithm as a slow relaxing component (the negative peaks around \(10^{4}\)).

  • This artefact manifests in the obtained relaxation distribution as a negative peak with a long relaxation time constant, as seen above.

  • The artefact may be separated from the features of interest by choosing a tau with a higher maximum value.