Skip to content

Commit

Permalink
Remove old function not updated to current standards
Browse files Browse the repository at this point in the history
  • Loading branch information
@CPrescher
CPrescher committed Jul 23, 2024
1 parent 9ff17a2 commit f4eda33
Show file tree
Hide file tree
Showing 3 changed files with 126 additions and 137 deletions.
211 changes: 106 additions & 105 deletions glassure/optimization.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,6 @@

__all__ = [
"optimize_sq",
"optimize_density",
]


Expand Down Expand Up @@ -91,107 +90,109 @@ def optimize_sq(
return sq_pattern


def optimize_density(
data_pattern,
background_pattern,
initial_background_scaling,
composition,
initial_density,
background_min,
background_max,
density_min,
density_max,
iterations,
r_cutoff,
use_modification_fcn=False,
extrapolation_cutoff=None,
r_step=0.01,
fcn_callback=None,
):
"""
Performs an optimization of the background scaling and density using a figure of merit function defined by the low
r region in F(r) as described in Eggert et al. (2002) PRB, 65, 174105.
:param data_pattern: raw data pattern in Q space (A^-1)
:param background_pattern: raw background pattern in Q space (A^-1)
:param initial_background_scaling:
start value for the background scaling optimization
:param composition: composition of the sample as a dictionary with elements as keys and abundances as values
:param initial_density: start value for the density optimization in g/cm^3
:param background_min: minimum value for the background scaling
:param background_max: maximum value for the background scaling
:param density_min: minimum value for the density
:param density_max: maximum value for the density
:param iterations: number of iterations of S(Q) (see optimize_sq(...) prior to calculating chi2
:param r_cutoff: cutoff value below which there is no signal expected (below the first peak in g(r))
:param use_modification_fcn:
Whether to use the Lorch modification function during the Fourier transform.
Warning: When using the Lorch modification function, more iterations are needed
to get to the wanted result.
:param extrapolation_cutoff:
Determines up to which q value the S(Q) will be extrapolated to zero. The default
(None) will use the minimum q value plus 0.2 A^-1
:param r_step: Step size for the r-space for calculating f(r) during each iteration.
:param fcn_callback: Function which will be called after each iteration. The function should take four
arguments: iteration number, chi2, density, and background scaling. Additionally, the
function should return a boolean value, where True continues the optimization and False
will stop the optimization procedure
:return: (tuple) - density, density standard error, background scaling, background scaling standard error
"""
params = lmfit.Parameters()
params.add("density", value=initial_density, min=density_min, max=density_max)
params.add(
"background_scaling",
value=initial_background_scaling,
min=background_min,
max=background_max,
)

r = np.arange(0, r_cutoff + r_step / 2.0, r_step)

def optimization_fcn(params, extrapolation_max, r, r_cutoff, use_modification_fcn):
density = params["density"].value
atomic_density = convert_density_to_atoms_per_cubic_angstrom(
composition, density
)
background_pattern.scaling = params["background_scaling"].value

sq = calculate_sq(data_pattern - background_pattern, density, composition)
extrapolation_max = extrapolation_max or np.min(sq._x[0]) + 0.2
sq = extrapolate_to_zero_poly(sq, extrapolation_max)
sq_optimized = optimize_sq(
sq, r_cutoff, iterations, atomic_density, use_modification_fcn
)
fr = calculate_fr(sq_optimized, r=r, use_modification_fcn=use_modification_fcn)

min_r, min_fr = fr.data

output = (min_fr + 4 * np.pi * atomic_density * min_r) ** 2 * r_step

if fcn_callback is not None:
if not fcn_callback(
optimization_fcn.iteration,
np.sum(output),
density,
params["background_scaling"].value,
):
return None
optimization_fcn.iteration += 1
return output

optimization_fcn.iteration = 1

lmfit.minimize(
optimization_fcn,
params,
args=(extrapolation_cutoff, r, r_cutoff, use_modification_fcn),
)
lmfit.report_fit(params)

return (
params["density"].value,
params["density"].stderr,
params["background_scaling"].value,
params["background_scaling"].stderr,
)
# NEEDS Update to new configuration way of calculating density - this will also ensure, that every possible option is
# properly used in the density optimization (e.g. Klein Nishima correction or different S(Q) extrapolation methods)
# def optimize_density(
# data_pattern,
# background_pattern,
# initial_background_scaling,
# composition,
# initial_density,
# background_min,
# background_max,
# density_min,
# density_max,
# iterations,
# r_cutoff,
# use_modification_fcn=False,
# extrapolation_cutoff=None,
# r_step=0.01,
# fcn_callback=None,
# ):
# """
# Performs an optimization of the background scaling and density using a figure of merit function defined by the low
# r region in F(r) as described in Eggert et al. (2002) PRB, 65, 174105.

# :param data_pattern: raw data pattern in Q space (A^-1)
# :param background_pattern: raw background pattern in Q space (A^-1)
# :param initial_background_scaling:
# start value for the background scaling optimization
# :param composition: composition of the sample as a dictionary with elements as keys and abundances as values
# :param initial_density: start value for the density optimization in g/cm^3
# :param background_min: minimum value for the background scaling
# :param background_max: maximum value for the background scaling
# :param density_min: minimum value for the density
# :param density_max: maximum value for the density
# :param iterations: number of iterations of S(Q) (see optimize_sq(...) prior to calculating chi2
# :param r_cutoff: cutoff value below which there is no signal expected (below the first peak in g(r))
# :param use_modification_fcn:
# Whether to use the Lorch modification function during the Fourier transform.
# Warning: When using the Lorch modification function, more iterations are needed
# to get to the wanted result.
# :param extrapolation_cutoff:
# Determines up to which q value the S(Q) will be extrapolated to zero. The default
# (None) will use the minimum q value plus 0.2 A^-1
# :param r_step: Step size for the r-space for calculating f(r) during each iteration.
# :param fcn_callback: Function which will be called after each iteration. The function should take four
# arguments: iteration number, chi2, density, and background scaling. Additionally, the
# function should return a boolean value, where True continues the optimization and False
# will stop the optimization procedure

# :return: (tuple) - density, density standard error, background scaling, background scaling standard error
# """
# params = lmfit.Parameters()
# params.add("density", value=initial_density, min=density_min, max=density_max)
# params.add(
# "background_scaling",
# value=initial_background_scaling,
# min=background_min,
# max=background_max,
# )

# r = np.arange(0, r_cutoff + r_step / 2.0, r_step)

# def optimization_fcn(params, extrapolation_max, r, r_cutoff, use_modification_fcn):
# density = params["density"].value
# atomic_density = convert_density_to_atoms_per_cubic_angstrom(
# composition, density
# )
# background_pattern.scaling = params["background_scaling"].value

# sq = calculate_sq(data_pattern - background_pattern, density, composition)
# extrapolation_max = extrapolation_max or np.min(sq._x[0]) + 0.2
# sq = extrapolate_to_zero_poly(sq, extrapolation_max)
# sq_optimized = optimize_sq(
# sq, r_cutoff, iterations, atomic_density, use_modification_fcn
# )
# fr = calculate_fr(sq_optimized, r=r, use_modification_fcn=use_modification_fcn)

# min_r, min_fr = fr.data

# output = (min_fr + 4 * np.pi * atomic_density * min_r) ** 2 * r_step

# if fcn_callback is not None:
# if not fcn_callback(
# optimization_fcn.iteration,
# np.sum(output),
# density,
# params["background_scaling"].value,
# ):
# return None
# optimization_fcn.iteration += 1
# return output

# optimization_fcn.iteration = 1

# lmfit.minimize(
# optimization_fcn,
# params,
# args=(extrapolation_cutoff, r, r_cutoff, use_modification_fcn),
# )
# lmfit.report_fit(params)

# return (
# params["density"].value,
# params["density"].stderr,
# params["background_scaling"].value,
# params["background_scaling"].stderr,
# )
50 changes: 19 additions & 31 deletions glassure/transform.py
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
from typing import Optional

import numpy as np
from scipy.integrate import simpson

from .pattern import Pattern
from .methods import FourierTransformMethod
Expand Down Expand Up @@ -65,8 +66,8 @@ def calculate_fr(
fr = (
2.0
/ np.pi
* np.trapz(
modification * q * (sq - 1) * np.array(np.sin(np.outer(q.T, r))).T, q
* simpson(
modification * q * (sq - 1) * np.array(np.sin(np.outer(q.T, r))).T, x=q
)
)
elif method == "fft" or method == FourierTransformMethod.FFT:
Expand All @@ -92,45 +93,35 @@ def calculate_fr(
return Pattern(r, fr)


def calculate_sq_from_fr(
fr_pattern: Pattern, q: np.ndarray, method: str = "integral"
) -> Pattern:
def calculate_sq_from_fr(fr_pattern: Pattern, q: np.ndarray) -> Pattern:
"""
Calculates S(Q) from an F(r) pattern for given q values.
:param fr_pattern: input F(r) pattern
:param q: numpy array giving the q-values for which S(q) will be calculated,
:param method: determines the method use for calculating fr, possible values are:
- 'integral' solves the Fourier integral, by calculating the integral
- 'fft' solves the Fourier integral by using fast fourier transformation
:return: F(r) pattern
"""
r, fr = fr_pattern.data
sq = simpson(fr * np.array(np.sin(np.outer(r.T, q))).T, x=r) / q + 1

if method == "integral":
sq = np.trapz(fr * np.array(np.sin(np.outer(r.T, q))).T, r) / q + 1

elif method == "fft":
q_step = q[1] - q[0]
r_step = r[1] - r[0]

n_out = int(np.pi / (r_step * q_step))
# there should be a way of doing this with fft, but the current implementation is not working
# correctly (it will not reproduce the S(q) pattern after transforming to fr and back accurately)
# elif method == "fft":
# q_step = q[1] - q[0]
# r_step = r[1] - r[0]

r_max_for_ifft = 2 * n_out * r_step
ifft_x_step = 2 * np.pi / r_max_for_ifft
ifft_x = np.arange(n_out) * ifft_x_step
# n_out = int(np.pi / (r_step * q_step))

y_for_ifft = np.concatenate((fr, np.zeros(2 * n_out - len(r))))
ifft_result = np.fft.ifft(y_for_ifft) * r_max_for_ifft
ifft_imag = np.imag(ifft_result)[:n_out]
# r_max_for_ifft = 2 * n_out * r_step
# ifft_x_step = 2 * np.pi / r_max_for_ifft
# ifft_x = np.arange(n_out) * ifft_x_step

sq = np.interp(q, ifft_x, ifft_imag) / q + 1
else:
raise NotImplementedError(
"{} is not an allowed method for calculate_sq_from_fr".format(method)
)
# y_for_ifft = np.concatenate((fr, np.zeros(2 * n_out - len(r))))
# ifft_result = np.fft.ifft(y_for_ifft) * r_max_for_ifft
# ifft_imag = np.imag(ifft_result)[:n_out]

# sq = np.interp(q, ifft_x, ifft_imag) / q + 1
return Pattern(q, sq)


Expand All @@ -146,9 +137,6 @@ def calculate_sq_from_gr(
:param gr_pattern: g(r) pattern
:param q: numpy array of q values for which S(Q) should be calculated
:param atomic_density: number_density in atoms/A^3
:param method: determines the method used for calculating fr, possible values are:
- 'integral' solves the Fourier integral, by calculating the integral
- 'fft' solves the Fourier integral by using fast fourier transformation
:return: S(Q) pattern
"""
Expand All @@ -158,7 +146,7 @@ def calculate_sq_from_gr(
if np.isnan(gr[0]):
gr[0] = 0
fr_pattern = Pattern(r, (gr - 1) * (4.0 * np.pi * r * atomic_density))
return calculate_sq_from_fr(fr_pattern, q, method)
return calculate_sq_from_fr(fr_pattern, q)


def calculate_gr(fr_pattern: Pattern, atomic_density: float) -> Pattern:
Expand Down
2 changes: 1 addition & 1 deletion tests/test_transform.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -114,4 +114,4 @@ def test_calculate_sq_from_fr(self):
fr = calculate_fr(sq)

sq_from_fr = calculate_sq_from_fr(fr, sq.x)
assert np.allclose(sq.y, sq_from_fr.y, atol=0.3)
assert np.allclose(sq.y, sq_from_fr.y, atol=0.3)

0 comments on commit f4eda33

Please sign in to comment.