.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/ILT/2DILT.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_ILT_2DILT.py: 2D ILT ====== Here, we're going to provide a few demonstrations of the ILT functionality. Let's start with Fig 1.10 in A. Beaton's thesis, which is based off the figures in Venkataramanan. .. GENERATED FROM PYTHON SOURCE LINES 8-56 .. code-block:: Python from pylab import ( figure, title, show, linspace, logspace, log10, exp, sqrt, rcParams, ) import numpy as np from pyspecdata import nddata, image rcParams["image.aspect"] = "auto" # needed for sphinx gallery # sphinx_gallery_thumbnail_number = 2 NT1 = 300 # Number of T1 values NT2 = 300 # Number of T2 values LT1_name = r"$\log(T_1)$" LT1 = nddata(linspace(-2.5, 0.5, NT1), LT1_name) LT2_name = r"$\log(T_2)$" LT2 = nddata(linspace(-2.5, 0.3, NT2), LT2_name) mu = [-1.25, -1.75] sigma = [0.1, 0.1] exact_data = exp( -((LT1 - mu[0]) ** 2) / 2 / sigma[0] ** 2 - (LT2 - mu[1]) ** 2 / 2 / sigma[1] ** 2 ) slanted_coord1 = (LT1 + LT2) / sqrt(2) slanted_coord2 = (LT2 - LT1) / sqrt(2) mu = [-1.0, -0.4] mu = [ # convert to slanted coords (mu[0] + mu[1]) / sqrt(2), (mu[1] - mu[0]) / sqrt(2), ] sigma = [0.5, 0.05] # in slanted exact_data += exp( -((slanted_coord1 - mu[0]) ** 2) / 2 / sigma[0] ** 2 - (slanted_coord2 - mu[1]) ** 2 / 2 / sigma[1] ** 2 ) exact_data.reorder(LT2_name) # T₂ along y axis figure(1) title("exact data") image(exact_data) .. image-sg:: /auto_examples/ILT/images/sphx_glr_2DILT_001.png :alt: exact data :srcset: /auto_examples/ILT/images/sphx_glr_2DILT_001.png, /auto_examples/ILT/images/sphx_glr_2DILT_001_2_00x.png 2.00x :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 57-59 Now add the experimental decay dimensions (:math:`\tau_1` and :math:`\tau_2`) .. GENERATED FROM PYTHON SOURCE LINES 59-63 .. code-block:: Python tau1 = nddata(logspace(log10(5.0e-4), log10(4), 30), "tau1") tau2 = nddata(linspace(5.0e-4, 3.8, 1000), "tau2") .. GENERATED FROM PYTHON SOURCE LINES 64-66 Pre-allocate the :math:`\tau_1\times\tau_2` result via ndshape’s `alloc`, inline .. GENERATED FROM PYTHON SOURCE LINES 66-70 .. code-block:: Python simulated_data = (tau1.shape | tau2.shape).alloc(dtype=np.float64) simulated_data.reorder("tau2") # T₂ along y axis .. rst-class:: sphx-glr-script-out .. code-block:: none array([[0., 0., 0., ..., 0., 0., 0.], [0., 0., 0., ..., 0., 0., 0.], [0., 0., 0., ..., 0., 0., 0.], ..., [0., 0., 0., ..., 0., 0., 0.], [0., 0., 0., ..., 0., 0., 0.], [0., 0., 0., ..., 0., 0., 0.]]) dimlabels=['tau2', 'tau1'] axes={`tau2':None +/-None, `tau1':None +/-None} .. GENERATED FROM PYTHON SOURCE LINES 71-78 pySpecData makes it easy to construct fake data like this. Typically this is very easy, but here, we must contend with the fact that we are memory-limited, so if we want a highly resolved fit basis, we need to chunk up the calculation. Nonetheless, pySpecData still makes that easy: let's see how! Block sizes (tune to available RAM) .. GENERATED FROM PYTHON SOURCE LINES 78-82 .. code-block:: Python bLT1 = 20 bLT2 = 20 .. GENERATED FROM PYTHON SOURCE LINES 83-93 Loop over `LT1` and `LT2` in blocks, vectorized over :math:`\tau_{1}` and :math:`\tau_{2}` dims each time .. math:: T_1 = 10^{\log(T_1)} R_1 = 10^{-\log(T_1)} \ln(R_1) = -\log(T_1) \ln(10) .. GENERATED FROM PYTHON SOURCE LINES 93-114 .. code-block:: Python print( "Generating the fake data can take some time. I need to loop a" f" calculation in chunks over a {LT1.shape[LT1_name]} ×" f" {LT2.shape[LT2_name]} grid" ) for i in range(0, LT1.shape[LT1_name], bLT1): LT1_blk = LT1[LT1_name, slice(i, i + bLT1)] B1 = 1 - 2 * exp(-tau1 / 10**LT1_blk) # dims: (tau1, LT1_blk) for j in range(0, LT2.shape[LT2_name], bLT2): print(i, j) LT2_blk = LT2[LT2_name, slice(j, j + bLT2)] B2 = exp(-tau2 / 10**LT2_blk) # dims: (tau2, LT2_blk) # Extract matching block of exact_data data_blk = exact_data[LT1_name, slice(i, i + bLT1)][ LT2_name, slice(j, j + bLT2) ] # dims: (tau1, tau2, LT1_blk, LT2_blk) # Multiply, sum out both LT axes, and accumulate simulated_data += (B2 * B1 * data_blk).real.sum(LT1_name).sum(LT2_name) print("done generating") .. rst-class:: sphx-glr-script-out .. code-block:: none Generating the fake data can take some time. I need to loop a calculation in chunks over a 300 × 300 grid 0 0 0 20 0 40 0 60 0 80 0 100 0 120 0 140 0 160 0 180 0 200 0 220 0 240 0 260 0 280 20 0 20 20 20 40 20 60 20 80 20 100 20 120 20 140 20 160 20 180 20 200 20 220 20 240 20 260 20 280 40 0 40 20 40 40 40 60 40 80 40 100 40 120 40 140 40 160 40 180 40 200 40 220 40 240 40 260 40 280 60 0 60 20 60 40 60 60 60 80 60 100 60 120 60 140 60 160 60 180 60 200 60 220 60 240 60 260 60 280 80 0 80 20 80 40 80 60 80 80 80 100 80 120 80 140 80 160 80 180 80 200 80 220 80 240 80 260 80 280 100 0 100 20 100 40 100 60 100 80 100 100 100 120 100 140 100 160 100 180 100 200 100 220 100 240 100 260 100 280 120 0 120 20 120 40 120 60 120 80 120 100 120 120 120 140 120 160 120 180 120 200 120 220 120 240 120 260 120 280 140 0 140 20 140 40 140 60 140 80 140 100 140 120 140 140 140 160 140 180 140 200 140 220 140 240 140 260 140 280 160 0 160 20 160 40 160 60 160 80 160 100 160 120 160 140 160 160 160 180 160 200 160 220 160 240 160 260 160 280 180 0 180 20 180 40 180 60 180 80 180 100 180 120 180 140 180 160 180 180 180 200 180 220 180 240 180 260 180 280 200 0 200 20 200 40 200 60 200 80 200 100 200 120 200 140 200 160 200 180 200 200 200 220 200 240 200 260 200 280 220 0 220 20 220 40 220 60 220 80 220 100 220 120 220 140 220 160 220 180 220 200 220 220 220 240 220 260 220 280 240 0 240 20 240 40 240 60 240 80 240 100 240 120 240 140 240 160 240 180 240 200 240 220 240 240 240 260 240 280 260 0 260 20 260 40 260 60 260 80 260 100 260 120 260 140 260 160 260 180 260 200 260 220 260 240 260 260 260 280 280 0 280 20 280 40 280 60 280 80 280 100 280 120 280 140 280 160 280 180 280 200 280 220 280 240 280 260 280 280 done generating .. GENERATED FROM PYTHON SOURCE LINES 115-118 `simulated_data` now holds the :math:`\tau_1\times\tau_2` synthetic data. So, add noise, and scale data so that noise has norm of 1 .. GENERATED FROM PYTHON SOURCE LINES 118-123 .. code-block:: Python simulated_data.add_noise(0.1) simulated_data /= 0.1 .. GENERATED FROM PYTHON SOURCE LINES 124-127 Use BRD to find the value of $\lambda$ ($\alpha$). Note that BRD assumes that you have scaled your data so that the stdev of the noise is 1.0. .. GENERATED FROM PYTHON SOURCE LINES 127-143 .. code-block:: Python simulated_data.nnls( ("tau1", "tau2"), (LT1, LT2), ( lambda tau1, LT1: 1 - 2 * exp(-tau1 * 10**-LT1), lambda tau2, LT2: exp(-tau2 * 10**-LT2), ), l="BRD", ) figure(2) title("BRD") image(simulated_data) show() .. image-sg:: /auto_examples/ILT/images/sphx_glr_2DILT_002.png :alt: BRD :srcset: /auto_examples/ILT/images/sphx_glr_2DILT_002.png, /auto_examples/ILT/images/sphx_glr_2DILT_002_2_00x.png 2.00x :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (2 minutes 14.289 seconds) .. _sphx_glr_download_auto_examples_ILT_2DILT.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 2DILT.ipynb <2DILT.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 2DILT.py <2DILT.py>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_