diff --git a/cleedpy/preprocessing.py b/cleedpy/preprocessing.py index 16f1d56..efc1777 100644 --- a/cleedpy/preprocessing.py +++ b/cleedpy/preprocessing.py @@ -1,20 +1,21 @@ -import numpy as np -from cleedpy.rfactor import r2_factor, rp_factor +import numpy as np from scipy.interpolate import CubicHermiteSpline +from cleedpy.rfactor import r2_factor, rp_factor + def lorentzian_smoothing(curve, vi): """ - Draws a Lorentzian curve in the same grid as the given curve. - Performs the convolution of the two curves. - The given curve is represented by an array of (e,i) points, where e = energy and i = intensity. + Draws a Lorentzian curve in the same grid as the given curve. + Performs the convolution of the two curves. + The given curve is represented by an array of (e,i) points, where e = energy and i = intensity. - Inputs: - - curve: numpy array of [e,i] arrays - - vi: predefined constant representing the width + Inputs: + - curve: numpy array of [e,i] arrays + - vi: predefined constant representing the width - Output: - The result of the convolution: a new array of [e,i] arrays + Output: + The result of the convolution: a new array of [e,i] arrays """ E = curve[:, 0] @@ -25,8 +26,8 @@ def lorentzian_smoothing(curve, vi): c = 0 L_value = 0 for i in range(E.size): - c += vi**2 / ((E[i] - E[j])**2 + vi**2) - L_value += (I[i] * vi**2 / ((E[i] - E[j])**2 + vi**2)) + c += vi**2 / ((E[i] - E[j]) ** 2 + vi**2) + L_value += I[i] * vi**2 / ((E[i] - E[j]) ** 2 + vi**2) L.append(L_value / c) return np.array(list(zip(E, L))) @@ -34,31 +35,31 @@ def lorentzian_smoothing(curve, vi): def preprocessing_loop(theoretical_curves, experimental_curves, shift, r_factor, vi): """ - Performs all the preprocessing steps and R-factor calculations, required before the Search. - A curve is represented by an array of (e,i) points, where e = energy and i = intensity. - - Inputs: - - theoretical_curves: numpy array of arrays of [e,i] arrays - - experimental_curves: numpy array of arrays of [e,i] arrays - where: - - number of theoretical curves = number of experimental curves - - shift: provided by user, used to edit the energy value per step - - r_factor: provided by user, name of the preferred R-factor method - - Design: - 1. Lorentzian smoothing of all curves - 2. For each energy point: - a. Shift theoretical curve - b. Find boundaries of overlap - c. Filter both curves to keep only the elements in the overlap - d. For each experimental curve: - i. Spline = interpolate it to the theoretical grid - ii. Calculate R-factor - e. Calculate average R-factor of all curves in this shift - 3. Find the min of the average values - - Output: - The optimal R-factor value per curve. + Performs all the preprocessing steps and R-factor calculations, required before the Search. + A curve is represented by an array of (e,i) points, where e = energy and i = intensity. + + Inputs: + - theoretical_curves: numpy array of arrays of [e,i] arrays + - experimental_curves: numpy array of arrays of [e,i] arrays + where: + - number of theoretical curves = number of experimental curves + - shift: provided by user, used to edit the energy value per step + - r_factor: provided by user, name of the preferred R-factor method + + Design: + 1. Lorentzian smoothing of all curves + 2. For each energy point: + a. Shift theoretical curve + b. Find boundaries of overlap + c. Filter both curves to keep only the elements in the overlap + d. For each experimental curve: + i. Spline = interpolate it to the theoretical grid + ii. Calculate R-factor + e. Calculate average R-factor of all curves in this shift + 3. Find the min of the average values + + Output: + The optimal R-factor value per curve. """ for i in range(len(theoretical_curves)): @@ -70,16 +71,28 @@ def preprocessing_loop(theoretical_curves, experimental_curves, shift, r_factor, print(shifted_theoretical_curves) # Find the boundaries of the overlap after the shift - min_bound = np.maximum(shifted_theoretical_curves[0, 0], experimental_curves[0, 0])[0] - max_bound = np.minimum(shifted_theoretical_curves[0, -1], experimental_curves[0, -1])[0] + min_bound = np.maximum(shifted_theoretical_curves[0, 0], experimental_curves[0, 0])[ + 0 + ] + max_bound = np.minimum( + shifted_theoretical_curves[0, -1], experimental_curves[0, -1] + )[0] print(min_bound, max_bound) # Filter both curves based on the boundaries - filtered_theoretical_curves = shifted_theoretical_curves[shifted_theoretical_curves[:, 0] >= min_bound] - filtered_theoretical_curves = filtered_theoretical_curves[filtered_theoretical_curves[:, 0] <= max_bound] - - filtered_experimental_curves = experimental_curves[experimental_curves[:, 0] >= min_bound] - filtered_experimental_curves = filtered_experimental_curves[filtered_experimental_curves[:, 0] <= max_bound] + filtered_theoretical_curves = shifted_theoretical_curves[ + shifted_theoretical_curves[:, 0] >= min_bound + ] + filtered_theoretical_curves = filtered_theoretical_curves[ + filtered_theoretical_curves[:, 0] <= max_bound + ] + + filtered_experimental_curves = experimental_curves[ + experimental_curves[:, 0] >= min_bound + ] + filtered_experimental_curves = filtered_experimental_curves[ + filtered_experimental_curves[:, 0] <= max_bound + ] print(filtered_theoretical_curves, filtered_experimental_curves) # Spline @@ -89,6 +102,5 @@ def preprocessing_loop(theoretical_curves, experimental_curves, shift, r_factor, theoretical_E_grid = filtered_theoretical_curves[i][:, 0] CubicHermiteSpline(E_values, I_values, theoretical_E_grid) - print(globals().get(r_factor)) return diff --git a/cleedpy/rfactor.py b/cleedpy/rfactor.py index bf9a9a1..0caa424 100644 --- a/cleedpy/rfactor.py +++ b/cleedpy/rfactor.py @@ -3,19 +3,19 @@ def r2_factor(theoretical_curve, experimental_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. + 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. - Inputs: - - theoretical curve: numpy array of [e,i] arrays - - experimental curve: numpy array of [e,i] arrays + Inputs: + - theoretical curve: numpy array of [e,i] arrays + - experimental curve: numpy array of [e,i] arrays - Design: - R2 = sqrt{ S(It - c*Ie)^2 / S(It - It_avg)^2 } + Design: + R2 = sqrt{ S(It - c*Ie)^2 / S(It - It_avg)^2 } - where: - c = sqrt( S|It|^2 / S|Ie|^ 2) - It_avg = (S It)/ dE + where: + c = sqrt( S|It|^2 / S|Ie|^ 2) + It_avg = (S It)/ dE """ # [TODO] (not sure) Use the length of the shortest curve @@ -29,11 +29,11 @@ def r2_factor(theoretical_curve, experimental_curve): # Calculate normalization factor c and It_avg c = np.sqrt(np.sum(It**2) / np.sum(Ie**2)) - It_avg = np.sum(It) / It.size # dE = number of energy steps + It_avg = np.sum(It) / It.size # dE = number of energy steps # Calculate the numerator and denominator of R2 - numerator = np.sum((It - c*Ie)**2) - denominator = np.sum((It - It_avg)**2) + numerator = np.sum((It - c * Ie) ** 2) + denominator = np.sum((It - It_avg) ** 2) # Calculate R2 R2 = np.sqrt(numerator / denominator) @@ -44,15 +44,15 @@ def r2_factor(theoretical_curve, experimental_curve): def rp_factor(theoretical_curve, experimental_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. + 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. - Inputs: - - theoretical curve: numpy array of [e,i] arrays - - experimental curve: numpy array of [e,i] arrays + Inputs: + - theoretical curve: numpy array of [e,i] arrays + - experimental curve: numpy array of [e,i] arrays - Design: - Rp = S(Ye - Yt)^2 / S(Ye^2 + Yt^2) + Design: + Rp = S(Ye - Yt)^2 / S(Ye^2 + Yt^2) """ # Extract the It, Ie values for theoretical and experimental intensity @@ -61,7 +61,7 @@ def rp_factor(theoretical_curve, experimental_curve): # Extract the Et, Ee values for theoretical and experimental energy Et = theoretical_curve[:, 0] - Ee = experimental_curve[:,0] + Ee = experimental_curve[:, 0] # Calculate theoretical and experimental energy steps step_Et = energy_step(Et) @@ -72,7 +72,7 @@ def rp_factor(theoretical_curve, experimental_curve): Ye = Y_function(Ie, step_Ee) # Calculate the numerator and denominator of Rp - numerator = np.sum((Yt - Ye)**2 * step_Et) + numerator = np.sum((Yt - Ye) ** 2 * step_Et) denominator = np.sum(Ye**2 * step_Ee) + np.sum(Yt**2 * step_Et) # Calculate Rp @@ -84,17 +84,17 @@ def rp_factor(theoretical_curve, experimental_curve): def Y_function(I, energy_step): """ - Calculates Y for a given curve. + Calculates Y for a given curve. - Inputs: - - I: array of intensity values of the curve - - E: array of energy values of the curve + Inputs: + - I: array of intensity values of the curve + - E: array of energy values of the curve - Design: - Y = L / (1 + L^2 * vi^2) + Design: + Y = L / (1 + L^2 * vi^2) - where: - L = (I[i] - I[i-1]) / (energy_step * 0.5 * (I[i] + I[i-1])) + where: + L = (I[i] - I[i-1]) / (energy_step * 0.5 * (I[i] + I[i-1])) """ # [TODO] constants should be defined elsewhere @@ -108,13 +108,13 @@ def Y_function(I, energy_step): def energy_step(E): """ - Calculates the pairwise energy step. Returns an array. + Calculates the pairwise energy step. Returns an array. - Inputs: - - E: array of energy values of the curve + Inputs: + - E: array of energy values of the curve - Design: - step = E[i] - E[i-1] + Design: + step = E[i] - E[i-1] """ return E[1:] - E[:-1] diff --git a/tests/curves_helper.py b/tests/curves_helper.py index 55cb4f3..5ded4f6 100644 --- a/tests/curves_helper.py +++ b/tests/curves_helper.py @@ -1,78 +1,78 @@ def curve_A(): return [ - [70, 0.0037557780], - [74, 0.0013476560], - [78, 0.0010809600], - [82, 0.0009844219], - [86, 0.0009538341], - [90, 0.0011735780], - [94, 0.0020866750], - [98, 0.0040076710], - [102, 0.0049233240], - [106, 0.0051959950], - [110, 0.0044225620], - [114, 0.0020602790], - [118, 0.0005740963], - [122, 0.0005428890], - [126, 0.0002391057], - ] + [70, 0.0037557780], + [74, 0.0013476560], + [78, 0.0010809600], + [82, 0.0009844219], + [86, 0.0009538341], + [90, 0.0011735780], + [94, 0.0020866750], + [98, 0.0040076710], + [102, 0.0049233240], + [106, 0.0051959950], + [110, 0.0044225620], + [114, 0.0020602790], + [118, 0.0005740963], + [122, 0.0005428890], + [126, 0.0002391057], + ] def curve_B(): return [ - [70, 0.0101207100], - [74, 0.0100794400], - [78, 0.0056215730], - [82, 0.0028791560], - [86, 0.0013405740], - [90, 0.0005220333], - [94, 0.0004865836], - [98, 0.0002609056], - [102, 0.0000723265], - [106, 0.0002626837], - [110, 0.0016541830], - [114, 0.0041199060], - [118, 0.0062692520], - [122, 0.0079176060], - [126, 0.0138149000], - ] + [70, 0.0101207100], + [74, 0.0100794400], + [78, 0.0056215730], + [82, 0.0028791560], + [86, 0.0013405740], + [90, 0.0005220333], + [94, 0.0004865836], + [98, 0.0002609056], + [102, 0.0000723265], + [106, 0.0002626837], + [110, 0.0016541830], + [114, 0.0041199060], + [118, 0.0062692520], + [122, 0.0079176060], + [126, 0.0138149000], + ] def curve_C(): return [ - [70, 0.0101222800], - [74, 0.0100787500], - [78, 0.0056227330], - [82, 0.0028790560], - [86, 0.0013404120], - [90, 0.0005221121], - [94, 0.0004865952], - [98, 0.0002609447], - [102, 0.0000722322], - [106, 0.0002627033], - [110, 0.0016542040], - [114, 0.0041198820], - [118, 0.0062694210], - [122, 0.0079189570], - [126, 0.0138144800], - ] + [70, 0.0101222800], + [74, 0.0100787500], + [78, 0.0056227330], + [82, 0.0028790560], + [86, 0.0013404120], + [90, 0.0005221121], + [94, 0.0004865952], + [98, 0.0002609447], + [102, 0.0000722322], + [106, 0.0002627033], + [110, 0.0016542040], + [114, 0.0041198820], + [118, 0.0062694210], + [122, 0.0079189570], + [126, 0.0138144800], + ] def curve_A_smoothed(): return [ - [70, 0.00258495], - [74, 0.00182937], - [78, 0.00145542], - [82, 0.00135855], - [86, 0.00144012], - [90, 0.00176726], - [94, 0.00244699], - [98, 0.00335730], - [102, 0.00395957], - [106, 0.00404653], - [110, 0.00348618], - [114, 0.00238778], - [118, 0.00145559], - [122, 0.00101030], - [126, 0.00079836], - ] + [70, 0.00258495], + [74, 0.00182937], + [78, 0.00145542], + [82, 0.00135855], + [86, 0.00144012], + [90, 0.00176726], + [94, 0.00244699], + [98, 0.00335730], + [102, 0.00395957], + [106, 0.00404653], + [110, 0.00348618], + [114, 0.00238778], + [118, 0.00145559], + [122, 0.00101030], + [126, 0.00079836], + ] diff --git a/tests/test_preprocessing.py b/tests/test_preprocessing.py index f00927f..ad58b42 100644 --- a/tests/test_preprocessing.py +++ b/tests/test_preprocessing.py @@ -1,27 +1,29 @@ -import pytest import math -import numpy as np + import matplotlib.pyplot as plt +import numpy as np +import pytest + from cleedpy.preprocessing import lorentzian_smoothing, preprocessing_loop -from tests.curves_helper import curve_A, curve_B, curve_A_smoothed +from tests.curves_helper import curve_A, curve_A_smoothed, curve_B @pytest.mark.parametrize( "curve, vi, expected", [ (np.array(curve_A()), 4, np.array(curve_A_smoothed())), - ] + ], ) def test_lorentzian_smoothing(curve, vi, expected): L_curve = lorentzian_smoothing(curve, vi) - + x_values = curve[:, 0] y_values = curve[:, 1] - plt.plot(x_values, y_values, marker='o', color='g', label='Input: Initial curve') + plt.plot(x_values, y_values, marker="o", color="g", label="Input: Initial curve") x_values = L_curve[:, 0] y_values = L_curve[:, 1] - plt.plot(x_values, y_values, marker='x', color='r', label='Output: Smoothed curve') + plt.plot(x_values, y_values, marker="x", color="r", label="Output: Smoothed curve") y_min = 0 y_max = 0.0110 @@ -32,18 +34,20 @@ def test_lorentzian_smoothing(curve, vi, expected): plt.legend() plt.savefig("test_lorentzian_smoothing_output.png") - assert np.allclose( - expected, L_curve - ) + assert np.allclose(expected, L_curve) @pytest.mark.parametrize( "the_curve, exp_curve, shift, r_factor, vi, expected", [ ([curve_A(), curve_B()], [curve_A(), curve_B()], 1, "r2_factor", 4, 0), - ] + ], ) def test_preprocessing_loop(the_curve, exp_curve, shift, r_factor, vi, expected): assert math.isclose( - expected, preprocessing_loop(np.array(the_curve), np.array(exp_curve), shift, r_factor, vi), abs_tol=5 + expected, + preprocessing_loop( + np.array(the_curve), np.array(exp_curve), shift, r_factor, vi + ), + abs_tol=5, ) diff --git a/tests/test_rfactor.py b/tests/test_rfactor.py index b7a560a..4eb005f 100644 --- a/tests/test_rfactor.py +++ b/tests/test_rfactor.py @@ -1,6 +1,8 @@ -import pytest import math + import numpy as np +import pytest + from cleedpy.rfactor import r2_factor, rp_factor from tests.curves_helper import curve_A, curve_B, curve_C @@ -10,7 +12,7 @@ [ (curve_A(), curve_B(), 1.8688830931), (curve_B(), curve_C(), 0.0001429368), - ] + ], ) def test_r2_factor(the_curve, exp_curve, expected_r): assert math.isclose( @@ -23,7 +25,7 @@ def test_r2_factor(the_curve, exp_curve, expected_r): [ (curve_A(), curve_B(), 1.4898282448), (curve_B(), curve_C(), 0), - ] + ], ) def test_rp_factor(the_curve, exp_curve, expected_r): assert math.isclose(