Very quick and dirty implementation of excluded regions

Author: AndrewSazonov

src/easydiffraction/analysis/minimization.py

@@ -121,6 +121,21 @@ def _residual_function(self,
         # Sync parameters back to objects
         self.minimizer._sync_result_to_parameters(parameters, engine_params)
 
+        # Update the excluded points in experiments
+        # TODO: This should not be handled here every time
+        #  This implementation is just very quick and dirty
+        for experiment in experiments._items.values():
+            experiment.datastore.pattern.excluded = np.full(experiment.datastore.pattern.x.shape,
+                                                            fill_value=False,
+                                                            dtype=bool)  # Reset excluded points
+            excluded_regions = experiment.excluded_regions._items  # List of excluded regions
+            if excluded_regions:  # If there are any excluded regions
+                for idx, point in enumerate(experiment.datastore.pattern.x):  # Set excluded points
+                    for region in excluded_regions.values():
+                        if region.minimum.value <= point <= region.maximum.value:
+                            experiment.datastore.pattern.excluded[idx] = True
+                            break
+
         # Prepare weights for joint fitting
         num_expts: int = len(experiments.ids)
         if weights is None:
@@ -145,7 +160,14 @@ def _residual_function(self,
                                                              called_by_minimizer=True)  # True False
             y_meas: np.ndarray = experiment.datastore.pattern.meas
             y_meas_su: np.ndarray = experiment.datastore.pattern.meas_su
-            diff: np.ndarray = (y_meas - y_calc) / y_meas_su
+            excluded: np.ndarray = experiment.datastore.pattern.excluded
+            # TODO: Excluded points must be handled differently.
+            #  They should not contribute to the residuals.
+            diff = np.where(
+                excluded,
+                0.0,  # Excluded points contribute zero to the residuals
+                (y_meas - y_calc) / y_meas_su  # Normalized residuals
+            )
             diff *= np.sqrt(weight)  # Residuals are squared before going into reduced chi-squared
             residuals.extend(diff)
 

src/easydiffraction/experiments/collections/datastore.py

@@ -16,6 +16,7 @@ def __init__(self, experiment: Experiment) -> None:
         self.meas: Optional[np.ndarray] = None
         self.meas_su: Optional[np.ndarray] = None
         self.bkg: Optional[np.ndarray] = None
+        self.excluded: Optional[np.ndarray] = None  # Excluded points
         self._calc: Optional[np.ndarray] = None  # Cached calculated intensities
 
     @property
@@ -33,6 +34,7 @@ class PowderPattern(Pattern):
     """
     Specialized pattern for powder diffraction (can be extended in the future).
     """
+    # TODO: Check if this class is needed or if it can be merged with Pattern
     def __init__(self, experiment: Experiment) -> None:
         super().__init__(experiment)
         # Additional powder-specific initialization if needed
@@ -49,7 +51,7 @@ def __init__(self, sample_form: str, experiment: Experiment) -> None:
         if sample_form == "powder":
             self.pattern: Pattern = PowderPattern(experiment)
         elif sample_form == "single_crystal":
-            self.pattern: Pattern = Pattern(experiment)
+            self.pattern: Pattern = Pattern(experiment)  # TODO: Find better name for single crystal pattern
         else:
             raise ValueError(f"Unknown sample form '{sample_form}'")
 

src/easydiffraction/experiments/collections/excluded_regions.py

@@ -0,0 +1,55 @@
+from typing import Type
+
+from easydiffraction.core.objects import (
+    Parameter,
+    Descriptor,
+    Component,
+    Collection
+)
+
+
+class ExcludedRegion(Component):
+    @property
+    def category_key(self) -> str:
+        return "excluded_region"
+
+    @property
+    def cif_category_key(self) -> str:
+        return "excluded_region"
+
+    def __init__(self,
+                 minimum: float,
+                 maximum: float):
+        super().__init__()
+
+        self.minimum = Descriptor(
+            value=minimum,
+            name="minimum",
+            cif_name="minimum"
+        )
+        self.maximum = Parameter(
+            value=maximum,
+            name="maximum",
+            cif_name="maximum"
+        )
+
+        # Select which of the input parameters is used for the
+        # as ID for the whole object
+        self._entry_id = f'{minimum}-{maximum}'
+
+        # Lock further attribute additions to prevent
+        # accidental modifications by users
+        self._locked = True
+
+
+class ExcludedRegions(Collection):
+    """
+    Collection of LinkedPhase instances.
+    """
+    @property
+    def _type(self) -> str:
+        return "category"  # datablock or category
+
+    @property
+    def _child_class(self) -> Type[ExcludedRegion]:
+        return ExcludedRegion

src/easydiffraction/experiments/experiment.py

@@ -15,6 +15,7 @@
 
 from easydiffraction.experiments.collections.linked_phases import LinkedPhases
 from easydiffraction.experiments.collections.background import BackgroundFactory
+from easydiffraction.experiments.collections.excluded_regions import ExcludedRegions
 from easydiffraction.experiments.collections.datastore import DatastoreFactory
 
 from easydiffraction.utils.formatting import paragraph, warning
@@ -129,6 +130,10 @@ def as_cif(self, max_points: Optional[int] = None) -> str:
         if hasattr(self, "background") and self.background:
             cif_lines += ["", self.background.as_cif()]
 
+        # Excluded regions
+        if hasattr(self, "excluded_regions") and self.excluded_regions:
+            cif_lines += ["", self.excluded_regions.as_cif()]
+
         # Measured data
         if hasattr(self, "datastore") and hasattr(self.datastore, "pattern"):
             cif_lines.append("")
@@ -191,7 +196,8 @@ def __init__(self,
             beam_mode=self.type.beam_mode.value,
             profile_type=self._peak_profile_type)
 
-        self.linked_phases = LinkedPhases()
+        self.linked_phases: LinkedPhases = LinkedPhases()
+        self.excluded_regions: ExcludedRegions = ExcludedRegions()
 
     @abstractmethod
     def _load_ascii_data_to_experiment(self, data_path: str) -> None:

src/easydiffraction/plotting/plotting.py

@@ -140,6 +140,11 @@ def plot_meas(self,
                                         x_min=x_min,
                                         x_max=x_max)
 
+        # Exclude points based on the pattern's excluded mask
+        excluded = pattern.excluded
+        x = x[~excluded]
+        y_meas = y_meas[~excluded]
+
         y_series = [y_meas]
         y_labels = ['meas']
 
@@ -176,6 +181,11 @@ def plot_calc(self,
                                         x_min=x_min,
                                         x_max=x_max)
 
+        # Exclude points based on the pattern's excluded mask
+        excluded = pattern.excluded
+        x = x[~excluded]
+        y_calc = y_calc[~excluded]
+
         y_series = [y_calc]
         y_labels = ['calc']
 
@@ -220,6 +230,12 @@ def plot_meas_vs_calc(self,
                                         x_min=x_min,
                                         x_max=x_max)
 
+        # Exclude points based on the pattern's excluded mask
+        excluded = pattern.excluded
+        x = x[~excluded]
+        y_meas = y_meas[~excluded]
+        y_calc = y_calc[~excluded]
+
         y_series = [y_meas, y_calc]
         y_labels = ['meas', 'calc']
 

tests/functional_tests/fitting/test_powder-diffraction_multiphase.py

@@ -84,7 +84,7 @@ def test_single_fit_neutron_pd_tof_mcstas_lbco_si() -> None:
 
     # Compare fit quality
     assert_almost_equal(project.analysis.fit_results.reduced_chi_square,
-                        desired=2.87,
+                        desired=1.79,  # 2.87
                         decimal=1)
 
 

tutorials/cryst-struct_pd-neut-tof_multphase-LBCO-Si_McStas.py

@@ -211,6 +211,12 @@
 # %%
 project.analysis.current_minimizer = 'lmfit (leastsq)'
 
+# %% [markdown]
+# ### Set Excluded Regions
+
+# %%
+project.experiments['mcstas'].excluded_regions.add(minimum=108000, maximum=200000)
+
 # %% [markdown]
 # ### Set Fitting Parameters
 #