Skip to content

Commit

Permalink
[pre-commit.ci] auto fixes from pre-commit.com hooks
Browse files Browse the repository at this point in the history 
for more information, see https://pre-commit.ci
  • Loading branch information
@pre-commit-ci
pre-commit-ci[bot] committed Feb 2, 2024
1 parent 2a7ac89 commit e2fa2ce
Show file tree
Hide file tree
Showing 5 changed files with 177 additions and 159 deletions.
102 changes: 57 additions & 45 deletions cleedpy/preprocessing.py
View file
Original file line number Diff line number Diff line change
@@ -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]
Expand All @@ -25,40 +26,40 @@ 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)))


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)):
Expand All @@ -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
Expand All @@ -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
70 changes: 35 additions & 35 deletions cleedpy/rfactor.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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]
Loading

0 comments on commit e2fa2ce

Please sign in to comment.