Improve missing/extra peak detection; add residual and intensity-mism…

notebooks/automated_refinement.ipynb

@@ -135513,12 +135513,7 @@
    "cell_type": "markdown",
    "id": "1b2c0fdd",
    "metadata": {},
-   "source": [
-    "After finishing refinement, you can read information from the `RefinementResult` object. The object contains the following attributes:\n",
-    "- `lst_data`: information about phases, metrics of refinement (from the .lst file in BGMN)\n",
-    "- `peak_data`: the simulated peaks in the calculated pattern\n",
-    "- `plot_data`: `x` (two-theta), `y_obs`, `y_calc`, `y_bkg`, contribution from each phase. This is mainly used for visualization."
-   ]
+   "source": "After finishing refinement, you can read information from the `RefinementResult` object. The object contains the following attributes:\n- `lst_data`: information about phases, metrics of refinement (from the .lst file in BGMN)\n- `peak_data`: the simulated peaks in the calculated pattern\n- `plot_data`: `x` (two-theta), `y_obs`, `y_calc`, `y_bkg`, contribution from each phase. This is mainly used for visualization.\n- `refinement_metrics`: peak-matching diagnostics (`missing_peaks`, `extra_peaks`, `intensity_mismatch_peaks`) plus `rwp`."
   },
   {
    "cell_type": "markdown",
@@ -316947,4 +316942,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 5
-}
+}

src/dara/plot.py

@@ -14,8 +14,16 @@ def visualize(
     diff_offset: bool = False,
     missing_peaks: list[list[float]] | np.ndarray | None = None,
     extra_peaks: list[list[float]] | np.ndarray | None = None,
+    intensity_mismatch_peaks: list[list[float]] | np.ndarray | None = None,
 ):
-    """Visualize the result from the refinement. It uses plotly as the backend engine."""
+    """Visualize the result from the refinement. It uses plotly as the backend engine.
+
+    Args:
+        intensity_mismatch_peaks: optional (N, 2) array of [2theta, 0.0] markers
+            for matched peaks whose calculated height/area deviates from
+            observed; drawn the same way as ``missing_peaks``/``extra_peaks``
+            when provided.
+    """
     colormap = [
         "#1f77b4",
         "#ff7f0e",
@@ -160,6 +168,7 @@ def visualize(
                 mode="markers",
                 marker=dict(color="#f9726a", symbol=53, size=10, opacity=0.8),
                 name="Missing peaks",
+                showlegend=True,
                 visible="legendonly",
                 text=[f"{x:.2f}, {y:.2f}" for x, y in missing_peaks],
             )
@@ -174,12 +183,29 @@ def visualize(
                 mode="markers",
                 marker=dict(color="#335da0", symbol=53, size=10, opacity=0.8),
                 name="Extra peaks",
+                showlegend=True,
                 visible="legendonly",
                 text=[f"{x:.2f}, {y:.2f}" for x, y in extra_peaks],
                 hovertemplate="%{text}",
             )
         )
 
+    if intensity_mismatch_peaks is not None:
+        intensity_mismatch_peaks = np.array(intensity_mismatch_peaks).reshape(-1, 2)
+        fig.add_trace(
+            go.Scatter(
+                x=intensity_mismatch_peaks[:, 0],
+                y=np.zeros_like(intensity_mismatch_peaks[:, 0]),
+                mode="markers",
+                marker=dict(color="#FFD700", symbol=53, size=10, opacity=0.8),
+                name="Intensity mismatch",
+                showlegend=True,
+                visible="legendonly",
+                text=[f"{x:.2f}, {y:.2f}" for x, y in intensity_mismatch_peaks],
+                hovertemplate="%{text}",
+            )
+        )
+
     title = f"{result.lst_data.pattern_name} (Rwp={result.lst_data.rwp:.2f}%)"
 
     # Updating layout with titles and labels

src/dara/refine.py

@@ -7,12 +7,13 @@
 from pathlib import Path
 from typing import Any, Literal
 
+import numpy as np
 from pydantic import BaseModel, ConfigDict, Field, field_validator
 
 from dara.bgmn_worker import BGMNWorker
 from dara.cif2str import cif2str
 from dara.generate_control_file import generate_control_file
-from dara.result import RefinementResult, get_result
+from dara.result import RefinementMetrics, RefinementResult, get_result
 from dara.xrd import convert_pattern_to_xy
 
 
@@ -65,6 +66,101 @@ def make(cls, path_obj: RefinementPhase | Path | str) -> RefinementPhase:
         )
 
 
+def _attach_peak_markers(
+    result: RefinementResult,
+    pattern_path: Path,
+    wavelength: Literal["Cu", "Co", "Cr", "Fe", "Mo"] | float,
+    instrument_profile: str | Path,
+    use_residual: bool = True,
+    residual_integral_fraction: float = 0.010,
+    residual_calc_coverage_ratio: float = 0.35,
+    residual_window_detect_fraction: float = 0.003,
+    missing_intensity_ratio: float = 0.005,
+    extra_intensity_ratio: float = 0.03,
+    intensity_mismatch_height_tolerance: float = 0.40,
+    intensity_mismatch_area_tolerance: float = 0.60,
+) -> None:
+    """Run the full peak-matching pipeline and store results on *result* in-place.
+
+    Detects observed peaks, matches them against the refined calc pattern, and
+    stores the resulting missing/extra/intensity-mismatch markers (plus rwp) on
+    ``result.refinement_metrics``.
+
+    Args:
+        use_residual: whether to additionally flag broad unfit regions via the
+            integrated-residual scan (see ``find_residual_regions``).
+        residual_integral_fraction: passed through to ``find_residual_regions``
+            as ``integral_fraction``.
+        residual_calc_coverage_ratio: passed through to ``find_residual_regions``
+            as ``calc_coverage_ratio``.
+        residual_window_detect_fraction: passed through to
+            ``find_residual_regions`` as ``window_detect_fraction``.
+        missing_intensity_ratio: minimum intensity ratio for isolated missing
+            peaks (see ``PeakMatcher.get_isolated_peaks``).
+        extra_intensity_ratio: minimum intensity ratio for isolated extra peaks
+            (see ``PeakMatcher.get_isolated_peaks``).
+        intensity_mismatch_height_tolerance: passed through to
+            ``find_intensity_mismatch_peaks`` as ``height_tolerance``.
+        intensity_mismatch_area_tolerance: passed through to
+            ``find_intensity_mismatch_peaks`` as ``area_tolerance``.
+    """
+    from dara.peak_detection import detect_peaks
+    from dara.search.peak_matcher import (
+        PeakMatcher,
+        find_intensity_mismatch_peaks,
+        find_residual_regions,
+        suppress_coincident_marker_pairs,
+    )
+
+    edf      = detect_peaks(str(pattern_path), wavelength=wavelength,
+                            instrument_profile=str(instrument_profile))
+    obs_raw  = edf[["2theta", "intensity"]].values
+    calc_raw = result.peak_data[["2theta", "intensity"]].values
+
+    px    = np.asarray(result.plot_data.x)
+    yobs  = np.asarray(result.plot_data.y_obs)
+    ycalc = np.asarray(result.plot_data.y_calc)
+    ybkg  = np.asarray(result.plot_data.y_bkg)
+
+    pm = PeakMatcher(calc_raw, obs_raw, intensity_resolution=0.005,
+                     profile_x=px, profile_y_calc=ycalc,
+                     profile_y_obs=yobs, profile_y_bkg=ybkg)
+    miss_f, extra_f = suppress_coincident_marker_pairs(
+        pm.get_isolated_peaks("missing", min_intensity_ratio=missing_intensity_ratio),
+        pm.get_isolated_peaks("extra",   min_intensity_ratio=extra_intensity_ratio),
+    )
+
+    parts = [a[:, 0] for a in (miss_f, extra_f) if len(a)]
+    known = np.concatenate(parts) if parts else None
+
+    residual = find_residual_regions(
+        px, yobs, ycalc,
+        profile_y_bkg=ybkg,
+        matched_peak_positions=known,
+        enabled=use_residual,
+        window_detect_fraction=residual_window_detect_fraction,
+        integral_fraction=residual_integral_fraction,
+        calc_coverage_ratio=residual_calc_coverage_ratio,
+    )
+
+    miss_combined = np.vstack([miss_f, residual]) if len(residual) and len(miss_f) else (
+        residual if len(residual) else miss_f
+    )
+
+    mismatch = find_intensity_mismatch_peaks(
+        pm, px, yobs, ycalc, profile_y_bkg=ybkg,
+        height_tolerance=intensity_mismatch_height_tolerance,
+        area_tolerance=intensity_mismatch_area_tolerance,
+    )
+
+    result.refinement_metrics = RefinementMetrics(
+        missing_peaks=miss_combined if len(miss_combined) else None,
+        extra_peaks=extra_f if len(extra_f) else None,
+        intensity_mismatch_peaks=mismatch if len(mismatch) else None,
+        rwp=result.lst_data.rwp,
+    )
+
+
 def do_refinement(
     pattern_path: Path | str,
     phases: list[RefinementPhase | Path | str],
@@ -74,8 +170,27 @@ def do_refinement(
     phase_params: dict | None = None,
     refinement_params: dict | None = None,
     show_progress: bool = False,
+    use_residual: bool = True,
+    residual_integral_fraction: float = 0.010,
+    residual_calc_coverage_ratio: float = 0.35,
+    residual_window_detect_fraction: float = 0.003,
+    missing_intensity_ratio: float = 0.005,
+    extra_intensity_ratio: float = 0.03,
+    intensity_mismatch_height_tolerance: float = 0.40,
+    intensity_mismatch_area_tolerance: float = 0.60,
 ) -> RefinementResult:
-    """Refine the structure using BGMN."""
+    """Refine the structure using BGMN.
+
+    Args:
+        use_residual: see ``_attach_peak_markers``.
+        residual_integral_fraction: see ``_attach_peak_markers``.
+        residual_calc_coverage_ratio: see ``_attach_peak_markers``.
+        residual_window_detect_fraction: see ``_attach_peak_markers``.
+        missing_intensity_ratio: see ``_attach_peak_markers``.
+        extra_intensity_ratio: see ``_attach_peak_markers``.
+        intensity_mismatch_height_tolerance: see ``_attach_peak_markers``.
+        intensity_mismatch_area_tolerance: see ``_attach_peak_markers``.
+    """
     pattern_path = Path(pattern_path)
     working_dir = (
         Path(working_dir)
@@ -100,7 +215,6 @@ def do_refinement(
         phase = RefinementPhase.make(phase_path)
         phase_path_ = phase.path
         phase_params_ = phase_params.copy()
-        # Update the default phase parameters with the specific parameters for the phase
         phase_params_.update(phase.params)
         if phase_path_.suffix == ".cif":
             str_path = cif2str(phase_path_, "", working_dir, **phase_params_)
@@ -121,7 +235,19 @@ def do_refinement(
 
     bgmn_worker = BGMNWorker()
     bgmn_worker.run_refinement_cmd(control_file_path, show_progress=show_progress)
-    return get_result(control_file_path)
+    result = get_result(control_file_path)
+    _attach_peak_markers(
+        result, pattern_path, wavelength, instrument_profile,
+        use_residual=use_residual,
+        residual_integral_fraction=residual_integral_fraction,
+        residual_calc_coverage_ratio=residual_calc_coverage_ratio,
+        residual_window_detect_fraction=residual_window_detect_fraction,
+        missing_intensity_ratio=missing_intensity_ratio,
+        extra_intensity_ratio=extra_intensity_ratio,
+        intensity_mismatch_height_tolerance=intensity_mismatch_height_tolerance,
+        intensity_mismatch_area_tolerance=intensity_mismatch_area_tolerance,
+    )
+    return result
 
 
 def do_refinement_no_saving(
@@ -132,18 +258,43 @@ def do_refinement_no_saving(
     phase_params: dict | None = None,
     refinement_params: dict | None = None,
     show_progress: bool = False,
+    use_residual: bool = True,
+    residual_integral_fraction: float = 0.010,
+    residual_calc_coverage_ratio: float = 0.35,
+    residual_window_detect_fraction: float = 0.003,
+    missing_intensity_ratio: float = 0.005,
+    extra_intensity_ratio: float = 0.03,
+    intensity_mismatch_height_tolerance: float = 0.40,
+    intensity_mismatch_area_tolerance: float = 0.60,
 ) -> RefinementResult:
-    """Refine the structure using BGMN in a temporary directory without saving."""
-    with tempfile.TemporaryDirectory() as tmpdir:
-        working_dir = Path(tmpdir)
+    """Refine the structure using BGMN in a temporary directory without saving.
 
+    Args:
+        use_residual: see ``_attach_peak_markers``.
+        residual_integral_fraction: see ``_attach_peak_markers``.
+        residual_calc_coverage_ratio: see ``_attach_peak_markers``.
+        residual_window_detect_fraction: see ``_attach_peak_markers``.
+        missing_intensity_ratio: see ``_attach_peak_markers``.
+        extra_intensity_ratio: see ``_attach_peak_markers``.
+        intensity_mismatch_height_tolerance: see ``_attach_peak_markers``.
+        intensity_mismatch_area_tolerance: see ``_attach_peak_markers``.
+    """
+    with tempfile.TemporaryDirectory() as tmpdir:
         return do_refinement(
             pattern_path=pattern_path,
             phases=phases,
             wavelength=wavelength,
             instrument_profile=instrument_profile,
-            working_dir=working_dir,
+            working_dir=Path(tmpdir),
             phase_params=phase_params,
             refinement_params=refinement_params,
             show_progress=show_progress,
+            use_residual=use_residual,
+            residual_integral_fraction=residual_integral_fraction,
+            residual_calc_coverage_ratio=residual_calc_coverage_ratio,
+            residual_window_detect_fraction=residual_window_detect_fraction,
+            missing_intensity_ratio=missing_intensity_ratio,
+            extra_intensity_ratio=extra_intensity_ratio,
+            intensity_mismatch_height_tolerance=intensity_mismatch_height_tolerance,
+            intensity_mismatch_area_tolerance=intensity_mismatch_area_tolerance,
         )

src/dara/result.py

@@ -160,22 +160,66 @@ class DiaResult(BaseModel):
     structs: dict[str, list[float]]
 
 
+class RefinementMetrics(BaseModel):
+    """Derived peak-matching diagnostics and refinement quality metrics.
+
+    Holds the missing/extra/intensity-mismatch peak markers produced by peak
+    matching after a refinement, plus the refinement's Rwp.
+    """
+
+    model_config = ConfigDict(arbitrary_types_allowed=True, populate_by_name=True)
+
+    missing_peaks: np.ndarray | None = Field(default=None, repr=False)
+    extra_peaks: np.ndarray | None = Field(default=None, repr=False)
+    intensity_mismatch_peaks: np.ndarray | None = Field(default=None, repr=False)
+    rwp: float
+
+
 class RefinementResult(BaseModel):
-    """The result from the refinement, which is parsed from the .lst and .dia files."""
+    """The result from the refinement, which is parsed from the .lst and .dia files.
+
+    ``refinement_metrics`` holds the peak-matching diagnostics attached after
+    refinement (see ``RefinementMetrics``).
+    """
 
     model_config = ConfigDict(arbitrary_types_allowed=True, populate_by_name=True)
 
     lst_data: LstResult
     plot_data: DiaResult = Field(repr=False)
     peak_data: pd.DataFrame = Field(repr=False)
+    refinement_metrics: RefinementMetrics
 
     @field_validator("peak_data", mode="before")
     @classmethod
     def transform(cls, data: dict) -> pd.DataFrame:
         """Create pandas dataframe from peak data dict."""
         return pd.DataFrame(data)
 
-    def visualize(self, diff_offset=False):
+    def visualize(self, diff_offset=False, plot_refinement_metrics: bool = True):
+        """Visualize the refinement result with plotly.
+
+        Args:
+            diff_offset: passed through to ``dara.plot.visualize``.
+            plot_refinement_metrics: when True (default) and ``refinement_metrics``
+                is present, draw its missing/extra/intensity-mismatch peak markers.
+                When False, or when ``refinement_metrics`` is unavailable, no peak
+                markers are drawn.
+
+        Returns
+        -------
+            the plotly ``Figure`` for the refinement plot
+        """
+        metrics = getattr(self, "refinement_metrics", None)
+        if metrics is not None and plot_refinement_metrics:
+            # A trace with zero points never appears in the legend even with
+            # showlegend=True, so an empty category gets a single non-rendering
+            # NaN point instead, just to keep its legend entry toggleable.
+            placeholder = np.array([[np.nan, np.nan]])
+            return visualize(self, diff_offset=diff_offset,
+                             missing_peaks=metrics.missing_peaks if metrics.missing_peaks is not None else placeholder,
+                             extra_peaks=metrics.extra_peaks if metrics.extra_peaks is not None else placeholder,
+                             intensity_mismatch_peaks=metrics.intensity_mismatch_peaks
+                             if metrics.intensity_mismatch_peaks is not None else placeholder)
         return visualize(self, diff_offset=diff_offset)
 
     def get_phase_weights(self, normalize=True) -> dict[str, float]:
@@ -251,10 +295,12 @@ def get_result(control_file: Path) -> RefinementResult:
         dia_path = control_file.parent / f"{control_file.stem}.dia"
         par_path = control_file.parent / f"{control_file.stem}.par"
 
+        lst_data = parse_lst(lst_path, phase_names=phase_names)
         result = {
-            "lst_data": parse_lst(lst_path, phase_names=phase_names),
+            "lst_data": lst_data,
             "plot_data": parse_dia(dia_path, phase_names=phase_names),
             "peak_data": parse_par(par_path, phase_names=phase_names),
+            "refinement_metrics": RefinementMetrics(rwp=lst_data.rwp),
         }
 
         return RefinementResult(**result)