.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/tutorial_compart_modeling.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_tutorial_compart_modeling.py: Comparment Modelling ==================== Here in this tutorial, we will use skultrafast to do global analysis. We are using the approach as presented by van Stokkum (doi:10.1016/j.bbabio.2004.04.01). The general outline is as follows: 1. Generate a decay associated spectrum 2. Build the transfer matrix K, choosing starting concentrations j 3. Generate transformation matricies from K and j 4. Apply them to the basis vectors of the DAS Here we assume step 1 is already done. skultrafast has a module to help with step 2. For that we need the `Model` class. Please note that the module is quite barebones and may still have bugs. .. GENERATED FROM PYTHON SOURCE LINES 20-20 .. code-block:: Python :dedent: 1 .. GENERATED FROM PYTHON SOURCE LINES 22-29 .. code-block:: Python from skultrafast import dataset, data_io, plot_helpers import numpy as np from matplotlib import pyplot as plt from skultrafast.kinetic_model import Model model = Model() .. GENERATED FROM PYTHON SOURCE LINES 30-34 To build the transfer matrix, we have to add transitions. The if one the compartments is not part of the model it will be added automatically. The following line means that the model has an transition from S2 to S1 with a rate k2. Since no yield is given, it defaults to 1. .. GENERATED FROM PYTHON SOURCE LINES 34-37 .. code-block:: Python model.add_transition('S2', 'S1', 'k2') .. GENERATED FROM PYTHON SOURCE LINES 38-40 Lets also assume, that S1 decays to zero. The comparment zero is special, since it is not modelled explicty and the name `zero` is reserved. .. GENERATED FROM PYTHON SOURCE LINES 40-43 .. code-block:: Python model.add_transition('S1', 'zero', 'k1') .. GENERATED FROM PYTHON SOURCE LINES 44-46 To get the corresponding matrix from the model we have to call the `build_matrix`-function .. GENERATED FROM PYTHON SOURCE LINES 46-50 .. code-block:: Python mat = model.build_matrix() mat .. rst-class:: sphx-glr-script-out .. code-block:: none [[[(-1, 'k2', 1)], []], [[(1, 'k2', 1)], [(-1, 'k1', 1)]]] .. GENERATED FROM PYTHON SOURCE LINES 51-54 To fill the matrix with numeric values we can use the `build_matfunc` method. This gives us a function which takes all free parameters, first the rates and then the yields. .. GENERATED FROM PYTHON SOURCE LINES 54-58 .. code-block:: Python f = model.build_mat_func() f .. rst-class:: sphx-glr-script-out .. code-block:: none .K_func at 0x77f8a287be20> .. GENERATED FROM PYTHON SOURCE LINES 59-62 Now, what can we do with that matrix? We can use it to project the DAS to the described model. For that, we first evaluate the matrix nummerically by supplying values for the parameters .. GENERATED FROM PYTHON SOURCE LINES 62-66 .. code-block:: Python num_mat = f(k2=2, k1=1) num_mat .. rst-class:: sphx-glr-script-out .. code-block:: none array([[-2., 0.], [ 2., -1.]]) .. GENERATED FROM PYTHON SOURCE LINES 67-74 Next, we need the eigenvectors of the matrix. Note, that the presented approach assumes that the eigenvalues are simple eigenvalues. If this is not the case, one has to use the jordan normal form. As we see, the eigenvalues of the matrix are the negative rates. That basically means, that the eigenbasis of the problem is given by a diagonal transfer matrix, which is the parallel model described by an DAS. Hence, the eigenvectors allow us to transform the DAS to a SAS and vice versa by using the inverse. .. GENERATED FROM PYTHON SOURCE LINES 74-78 .. code-block:: Python vals, vecs = np.linalg.eig(num_mat) vals .. rst-class:: sphx-glr-script-out .. code-block:: none array([-1., -2.]) .. GENERATED FROM PYTHON SOURCE LINES 79-80 To continue, we also need to choose the starting concentrations. .. GENERATED FROM PYTHON SOURCE LINES 80-85 .. code-block:: Python j = np.zeros(len(vals)) j[0] = 1 j .. rst-class:: sphx-glr-script-out .. code-block:: none array([1., 0.]) .. GENERATED FROM PYTHON SOURCE LINES 86-87 The transformation matrix is then given by .. GENERATED FROM PYTHON SOURCE LINES 87-91 .. code-block:: Python A = vecs @ np.diag(np.linalg.inv(vecs) @ j) A_inv = np.linalg.inv(A) .. GENERATED FROM PYTHON SOURCE LINES 92-95 Now `DAS @ A_inv` should give the SAS, while `A_inv @ basis_vecs` should return the time-depedence of the concentrations. Let's test that on some data. Load test data and correct the dispersion. .. GENERATED FROM PYTHON SOURCE LINES 95-103 .. code-block:: Python plot_helpers.enable_style() wl, t, d = data_io.load_example() ds = dataset.TimeResSpec(wl, t, d) dsb = ds.bin_freqs(50) # Bin dataset to save computation while building docs res = dsb.estimate_dispersion(heuristic_args=(1.5, ), deg=3, shift_result=.15) .. image-sg:: /auto_examples/images/sphx_glr_tutorial_compart_modeling_001.png :alt: tutorial compart modeling :srcset: /auto_examples/images/sphx_glr_tutorial_compart_modeling_001.png, /auto_examples/images/sphx_glr_tutorial_compart_modeling_001_1_50x.png 1.50x, /auto_examples/images/sphx_glr_tutorial_compart_modeling_001_2_00x.png 2.00x :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/skultrafast/envs/latest/lib/python3.11/site-packages/skultrafast/data_io.py:75: UserWarning: Reading `.npy` or `.npz` file required additional header parsing as it was created on Python 2. Save the file again to speed up loading and avoid this warning. wl, data, t = a["wl"], a["data"], a["t"] .. GENERATED FROM PYTHON SOURCE LINES 104-105 Fit the DAS first. .. GENERATED FROM PYTHON SOURCE LINES 105-111 .. code-block:: Python ids = res.correct_ds fr = ids.fit_exp([0.0, 0.08, 1, 500000], model_coh=True, fix_sigma=False, fix_t0=False) ids.plot.das() .. image-sg:: /auto_examples/images/sphx_glr_tutorial_compart_modeling_002.png :alt: tutorial compart modeling :srcset: /auto_examples/images/sphx_glr_tutorial_compart_modeling_002.png, /auto_examples/images/sphx_glr_tutorial_compart_modeling_002_1_50x.png 1.50x, /auto_examples/images/sphx_glr_tutorial_compart_modeling_002_2_00x.png 2.00x :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [, ] .. GENERATED FROM PYTHON SOURCE LINES 112-113 Make sequential Model .. GENERATED FROM PYTHON SOURCE LINES 113-119 .. code-block:: Python m = Model() m.add_transition('S1hot', 'S1', 'k1') m.add_transition('S1', 'zero', 'const') m.build_matrix() .. rst-class:: sphx-glr-script-out .. code-block:: none [[[(-1, 'k1', 1)], []], [[(1, 'k1', 1)], [(-1, 'const', 1)]]] .. GENERATED FROM PYTHON SOURCE LINES 120-125 Make transformation matricies. Notice that we reverse the order of the eigenvectors. Initally, they are sorted by they eigenvalue, therefore the slowest componet comes first. Since this order is inverse of the skultrafast order, we reserve it so we don't have to reverse our spectra and basis vectors later. .. GENERATED FROM PYTHON SOURCE LINES 125-140 .. code-block:: Python func = m.build_mat_func() num_mat = func(1 / fr.lmfit_res.params['t0'], 1/fr.lmfit_res.params['t1']) vals, vecs = np.linalg.eig(num_mat) # Reverse order vecs = vecs[:, ::-1] j = np.zeros(len(vals)) j[0] = 1 A = (vecs @ np.diag(np.linalg.inv(vecs) @ j)).T A_inv = np.linalg.inv(A) .. GENERATED FROM PYTHON SOURCE LINES 141-145 The DAS are members of the fitter object. Since we are also modeling coherent contributions, we only take the first two components. Also, the eigenvalue routine sorts the eigenvectors after their eigenvalue, which is the opposite of our usal sorting, where the fast component comes first. T .. GENERATED FROM PYTHON SOURCE LINES 145-158 .. code-block:: Python fig, ax = plt.subplots(2, figsize=(3, 4)) das = fr.fitter.c[:, :2] ax[0].plot(dsb.wn, das) ax[0].set_title('DAS') plot_helpers.lbl_spec(ax[0]) sas = A_inv @ das.T edas = np.cumsum(das, axis=1) ax[1].plot(dsb.wn, sas.T) ax[1].set_title('SAS') plot_helpers.lbl_spec(ax[1]) .. image-sg:: /auto_examples/images/sphx_glr_tutorial_compart_modeling_003.png :alt: DAS, SAS :srcset: /auto_examples/images/sphx_glr_tutorial_compart_modeling_003.png, /auto_examples/images/sphx_glr_tutorial_compart_modeling_003_1_50x.png 1.50x, /auto_examples/images/sphx_glr_tutorial_compart_modeling_003_2_00x.png 2.00x :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 159-161 As we can see, we sucessfully get the SAS, which in this case are just EDAS. Let's also look at the concentrations. .. GENERATED FROM PYTHON SOURCE LINES 161-170 .. code-block:: Python fig, ax = plt.subplots(2, figsize=(3, 4)) ax[0].plot(dsb.t, fr.fitter.x_vec[:, :2]) ax[0].set_title('DAS') plot_helpers.lbl_trans(ax[0], use_symlog=False) ct = fr.fitter.x_vec[:, :2] @ A ax[1].plot(dsb.t, ct) ax[1].set_title('SAS') plot_helpers.lbl_trans(ax[1], use_symlog=False) .. image-sg:: /auto_examples/images/sphx_glr_tutorial_compart_modeling_004.png :alt: DAS, SAS :srcset: /auto_examples/images/sphx_glr_tutorial_compart_modeling_004.png, /auto_examples/images/sphx_glr_tutorial_compart_modeling_004_1_50x.png 1.50x, /auto_examples/images/sphx_glr_tutorial_compart_modeling_004_2_00x.png 2.00x :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 171-173 So why does it work? The dataset is given by the outer product of the concentrations and the spectrum .. GENERATED FROM PYTHON SOURCE LINES 173-180 .. code-block:: Python C = fr.fitter.x_vec[:, :2] S = das[:, :2] fit = C @ S.T fit - fr.fitter.model .. rst-class:: sphx-glr-script-out .. code-block:: none array([[-1.35052986e-251, -1.54721491e-251, -1.18337669e-251, ..., -3.66060512e-252, -1.35137006e-251, -5.72057605e-252], [-5.58819129e-245, -6.40153469e-245, -4.89550858e-245, ..., -1.51465272e-245, -5.59255679e-245, -2.36722887e-245], [-1.88998042e-238, -2.16488894e-238, -1.65534791e-238, ..., -5.12260691e-239, -1.89176359e-238, -8.00682388e-239], ..., [ 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, ..., 0.00000000e+000, 0.00000000e+000, 0.00000000e+000], [ 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, ..., 0.00000000e+000, 0.00000000e+000, 0.00000000e+000], [ 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, ..., 0.00000000e+000, 0.00000000e+000, 0.00000000e+000]], shape=(409, 50)) .. GENERATED FROM PYTHON SOURCE LINES 181-183 Now we can insert `1 = A @ A_inv`. As expected this does not change the product. .. GENERATED FROM PYTHON SOURCE LINES 183-186 .. code-block:: Python ((C@A) @ (A_inv @ S.T)) - (C @ S.T) .. rst-class:: sphx-glr-script-out .. code-block:: none array([[ 0.00000000e+000, 0.00000000e+000, 1.18305219e-271, ..., -3.69703808e-273, 0.00000000e+000, 0.00000000e+000], [ 4.96208052e-265, 0.00000000e+000, 0.00000000e+000, ..., 0.00000000e+000, 0.00000000e+000, -2.48104026e-265], [ 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, ..., -3.25194909e-260, 0.00000000e+000, -5.20311854e-259], ..., [-3.55271368e-015, 0.00000000e+000, 0.00000000e+000, ..., 0.00000000e+000, 0.00000000e+000, 0.00000000e+000], [ 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, ..., 0.00000000e+000, 0.00000000e+000, -1.77635684e-015], [ 0.00000000e+000, 0.00000000e+000, 0.00000000e+000, ..., 0.00000000e+000, 0.00000000e+000, 0.00000000e+000]], shape=(409, 50)) .. GENERATED FROM PYTHON SOURCE LINES 187-189 All steps above are done in the `make_sas` method of the the fitresult, taking a `Model` as a parameter. If yields are required, they also have to be supplied. .. GENERATED FROM PYTHON SOURCE LINES 189-200 .. code-block:: Python sas, ct = fr.make_sas(model, {}) fig, ax = plt.subplots(2, figsize=(3, 4)) ax[0].plot(fr.fitter.wl, sas.T) ax[0].set_title('SAS') plot_helpers.lbl_trans(ax[0], use_symlog=False) ct = fr.fitter.x_vec[:, :2] @ A ax[1].plot(dsb.t, ct) ax[1].set_title('Convoluted Concentrations') plot_helpers.lbl_trans(ax[1], use_symlog=False) .. image-sg:: /auto_examples/images/sphx_glr_tutorial_compart_modeling_005.png :alt: SAS, Convoluted Concentrations :srcset: /auto_examples/images/sphx_glr_tutorial_compart_modeling_005.png, /auto_examples/images/sphx_glr_tutorial_compart_modeling_005_1_50x.png 1.50x, /auto_examples/images/sphx_glr_tutorial_compart_modeling_005_2_00x.png 2.00x :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 2.640 seconds) .. _sphx_glr_download_auto_examples_tutorial_compart_modeling.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: tutorial_compart_modeling.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: tutorial_compart_modeling.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: tutorial_compart_modeling.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_