Commit 04af07df authored by Henri Payno's avatar Henri Payno
Browse files

continue Spectrum + FT clean up

parent b3bd1b3a
Pipeline #78820 passed with stages
in 6 minutes and 55 seconds
......@@ -200,14 +200,6 @@ class Spectra:
else:
return array.reshape(shape)
def keys(self) -> tuple:
"""
return the tuple of keys contained in the spectra
"""
if len(self.data.flat) > 0:
assert isinstance(self.data.flat[0], Spectrum)
return tuple(self.data.flat[0].keys())
def __eq__(self, other):
if not isinstance(other, Spectra):
return False
......
import numpy
from typing import Union
from typing import Iterable
import logging
from est.units import units
from est.units import ur
......@@ -87,8 +86,7 @@ class Spectrum:
self.__e0 = None
self.__noise_savgol = None
self.__norm_noise_savgol = None
self.__other_parameters = {}
self.ft = {}
self.ft = _FT()
self.__pymca_dict = {}
# pymca dict use to store some processing information specific to pymca
self.__larch_dict = {}
......@@ -239,8 +237,10 @@ class Spectrum:
def ft(self, ft):
if isinstance(ft, _FT):
self.__ft = ft
elif isinstance(ft, dict):
self.__ft = _FT.from_dict(ft)
else:
self.__ft = _FT(ddict=ft)
raise TypeError
@property
def pymca_dict(self):
......@@ -271,9 +271,6 @@ class Spectrum:
return (_energy_len, _mu_len)
def extra_keys(self) -> tuple:
return tuple(self.__other_parameters.keys())
def load_frm_dict(self, ddict: dict):
assert isinstance(ddict, dict)
......@@ -334,7 +331,6 @@ class Spectrum:
self._PYMCA_DICT_KEY: self.pymca_dict,
self._LARCH_DICT_KEY: self.larch_dict,
}
res.update(self.__other_parameters)
return res
def __str__(self):
......@@ -352,46 +348,8 @@ class Spectrum:
"normalized_energy",
):
main_info = add_info(str_=main_info, attr=info)
def add_third_info(str_, key):
sub_str = ("- " + key + ": " + str(self[key])) + "\n"
return str_ + sub_str
for key in self.__other_parameters:
main_info = add_third_info(str_=main_info, key=key)
return main_info
def update(self, obj):
"""
Update the contained values from the given obj.
:param obj:
:type obj: Union[XASObject, dict]
"""
if isinstance(obj, Spectrum):
_obj = obj.to_dict()
else:
_obj = obj
assert isinstance(_obj, dict)
for key, value in _obj.items():
if isinstance(value, str) and value == "None":
self[key] = None
else:
self[key] = value
def get_missing_keys(self, keys: Iterable) -> tuple:
"""Return missing keys on the spectrum"""
missing = []
for key in keys:
if key not in self or self[key] is None:
missing.append(key)
return tuple(missing)
def keys(self) -> list:
keys = list(self.__other_parameters.keys())
keys += list(self.__key_mapper.keys())
return keys
def copy(self):
return Spectrum.from_dict(self.to_dict())
......@@ -411,21 +369,10 @@ class _FT:
_IMAGINARY_KEY = "FTImaginary"
def __init__(self, ddict):
self.__radius = numpy.nan
self.__intensity = numpy.nan
self.__imaginary = numpy.nan
self.__other_parameters = {}
self.__key_mapper = {
self._RADIUS_KEY: self.__class__.radius,
self._INTENSITY_KEY: self.__class__.intensity,
self._IMAGINARY_KEY: self.__class__.imaginary,
}
if ddict is not None:
for key, values in ddict.items():
self[key] = values
def __init__(self):
self.__radius = None
self.__intensity = None
self.__imaginary = None
@property
def radius(self):
......@@ -451,36 +398,24 @@ class _FT:
def imaginary(self, imaginery):
self.__imaginary = imaginery
def __getitem__(self, key):
"""Need for pymca compatibility"""
if key in self.__key_mapper:
return self.__key_mapper[key].fget(self)
else:
return self.__other_parameters[key]
def __setitem__(self, key, value):
"""Need for pymca compatibility"""
if key in self.__key_mapper:
self.__key_mapper[key].fset(self, value)
else:
self.__other_parameters[key] = value
def __contains__(self, item):
return item in self.__key_mapper or item in self.__other_parameters
def to_dict(self) -> dict:
res = {
self._RADIUS_KEY: self.radius,
self._INTENSITY_KEY: self.intensity,
self._IMAGINARY_KEY: self.imaginary,
}
res.update(self.__other_parameters)
return res
def get_missing_keys(self, keys: Iterable) -> tuple:
"""Return missing keys on the spectrum"""
missing = []
for key in keys:
if key not in self:
missing.append(key)
return tuple(missing)
def load_from_dict(self, ddict: dict):
if self._RADIUS_KEY in ddict:
self.radius = ddict[self._RADIUS_KEY]
if self._INTENSITY_KEY in ddict:
self.intensity = ddict[self._INTENSITY_KEY]
if self._IMAGINARY_KEY in ddict:
self.imaginary = ddict[self._IMAGINARY_KEY]
@staticmethod
def from_dict(ddict: dict):
res = _FT()
res.load_from_dict(ddict=ddict)
return res
%% Cell type:markdown id: tags:
# PyMca XAS data processing
PyMca is a spectroscopy library from python: https://github.com/vasole/pymca
%% Cell type:markdown id: tags:
## Example EXAFS spectrum
%% Cell type:code id: tags:
``` python
%matplotlib inline
import matplotlib.pyplot as plt
from est.tests.data import example_spectrum
from est.resources import resource_path
filename = str(resource_path("exafs/EXAFS_Cu.dat"))
energy, mu = example_spectrum("exafs/EXAFS_Cu.dat")
plt.figure(figsize=(20, 10))
plt.plot(energy, mu)
plt.xlabel("Energy (eV)")
plt.ylabel("mu")
print("Example EXFAS spectrum")
```
%% Output
Example EXFAS spectrum
%% Cell type:markdown id: tags:
## Import data
The main class dealing with XAS data in *pymca* is the `XASClass`
%% Cell type:code id: tags:
``` python
from PyMca5.PyMcaPhysics.xas.XASClass import XASClass
pymca_xas = XASClass()
pymca_xas.setSpectrum(energy, mu)
```
%% Cell type:markdown id: tags:
In *est* the main class is `XASObject`
%% Cell type:code id: tags:
``` python
from est.core.io import read as read_pymca_xas
from est.io.utils import url as url_utils
xas_data = read_pymca_xas(
spectra_url=url_utils.build_spec_data_url(
file_path=filename,
col_name="Column 2",
),
channel_url=url_utils.build_spec_data_url(
file_path=filename,
col_name="Column 1",
),
dimensions=(2, 1, 0),
)
```
%% Cell type:markdown id: tags:
### Normalization
Pre-edge and post-edge fitting and normalization
$$
\mu_{\mathrm{norm}}(E) = \frac{\mu(E)-\mathrm{pre}(E)}{\mathrm{post}(E)-\mathrm{pre}(E)}
$$
%% Cell type:code id: tags:
``` python
norm_pymca = pymca_xas.normalize()
ax1, ax2 = plt.subplots(1, 2, figsize=(20, 10))[-1]
ax1.plot(energy, mu, label="ln(I0/I)")
ax1.plot(
norm_pymca["NormalizedEnergy"], norm_pymca["NormalizedBackground"], label="signal"
)
ax1.plot(
norm_pymca["NormalizedEnergy"], norm_pymca["NormalizedSignal"], label="background"
)
ax1.set_xlabel("Energy (eV)")
ax1.set_ylabel("mu")
ax1.legend()
ax2.plot(norm_pymca["NormalizedEnergy"], norm_pymca["NormalizedMu"])
ax2.set_xlabel("Energy (eV)")
ax2.set_ylabel("mu_norm")
print(list(norm_pymca))
```
%% Output
['Jump', 'JumpNormalizationMethod', 'Edge', 'NormalizedEnergy', 'NormalizedMu', 'NormalizedBackground', 'NormalizedSignal', 'NormalizedPlotMin', 'NormalizedPlotMax']
%% Cell type:code id: tags:
``` python
from est.core.process.pymca.normalization import pymca_normalization
xas_norm = pymca_normalization(xas_data.copy())
norm_est = xas_norm.spectra.data.flatten()[0]
ax1, ax2 = plt.subplots(1, 2, figsize=(20, 10))[-1]
ax1.plot(energy, mu)
ax1.plot(norm_est.normalized_energy, norm_est.pre_edge)
ax1.plot(norm_est.normalized_energy, norm_est.post_edge)
ax1.set_xlabel("Energy (eV)")
ax1.set_ylabel("mu")
ax1.legend()
ax2.plot(norm_est.normalized_energy, norm_est.normalized_mu)
ax2.set_xlabel("Energy (eV)")
ax2.set_ylabel("mu_norm")
```
%% Output
[est.core.types.spectra] INFO : fail to access to NormalizedEnergy
[est.core.types.spectra] INFO : fail to access to NormalizedEnergy
[est.core.types.spectra] INFO : fail to access to NormalizedMu
[est.core.types.spectra] INFO : fail to access to NormalizedMu
[est.core.types.spectra] INFO : fail to access to NormalizedSignal
[est.core.types.spectra] INFO : fail to access to NormalizedSignal
[est.io.io] WARNING : Unable to write at results/NormalizedEnergy reason is One of data, shape or dtype must be specified
[est.io.io] WARNING : Unable to write at results/NormalizedMu reason is One of data, shape or dtype must be specified
[est.io.io] WARNING : Unable to write at results/NormalizedSignal reason is One of data, shape or dtype must be specified
[matplotlib.legend] WARNING : No artists with labels found to put in legend. Note that artists whose label start with an underscore are ignored when legend() is called with no argument.
Text(0, 0.5, 'mu_norm')
%% Cell type:markdown id: tags:
## EXAFS signal
$$
\chi(k) = \frac{\mu_{\mathrm{norm}}(k) - \mu_0(k)}{\mu_0(k)}
$$
$$
k = \frac{\sqrt{2m(E-E_0)}}{\hbar^2}
$$
%% Cell type:code id: tags:
``` python
from PyMca5.PyMcaPhysics.xas.XASClass import e2k
k = e2k(energy - norm_pymca["Edge"])
# norm_mu = norm_est["NormalizedMu"]
norm_mu = mu - norm_est.pre_edge
exafs_pymca = pymca_xas.postEdge(k=k, mu=norm_mu)
ax1, ax2 = plt.subplots(1, 2, figsize=(20, 10))[-1]
ax1.plot(k, norm_mu, label="mu_norm")
ax1.plot(k, exafs_pymca["PostEdgeB"], label="mu_0")
ax1.set_xlim(0, None)
ax1.set_ylim(0, 5)
ax1.set_xlabel("Wavenumber (Å^-1)")
ax1.set_ylabel("χ(k)")
ax1.legend()
chi = (norm_mu - exafs_pymca["PostEdgeB"]) / exafs_pymca["PostEdgeB"]
exafs = chi * k**2
ax2.plot(k, exafs)
ax2.set_xlabel("Wavenumber (Å^-1)")
ax2.set_ylabel("k^2 χ(k) (Å^-2)")
ax2.set_xlim(0, None)
ax2.set_ylim(-5, 4)
print(list(exafs_pymca))
```
%% Output
[PyMca5.PyMcaPhysics.xas.XASClass] WARNING : Error: dimension of knots must be dimension of polDegree+1
[PyMca5.PyMcaPhysics.xas.XASClass] WARNING : Forced automatic (equidistant) knot definition.
['PostEdgeK', 'PostEdgeB', 'KnotsX', 'KnotsY', 'KMin', 'KMax', 'KWeight']
%% Cell type:code id: tags:
``` python
from est.core.process.pymca.exafs import pymca_exafs
xas_exafs = pymca_exafs(xas_norm.copy(), exafs={"KWeight": 2})
exafs_est = xas_exafs.spectra.data.flatten()[0]
ax1, ax2 = plt.subplots(1, 2, figsize=(20, 10))[-1]
ax1.plot(exafs_est.k, exafs_est.chi, label="mu_norm")
ax1.plot(exafs_est.k, exafs_est.pymca_dict["PostEdgeB"], label="mu_0")
ax1.set_xlim(0, None)
ax1.set_ylim(0, 5)
ax1.set_xlabel("Wavenumber (Å^-1)")
ax1.set_ylabel("χ(k)")
ax1.legend()
ax2.plot(exafs_est.k, exafs_est.pymca_dict["EXAFSNormalized"])
ax2.set_xlabel("Wavenumber (Å^-1)")
ax2.set_ylabel("k^2 χ(k) (Å^-2)")
ax2.set_xlim(0, None)
ax2.set_ylim(-5, 4)
```
%% Output
[PyMca5.PyMcaPhysics.xas.XASClass] WARNING : Error: dimension of knots must be dimension of polDegree+1
[PyMca5.PyMcaPhysics.xas.XASClass] WARNING : Forced automatic (equidistant) knot definition.
exafs: [####################] 100% DONE
[est.core.types.spectra] INFO : fail to access to EXAFSKValues
[est.core.types.spectra] INFO : fail to access to EXAFSKValues
[est.core.types.spectra] INFO : fail to access to EXAFSSignal
[est.core.types.spectra] INFO : fail to access to EXAFSSignal
[est.io.io] WARNING : Unable to write at results/EXAFSKValues reason is One of data, shape or dtype must be specified
[est.io.io] WARNING : Unable to write at results/EXAFSSignal reason is One of data, shape or dtype must be specified
(-5.0, 4.0)
%% Cell type:markdown id: tags:
### Fourier transform
$$
\mathrm{ft}(R) = \mathrm{FFT}(k^2\chi(k)w(k))
$$
%% Cell type:code id: tags:
``` python
spectrum = xas_exafs.spectra.data.flatten()[0]
ft_pymca = pymca_xas.fourierTransform(
k=k, mu=exafs, kMin=spectrum.pymca_dict["KMin"], kMax=spectrum.pymca_dict["KMax"]
)
plt.figure(figsize=(20, 10))
plt.plot(ft_pymca["FTRadius"], ft_pymca["FTIntensity"])
plt.xlabel("Radius (Å)")
plt.ylabel("|χ(k)| (Å^-3)")
print(list(ft_pymca))
```
%% Output
['Set', 'InterpolatedK', 'InterpolatedSignal', 'KWeight', 'K', 'WindowWeight', 'FTRadius', 'FTIntensity', 'FTReal', 'FTImaginary']
%% Cell type:code id: tags:
``` python
from est.core.process.pymca.ft import pymca_ft
xas_ft = pymca_ft(xas_exafs.copy())
ft_est = xas_ft.spectra.data.flatten()[0]
plt.figure(figsize=(20, 10))
plt.plot(ft_est.ft["FTRadius"], ft_est.ft["FTIntensity"])
plt.plot(ft_est.ft.radius, ft_est.ft.intensity)
plt.xlabel("Radius (Å)")
plt.ylabel("|χ(k)| (Å^-3)")
```
%% Output
ft: [####################] 100% DONE
Text(0, 0.5, '|χ(k)| (Å^-3)')
%% Cell type:code id: tags:
``` python
```
......
......@@ -27,9 +27,9 @@ def test_single_spectrum(spectrum_cu_from_pymca):
xas_obj = pymca_k_weight(xas_obj=xas_obj)
res = pymca_ft(xas_obj=xas_obj)
assert isinstance(res, XASObject)
assert "FTRadius" in res.spectra.data.flat[0].ft
assert "FTImaginary" in res.spectra.data.flat[0].ft
assert "FTIntensity" in res.spectra.data.flat[0].ft
assert res.spectra.data.flat[0].ft.radius is not None
assert res.spectra.data.flat[0].ft.imaginary is not None
assert res.spectra.data.flat[0].ft.intensity is not None
@pytest.mark.skipif(PyMca5 is None, reason="PyMca5 is not installed")
......@@ -47,6 +47,6 @@ def test_single_spectrum_asdict(spectrum_cu_from_pymca):
xas_obj = pymca_k_weight(xas_obj=xas_obj)
res = pymca_ft(xas_obj=xas_obj.to_dict())
assert isinstance(res, XASObject)
assert "FTRadius" in res.spectra.data.flat[0].ft
assert "FTImaginary" in res.spectra.data.flat[0].ft
assert "FTIntensity" in res.spectra.data.flat[0].ft
assert res.spectra.data.flat[0].ft.radius is not None
assert res.spectra.data.flat[0].ft.imaginary is not None
assert res.spectra.data.flat[0].ft.intensity is not None
......@@ -17,6 +17,5 @@ def test_from_mock():
spectra = xas_obj.spectra
assert xas_obj.n_spectrum == 20 * 10
assert xas_obj.n_spectrum == 20 * 10
spectra.keys()
assert spectra.data.flat[0] == spectra[0, 0]
spectra.map_to("mu")
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment