From 724882dddb180a09c473be22e451e421e8ec0790 Mon Sep 17 00:00:00 2001 From: Aliaksandr Yakutovich Date: Thu, 24 Oct 2024 11:41:51 +0200 Subject: [PATCH] First working version of the rfactor module. --- cleedpy/cli/rfactor.py | 41 ++++++++++++-------- cleedpy/preprocessing.py | 2 +- cleedpy/rfactor.py | 80 ++++++++++++++++++++++++++++++++++++---- 3 files changed, 100 insertions(+), 23 deletions(-) diff --git a/cleedpy/cli/rfactor.py b/cleedpy/cli/rfactor.py index 849c0a8..371320c 100644 --- a/cleedpy/cli/rfactor.py +++ b/cleedpy/cli/rfactor.py @@ -1,23 +1,34 @@ -import click import numpy as np -import yaml - -from .. import rfactor - - -@click.command("cli") -@click.option("--config", "-c", help="Config file", required=True) -def cli(config): +import typer + +from ..rfactor import compute_rfactor + + +def rfactor( + experimental_iv: str = typer.Option( # noqa: B008 + "experimental.txt", "--experimental", "-e", help="Experimental IV curves" + ), + theoretical_iv: str = typer.Option( # noqa: B008 + "theoretical.txt", "--theoretical", "-t", help="Theoretical IV curves" + ), + rfactor_type: str = typer.Option( # noqa: B008 + "r2", "--rfactor", "-r", help="Rfactor type: r2 or pendry" + ), +): """Command line interface for the rfactor tool""" - with open(config) as f: - data = yaml.safe_load(f) + exper_iv = np.loadtxt(experimental_iv) + theor_iv = np.loadtxt(theoretical_iv) + + r = compute_rfactor( + theoretical_iv=theor_iv, experimental_iv=exper_iv, rfactor_type=rfactor_type + ) - y_true = np.array(data["y_true"]) - y_pred = np.array(data["y_pred"]) + print(f"Rfactor={r}") - r_factor = rfactor.mean_square_error(y_true=y_true, y_pred=y_pred) - print(f"R factor: {r_factor}") +def cli(): + """Leed CLI.""" + typer.run(rfactor) if __name__ == "__main__": diff --git a/cleedpy/preprocessing.py b/cleedpy/preprocessing.py index 0501df6..b9cef89 100644 --- a/cleedpy/preprocessing.py +++ b/cleedpy/preprocessing.py @@ -2,7 +2,7 @@ from scipy.interpolate import CubicHermiteSpline -def lorentzian_smoothing(curve, vi): +def lorentzian_smoothing(curve, vi=4.0): """ Draws a Lorentzian curve in the same grid as the given curve. Performs the convolution of the two curves. diff --git a/cleedpy/rfactor.py b/cleedpy/rfactor.py index e201b1b..c128362 100644 --- a/cleedpy/rfactor.py +++ b/cleedpy/rfactor.py @@ -1,7 +1,10 @@ import numpy as np +from scipy.interpolate import CubicSpline +from .preprocessing import lorentzian_smoothing -def r2_factor(theoretical_curve, experimental_curve): + +def r2_factor(experimental_curve, theoretical_curve): """ Calculates R2-factor of the curve of a specific spot (x,y) of the experiment. A curve is represented by an array of (e,i) points, where e = energy and i = intensity. @@ -18,11 +21,6 @@ def r2_factor(theoretical_curve, experimental_curve): it_avg = (S it)/ dE """ - # [TODO] (not sure) Use the length of the shortest curve - min_length = min(len(theoretical_curve), len(experimental_curve)) - theoretical_curve = theoretical_curve[:min_length] - experimental_curve = experimental_curve[:min_length] - # Extract the it, ie values for theoretical and experimental intensity it = theoretical_curve[:, 1] ie = experimental_curve[:, 1] @@ -42,7 +40,7 @@ def r2_factor(theoretical_curve, experimental_curve): return r2 -def rp_factor(theoretical_curve, experimental_curve): +def rp_factor(experimental_curve, theoretical_curve): """ Calculates Pendry's R-factor of the curve of a specific spot (x,y) of the experiment. A curve is represented by an array of (e,i) points, where e = energy and i = intensity. @@ -118,3 +116,71 @@ def energy_step(energies): """ return energies[1:] - energies[:-1] + + +def split_in_pairs(exp_iv, theo_iv): + """This function takes a complete datasets with a experimental and theoretical iv curves splits them them in pairs corresponding to the same index.""" + + beam_indices = np.unique(np.vstack([exp_iv, theo_iv])[:, [0, 1]], axis=0) + + for indx1, indx2 in beam_indices: + exp_curve = exp_iv[(exp_iv[:, 0] == indx1) & (exp_iv[:, 1] == indx2)][ + :, [-2, -1] + ] + theo_curve = theo_iv[(theo_iv[:, 0] == indx1) & (theo_iv[:, 1] == indx2)][ + :, [-2, -1] + ] + yield exp_curve, theo_curve + + +def find_common_x_axis(x1, x2): + min_x = max(np.min(x1), np.min(x2)) + max_x = min(np.max(x1), np.max(x2)) + + x1 = x1[(x1 >= min_x) & (x1 <= max_x)] + x2 = x2[(x2 >= min_x) & (x2 <= max_x)] + + return np.unique(np.sort(np.concatenate([x2, x1]))) + + +def compute_rfactor(experimental_iv, theoretical_iv, shift=0.0, rfactor_type="r2"): + + r_tot = 0.0 + rfactor = { + "r2": r2_factor, + "pendry": rp_factor, + } + + # Looping over the pairs of experimental and theoretical curves corresponding to the same index. + for exp_curve, theo_curve in split_in_pairs(experimental_iv, theoretical_iv): + if len(exp_curve) < 1 or len(theo_curve) < 1: + continue + # Perform Lorentzian smoothing. + exp_curve = lorentzian_smoothing(exp_curve) + theo_curve = lorentzian_smoothing(exp_curve) + + # Shifting the theoretical curve. + theo_curve[:, 0] += shift + + # Finding common x axis. + common_x = find_common_x_axis(theo_curve[:, 0], exp_curve[:, 0]) + + exp_spline = CubicSpline(x=exp_curve[:, 0], y=exp_curve[:, 1])(common_x) + theo_spline = CubicSpline(x=theo_curve[:, 0], y=theo_curve[:, 1])(common_x) + + # Uncomment for testing. + # import matplotlib.pyplot as plt # noqa: E800 + # plt.plot(theo_curve[:,0], theo_curve[:,1], 'o', label='Data Points') # noqa: E800 + # plt.plot(common_x, theo_spline, label='Cubic Spline') # noqa: E800 + # #plt.plot(exp_curve_before_smooth[:,0], exp_curve_before_smooth[:,1], label='Before smooth') # noqa: E800 + # plt.legend() # noqa: E800 + # plt.show() # noqa: E800 + + r = rfactor[rfactor_type]( + np.column_stack([common_x, exp_spline]), + np.column_stack([common_x, theo_spline]), + ) + + r_tot += r + + return r_tot