MCRALS example using hardConc #609
-
|
Hello, is it possible to have a working example using the hardConc constraint? I have three concentration profiles as initial guess here (c_initialguess.txt) for this data set (cell4pos1.txt). I would like to fix only the values at the very start so that one of them is set to 1 (100% of that species) and the others to 0. I am afraid I might be using the hardConc functionality wrong.
It only fixes the first point of the concentration profiles but there is a sharp change right after it so I guess it is not working very well... is there a better way to do it? Thank you. |
Beta Was this translation helpful? Give feedback.
Replies: 5 comments 1 reply
-
Beta Was this translation helpful? Give feedback.
-
|
Hi, that's a good idea but I did not know you could enforce a hard constraint on the spectra. Is it also with hardConc? |
Beta Was this translation helpful? Give feedback.
-
Beta Was this translation helpful? Give feedback.
-
|
Hello @GhostDeini , We have added the possibility of adding hard constraints on spectra (to be added soon in the dev version of spectrochempy). Here are the main results, for C_hard:
and for C (obtained after calculation of St from C_hard and data)
The behaviour is probably closer from what you expected but note that there are some negative values on concentrations (also present when fixing the initial concentrations). As fixing a spectral profile adds more constraints than fixing initial concentrations, there are more residuals, but the magnitude of the error is comparable in both cases. Below are the squared rtesiduals (top: when fixing the spectrum, bottom when fixing the initial concentrations)
and the sum of residuals vs time or energy (line: fixed inital concentration, dotted line: fixed spectral profile):
Does this help you ? or do you need to fix more constraints ? We will release the next dev version with these changes within few days -- it will come with some changes of the synthax (this will be in the documentation) and I will post here the script to get the figures above. I also notice that you are importing data from from textfiles, probably obtained from another source. You might consider writing a specific reader -- this would allow you to get directly the correct labelling of energy and time axes, units, etc..., We can help if you are interested in this (see here) Best regards |
Beta Was this translation helpful? Give feedback.
-
|
Dear @atravert , Thanks! This looks much better. I am looking forward to the new version. No need to use the reader for the moment. Best, |
Beta Was this translation helpful? Give feedback.
-
|
Dear Lucia, I apologize for the late reply. Here is the script, which works with the current stable and dev version of Scpy. Note that there are some changes in the synthax (first initiate a MCRALS instance with options, and then use fit() method), but the internal algorithm is the same as previously. Best, # %%
import spectrochempy as scp
import numpy as np
# %%
#energy = np.loadtxt('energy.txt')
arr_pos1 = np.loadtxt('cell4pos1.txt')
c_init = np.loadtxt('c_initialguess.txt')
data = scp.NDDataset(
#np.transpose(arr_pos1),
arr_pos1,
name="Dataset Cell 4 Pos 1",
author="Lucia P.",
description="XAS dataset from a battery cycling operando experiment"
)
c_guess = scp.NDDataset(
c_init,
name="Dataset initial guesses concentration",
author="Lucia P.",
description="Initial concentration guesses",
)
# %%
# fix initial concentration
def get_C2(C):
C[0,0]=0
C[0,1]=1
C[0,2]=0
return C
mcr = scp.MCRALS(unimodConc=[],
hardConc=[0,1,2], getConc=get_C2, tol=30.0)
mcr.fit(data, c_guess)
_ = mcr.C_hard.T.plot() ; _.set_title('Chard, mccr1') ; scp.show()
_ = mcr.St.plot() ; _.set_title('St, mccr1') ; scp.show()
_ = mcr.C.T.plot() ; _.set_title('C, mcr1') ; scp.show()
data_hat1 = mcr.inverse_transform()
# %%
# %%
# fix spectral profile #1 equal to first spectrum
def get_St(St): #St must be passed, even if not actuially used
return data[0,:]
mcr2 = scp.MCRALS(log_level="INFO", max_iter=350, unimodConc=[],
hardSpec = [1], hardSt_to_St_idx = [0], getSpec=get_St, tol=3.0)
mcr2.fit(data, c_guess)
# plot concentrationand spactra profiles
_ = mcr2.C.T.plot() ; _.set_title('C, mcr2') ; scp.show()
_ = mcr2.St.plot() ; _.set_title('St, mcr2') ; scp.show()
data_hat2 = mcr2.inverse_transform()
_ = mcr2.plotmerit(data, data_hat2) ; _.set_title('merit, mcr2') ; scp.show()
# plot squared residuals
_ = ((data - data_hat1)**2).plot()
_ = ((data - data_hat2)**2 + 0.0002).plot(clear=False) ; _.set_title('squared residuals') ; _.set_ylim(0, 0.0004) ; scp.show()
# plot sum of squared residuals vs time
_ = scp.sum((data - data_hat1)**2, dim=1).plot()
_ = scp.sum((data - data_hat2)**2, dim=1).plot(clear=False) ; _.set_title('squared residuals vs time') ; scp.show()
# plot sum of squared residuals vs time
_ = scp.sum((data - data_hat1)**2, dim=0).plot()
_ = scp.sum((data - data_hat2)**2, dim=0).plot(clear=False) ; _.set_title('squared residuals vs energy') ; scp.show() |
Beta Was this translation helpful? Give feedback.
Dear Lucia,
I apologize for the late reply. Here is the script, which works with the current stable and dev version of Scpy. Note that there are some changes in the synthax (first initiate a MCRALS instance with options, and then use fit() method), but the internal algorithm is the same as previously.
Currently, the profiles are computed with lstsq, and non-negativity constraints applied by brute force (all values are set to 0). We are currently working on adding nnls, I will lety you know when done.
Do not hesitate ti ask if you have any question.
Best,
Arnaud