diff --git a/notebooks/automated_refinement.ipynb b/notebooks/automated_refinement.ipynb index 31f6567..2a6e8b8 100644 --- a/notebooks/automated_refinement.ipynb +++ b/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 -} +} \ No newline at end of file diff --git a/sample_1_check.html b/sample_1_check.html new file mode 100644 index 0000000..224cdec --- /dev/null +++ b/sample_1_check.html @@ -0,0 +1,14 @@ + + + +
+
+ + \ No newline at end of file diff --git a/sample_2_check.html b/sample_2_check.html new file mode 100644 index 0000000..1eab0f8 --- /dev/null +++ b/sample_2_check.html @@ -0,0 +1,14 @@ + + + +
+
+ + \ No newline at end of file diff --git a/sample_960C_check.html b/sample_960C_check.html new file mode 100644 index 0000000..40ac6f3 --- /dev/null +++ b/sample_960C_check.html @@ -0,0 +1,14 @@ + + + +
+
+ + \ No newline at end of file diff --git a/src/dara/plot.py b/src/dara/plot.py index 9f3cd86..5471198 100644 --- a/src/dara/plot.py +++ b/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 diff --git a/src/dara/refine.py b/src/dara/refine.py index dd5d9fb..9710f12 100644 --- a/src/dara/refine.py +++ b/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, ) diff --git a/src/dara/result.py b/src/dara/result.py index 7227a24..8698038 100644 --- a/src/dara/result.py +++ b/src/dara/result.py @@ -160,14 +160,34 @@ 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 @@ -175,7 +195,31 @@ 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) diff --git a/src/dara/search/peak_matcher.py b/src/dara/search/peak_matcher.py index d154dc9..0fa688d 100644 --- a/src/dara/search/peak_matcher.py +++ b/src/dara/search/peak_matcher.py @@ -1,6 +1,9 @@ +from __future__ import annotations + from typing import Any, Literal import numpy as np +from scipy.signal import find_peaks from scipy.spatial.distance import cdist DEFAULT_ANGLE_TOLERANCE = 0.2 # maximum difference in angle @@ -8,6 +11,21 @@ # maximum ratio of the intensities to be considered as missing instead of wrong intensity DEFAULT_MAX_INTENSITY_TOLERANCE = 5 +WRONG_INTENSITY_HEIGHT_RATIO_THRESHOLD: float = 1 / 3 +COINCIDENCE_SUPPRESSION_TOLERANCE: float = 0.06 +DEFAULT_MISSING_MIN_INTENSITY_RATIO: float = 0.005 +DEFAULT_EXTRA_MIN_INTENSITY_RATIO: float = 0.03 +INTENSITY_MISMATCH_HEIGHT_TOLERANCE: float = 0.40 # height band [0.60, 1.40]; catches ~2x height mismatches +INTENSITY_MISMATCH_AREA_TOLERANCE: float = 0.60 # area band [0.40, 1.60] +INTENSITY_MISMATCH_WINDOW: float = 0.15 # degrees either side of peak apex for profile integration +RESIDUAL_WINDOW_WIDTH: float = 0.5 # degrees 2θ; sliding window width for integrated residual scan +RESIDUAL_WINDOW_STEP: float = 0.1 # degrees 2θ; step size between window centres +RESIDUAL_WINDOW_DETECT_FRACTION: float = 0.003 # per-window detection floor for forming/merging candidate regions +# merged-region threshold: flag only if total >= this fraction of total integrated obs +RESIDUAL_INTEGRAL_FRACTION: float = 0.010 +# suppress merged region if integrated(calc-bkg)/integrated(obs-bkg) >= this +RESIDUAL_CALC_COVERAGE_RATIO: float = 0.35 + def absolute_log_error(x: np.ndarray, y: np.ndarray) -> np.ndarray: """ @@ -133,7 +151,6 @@ def find_best_match( all_assigned = {m[1] for m in matched} missing = [i for i in range(len(peak_obs)) if i not in all_assigned] - # tell if a peak has wrong intensity by the sum of the intensities of the matched peaks to_be_deleted = set() for i in range(len(matched)): peak_idx = matched[i][1] @@ -175,7 +192,6 @@ def merge_peaks(peaks: np.ndarray, resolution: float = 0.0) -> np.ndarray: if len(peaks) <= 1 or resolution == 0.0: return peaks - # sorted by 0th column peaks = peaks[np.argsort(peaks[:, 0])] merge_to = np.arange(len(peaks)) @@ -218,6 +234,16 @@ class PeakMatcher: intensity_tolerance: the maximum ratio of the intensities, default to 2 max_intensity_tolerance: the maximum ratio of the intensities to be considered as missing or extra, default to 10 + profile_x: 2θ grid (degrees) for the full diffraction profile. When this and + the three profile_y_* args are all given, matched peaks are re-examined + against the profile (see ``_reclassify_wrong_intensity_by_height`` and + ``_relocate_extra_to_apex``); omit all four to skip this step. + profile_y_calc: calculated profile on ``profile_x``. + profile_y_obs: observed profile on ``profile_x``. + profile_y_bkg: background profile on ``profile_x``. + height_ratio_threshold: passed to ``_reclassify_wrong_intensity_by_height``; + wrong_intensity pairs whose calc/obs profile-height ratio falls below + this are reclassified as missing+extra instead. """ def __init__( @@ -229,6 +255,11 @@ def __init__( angle_tolerance: float = DEFAULT_ANGLE_TOLERANCE, intensity_tolerance: float = DEFAULT_INTENSITY_TOLERANCE, max_intensity_tolerance: float = DEFAULT_MAX_INTENSITY_TOLERANCE, + profile_x: np.ndarray | None = None, + profile_y_calc: np.ndarray | None = None, + profile_y_obs: np.ndarray | None = None, + profile_y_bkg: np.ndarray | None = None, + height_ratio_threshold: float = WRONG_INTENSITY_HEIGHT_RATIO_THRESHOLD, ): self.intensity_resolution = intensity_resolution self.angle_resolution = angle_resolution @@ -258,6 +289,77 @@ def __init__( max_intensity_tolerance=max_intensity_tolerance, ) + if all(v is not None for v in + [profile_x, profile_y_calc, profile_y_obs, profile_y_bkg]): + px = np.asarray(profile_x) + py_calc = np.asarray(profile_y_calc) + py_obs = np.asarray(profile_y_obs) + py_bkg = np.asarray(profile_y_bkg) + self._reclassify_wrong_intensity_by_height( + px, py_calc, py_obs, py_bkg, height_ratio_threshold, + ) + self._relocate_extra_to_apex(px, py_calc, py_bkg) + + def _reclassify_wrong_intensity_by_height( + self, + profile_x: np.ndarray, + profile_y_calc: np.ndarray, + profile_y_obs: np.ndarray, + profile_y_bkg: np.ndarray, + height_ratio_threshold: float, + ) -> None: + """Re-examine wrong_intensity pairs using background-subtracted profile heights. + + When the ratio calc_height / obs_height is below *height_ratio_threshold* + the calc peak is too weak relative to the observed feature to count as + an explanation. The pair is reclassified: obs peak -> missing, + calc peak -> extra. Pairs that pass the height check remain wrong_intensity. + """ + keep_wi: list = [] + for c_idx, o_idx in self._result["wrong_intensity"]: + tt_c = float(self.peak_calc[c_idx, 0]) + tt_o = float(self.peak_obs[o_idx, 0]) + calc_h = float(np.interp(tt_c, profile_x, profile_y_calc - profile_y_bkg)) + obs_h = float(np.interp(tt_o, profile_x, profile_y_obs - profile_y_bkg)) + if obs_h > 0 and calc_h / obs_h < height_ratio_threshold: + if o_idx not in self._result["missing"]: + self._result["missing"].append(o_idx) + if c_idx not in self._result["extra"]: + self._result["extra"].append(c_idx) + else: + keep_wi.append((c_idx, o_idx)) + self._result["wrong_intensity"] = keep_wi + + def _relocate_extra_to_apex( + self, + profile_x: np.ndarray, + profile_y_calc: np.ndarray, + profile_y_bkg: np.ndarray, + apex_window: float = 0.3, + ) -> None: + """Replace centroid 2θ of extra-calc peaks with the nearest profile-apex 2θ. + + merge_peaks() uses an intensity-weighted centroid for merged sub-peaks, which + can land between two profile maxima and corrupt coincidence-suppression distances. + Each extra peak is reassigned to the nearest local maximum of the + background-subtracted calc profile within *apex_window* degrees. + Peaks with no local maximum in their window keep their centroid position. + """ + net = profile_y_calc - profile_y_bkg + dx = float(profile_x[1] - profile_x[0]) if len(profile_x) > 1 else 0.01 + min_dist_pts = max(1, int(0.05 / dx)) + local_max_pts, _ = find_peaks(net, distance=min_dist_pts) + if len(local_max_pts) == 0: + return + local_max_x = profile_x[local_max_pts] + for c_idx in self._result["extra"]: + tt_old = float(self.peak_calc[c_idx, 0]) + in_window = (local_max_x >= tt_old - apex_window) & (local_max_x <= tt_old + apex_window) + candidates = local_max_x[in_window] + if len(candidates) == 0: + continue + self.peak_calc[c_idx, 0] = float(candidates[np.argmin(np.abs(candidates - tt_old))]) + @property def missing(self) -> np.ndarray: """Get the missing peaks in the `observed peaks`. The shape should be (N, 2) with [position, intensity].""" @@ -383,7 +485,7 @@ def get_isolated_peaks( self, peak_type: Literal["missing", "extra"], min_angle_difference: float = 0.3, - min_intensity_ratio: float = 0.03, + min_intensity_ratio: float | None = None, ) -> np.ndarray: """ Get the isolated missing peaks in the `observed peaks`. @@ -394,12 +496,21 @@ def get_isolated_peaks( Args: peak_type: the type of the peaks to consider, either "missing" or "extra" min_angle_difference: the tolerance to consider a peak as close to another peak, default to 0.3 degree - min_intensity_ratio: the minimum ratio of the intensity to be considered as a peak, default to 0.01 + min_intensity_ratio: minimum intensity relative to the observed maximum. + Defaults to DEFAULT_MISSING_MIN_INTENSITY_RATIO for missing peaks and + DEFAULT_EXTRA_MIN_INTENSITY_RATIO for extra peaks when None. Returns ------- the isolated missing peaks with [position, intensity] """ + if min_intensity_ratio is None: + min_intensity_ratio = ( + DEFAULT_MISSING_MIN_INTENSITY_RATIO + if peak_type == "missing" + else DEFAULT_EXTRA_MIN_INTENSITY_RATIO + ) + if peak_type == "missing": peaks = self.missing matched = self.matched[1] @@ -483,8 +594,260 @@ def visualize(self): label="wrong intens", ) - # add a line y=0 ax.axhline(0, color="black", lw=0.5) ax.set_xlabel("2theta") ax.set_ylabel("Intensity") ax.legend() + + +def suppress_coincident_marker_pairs( + missing_peaks: np.ndarray, + extra_peaks: np.ndarray, + tolerance: float = COINCIDENCE_SUPPRESSION_TOLERANCE, +) -> tuple[np.ndarray, np.ndarray]: + """Suppress missing/extra marker pairs within *tolerance* degrees of each other. + + Near-coincident missing+extra pairs indicate refinement position wobble: the + model placed a reflection slightly off-angle rather than genuinely missing it. + Both markers are dropped. + + Returns + ------- + (missing_peaks, extra_peaks) with coincident pairs removed, same shapes + """ + missing_peaks = np.asarray(missing_peaks).reshape(-1, 2) + extra_peaks = np.asarray(extra_peaks).reshape(-1, 2) + + if len(missing_peaks) == 0 or len(extra_peaks) == 0: + return missing_peaks, extra_peaks + + miss_suppress: set[int] = set() + extra_suppress: set[int] = set() + for m_i, (mt, _) in enumerate(missing_peaks): + for e_i, (et, _) in enumerate(extra_peaks): + if abs(float(mt) - float(et)) <= tolerance: + miss_suppress.add(m_i) + extra_suppress.add(e_i) + + keep_m = [i for i in range(len(missing_peaks)) if i not in miss_suppress] + keep_e = [i for i in range(len(extra_peaks)) if i not in extra_suppress] + + out_m = missing_peaks[keep_m] if keep_m else np.empty((0, 2)) + out_e = extra_peaks[keep_e] if keep_e else np.empty((0, 2)) + return out_m, out_e + + +def find_residual_regions( + profile_x: np.ndarray, + profile_y_obs: np.ndarray, + profile_y_calc: np.ndarray, + profile_y_bkg: np.ndarray | None = None, + window_width: float = RESIDUAL_WINDOW_WIDTH, + window_step: float = RESIDUAL_WINDOW_STEP, + window_detect_fraction: float = RESIDUAL_WINDOW_DETECT_FRACTION, + integral_fraction: float = RESIDUAL_INTEGRAL_FRACTION, + calc_profile_ratio: float = DEFAULT_EXTRA_MIN_INTENSITY_RATIO, + calc_coverage_ratio: float = RESIDUAL_CALC_COVERAGE_RATIO, + matched_peak_positions: np.ndarray | None = None, + enabled: bool = True, +) -> np.ndarray: + """Find regions where the integrated positive residual indicates an unfit feature. + + Scans the pattern with a fixed-width sliding window. Each window's integrated + positive residual (sum of max(y_obs - y_calc, 0) * dx) is compared against a + fraction of the total integrated observed intensity. Overlapping flagged windows + are merged into contiguous regions. + + Bragg-peak filter (corrected form): a merged region is suppressed only if the + calc *profile* is already adequate inside it — i.e. max(y_calc - y_bkg) within + the region exceeds ``calc_profile_ratio * max(y_obs)``. A tabulated Bragg peak + whose calc profile stays near background does NOT suppress the region. + + Args: + profile_x: 2θ grid (degrees), uniformly spaced. + profile_y_obs: observed counts on that grid. + profile_y_calc: calculated profile on that grid. + profile_y_bkg: background on that grid. When None, zeros are used (filter + compares calc profile directly against obs max). + window_width: sliding window width in degrees (default 0.5°). + window_step: step between window centres in degrees (default 0.1°). + window_detect_fraction: per-window detection floor for forming candidate + windows and merging them into regions (default 0.003 = 0.3%). This + is intentionally lower than ``integral_fraction`` so that broad + features are not penalised for being diffuse. + integral_fraction: merged-region threshold — a region is flagged only if + its total integrated positive residual exceeds this fraction of the + total integrated observed intensity (default 0.010 = 1.0%). + calc_profile_ratio: suppress a window if max(y_calc - y_bkg) inside it + exceeds this fraction of max(y_obs) (default 0.03 = 3%). + calc_coverage_ratio: suppress a merged region when + integrated(y_calc - y_bkg) / integrated(y_obs - y_bkg) over the + merged span exceeds this fraction — the calc already accounts for + this much of the observed area in that region (default 0.35). + matched_peak_positions: 1-D array of 2θ positions of already-identified + missing or extra peaks (after all normal suppression). Any merged + region whose [start, end] span contains one of these positions is + dropped — the feature has already been flagged by the matcher. + enabled: when False, skip all computation and return an empty (0, 4) + array. Default True reproduces normal behavior. + + Returns + ------- + Array of shape (N, 2) with [center_2theta, 0.0] for each flagged region, + formatted identically to missing-peak arrays for direct concatenation. + """ + if not enabled: + return np.empty((0, 2)) + + profile_x = np.asarray(profile_x, dtype=float) + profile_y_obs = np.asarray(profile_y_obs, dtype=float) + profile_y_calc = np.asarray(profile_y_calc, dtype=float) + if profile_y_bkg is not None: + profile_y_bkg = np.asarray(profile_y_bkg, dtype=float) + else: + profile_y_bkg = np.zeros_like(profile_y_obs) + + dx = float(profile_x[1] - profile_x[0]) if len(profile_x) > 1 else 0.01 + obs_max = float(profile_y_obs.max()) + half_w = window_width / 2.0 + + total_obs_integral = float(np.sum(np.maximum(profile_y_obs - profile_y_bkg, 0.0))) * dx + detect_threshold = window_detect_fraction * total_obs_integral # per-window floor + region_threshold = integral_fraction * total_obs_integral # merged-region gate + + positive_residual = np.maximum(profile_y_obs - profile_y_calc, 0.0) + + # Integrated-residual scan: flag broad features the matcher never detected. + calc_profile_threshold = calc_profile_ratio * obs_max + calc_net = profile_y_calc - profile_y_bkg + obs_net = np.maximum(profile_y_obs - profile_y_bkg, 0.0) + + flagged_centres: list[float] = [] + centre = float(profile_x[0]) + half_w + while centre <= float(profile_x[-1]) - half_w: + mask = (profile_x >= centre - half_w) & (profile_x <= centre + half_w) + win_integral = float(np.sum(positive_residual[mask])) * dx + calc_net_max = float(calc_net[mask].max()) + if win_integral >= detect_threshold and calc_net_max < calc_profile_threshold: + flagged_centres.append(centre) + centre += window_step + + if not flagged_centres: + return np.empty((0, 2)) + + merged: list[tuple[float, float]] = [] + span_start = flagged_centres[0] - half_w + span_end = flagged_centres[0] + half_w + for c in flagged_centres[1:]: + win_start = c - half_w + win_end = c + half_w + if win_start <= span_end + window_step / 2: + span_end = max(span_end, win_end) + else: + merged.append((span_start, span_end)) + span_start = win_start + span_end = win_end + merged.append((span_start, span_end)) + + known_positions: np.ndarray | None = None + if matched_peak_positions is not None: + known_positions = np.asarray(matched_peak_positions, dtype=float).reshape(-1) + + regions: list[list[float]] = [] + for s, e in merged: + mask = (profile_x >= s) & (profile_x <= e) + seg_res = positive_residual[mask] + merged_integ = float(np.sum(seg_res)) * dx + if merged_integ < region_threshold: + continue + region_obs_int = float(np.sum(obs_net[mask])) * dx + region_calc_int = float(np.sum(np.maximum(calc_net[mask], 0.0))) * dx + if region_obs_int > 0 and region_calc_int / region_obs_int >= calc_coverage_ratio: + continue + if known_positions is not None and np.any((known_positions >= s) & (known_positions <= e)): + continue + midpoint = (float(s) + float(e)) / 2.0 + regions.append([midpoint, 0.0]) + + return np.array(regions).reshape(-1, 2) if regions else np.empty((0, 2)) + + +def find_intensity_mismatch_peaks( + pm: PeakMatcher, + profile_x: np.ndarray, + profile_y_obs: np.ndarray, + profile_y_calc: np.ndarray, + profile_y_bkg: np.ndarray | None = None, + height_tolerance: float = INTENSITY_MISMATCH_HEIGHT_TOLERANCE, + area_tolerance: float = INTENSITY_MISMATCH_AREA_TOLERANCE, + window: float = INTENSITY_MISMATCH_WINDOW, +) -> np.ndarray: + """Flag matched peaks whose profile height or area deviates beyond the respective band. + + Operates only on matched pairs — missing and extra take precedence. + height band: [1 - height_tolerance, 1 + height_tolerance] + area band: [1 - area_tolerance, 1 + area_tolerance] + Flags if either band is violated. + + Args: + pm: the PeakMatcher holding the matched calc/obs peak pairs. + profile_x: 2θ grid (degrees), uniformly spaced. + profile_y_obs: observed counts on that grid. + profile_y_calc: calculated profile on that grid. + profile_y_bkg: background on that grid. When None, zeros are used. + height_tolerance: half-width of the height band (default 0.40). + area_tolerance: half-width of the area band (default 0.60). + window: degrees either side of the peak apex used for height/area + integration (default 0.15). + + Returns + ------- + (N, 2) array of [obs_2theta, 0.0], same shape as missing_peaks + """ + matched_pairs = pm._result["matched"] + if not matched_pairs: + return np.empty((0, 2)) + + profile_x = np.asarray(profile_x, dtype=float) + profile_y_obs = np.asarray(profile_y_obs, dtype=float) + profile_y_calc = np.asarray(profile_y_calc, dtype=float) + profile_y_bkg = np.zeros_like(profile_y_obs) if profile_y_bkg is None else np.asarray(profile_y_bkg, dtype=float) + + dx = float(profile_x[1] - profile_x[0]) if len(profile_x) > 1 else 0.01 + h_lo = 1.0 - height_tolerance + h_hi = 1.0 + height_tolerance + a_lo = 1.0 - area_tolerance + a_hi = 1.0 + area_tolerance + + seen_obs: set[int] = set() + flagged: list[list[float]] = [] + + for _c_idx, o_idx in matched_pairs: + if o_idx in seen_obs: + continue + seen_obs.add(o_idx) + + tt_o = float(pm.peak_obs[o_idx, 0]) + + mask = (profile_x >= tt_o - window) & (profile_x <= tt_o + window) + obs_net = np.maximum(profile_y_obs[mask] - profile_y_bkg[mask], 0.0) + calc_net = np.maximum(profile_y_calc[mask] - profile_y_bkg[mask], 0.0) + obs_area = float(np.sum(obs_net)) * dx + calc_area = float(np.sum(calc_net)) * dx + + if obs_area <= 0: + continue + + obs_h = float(np.interp(tt_o, profile_x, profile_y_obs - profile_y_bkg)) + calc_h = float(np.interp(tt_o, profile_x, profile_y_calc - profile_y_bkg)) + + if obs_h <= 0: + continue + + area_ratio = calc_area / obs_area + height_ratio = calc_h / obs_h + + if not (h_lo <= height_ratio <= h_hi) or not (a_lo <= area_ratio <= a_hi): + flagged.append([tt_o, 0.0]) + + return np.array(flagged).reshape(-1, 2) if flagged else np.empty((0, 2))