refactor: more clean up of test optimization

Author: CPrescher

glassure/calc.py

@@ -69,7 +69,7 @@ def create_calculate_pdf_configs(
 
 def calculate_pdf(
     data_config: DataConfig, calculation_config: CalculationConfig
-) -> Pattern:
+) -> Result:
     """
     Process the input configuration and return the result.
     """

glassure/optimization.py

@@ -8,6 +8,7 @@
 
 from . import Pattern
 from .transform import calculate_fr, calculate_gr, calculate_sq_from_fr
+from .utility import convert_density_to_atoms_per_cubic_angstrom
 
 __all__ = [
     "optimize_sq",
@@ -20,6 +21,7 @@ def optimize_sq(
     r_cutoff: float,
     iterations: int,
     atomic_density: float,
+    r_step: float = 0.01,
     use_modification_fcn: bool = False,
     attenuation_factor: float = 1,
     fcn_callback=None,
@@ -39,6 +41,9 @@ def optimize_sq(
         number of back and forward transforms
     :param atomic_density:
         density in atoms/A^3
+    :param r_step:
+        step size for the r-axis, default is 0.01. Use smaller values for better accuracy (especially if needed for
+         fft)
     :param use_modification_fcn:
         Whether to use the Lorch modification function during the Fourier transform.
         Warning: When using the Lorch modification function, usually more iterations are needed to get to the
@@ -59,7 +64,7 @@ def optimize_sq(
     :return:
         optimized S(Q) pattern
     """
-    r = np.arange(0, r_cutoff, 0.01)
+    r = np.arange(0, r_cutoff, r_step)
     sq_pattern = deepcopy(sq_pattern)
     for iteration in range(iterations):
         fr_pattern = calculate_fr(
@@ -71,13 +76,14 @@ def optimize_sq(
         delta_fr = fr_int + 4 * np.pi * r * atomic_density
 
         if fourier_transform_method == "fft":
-            sq_trans_fft = calculate_sq_from_fr(
-                Pattern(r, delta_fr), sq_pattern.x, method="fft"
-            ) - 1 
+            sq_trans_fft = (
+                calculate_sq_from_fr(Pattern(r, delta_fr), sq_pattern.x, method="fft")
+                - 1
+            )
             iq = sq_trans_fft.y
         else:
             in_integral = np.array(np.sin(np.outer(q.T, r))) * delta_fr
-            iq = simpson(in_integral, r) / q
+            iq = simpson(in_integral, x=r) / q
 
         sq_pattern = sq_pattern * (1 - iq / attenuation_factor)
 
@@ -177,6 +183,10 @@ def fcn(params):
         if type == "gr":
             r, gr = result.gr.limit(*min_range).data
             residual = np.trapz(gr**2, r)
+        if type == "fr":
+            atomic_density = optim_config.sample.atomic_density
+            r, fr = result.fr.limit(*min_range).data
+            residual = np.trapz(fr + 4 * np.pi * r * atomic_density)
         elif type == "sq":
             q, sq = result.sq.limit(*min_range).data
             sq_ref = reference_result.sq.limit(*min_range).y

tests/test_normalization.py

@@ -11,7 +11,6 @@
 )
 from glassure.normalization import normalize, normalize_fit, normalize_fit_lmfit
 from glassure.scattering_factors import calculate_coherent_scattering_factor
-from glassure.transform import calculate_sq
 
 from . import unittest_data_path
 

tests/test_optimization.py

@@ -1,7 +1,7 @@
 # -*- coding: utf-8 -*-
 
 import os
-import unittest
+import pytest
 import numpy as np
 
 from glassure import Pattern, convert_density_to_atoms_per_cubic_angstrom
@@ -14,7 +14,7 @@
     calculate_s0,
 )
 from glassure.transform import calculate_sq, calculate_fr
-from glassure.normalization import normalize_fit
+from glassure.normalization import normalize_fit, normalize
 from glassure.configuration import OptimizeConfig
 from glassure.optimization import optimize_sq, optimize_density
 from glassure.calc import calculate_pdf, create_calculate_pdf_configs
@@ -32,58 +32,98 @@
 background_path_SiO2 = os.path.join(unittest_data_path, "SiO2_bkg.xy")
 
 
-def test_optimize_sq():
-    data = Pattern.from_file(data_path_alloy)
-    background = Pattern.from_file(background_path_alloy)
-    composition = {"Fe": 0.81, "S": 0.19}
-    density = 7.9
-    atomic_density = convert_density_to_atoms_per_cubic_angstrom(composition, density)
-    f_squared_mean = calculate_f_squared_mean(composition, data.x)
-    f_mean_squared = calculate_f_mean_squared(composition, data.x)
-    incoherent_scattering = calculate_incoherent_scattering(composition, data.x)
-    background_scaling = 0.97
+@pytest.fixture
+def data_path():
+    """Path to the test data file."""
+    return os.path.join(unittest_data_path, "SiO2.xy")
 
-    sample_pattern = data - background_scaling * background
-
-    sq = calculate_sq(sample_pattern, f_squared_mean, f_mean_squared)
-    sq = extrapolate_to_zero_poly(sq, np.min(sq.x) + 0.3)
-    sq_optimized = optimize_sq(sq, 1.6, 5, atomic_density)
-    assert not np.allclose(sq.y, sq_optimized.y)
 
+@pytest.fixture
+def bkg_path():
+    """Path to the background data file."""
+    return os.path.join(unittest_data_path, "SiO2_bkg.xy")
 
-def test_optimize_sq_fft():
-    data = Pattern.from_file(data_path_SiO2)
-    background = Pattern.from_file(background_path_SiO2)
-    composition = {"Si": 1, "O": 2}
-    density = 2.2
-    atomic_density = convert_density_to_atoms_per_cubic_angstrom(composition, density)
-    background_scaling = 1.0
 
-    sample_pattern = data - background_scaling * background
-    sample_pattern = sample_pattern.limit(0, 17).rebin(0.05)
-    q = sample_pattern.x
-    f_squared_mean = calculate_f_squared_mean(composition, q)
-    f_mean_squared = calculate_f_mean_squared(composition, q)
-    incoherent_scattering = calculate_incoherent_scattering(composition, q)
+@pytest.fixture
+def sample(data_path, bkg_path):
+    """Create a sample pattern by subtracting background from data."""
+    data = Pattern.from_file(data_path)
+    bkg = Pattern.from_file(bkg_path)
+    sample = data - bkg
+    return sample.limit(1, 17)
 
-    _, norm_pattern = normalize_fit(
-        sample_pattern, f_squared_mean, incoherent_scattering
-    )
 
-    sq = calculate_sq(norm_pattern, f_squared_mean, f_mean_squared)
+@pytest.fixture
+def sq(normalized_pattern, f_squared_mean, f_mean_squared, composition):
+    """Create a sq pattern for testing."""
+    sq = calculate_sq(normalized_pattern, f_squared_mean, f_mean_squared)
     sq = extrapolate_to_zero_linear(sq, y0=calculate_s0(composition))
+    sq = sq.rebin(0.05)
     sq.x[0] = 1e-10
-    iterations = 5 
+    return sq
+
+
+@pytest.fixture
+def composition():
+    """Sample composition for testing."""
+    return {"Si": 1, "O": 2}
+
+
+@pytest.fixture
+def density():
+    """Sample density for testing."""
+    return 2.2
+
+
+@pytest.fixture
+def atomic_density(composition, density):
+    """Calculate atomic density from composition and density."""
+    return convert_density_to_atoms_per_cubic_angstrom(composition, density)
+
+
+@pytest.fixture
+def f_squared_mean(composition, sample):
+    """Calculate f squared mean for the sample."""
+    return calculate_f_squared_mean(composition, sample.x)
 
-    fr = calculate_fr(sq, method="fft")
 
+@pytest.fixture
+def f_mean_squared(composition, sample):
+    """Calculate f mean squared for the sample."""
+    return calculate_f_mean_squared(composition, sample.x)
+
+
+@pytest.fixture
+def incoherent_scattering(composition, sample):
+    """Calculate incoherent scattering for the sample."""
+    return calculate_incoherent_scattering(composition, sample.x)
+
+
+@pytest.fixture
+def normalized_pattern(sample, f_squared_mean, incoherent_scattering):
+    """Create normalized pattern for testing."""
+    _, normalized = normalize_fit(
+        sample, f_squared_mean, incoherent_scattering, q_cutoff=10
+    )
+    return normalized
+
+
+def test_optimize_sq(sq, atomic_density):
+    sq_optimized = optimize_sq(sq, 1.4, 5, atomic_density)
+    assert not np.allclose(sq.y, sq_optimized.y)
+
+
+def test_optimize_sq_fft(sq, atomic_density):
+    iterations = 5
+    r_step = 0.001 # need high value to be accurate for fft
     sq_optimized = optimize_sq(
-        sq, 1.3, iterations, atomic_density, fourier_transform_method="integral"
+        sq, 1.3, iterations, atomic_density, fourier_transform_method="integral", r_step=r_step
     )
     fr_optimized = calculate_fr(sq_optimized, method="fft")
 
+
     sq_optimized_fft = optimize_sq(
-        sq, 1.3, iterations, atomic_density, fourier_transform_method="fft"
+        sq, 1.3, iterations, atomic_density, fourier_transform_method="fft", r_step=r_step
     )
     fr_optimized_fft = calculate_fr(sq_optimized_fft, method="fft")