Note
Go to the end to download the full example code.
MCR-ALS with kinetic constraints
In this example, we perform MCR-ALS optimization on UV-Vis spectra from a
three-component reaction A -> B -> C investigated by UV-Vis spectroscopy.
Full details on the reaction and data acquisition conditions can be found in
Bijlsma et al. [2001].
The data can be downloaded from the Biosystems Data Analysis Group, University
of Amsterdam.
For convenience, this dataset is also available in the SpectroChemPy test-data
directory as matlabdata/METING9.MAT.
import numpy as np
import spectrochempy as scp
Loading a NDDataset
Load the data with the read function.
This file contains a pair of datasets. The first dataset contains the time in seconds since the start of the reaction (t=0). The second dataset contains the UV-VIS spectra of the reaction mixture, recorded at different time points. The first column of the matrix contains the wavelength axis and the remaining columns are the measured UV-VIS spectra (wavelengths x timepoints)
print("NDDataset names: " + str([d.name for d in ds]))
NDDataset names: ['RelTime', 'x9b']
We load the experimental spectra (in ds[1]), add the y (time) and x
(wavelength) coordinates, and keep one spectrum of out 4:
A first estimate of the concentrations can be obtained by EFA:
compute EFA...
We can get a better estimate of the concentration (C) and pure spectra profiles (St) by soft MCR-ALS:
mcr_1 = scp.MCRALS(log_level="INFO")
_ = mcr_1.fit(D, C0)
_ = mcr_1.C.T.plot()
_ = mcr_1.St.plot()
Concentration profile initialized with 3 components
Initial spectra profile computed
*** ALS optimisation log ***
#iter reconstruction_error residual_change profile_change trend
-----------------------------------------------------------------------
1 6.124937e-02 9.263378e-01 2.558410e-01 down
2 4.178693e-02 3.175432e-01 5.179759e-02 down
3 3.506613e-02 1.607696e-01 2.907913e-02 down
4 3.244504e-02 7.472313e-02 1.892498e-02 down
5 3.117520e-02 3.912723e-02 1.342159e-02 down
6 3.030733e-02 2.783116e-02 1.006982e-02 down
7 2.953519e-02 2.546995e-02 8.662148e-03 down
8 2.876214e-02 2.616711e-02 8.002595e-03 down
9 2.802142e-02 2.574785e-02 7.313129e-03 down
10 2.732705e-02 2.477523e-02 6.774127e-03 down
11 2.662667e-02 2.562446e-02 6.307022e-03 down
12 2.593077e-02 2.613071e-02 5.766931e-03 down
13 2.531636e-02 2.369021e-02 5.407822e-03 down
14 2.476338e-02 2.183938e-02 5.131193e-03 down
15 2.426729e-02 2.002995e-02 4.873074e-03 down
16 2.382325e-02 1.829488e-02 4.625059e-03 down
17 2.342662e-02 1.664638e-02 4.385941e-03 down
18 2.307295e-02 1.509472e-02 4.155979e-03 down
19 2.275811e-02 1.364342e-02 3.935322e-03 down
20 2.247830e-02 1.229318e-02 3.723783e-03 down
21 2.223001e-02 1.104400e-02 3.521169e-03 down
22 2.200979e-02 9.905328e-03 3.330119e-03 down
23 2.180784e-02 9.174134e-03 3.201974e-03 down
24 2.161986e-02 8.618922e-03 3.111721e-03 down
25 2.144041e-02 8.298906e-03 3.050595e-03 down
26 2.126779e-02 8.050080e-03 3.015103e-03 down
27 2.110246e-02 7.772713e-03 2.971999e-03 down
28 2.094262e-02 7.573616e-03 2.929829e-03 down
29 2.078812e-02 7.376385e-03 2.887971e-03 down
30 2.063882e-02 7.180903e-03 2.846299e-03 down
31 2.049460e-02 6.987324e-03 2.804773e-03 down
32 2.035530e-02 6.795783e-03 2.763376e-03 down
33 2.022081e-02 6.606411e-03 2.722099e-03 down
34 2.009099e-02 6.419330e-03 2.680940e-03 down
35 1.996572e-02 6.234656e-03 2.639901e-03 down
36 1.984486e-02 6.052502e-03 2.598989e-03 down
37 1.972830e-02 5.872971e-03 2.558215e-03 down
38 1.961591e-02 5.696164e-03 2.517590e-03 down
39 1.950690e-02 5.556998e-03 2.476993e-03 down
40 1.939950e-02 5.505160e-03 2.448477e-03 down
41 1.929338e-02 5.469423e-03 2.435079e-03 down
42 1.918859e-02 5.431212e-03 2.425148e-03 down
43 1.908511e-02 5.392239e-03 2.415839e-03 down
44 1.898293e-02 5.353123e-03 2.406447e-03 down
45 1.888204e-02 5.314011e-03 2.396803e-03 down
46 1.878243e-02 5.274921e-03 2.386873e-03 down
47 1.868408e-02 5.235836e-03 2.376656e-03 down
48 1.858698e-02 5.196735e-03 2.366161e-03 down
49 1.849110e-02 5.157595e-03 2.355399e-03 down
50 1.839645e-02 5.118396e-03 2.344380e-03 down
Convergence criterion not reached after 50 iterations.
Stop ALS optimization.
Kinetic constraints can be added, i.e., imposing that the concentration profiles obey a kinetic model. To do so we first define an ActionMAssKinetics object with roughly estimated rate constants:
reactions = ("A -> B", "B -> C")
species_concentrations = {"A": 5.0, "B": 0.0, "C": 0.0}
k0 = np.array((0.5, 0.05))
kin = scp.ActionMassKinetics(reactions, species_concentrations, k0)
The concentration profile obtained with this approximate model can be computed and compared with those of the soft MCR-ALS:
Ckin = kin.integrate(D.y.data)
_ = mcr_1.C.T.plot(linestyle="-", cmap=None)
_ = Ckin.T.plot(clear=False, cmap=None)
Even though very approximate, the same values can be used to run a hard-soft MCR-ALS,
using the public ModelProfile constraint:
import spectrochempy.analysis.constraints as ct
X = D[:, 300.0:500.0]
param_to_optimize = {"k[0]": 0.5, "k[1]": 0.05}
mcr_2 = scp.MCRALS(
constraints=[
ct.ModelProfile(
"C",
components=[0, 1, 2],
model=kin.fit_to_concentrations,
model_args=([0, 1, 2], [0, 1, 2], param_to_optimize),
model_kwargs={"ivp_solver_kwargs": {"return_NDDataset": False}},
)
]
)
_ = mcr_2.fit(X, Ckin)
Optimization terminated successfully.
Current function value: 5.293217
Iterations: 26
Function evaluations: 50
Optimization terminated successfully.
Current function value: 3.044501
Iterations: 24
Function evaluations: 48
Optimization terminated successfully.
Current function value: 2.061456
Iterations: 22
Function evaluations: 44
Optimization terminated successfully.
Current function value: 1.566531
Iterations: 19
Function evaluations: 38
Optimization terminated successfully.
Current function value: 1.293158
Iterations: 22
Function evaluations: 43
Optimization terminated successfully.
Current function value: 1.127349
Iterations: 18
Function evaluations: 36
Optimization terminated successfully.
Current function value: 1.020664
Iterations: 20
Function evaluations: 38
Optimization terminated successfully.
Current function value: 0.945658
Iterations: 20
Function evaluations: 38
Optimization terminated successfully.
Current function value: 0.890437
Iterations: 22
Function evaluations: 39
Optimization terminated successfully.
Current function value: 0.847462
Iterations: 21
Function evaluations: 40
Optimization terminated successfully.
Current function value: 0.810971
Iterations: 20
Function evaluations: 37
Optimization terminated successfully.
Current function value: 0.779836
Iterations: 17
Function evaluations: 33
Optimization terminated successfully.
Current function value: 0.751719
Iterations: 19
Function evaluations: 36
Optimization terminated successfully.
Current function value: 0.726579
Iterations: 22
Function evaluations: 41
Optimization terminated successfully.
Current function value: 0.703097
Iterations: 18
Function evaluations: 35
Optimization terminated successfully.
Current function value: 0.680920
Iterations: 18
Function evaluations: 36
Optimization terminated successfully.
Current function value: 0.660181
Iterations: 16
Function evaluations: 33
Optimization terminated successfully.
Current function value: 0.640675
Iterations: 16
Function evaluations: 32
Optimization terminated successfully.
Current function value: 0.622104
Iterations: 17
Function evaluations: 34
Optimization terminated successfully.
Current function value: 0.603947
Iterations: 18
Function evaluations: 37
Optimization terminated successfully.
Current function value: 0.586621
Iterations: 16
Function evaluations: 32
Optimization terminated successfully.
Current function value: 0.570040
Iterations: 19
Function evaluations: 37
Optimization terminated successfully.
Current function value: 0.554027
Iterations: 17
Function evaluations: 34
Optimization terminated successfully.
Current function value: 0.538478
Iterations: 16
Function evaluations: 33
Optimization terminated successfully.
Current function value: 0.523443
Iterations: 21
Function evaluations: 39
Optimization terminated successfully.
Current function value: 0.509274
Iterations: 21
Function evaluations: 40
Optimization terminated successfully.
Current function value: 0.495585
Iterations: 19
Function evaluations: 37
Optimization terminated successfully.
Current function value: 0.482329
Iterations: 19
Function evaluations: 37
Optimization terminated successfully.
Current function value: 0.469765
Iterations: 20
Function evaluations: 39
Optimization terminated successfully.
Current function value: 0.457282
Iterations: 20
Function evaluations: 38
Optimization terminated successfully.
Current function value: 0.445468
Iterations: 20
Function evaluations: 38
Optimization terminated successfully.
Current function value: 0.434075
Iterations: 20
Function evaluations: 39
Optimization terminated successfully.
Current function value: 0.423038
Iterations: 19
Function evaluations: 37
Optimization terminated successfully.
Current function value: 0.411655
Iterations: 18
Function evaluations: 35
Optimization terminated successfully.
Current function value: 0.401259
Iterations: 19
Function evaluations: 37
Optimization terminated successfully.
Current function value: 0.391298
Iterations: 20
Function evaluations: 39
Optimization terminated successfully.
Current function value: 0.381477
Iterations: 19
Function evaluations: 37
Optimization terminated successfully.
Current function value: 0.372801
Iterations: 18
Function evaluations: 35
Optimization terminated successfully.
Current function value: 0.363693
Iterations: 21
Function evaluations: 39
Optimization terminated successfully.
Current function value: 0.354968
Iterations: 19
Function evaluations: 38
Optimization terminated successfully.
Current function value: 0.346511
Iterations: 19
Function evaluations: 37
Optimization terminated successfully.
Current function value: 0.338322
Iterations: 18
Function evaluations: 35
Optimization terminated successfully.
Current function value: 0.330577
Iterations: 18
Function evaluations: 35
Optimization terminated successfully.
Current function value: 0.322953
Iterations: 17
Function evaluations: 34
Optimization terminated successfully.
Current function value: 0.315715
Iterations: 19
Function evaluations: 37
Optimization terminated successfully.
Current function value: 0.308588
Iterations: 19
Function evaluations: 38
Optimization terminated successfully.
Current function value: 0.301516
Iterations: 19
Function evaluations: 36
Optimization terminated successfully.
Current function value: 0.294675
Iterations: 20
Function evaluations: 39
Optimization terminated successfully.
Current function value: 0.288040
Iterations: 18
Function evaluations: 36
Optimization terminated successfully.
Current function value: 0.281669
Iterations: 16
Function evaluations: 32
Now, let’s compare the concentration profile of the hard-soft modeling with the previous one:
_ = mcr_2.C.T.plot()
_ = mcr_1.C_constrained.T.plot(clear=False, ls="--")
- Finally, let’s plot the pure spectra profiles St, and some the
reconstructed dataset (X_hat = C St) vs original dataset (X) and residuals.
_ = mcr_2.St.plot()
_ = mcr_2.plot_merit(nb_traces=10, offset=5)
This ends the example ! The following line can be uncommented if no plot shows when running the .py script with python
# scp.show()
Total running time of the script: (0 minutes 7.436 seconds)