Add Support for Excluded Regions in Fitting and Visualization (#58)

* Very quick and dirty implementation of excluded regions

* Adds parent handling and item addition callback to Collection

* Refactors excluded points handling in diffraction minimization

Moves the logic for updating excluded points from minimization to ExcludedRegions, ensuring better modularity.

* Extends unit test with excluded regions

* Refactors excluded region logic in experiments

Optimizes the process of updating the pattern's excluded points
by avoiding unnecessary resets with each new excluded region.

Initializes pattern's excluded points with default values on data load.

* Excludes specified points from calculations for improve accuracy

* Refines excluded region handling in experiments

* Refines data plotting and analysis settings

* Formats function signatures for readability

* Removes unused assignment in test_fit_with_params

* Adds 'excluded' flag to experiment mock data

* Improves handling of excluded data points during minimization

* Enhances exclusion handling in Plotter class

* Fixes empty background and excluded regions check

* Adds CIF display for experiment data

* Updates the multiphase tutorial

* Simplifies exclusion handling in data processing

* Renames tutorial files and updates documentation

* Updates tutorials and adds display for excluded regions

* Refines tutorial markdown headers for consistency

* Adjusts data preprocessing for stability in analysis

Rounds x-values to 4 decimal places to align data size
Replaces small uncertainty values with 1.0 to prevent
failures in minimization algorithms

Author: AndrewSazonov

docs/mkdocs.yml

@@ -67,11 +67,11 @@ nav:
           - LBCO quick: tutorials/quick_single-fit_pd-neut-cwl_LBCO-HRPT.ipynb
           - LBCO basic: tutorials/basic_single-fit_pd-neut-cwl_LBCO-HRPT.ipynb
           - PbSO4 advanced: tutorials/advanced_joint-fit_pd-neut-xray-cwl_PbSO4.ipynb
-      - Structure Refinement:
+      - Standard Diffraction:
           - Co2SiO4 pd-neut-cwl: tutorials/cryst-struct_pd-neut-cwl_CoSiO4-D20.ipynb
           - HS pd-neut-cwl: tutorials/cryst-struct_pd-neut-cwl_HS-HRPT.ipynb
           - Si pd-neut-tof: tutorials/cryst-struct_pd-neut-tof_Si-SEPD.ipynb
-          - NCAF pd-neut-tof: tutorials/cryst-struct_pd-neut-tof_NCAF-WISH.ipynb
+          - NCAF pd-neut-tof: tutorials/cryst-struct_pd-neut-tof_multidata_NCAF-WISH.ipynb
           - LBCO+Si McStas: tutorials/cryst-struct_pd-neut-tof_multphase-LBCO-Si_McStas.ipynb
       - Pair Distribution Function:
           - Ni pd-neut-cwl: tutorials/pdf_pd-neut-cwl_Ni.ipynb

docs/tutorials/index.md

@@ -10,7 +10,7 @@ as self-contained, step-by-step **guides** to help users grasp the workflow of d
 analysis using EasyDiffraction.
 
 Instructions on how to run the tutorials are provided in the
-[:material-cog-box: Installation & Setup](../installation-and-setup/index.md#running-tutorials)
+[:material-cog-box: Installation & Setup](../installation-and-setup/index.md#how-to-run-tutorials)
 section of the documentation.
 
 The tutorials are organized into the following categories.
@@ -46,9 +46,14 @@ The tutorials are organized into the following categories.
 - [Si `pd-neut-tof`](cryst-struct_pd-neut-tof_Si-SEPD.ipynb) –
   Demonstrates a Rietveld refinement of the Si crystal structure using
   time-of-flight neutron powder diffraction data from SEPD at Argonne.
-- [NCAF `pd-neut-tof`](cryst-struct_pd-neut-tof_NCAF-WISH.ipynb) –
+- [NCAF `pd-neut-tof`](cryst-struct_pd-neut-tof_multidata_NCAF-WISH.ipynb) –
   Demonstrates a Rietveld refinement of the Na2Ca3Al2F14 crystal structure 
-  using time-of-flight neutron powder diffraction data from WISH at ISIS.
+  using two time-of-flight neutron powder diffraction datasets (from two 
+  detector banks) of the WISH instrument at ISIS.
+- [LBCO+Si McStas](cryst-struct_pd-neut-tof_multiphase-LBCO-Si_McStas.ipynb) –
+  Demonstrates a Rietveld refinement of the La0.5Ba0.5CoO3 crystal structure
+  with a small amount of Si impurity as a secondary phase using time-of-flight 
+  neutron powder diffraction data simulated with McStas.
 
 ## Pair Distribution Function (PDF)
 

src/easydiffraction/analysis/calculation.py

@@ -31,7 +31,9 @@ def set_calculator(self, engine: str) -> None:
         """
         self._calculator = self.calculator_factory.create_calculator(engine)
 
-    def calculate_structure_factors(self, sample_models: SampleModels, experiments: Experiments) -> Optional[List[Any]]:
+    def calculate_structure_factors(self,
+                                    sample_models: SampleModels,
+                                    experiments: Experiments) -> Optional[List[Any]]:
         """
         Calculate HKL intensities (structure factors) for sample models and experiments.
 
@@ -44,15 +46,17 @@ def calculate_structure_factors(self, sample_models: SampleModels, experiments:
         """
         return self._calculator.calculate_structure_factors(sample_models, experiments)
 
-    def calculate_pattern(self, sample_models: SampleModels, experiment: Experiment) -> np.ndarray:
+    def calculate_pattern(self,
+                          sample_models: SampleModels,
+                          experiment: Experiment) -> np.ndarray:
         """
-        Generate diffraction pattern based on sample models and experiment.
+        Calculate diffraction pattern based on sample models and experiment.
 
         Args:
             sample_models: Collection of sample models.
             experiment: A single experiment object.
 
         Returns:
-            Diffraction pattern generated by the backend calculator.
+            Diffraction pattern calculated by the backend calculator.
         """
         return self._calculator.calculate_pattern(sample_models, experiment)

src/easydiffraction/analysis/calculators/calculator_cryspy.py

@@ -127,6 +127,7 @@ def _recreate_cryspy_dict(self,
 
         cryspy_model_id = f'crystal_{sample_model.name}'
         cryspy_model_dict = cryspy_dict[cryspy_model_id]
+
         # Cell
         cryspy_cell = cryspy_model_dict['unit_cell_parameters']
         cryspy_cell[0] = sample_model.cell.length_a.value
@@ -135,6 +136,7 @@ def _recreate_cryspy_dict(self,
         cryspy_cell[3] = np.deg2rad(sample_model.cell.angle_alpha.value)
         cryspy_cell[4] = np.deg2rad(sample_model.cell.angle_beta.value)
         cryspy_cell[5] = np.deg2rad(sample_model.cell.angle_gamma.value)
+
         # Atomic coordinates
         cryspy_xyz = cryspy_model_dict['atom_fract_xyz']
         for idx, atom_site in enumerate(sample_model.atom_sites):

src/easydiffraction/analysis/minimization.py

@@ -78,7 +78,11 @@ def _process_fit_results(self,
         f_obs, f_calc = None, None
 
         if self.results:
-            self.results.display_results(y_obs=y_obs, y_calc=y_calc, y_err=y_err, f_obs=f_obs, f_calc=f_calc)
+            self.results.display_results(y_obs=y_obs,
+                                         y_calc=y_calc,
+                                         y_err=y_err,
+                                         f_obs=f_obs,
+                                         f_calc=f_calc)
 
     def _collect_free_parameters(self,
                                  sample_models: SampleModels,
@@ -140,13 +144,19 @@ def _residual_function(self,
         residuals: List[float] = []
 
         for (expt_id, experiment), weight in zip(experiments._items.items(), _weights):
+
+            # Calculate the difference between measured and calculated patterns
             y_calc: np.ndarray = calculator.calculate_pattern(sample_models,
                                                              experiment,
-                                                             called_by_minimizer=True)  # True False
+                                                             called_by_minimizer=True)
             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
-            diff *= np.sqrt(weight)  # Residuals are squared before going into reduced chi-squared
+            diff = ((y_meas - y_calc) / y_meas_su)
+
+            # Residuals are squared before going into reduced chi-squared
+            diff *= np.sqrt(weight)
+
+            # Append the residuals for this experiment
             residuals.extend(diff)
 
         return self.minimizer.tracker.track(np.array(residuals), parameters)

src/easydiffraction/analysis/reliability_factors.py

@@ -4,7 +4,8 @@
 from easydiffraction.experiments.experiments import Experiments
 from easydiffraction.analysis.calculators.calculator_base import CalculatorBase
 
-def calculate_r_factor(y_obs: np.ndarray, y_calc: np.ndarray) -> float:
+def calculate_r_factor(y_obs: np.ndarray,
+                       y_calc: np.ndarray) -> float:
     """
     Calculate the R-factor (reliability factor) between observed and calculated data.
     
@@ -22,7 +23,9 @@ def calculate_r_factor(y_obs: np.ndarray, y_calc: np.ndarray) -> float:
     return numerator / denominator if denominator != 0 else np.nan
 
 
-def calculate_weighted_r_factor(y_obs: np.ndarray, y_calc: np.ndarray, weights: np.ndarray) -> float:
+def calculate_weighted_r_factor(y_obs: np.ndarray,
+                                y_calc: np.ndarray,
+                                weights: np.ndarray) -> float:
     """
     Calculate the weighted R-factor between observed and calculated data.
     
@@ -42,7 +45,8 @@ def calculate_weighted_r_factor(y_obs: np.ndarray, y_calc: np.ndarray, weights:
     return np.sqrt(numerator / denominator) if denominator != 0 else np.nan
 
 
-def calculate_rb_factor(y_obs: np.ndarray, y_calc: np.ndarray) -> float:
+def calculate_rb_factor(y_obs: np.ndarray,
+                        y_calc: np.ndarray) -> float:
     """
     Calculate the Bragg R-factor between observed and calculated data.
     
@@ -60,7 +64,8 @@ def calculate_rb_factor(y_obs: np.ndarray, y_calc: np.ndarray) -> float:
     return numerator / denominator if denominator != 0 else np.nan
 
 
-def calculate_r_factor_squared(y_obs: np.ndarray, y_calc: np.ndarray) -> float:
+def calculate_r_factor_squared(y_obs: np.ndarray,
+                               y_calc: np.ndarray) -> float:
     """
     Calculate the R-factor squared between observed and calculated data.
     
@@ -78,7 +83,8 @@ def calculate_r_factor_squared(y_obs: np.ndarray, y_calc: np.ndarray) -> float:
     return np.sqrt(numerator / denominator) if denominator != 0 else np.nan
 
 
-def calculate_reduced_chi_square(residuals: np.ndarray, num_parameters: int) -> float:
+def calculate_reduced_chi_square(residuals: np.ndarray,
+                                 num_parameters: int) -> float:
     """
     Calculate the reduced chi-square statistic.
     
@@ -99,7 +105,9 @@ def calculate_reduced_chi_square(residuals: np.ndarray, num_parameters: int) ->
         return np.nan
 
 
-def get_reliability_inputs(sample_models: SampleModels, experiments: Experiments, calculator: CalculatorBase) -> Tuple[np.ndarray, np.ndarray, Optional[np.ndarray]]:
+def get_reliability_inputs(sample_models: SampleModels,
+                           experiments: Experiments,
+                           calculator: CalculatorBase) -> Tuple[np.ndarray, np.ndarray, Optional[np.ndarray]]:
     """
     Collect observed and calculated data points for reliability calculations.
     
@@ -120,9 +128,13 @@ def get_reliability_inputs(sample_models: SampleModels, experiments: Experiments
         y_meas_su = experiment.datastore.pattern.meas_su
 
         if y_meas is not None and y_calc is not None:
+            # If standard uncertainty is not provided, use ones
+            if y_meas_su is None:
+                y_meas_su = np.ones_like(y_meas)
+
             y_obs_all.extend(y_meas)
             y_calc_all.extend(y_calc)
-            y_err_all.extend(y_meas_su if y_meas_su is not None else np.ones_like(y_meas))
+            y_err_all.extend(y_meas_su)
 
     return (
         np.array(y_obs_all),

src/easydiffraction/core/objects.py

@@ -341,7 +341,8 @@ class Collection(ABC):
     def _child_class(self):
         return None
 
-    def __init__(self):
+    def __init__(self, parent=None):
+        self._parent = parent  # Parent datablock
         self._datablock_id = None  # Parent datablock name to be set by the parent
         self._items = {}
 
@@ -373,6 +374,10 @@ def add(self, *args, **kwargs):
         child_obj.entry_id = child_obj.entry_id  # Forcing the entry_id to be reset to update its child parameters
         self._items[child_obj._entry_id] = child_obj
 
+        # Call on_item_added if it exists, i.e. defined in the derived class
+        if hasattr(self, "on_item_added"):
+            self.on_item_added(child_obj)
+
     def get_all_params(self):
         params = []
         for item in self._items.values():

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  # Flags for 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,12 +51,14 @@ 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}'")
 
     def load_measured_data(self, file_path: str) -> None:
         """Load measured data from an ASCII file."""
+        # TODO: Check if this method is used...
+        #  Looks like _load_ascii_data_to_experiment from experiments.py is used instead
         print(f"Loading measured data for {self.sample_form} diffraction from {file_path}")
 
         try:
@@ -73,6 +77,7 @@ def load_measured_data(self, file_path: str) -> None:
         self.pattern.x = x
         self.pattern.meas = y
         self.pattern.meas_su = sy
+        self.pattern.excluded = np.full(x.shape, fill_value=False, dtype=bool)  # No excluded points by default
 
         print(f"Loaded {len(x)} points for experiment '{self.pattern.experiment.name}'.")
 

src/easydiffraction/experiments/collections/excluded_regions.py

@@ -0,0 +1,93 @@
+from typing import List, Type
+
+from easydiffraction.utils.utils import render_table
+from easydiffraction.utils.formatting import paragraph
+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 ExcludedRegion instances.
+    """
+    @property
+    def _type(self) -> str:
+        return "category"  # datablock or category
+
+    @property
+    def _child_class(self) -> Type[ExcludedRegion]:
+        return ExcludedRegion
+
+    def on_item_added(self, item: ExcludedRegion) -> None:
+        """
+        Mark excluded points in the experiment pattern when a new region is added.
+        """
+        experiment = self._parent
+        pattern = experiment.datastore.pattern
+
+        # Boolean mask for points within the new excluded region
+        in_region = ((pattern.full_x >= item.minimum.value) &
+                     (pattern.full_x <= item.maximum.value))
+
+        # Update the exclusion mask
+        pattern.excluded[in_region] = True
+
+        # Update the excluded points in the datastore
+        pattern.x = pattern.full_x[~pattern.excluded]
+        pattern.meas = pattern.full_meas[~pattern.excluded]
+        pattern.meas_su = pattern.full_meas_su[~pattern.excluded]
+
+    def show(self) -> None:
+        # TODO: Consider moving this to the base class
+        #  to avoid code duplication with implementations in Background, etc.
+        #  Consider using parameter names as column headers
+        columns_headers: List[str] = ["minimum", "maximum"]
+        columns_alignment = ["left", "left"]
+        columns_data: List[List[float]] = []
+        for region in self._items.values():
+            minimum = region.minimum.value
+            maximum = region.maximum.value
+            columns_data.append([minimum, maximum])
+
+        print(paragraph("Excluded regions"))
+        render_table(columns_headers=columns_headers,
+                     columns_alignment=columns_alignment,
+                     columns_data=columns_data)