#!/usr/bin/env python3 from __future__ import annotations import argparse import itertools import importlib.util import io import math import sys from dataclasses import dataclass, asdict from pathlib import Path from typing import Iterable import matplotlib import numpy as np import pandas as pd import scipy.signal from openpyxl import Workbook from openpyxl.styles import Font import matplotlib.pyplot as plt """ 概要: 粉末X線回折データのためのピーク探索と格子定数推定を行うスクリプトです。 詳細説明: このスクリプトは、X線回折データからピークを検出し、その情報に基づいて 様々な結晶系の格子定数を推定する機能を提供します。 また、検出されたピークと理論的な回折線を比較し、結果をグラフやExcelファイルとして出力できます。 主な機能: 1. ピーク探索: Savitzky-Golayフィルターと3次導関数ゼロ交差法を用いて、XRDデータからピークを検出します。 2. Kα2割り当て: 検出されたピークの中からKα2ピークを識別し、対応するKα1ピークと関連付けます。 3. 格子定数推定: 検出ピークの逆d二乗値に基づいて、結晶系ごとの格子定数を推定します。 推定方法には、最小二乗法に基づく組み合わせ探索と、格子定数グリッド探索があります。 4. 結果比較とプロット: 推定された格子定数、またはユーザーが指定した格子定数と観測ピークを比較し、 理論的な回折線を重ね合わせたプロットを生成します。 5. データ入出力: 生のXRDデータファイル、ピーク情報ファイル、推定結果ファイル（すべてExcel形式を推奨）の読み書きをサポートします。 関連リンク: guess_lattice_parameters_usage """ CU_KA1 = 1.54056 CU_KA2 = 1.54439 ANGLE_EPS = 1.0e-12 SUPPORTED_CRYSTAL_SYSTEMS = ["hexagonal", "orthorhombic", "tetragonal", "cubic"] CRYSTAL_SYSTEM_ALIASES = { "auto": "auto", "all": "auto", "any": "auto", "hex": "hexagonal", "hexagonal": "hexagonal", "orth": "orthorhombic", "ortho": "orthorhombic", "orthorhombic": "orthorhombic", "tet": "tetragonal", "tetragonal": "tetragonal", "cub": "cubic", "cubic": "cubic", } KNOWN_MODES = {"search", "guess", "all", "compare", "calc", "index", "seaerch"} SUPPORTED_METHODS = ["combinatorial", "grid", "both"] @dataclass class Peak: """ 概要: 検出されたピークの情報を保持するデータクラスです。 属性: :param index: データ配列におけるピークのインデックス。 :type index: int :param two_theta: ピークの中心角度2θ（度）。 :type two_theta: float :param intensity: スムージング後のピーク強度。 :type intensity: float :param intensity_raw: 元データのピーク強度。 :type intensity_raw: float :param fwhm_deg: 半値全幅（度）。 :type fwhm_deg: float :param inv_d2: d間隔の二乗の逆数（1 / d^2）。 :type inv_d2: float :param d: d間隔（オングストローム）。 :type d: float :param q: 散乱ベクトルq（A^-1）。 :type q: float :param ka_role: Kα1またはKα2の役割を示す文字列（"ka1", "ka2", ""）。 :type ka_role: str :param ka_pair_index: Kαペアの相手のピークのインデックス（存在しない場合はNone）。 :type ka_pair_index: int | None :param source: ピークの検出元（"detected"またはファイルパス）。 :type source: str """ index: int two_theta: float intensity: float intensity_raw: float fwhm_deg: float inv_d2: float d: float q: float ka_role: str = "" ka_pair_index: int | None = None source: str = "detected" @dataclass class Candidate: """ 概要: 格子定数推定の候補結果を保持するデータクラスです。 属性: :param crystal_system: 候補の結晶系（例: "cubic", "tetragonal"）。 :type crystal_system: str :param ls_code: 結晶系に対応する最小二乗コード。 :type ls_code: int :param params: 格子定数と角度の辞書。 :type params: dict[str, float] :param score_matches: マッチしたピークの数。 :type score_matches: int :param score_rms_rel: 相対RMS誤差。 :type score_rms_rel: float :param selected_matches: 推定に使用されたピークのマッチング詳細のリスト。 :type selected_matches: list[dict] :param all_matches: 全ての観測ピークのマッチング詳細のリスト。 :type all_matches: list[dict] """ crystal_system: str ls_code: int params: dict[str, float] score_matches: int score_rms_rel: float selected_matches: list[dict] all_matches: list[dict] def bragg_d(two_theta_deg: float, wavelength: float = CU_KA1) -> float: """ 概要: ブラッグの法則に基づいてd間隔を計算します。 詳細説明: 指定された2θ角度とX線波長から、ブラッグの法則 (n * lambda = 2 * d * sin(theta)) を用いて 回折面間隔dを計算します。thetaはtwo_theta_degの半分です。 引数: :param two_theta_deg: 回折角度2θ（度）。 :type two_theta_deg: float :param wavelength: X線波長（オングストローム）。デフォルトはCU_KA1です。 :type wavelength: float 戻り値: :returns: 計算されたd間隔（オングストローム）。 :rtype: float """ theta = math.radians(two_theta_deg / 2.0) s = math.sin(theta) if s <= 0: return float("inf") return wavelength / (2.0 * s) def inv_d2_from_two_theta(two_theta_deg: float, wavelength: float = CU_KA1) -> float: """ 概要: 2θ角度からd間隔の二乗の逆数 (1 / d^2) を計算します。 詳細説明: ブラッグの法則を用いてd間隔を計算し、その二乗の逆数を返します。 引数: :param two_theta_deg: 回折角度2θ（度）。 :type two_theta_deg: float :param wavelength: X線波長（オングストローム）。デフォルトはCU_KA1です。 :type wavelength: float 戻り値: :returns: d間隔の二乗の逆数（inv_d2）。 :rtype: float """ d = bragg_d(two_theta_deg, wavelength) if not math.isfinite(d) or d <= 0: return float("nan") return 1.0 / (d * d) def two_theta_from_d(d: float, wavelength: float = CU_KA1) -> float | None: """ 概要: d間隔から2θ角度を計算します。 詳細説明: ブラッグの法則を逆算して、指定されたd間隔とX線波長に対応する回折角度2θを計算します。 2 * d * sin(theta) = wavelength の関係を使用します。 引数: :param d: d間隔（オングストローム）。 :type d: float :param wavelength: X線波長（オングストローム）。デフォルトはCU_KA1です。 :type wavelength: float 戻り値: :returns: 計算された2θ角度（度）、または回折が物理的に不可能な場合はNone。 :rtype: float | None """ x = wavelength / (2.0 * d) if not (0.0 < x < 1.0): return None return math.degrees(2.0 * math.asin(x)) def two_theta_from_inv_d2(inv_d2: float, wavelength: float = CU_KA1) -> float | None: """ 概要: d間隔の二乗の逆数 (1 / d^2) から2θ角度を計算します。 詳細説明: inv_d2値からd間隔を計算し、そのd間隔とX線波長に対応する2θ角度を返します。 引数: :param inv_d2: d間隔の二乗の逆数。 :type inv_d2: float :param wavelength: X線波長（オングストローム）。デフォルトはCU_KA1です。 :type wavelength: float 戻り値: :returns: 計算された2θ角度（度）、または計算が不可能な場合はNone。 :rtype: float | None """ if inv_d2 <= 0: return None return two_theta_from_d(1.0 / math.sqrt(inv_d2), wavelength) def sanitize_stem(text: str) -> str: """ 概要: ファイル名として安全なステム文字列を生成します。 詳細説明: 入力された文字列から、ファイル名に使用できない特殊文字を除去し、 英数字、ハイフン、アンダースコア、ピリオドのみを含む文字列を返します。 連続するアンダースコアは1つにまとめられます。 引数: :param text: 処理する文字列。 :type text: str 戻り値: :returns: サニタイズされたファイル名ステム。 :rtype: str """ return "".join(ch if ch.isalnum() or ch in "-_." else "_" for ch in text).strip("_") or "output" def read_xy_data(path: Path) -> tuple[np.ndarray, np.ndarray]: """ 概要: 指定されたファイルからX-Yデータ（2θ-強度データ）を読み込みます。 詳細説明: ファイルの種類（Excelまたはテキスト/CSV）を判別し、 最初の2列をそれぞれX（2θ）とY（強度）として読み込みます。 数値変換エラーや非有限値は除去され、X値でソートされます。 引数: :param path: 入力データファイルへのパス。 :type path: Path 戻り値: :returns: X値のnumpy.ndarray、Y値のnumpy.ndarrayのタプル。 :rtype: tuple 例外: :raises ValueError: スプレッドシートやテキストデータに必要な列数がない場合。 """ suffix = path.suffix.lower() if suffix in {".xlsx", ".xls"}: df = pd.read_excel(path) if df.shape[1] < 2: raise ValueError("スプレッドシートには少なくとも2列が必要です。") x = pd.to_numeric(df.iloc[:, 0], errors="coerce").to_numpy() y = pd.to_numeric(df.iloc[:, 1], errors="coerce").to_numpy() else: df = pd.read_csv(path, sep=None, engine="python", header=None, comment="#") if df.shape[1] < 2: raise ValueError("テキストデータには少なくとも2列が必要です。") x = pd.to_numeric(df.iloc[:, 0], errors="coerce").to_numpy() y = pd.to_numeric(df.iloc[:, 1], errors="coerce").to_numpy() mask = np.isfinite(x) & np.isfinite(y) x = x[mask].astype(float) y = y[mask].astype(float) order = np.argsort(x) x = x[order] y = y[order] return x, y def _interp_zero(x1: float, y1: float, x2: float, y2: float) -> float: """ 概要: 2点間の線形補間によってY値がゼロになるX値を推定します。 詳細説明: (x1, y1) と (x2, y2) の2点を通る直線上で、y = 0 となるxの値を計算します。 y2 - y1 が極めて小さい場合は、単純に2点の中間点を返します。 引数: :param x1: 最初の点のX座標。 :type x1: float :param y1: 最初の点のY座標。 :type y1: float :param x2: 2番目の点のX座標。 :type x2: float :param y2: 2番目の点のY座標。 :type y2: float 戻り値: :returns: Y値がゼロとなるX値。 :rtype: float """ if abs(y2 - y1) <= ANGLE_EPS: return 0.5 * (x1 + x2) return float(x1 - y1 * (x2 - x1) / (y2 - y1)) def _find_zero_crossing_bounds( x: np.ndarray, ydiff3: np.ndarray, idx: int, ) -> tuple[tuple[int, float] | None, tuple[int, float] | None]: """ 概要: 与えられたインデックスの周りで3次導関数のゼロ交差を見つけます。 詳細説明: ピーク位置のインデックス idx を中心に、ydiff3（3次導関数）の符号が変化する点を 左右方向に探索し、そのインデックスと補間されたX座標を返します。 引数: :param x: X座標のnumpy配列。 :type x: numpy.ndarray :param ydiff3: 3次導関数のnumpy配列。 :type ydiff3: numpy.ndarray :param idx: ピーク中心のインデックス。 :type idx: int 戻り値: :returns: 左側のゼロ交差情報 (インデックス, X座標) と右側のゼロ交差情報のタプル。 見つからない場合はNoneが含まれます。 :rtype: tuple[tuple[int, float] | None, tuple[int, float] | None] """ left_info = None right_info = None for i in range(idx, 0, -1): y1 = float(ydiff3[i - 1]) y2 = float(ydiff3[i]) if y1 == 0.0 or y1 * y2 <= 0.0: xz = _interp_zero(float(x[i - 1]), y1, float(x[i]), y2) left_info = (i, xz) break for i in range(idx + 1, len(x)): y1 = float(ydiff3[i - 1]) y2 = float(ydiff3[i]) if y2 == 0.0 or y1 * y2 <= 0.0: xz = _interp_zero(float(x[i - 1]), y1, float(x[i]), y2) right_info = (i, xz) break return left_info, right_info def estimate_fwhm( x: np.ndarray, ysmooth: np.ndarray, idx: int, ydiff3: np.ndarray | None = None, ) -> float | None: """ 概要: ピークの半値全幅（FWHM）を推定します。 詳細説明: ピークトップの周辺のデータを用いてFWHMを計算します。 ydiff3（3次導関数）が提供されている場合、3次導関数のゼロ交差を境界として 局所的なFWHMを優先的に推定します。これにより、隣接するピークの影響を軽減します。 ydiff3が提供されない場合や、局所的な推定が不可能な場合は、 ピークトップの周辺の最小強度をベースラインとしてFWHMを計算します。 引数: :param x: X座標のnumpy配列。 :type x: numpy.ndarray :param ysmooth: スムージングされたY座標のnumpy配列。 :type ysmooth: numpy.ndarray :param idx: FWHMを推定するピークトップのインデックス。 :type idx: int :param ydiff3: 3次導関数のnumpy配列。デフォルトはNone。 :type ydiff3: numpy.ndarray | None 戻り値: :returns: 推定されたFWHM（度）、または推定できなかった場合はNone。 :rtype: float | None """ ytop = float(ysmooth[idx]) if not math.isfinite(ytop) or ytop <= 0.0: return None left_info = None right_info = None if ydiff3 is not None: left_info, right_info = _find_zero_crossing_bounds(x, ydiff3, idx) if left_info is not None and right_info is not None: il, xl0 = left_info ir, xr0 = right_info if il < idx < ir: ybg_left = float(ysmooth[max(il - 1, 0)]) ybg_right = float(ysmooth[min(ir, len(ysmooth) - 1)]) ybg = min(ybg_left, ybg_right) half = ybg + 0.5 * (ytop - ybg) xl = None for i in range(idx, il - 1, -1): if i <= 0: break y1 = float(ysmooth[i - 1]) y2 = float(ysmooth[i]) if (y1 - half) == 0.0 or (y1 - half) * (y2 - half) <= 0.0: xl = _interp_zero(float(x[i - 1]), y1 - half, float(x[i]), y2 - half) break xr = None for i in range(idx + 1, ir + 1): if i >= len(x): break y1 = float(ysmooth[i - 1]) y2 = float(ysmooth[i]) if (y2 - half) == 0.0 or (y1 - half) * (y2 - half) <= 0.0: xr = _interp_zero(float(x[i - 1]), y1 - half, float(x[i]), y2 - half) break if xl is not None and xr is not None and xr > xl: return float(xr - xl) width_like = float(xr0 - xl0) if width_like > 0.0: return width_like lo = max(0, idx - 10) hi = min(len(x), idx + 11) local_bg = float(np.min(ysmooth[lo:hi])) half = local_bg + 0.5 * (ytop - local_bg) il = idx while il > 0 and ysmooth[il] > half: il -= 1 ir = idx while ir < len(x) - 1 and ysmooth[ir] > half: ir += 1 if il >= idx or ir <= idx: return None if abs(float(ysmooth[il + 1]) - float(ysmooth[il])) <= ANGLE_EPS: xl = float(x[il]) else: xl = float(np.interp(half, [ysmooth[il], ysmooth[il + 1]], [x[il], x[il + 1]])) if abs(float(ysmooth[ir]) - float(ysmooth[ir - 1])) <= ANGLE_EPS: xr = float(x[ir]) else: xr = float(np.interp(half, [ysmooth[ir - 1], ysmooth[ir]], [x[ir - 1], x[ir]])) return float(xr - xl) if xr > xl else None def peak_search_deriv3( x: np.ndarray, y: np.ndarray, nsmooth: int = 11, norder: int = 4, threshold: float = 1000.0, ydiff1_threshold: float = 1.0e-2, fwhm_min_deg: float = 0.06, fwhm_max_deg: float = 1.0, is_print: bool = True, ) -> tuple[list[Peak], dict[str, np.ndarray | float]]: """ 概要: Savitzky-Golayフィルターと3次導関数ゼロ交差法を用いてピークを探索します。 詳細説明: 入力X線回折データ (x, y) に対してSavitzky-Golayフィルターを適用し、 スムージングされたデータおよび1次、2次、3次導関数を計算します。 3次導関数のゼロ交差点をピーク候補点として抽出し、強度閾値、FWHM範囲、 2次導関数の符号などの基準でフィルタリングを行い、有効なピークを特定します。 重複するピークは強度に基づいてマージされます。 引数: :param x: 2θ角度のnumpy配列。 :type x: numpy.ndarray :param y: 強度のnumpy配列。 :type y: numpy.ndarray :param nsmooth: Savitzky-Golayフィルターのウィンドウサイズ（奇数）。 :type nsmooth: int :param norder: Savitzky-Golayフィルターの多項式の次数。 :type norder: int :param threshold: ピーク強度の最小閾値。 :type threshold: float :param ydiff1_threshold: 1次導関数の相対閾値（診断用）。 :type ydiff1_threshold: float :param fwhm_min_deg: 許容されるFWHMの最小値（度）。 :type fwhm_min_deg: float :param fwhm_max_deg: 許容されるFWHMの最大値（度）。 :type fwhm_max_deg: float :param is_print: 診断情報をコンソールに出力するかどうか。 :type is_print: bool 戻り値: :returns: 検出されたPeakオブジェクトのリストと、スムージングデータや導関数などの診断情報の辞書。 :rtype: tuple[list[Peak], dict[str, numpy.ndarray | float]] 例外: :raises ValueError: データポイントが不足している場合。 """ if len(x) < max(nsmooth + 2, 7): raise ValueError("データポイントが不足しています。") if nsmooth % 2 == 0: nsmooth += 1 if nsmooth <= norder: nsmooth = norder + 3 if (norder + 3) % 2 == 1 else norder + 4 h = float(np.median(np.diff(x))) ysmooth = scipy.signal.savgol_filter(y, nsmooth, norder, deriv=0) ydiff1 = scipy.signal.savgol_filter(y, nsmooth, norder, deriv=1) / h ydiff2 = scipy.signal.savgol_filter(ydiff1, nsmooth, norder, deriv=1) / h ydiff3 = scipy.signal.savgol_filter(ydiff2, nsmooth, norder, deriv=1) / h ydiff3 = scipy.signal.savgol_filter(ydiff3, nsmooth, norder, deriv=0) diff1_ratio = np.abs(ydiff1 / np.maximum(ysmooth, 1.0e-5)) max_diff1_ratio = float(np.max(diff1_ratio)) diff1_ratio_th = max_diff1_ratio * ydiff1_threshold peaks: list[Peak] = [] last_accept_idx = -10**9 visited_idx: set[int] = set() if is_print: print("") print("=== ピーク探索診断 ===") print(f"nsmooth : {nsmooth}") print(f"norder : {norder}") print(f"threshold : {threshold}") print(f"FWHM window : {fwhm_min_deg} - {fwhm_max_deg} deg") print(f"|dy/dx|/y thresh : {diff1_ratio_th:.6g} (診断のみ)") print("") for i in range(1, len(x)): if ydiff3[i - 1] == 0.0 or ydiff3[i - 1] * ydiff3[i] <= 0.0: lo = max(0, i - max(3, nsmooth // 2)) hi = min(len(x), i + max(4, nsmooth // 2 + 1)) idx = lo + int(np.argmax(ysmooth[lo:hi])) if idx in visited_idx: continue visited_idx.add(idx) if idx - last_accept_idx <= 1: continue ytop = float(ysmooth[idx]) d1ratio = float(abs(ydiff1[idx]) / max(ysmooth[idx], 1.0e-5)) d2 = float(ydiff2[idx]) reason = None fwhm = estimate_fwhm(x, ysmooth, idx, ydiff3=ydiff3) if ytop < threshold: reason = "強度が閾値未満" elif d2 >= 0.0 and (fwhm is None or fwhm < fwhm_min_deg): reason = "最小値のような形状" elif fwhm is None: reason = "FWHMが利用不可" elif fwhm < fwhm_min_deg: reason = f"FWHMが小さすぎる ({fwhm:.4f})" elif fwhm > fwhm_max_deg: reason = f"FWHMが大きすぎる ({fwhm:.4f})" if is_print: status = "採用" if reason is None else f"除外: {reason}" print( f"2theta={x[idx]:8.3f} I={ytop:10.3f} FWHM={fwhm if fwhm is not None else float('nan'):7.4f} " f"|dy/dx|/y={d1ratio:10.5g} d2={d2:10.5g} -> {status}" ) if reason is not None: continue inv_d2 = inv_d2_from_two_theta(float(x[idx]), CU_KA1) d = 1.0 / math.sqrt(inv_d2) q = 2.0 * math.pi / d peaks.append( Peak( index=int(idx), two_theta=float(x[idx]), intensity=float(ysmooth[idx]), intensity_raw=float(y[idx]), fwhm_deg=float(fwhm), inv_d2=float(inv_d2), d=float(d), q=float(q), ) ) last_accept_idx = idx # ごく近い重複をマージ merged: list[Peak] = [] for p in peaks: if merged and abs(p.two_theta - merged[-1].two_theta) < max(p.fwhm_deg, merged[-1].fwhm_deg) * 0.4: if p.intensity > merged[-1].intensity: merged[-1] = p else: merged.append(p) info = { "ysmooth": ysmooth, "ydiff1": ydiff1, "ydiff2": ydiff2, "ydiff3": ydiff3, "max_diff1_ratio": max_diff1_ratio, "diff1_ratio_th": diff1_ratio_th, } return merged, info def assign_ka2(peaks: list[Peak], tol_deg: float = 0.06, ratio_min: float = 0.18, ratio_max: float = 0.75) -> list[Peak]: """ 概要: 検出されたピークにKα2の役割を割り当てます。 詳細説明: Cu Kα1の波長に対応するピークに対して、Cu Kα2の波長で同じd間隔を持つピークを探索します。 2θ角度の許容範囲 (tol_deg) とKα2/Kα1強度比 (ratio_min, ratio_max) に基づいてペアを特定し、 それぞれのピークに "ka1" または "ka2" の役割を割り当てます。 最も強いピークから順に処理し、一度Kα2として使用されたピークは再利用されません。 引数: :param peaks: Peakオブジェクトのリスト。 :type peaks: list[Peak] :param tol_deg: Kα2ペアを識別するための2θ角度の許容誤差（度）。 :type tol_deg: float :param ratio_min: Kα2ピークとKα1ピークの強度比の最小値。 :type ratio_min: float :param ratio_max: Kα2ピークとKα1ピークの強度比の最大値。 :type ratio_max: float 戻り値: :returns: Kα役割が割り当てられたPeakオブジェクトのリスト。 :rtype: list[Peak] """ for p in peaks: p.ka_role = "" p.ka_pair_index = None # 強度/幅が強い/狭い順 order = sorted(range(len(peaks)), key=lambda i: peaks[i].intensity, reverse=True) used_secondary = set() for i in order: p1 = peaks[i] d = bragg_d(p1.two_theta, CU_KA1) tt2 = two_theta_from_d(d, CU_KA2) if tt2 is None: continue best_j = None best_delta = None for j, p2 in enumerate(peaks): if j == i or j in used_secondary: continue if p2.two_theta <= p1.two_theta: continue delta = abs(p2.two_theta - tt2) ratio = p2.intensity / max(p1.intensity, 1e-9) if delta <= tol_deg and ratio_min <= ratio <= ratio_max: if best_delta is None or delta < best_delta: best_delta = delta best_j = j if best_j is not None: peaks[i].ka_role = "ka1" peaks[i].ka_pair_index = best_j peaks[best_j].ka_role = "ka2" peaks[best_j].ka_pair_index = i used_secondary.add(best_j) return peaks def system_rows(system: str, hmax: int = 12) -> list[tuple[tuple[int, ...], tuple[int, int, int]]]: """ 概要: 指定された結晶系とhkl指数に対して、格子定数計算のための特徴量とhklの組み合わせを生成します。 詳細説明: 立方晶、正方晶、六方晶、斜方晶の各結晶系について、hmaxまでのhkl指数を生成します。 各hkl組み合わせに対して、inv_d2 (1 / d^2) の計算に必要な特徴量タプル (feat) と、 実際のhkl指数タプルをペアとしてリストに格納します。 重複する特徴量を持つ組み合わせは除外され、特徴量の合計値でソートされます。 引数: :param system: 結晶系名（例: "cubic", "hexagonal"）。 :type system: str :param hmax: h, k, l指数の最大値。 :type hmax: int 戻り値: :returns: (特徴量タプル, hklタプル) のリスト。 :rtype: list[tuple[tuple[int, ...], tuple[int, int, int]]] 例外: :raises ValueError: サポートされていない結晶系が指定された場合。 """ rows: list[tuple[tuple[int, ...], tuple[int, int, int]]] = [] seen: set[tuple[int, ...]] = set() for h in range(hmax + 1): for k in range(hmax + 1): for l in range(hmax + 1): if h == 0 and k == 0 and l == 0: continue hkl = (h, k, l) if system == "cubic": feat = (h * h + k * k + l * l,) hkl = tuple(sorted((h, k, l), reverse=True)) elif system == "tetragonal": feat = (h * h + k * k, l * l) hk = tuple(sorted((h, k), reverse=True)) hkl = (hk[0], hk[1], l) elif system == "hexagonal": feat = (h * h + h * k + k * k, l * l) elif system == "orthorhombic": feat = (h * h, k * k, l * l) else: raise ValueError(system) if all(v == 0 for v in feat): continue if feat in seen: continue seen.add(feat) rows.append((feat, hkl)) rows.sort(key=lambda t: (sum(t[0]), t[0])) return rows def params_from_beta(system: str, beta: np.ndarray) -> dict[str, float]: """ 概要: 最小二乗法の係数 (beta) から格子定数を計算します。 詳細説明: inv_d2 = beta[0]*feat[0] + beta[1]*feat[1] + ... の形式で得られた 最小二乗法の係数 beta から、各結晶系に対応するa, b, c, alpha, beta, gammaの 格子定数を導出します。 引数: :param system: 結晶系名。 :type system: str :param beta: 最小二乗法で得られた係数のnumpy配列。 :type beta: numpy.ndarray 戻り値: :returns: 格子定数を含む辞書。 :rtype: dict[str, float] 例外: :raises ValueError: サポートされていない結晶系が指定された場合。 """ if system == "cubic": a = 1.0 / math.sqrt(beta[0]) return {"a": a, "b": a, "c": a, "alpha": 90.0, "beta": 90.0, "gamma": 90.0} if system == "tetragonal": a = 1.0 / math.sqrt(beta[0]) c = 1.0 / math.sqrt(beta[1]) return {"a": a, "b": a, "c": c, "alpha": 90.0, "beta": 90.0, "gamma": 90.0} if system == "hexagonal": a = math.sqrt(4.0 / (3.0 * beta[0])) c = 1.0 / math.sqrt(beta[1]) return {"a": a, "b": a, "c": c, "alpha": 90.0, "beta": 90.0, "gamma": 120.0} if system == "orthorhombic": a = 1.0 / math.sqrt(beta[0]) b = 1.0 / math.sqrt(beta[1]) c = 1.0 / math.sqrt(beta[2]) return {"a": a, "b": b, "c": c, "alpha": 90.0, "beta": 90.0, "gamma": 90.0} raise ValueError(system) def predicted_inv_d2(system: str, feat: Iterable[int], beta: np.ndarray) -> float: """ 概要: 結晶系、特徴量、最小二乗係数から予測されるd間隔の二乗の逆数 (inv_d2) を計算します。 詳細説明: 与えられた結晶系、hkl指数から導出される特徴量 feat、 および格子定数に相当する最小二乗係数 beta を用いて、 理論的な inv_d2 値を計算します。 引数: :param system: 結晶系名。 :type system: str :param feat: hkl指数から導出される特徴量のイテラブル。 :type feat: Iterable[int] :param beta: 最小二乗法で得られた係数のnumpy配列。 :type beta: numpy.ndarray 戻り値: :returns: 予測されるinv_d2値。 :rtype: float 例外: :raises ValueError: サポートされていない結晶系が指定された場合。 """ if system in {"cubic", "tetragonal", "hexagonal", "orthorhombic"}: return float(np.dot(np.asarray(tuple(feat), dtype=float), beta)) raise ValueError(system) def score_system(system: str, selected: list[Peak], all_peaks: list[Peak], hmax: int = 12) -> Candidate | None: """ 概要: 特定の結晶系に対して、観測ピークと理論的回折線を比較し、格子定数候補を評価します。 詳細説明: この関数は、与えられた結晶系に対し、選択された観測ピーク (selected) を用いて 格子定数 (beta) を最小二乗法で推定します。 システム行 (system_rows) から特徴量行列を構築し、観測された inv_d2 値との 組み合わせでbetaを求めます。 beta値が正である組み合わせの中から、マッチするピークの数と相対RMS誤差を スコアとして最適な候補を選定します。 全ての観測ピーク (all_peaks) に対するマッチング情報も生成します。 引数: :param system: 結晶系名。 :type system: str :param selected: 格子定数推定に使用する、フィルタリングされたPeakオブジェクトのリスト。 :type selected: list[Peak] :param all_peaks: 全ての観測されたPeakオブジェクトのリスト。 :type all_peaks: list[Peak] :param hmax: h, k, l指数の最大値。 :type hmax: int 戻り値: :returns: 最適な格子定数候補を表すCandidateオブジェクト、または見つからない場合はNone。 :rtype: Candidate | None """ rows = system_rows(system, hmax=hmax) p = len(rows[0][0]) pool = min(len(rows), 12) first_n = min(len(selected), max(p + 2, 5)) if len(selected) < p: return None ls_map = {"orthorhombic": 4, "tetragonal": 5, "cubic": 6, "hexagonal": 8} best = None obs = np.array([pk.inv_d2 for pk in selected], dtype=float) for obs_idx in __import__("itertools").combinations(range(first_n), p): y = np.array([obs[i] for i in obs_idx], dtype=float) for cand_idx in __import__("itertools").combinations(range(pool), p): X = np.array([rows[i][0] for i in cand_idx], dtype=float) try: beta = np.linalg.solve(X, y) except np.linalg.LinAlgError: continue if np.any(beta <= 0): continue pred_lines = np.array([predicted_inv_d2(system, feat, beta) for feat, _hkl in rows], dtype=float) selected_matches = [] used = set() rels = [] for pk in selected: rel = np.abs(pred_lines - pk.inv_d2) / np.maximum(pk.inv_d2, 1e-12) j = int(np.argmin(rel)) if rel[j] <= 0.08 and j not in used: used.add(j) selected_matches.append({ "two_theta_obs_deg": pk.two_theta, "inv_d2_obs": pk.inv_d2, "h": rows[j][1][0], "k": rows[j][1][1], "l": rows[j][1][2], "inv_d2_calc": float(pred_lines[j]), "rel_error_inv_d2": float(rel[j]), }) rels.append(float(rel[j])) if not selected_matches: continue rms = float(np.sqrt(np.mean(np.square(rels)))) score = (len(selected_matches), -rms) if best is None or score > best["score"]: # 全てのピークを割り当て all_matches = [] for pk in all_peaks: rel = np.abs(pred_lines - pk.inv_d2) / np.maximum(pk.inv_d2, 1e-12) j = int(np.argmin(rel)) matched = rel[j] <= 0.10 tt_calc = two_theta_from_inv_d2(float(pred_lines[j]), CU_KA1) if matched else None rel_t = abs((tt_calc - pk.two_theta) / pk.two_theta) if (matched and tt_calc and pk.two_theta != 0) else None all_matches.append({ "two_theta_obs_deg": pk.two_theta, "intensity": pk.intensity, "ka_role": pk.ka_role, "fwhm_deg": pk.fwhm_deg, "inv_d2_obs": pk.inv_d2, "h": rows[j][1][0] if matched else None, "k": rows[j][1][1] if matched else None, "l": rows[j][1][2] if matched else None, "inv_d2_calc": float(pred_lines[j]) if matched else None, "rel_error_inv_d2": float(rel[j]) if matched else None, "two_theta_calc_deg": tt_calc if matched else None, "rel_error_2theta": rel_t if matched else None, "matched": bool(matched), }) best = { "score": score, "beta": beta, "params": params_from_beta(system, beta), "selected_matches": selected_matches, "all_matches": all_matches, } if best is None: return None return Candidate( crystal_system=system, ls_code=ls_map[system], params=best["params"], score_matches=best["score"][0], score_rms_rel=float(-best["score"][1]), selected_matches=best["selected_matches"], all_matches=best["all_matches"], ) def import_lsq_module(path: Path): """ 概要: 指定されたパスからlsq_latt2モジュールを動的にインポートします。 詳細説明: lsq_latt2.py スクリプトを動的に読み込み、その中の関数やクラスを 現在の実行コンテキストで利用できるようにします。 引数: :param path: lsq_latt2.py スクリプトへのパス。 :type path: Path 戻り値: :returns: インポートされたモジュールオブジェクト。 :rtype: module 例外: :raises ImportError: モジュールをロードできない場合。 """ spec = importlib.util.spec_from_file_location("lsq_latt2_imported", path) if spec is None or spec.loader is None: raise ImportError(f"モジュール {path} をロードできません。") mod = importlib.util.module_from_spec(spec) sys.modules[spec.name] = mod spec.loader.exec_module(mod) return mod def refine_with_lsq(lsq_script: Path, candidate: Candidate, matched_rows: list[dict], wavelength: float = CU_KA1): """ 概要: 最小二乗法 (lsq_latt2) を用いて格子定数を精密化します。 詳細説明: 外部スクリプト lsq_latt2.py をインポートし、与えられた格子定数候補 (candidate) と マッチングされた回折ピーク情報 (matched_rows) を使って、格子定数の最小二乗精密化を実行します。 精密化された格子定数と標準偏差、および各ピークの精密化結果を返します。 引数: :param lsq_script: lsq_latt2.py スクリプトへのパス。 :type lsq_script: Path :param candidate: 精密化のベースとなる格子定数候補。 :type candidate: Candidate :param matched_rows: マッチングされた回折ピークの情報のリスト。 :type matched_rows: list[dict] :param wavelength: X線波長（オングストローム）。デフォルトはCU_KA1です。 :type wavelength: float 戻り値: :returns: 精密化された格子定数の要約辞書と、精密化されたマッチング結果のリストのタプル。 マッチングされた観測がない場合はNoneを返します。 :rtype: tuple[dict, list[dict]] | None """ mod = import_lsq_module(lsq_script) obs_list = [] serial = 0 for row in matched_rows: if not row.get("matched"): continue serial += 1 h, k, l = int(row["h"]), int(row["k"]), int(row["l"]) design = mod.lattice_design_row(candidate.ls_code, h, k, l) obs_list.append( mod.Observation( serial_index=serial, h=h, k=k, l=l, raw_position=float(row["two_theta_obs_deg"]), transformed_position=float(row["inv_d2_obs"]), weight=1.0, design_row=design, wavelength=wavelength, ) ) if not obs_list: return None fit = mod.weighted_least_squares(obs_list, io.StringIO()) cell = mod.derive_cell_constants(candidate.ls_code, fit.parameters, fit.parameter_sigma) refined_rows = [] for obs, fitted, resid in zip(fit.kept_observations, fit.fitted_y, fit.residuals): tt_calc = two_theta_from_inv_d2(float(fitted), wavelength) rel_inv = abs(resid) / max(obs.transformed_position, 1e-12) rel_tt = abs((tt_calc - obs.raw_position) / obs.raw_position) if tt_calc is not None and obs.raw_position != 0 else None refined_rows.append({ "serial": obs.serial_index, "h": obs.h, "k": obs.k, "l": obs.l, "two_theta_obs_deg": obs.raw_position, "two_theta_calc_deg": tt_calc, "rel_error_2theta": rel_tt, "inv_d2_obs": obs.transformed_position, "inv_d2_calc": float(fitted), "rel_error_inv_d2": float(rel_inv), "residual_inv_d2": float(resid), }) rms_inv = float(np.sqrt(np.mean([r["rel_error_inv_d2"] ** 2 for r in refined_rows]))) if refined_rows else None rms_tt = float(np.sqrt(np.mean([r["rel_error_2theta"] ** 2 for r in refined_rows if r["rel_error_2theta"] is not None]))) if refined_rows else None summary = { "crystal_system": candidate.crystal_system, "ls_code": candidate.ls_code, "a": float(cell.direct_lengths[0]), "b": float(cell.direct_lengths[1]), "c": float(cell.direct_lengths[2]), "alpha": float(cell.direct_angles_deg[0]), "beta": float(cell.direct_angles_deg[1]), "gamma": float(cell.direct_angles_deg[2]), "sigma_a": float(cell.direct_length_sigma[0]), "sigma_b": float(cell.direct_length_sigma[1]), "sigma_c": float(cell.direct_length_sigma[2]), "sigma_alpha": float(cell.direct_angle_sigma_deg[0]), "sigma_beta": float(cell.direct_angle_sigma_deg[1]), "sigma_gamma": float(cell.direct_angle_sigma_deg[2]), "n_used": len(refined_rows), "rms_rel_error_inv_d2": rms_inv, "rms_rel_error_2theta": rms_tt, } return summary, refined_rows def _canonical_column_name(name: object) -> str: """ 概要: データフレームの列名を標準的な形式に変換します。 詳細説明: 列名文字列から特殊文字を除去し、小文字に変換します。 例えば、"2Theta (deg)" は "two_theta_deg" に変換されます。 引数: :param name: 列名オブジェクト。 :type name: object 戻り値: :returns: 標準化された列名文字列。 :rtype: str """ s = str(name).strip().lower() s = s.replace("θ", "theta") s = s.replace("°", "deg") for ch in " -/.()[]{}": s = s.replace(ch, "_") while "__" in s: s = s.replace("__", "_") return s.strip("_") def _find_column(df: pd.DataFrame, candidates: set[str]) -> str | None: """ 概要: データフレーム内で候補リストに一致する列名を検索します。 詳細説明: データフレーム df の列名を標準化 (canonical_column_name) した後、 与えられた候補セット candidates のいずれかに一致するかどうかを調べます。 最初に見つかった一致する元の列名を返します。 引数: :param df: 検索対象のPandasデータフレーム。 :type df: pandas.DataFrame :param candidates: 検索する列名の候補セット。 :type candidates: set[str] 戻り値: :returns: 一致する元の列名、または見つからない場合はNone。 :rtype: str | None """ for col in df.columns: if _canonical_column_name(col) in candidates: return str(col) return None def _read_peak_dataframe(path: Path) -> pd.DataFrame: """ 概要: ピーク情報ファイルからPandasデータフレームを読み込みます。 詳細説明: Excelファイルの場合、"peaks", "all_detected_peaks", "selected_for_indexing" のいずれかのシートを優先して読み込みます。これらのシートがない場合は最初のシートを読み込みます。 CSV/テキストファイルの場合は、自動的に区切り文字を判別して読み込みます。 引数: :param path: ピーク情報ファイルへのパス。 :type path: Path 戻り値: :returns: 読み込まれたピーク情報を含むPandasデータフレーム。 :rtype: pandas.DataFrame """ suffix = path.suffix.lower() if suffix in {".xlsx", ".xls"}: xls = pd.ExcelFile(path) preferred = ["peaks", "all_detected_peaks", "selected_for_indexing"] sheet = None lower_map = {str(s).lower(): s for s in xls.sheet_names} for name in preferred: if name.lower() in lower_map: sheet = lower_map[name.lower()] break if sheet is None: sheet = xls.sheet_names[0] return pd.read_excel(path, sheet_name=sheet) return pd.read_csv(path, sep=None, engine="python", comment="#") def read_peaks_file(path: Path, wavelength: float = CU_KA1) -> list[Peak]: """ 概要: ピーク情報ファイルからPeakオブジェクトのリストを読み込みます。 詳細説明: ExcelまたはCSV形式のピーク情報ファイルを読み込み、Pandasデータフレームとして扱います。 "two_theta_deg", "2theta", "inv_d2"などの一般的な列名から2θ角度またはinv_d2値を抽出し、 強度、FWHM、Kα役割などの情報をPeakオブジェクトに変換します。 ファイルにKα役割の割り当てがない場合は、assign_ka2関数を呼び出して割り当てを試みます。 引数: :param path: ピーク情報ファイルへのパス。 :type path: Path :param wavelength: X線波長（オングストローム）。デフォルトはCU_KA1です。 :type wavelength: float 戻り値: :returns: 読み込まれたPeakオブジェクトのリスト。 :rtype: list[Peak] 例外: :raises ValueError: ファイルが空である場合、または回折角度の列が見つからない場合、 または利用可能なピークが見つからない場合。 """ df = _read_peak_dataframe(path) if df.empty: raise ValueError(f"ピークファイルが空です: {path}") two_theta_col = _find_column(df, { "two_theta_deg", "two_theta", "twotheta", "2theta", "2theta_deg", "theta2", "angle", "angle_deg", "position", "position_deg", "two_theta_obs_deg", }) # "inv_d2_A-2" の標準形式は "inv_d2_a_2" となります。 inv_d2_col = _find_column(df, { "inv_d2", "inv_d2_obs", "inv_d2_calc", "inv_d2_a_2", "inv_d2_a2", "one_over_d2", "one_over_d_2", "d_star2", }) intensity_col = _find_column(df, {"intensity", "i", "counts", "count", "y"}) intensity_raw_col = _find_column(df, {"intensity_raw", "raw_intensity", "i_raw", "counts_raw"}) fwhm_col = _find_column(df, {"fwhm_deg", "fwhm", "fwhm_degree"}) role_col = _find_column(df, {"ka_role", "role", "kalpha_role"}) pair_col = _find_column(df, {"ka_pair_peak_id", "pair_peak_id", "ka_pair_index"}) if two_theta_col is None and inv_d2_col is None: # 最終的なフォールバック: 最初の数値列を2θとして使用します。 for col in df.columns: values = pd.to_numeric(df[col], errors="coerce") if values.notna().sum() > 0: two_theta_col = str(col) break if two_theta_col is None and inv_d2_col is None: raise ValueError( "ピークファイルに回折角度の列が見つかりません。 " "two_theta_deg, 2theta, angle, または inv_d2_A-2 を使用してください。" ) peaks: list[Peak] = [] for row_index, row in df.iterrows(): two_theta = float("nan") inv_d2 = float("nan") if two_theta_col is not None: two_theta = float(pd.to_numeric(row.get(two_theta_col), errors="coerce")) if inv_d2_col is not None: inv_d2 = float(pd.to_numeric(row.get(inv_d2_col), errors="coerce")) if not math.isfinite(two_theta): if math.isfinite(inv_d2) and inv_d2 > 0: tt = two_theta_from_inv_d2(inv_d2, wavelength) two_theta = float(tt) if tt is not None else float("nan") if not math.isfinite(inv_d2) and math.isfinite(two_theta): inv_d2 = inv_d2_from_two_theta(two_theta, wavelength) if not (math.isfinite(two_theta) and math.isfinite(inv_d2) and inv_d2 > 0): continue intensity = float(pd.to_numeric(row.get(intensity_col), errors="coerce")) if intensity_col else 1.0 if not math.isfinite(intensity): intensity = 1.0 intensity_raw = float(pd.to_numeric(row.get(intensity_raw_col), errors="coerce")) if intensity_raw_col else intensity if not math.isfinite(intensity_raw): intensity_raw = intensity fwhm = float(pd.to_numeric(row.get(fwhm_col), errors="coerce")) if fwhm_col else float("nan") if not math.isfinite(fwhm): fwhm = 0.0 d = 1.0 / math.sqrt(inv_d2) q = 2.0 * math.pi / d role = "" if role_col is not None and not pd.isna(row.get(role_col)): role = str(row.get(role_col)).strip().lower() pair_index = None if pair_col is not None and not pd.isna(row.get(pair_col)): try: pair_id = int(row.get(pair_col)) if pair_id > 0: pair_index = pair_id - 1 except Exception: pair_index = None peaks.append( Peak( index=int(row_index), two_theta=float(two_theta), intensity=float(intensity), intensity_raw=float(intensity_raw), fwhm_deg=float(fwhm), inv_d2=float(inv_d2), d=float(d), q=float(q), ka_role=role, ka_pair_index=pair_index, source=str(path), ) ) peaks.sort(key=lambda p: p.two_theta) if not peaks: raise ValueError(f"'{path}'から利用可能なピークが見つかりませんでした。") if not any(p.ka_role for p in peaks): assign_ka2(peaks) return peaks def normalize_mode(mode: str) -> str: """ 概要: 入力されたモード文字列を標準的なモード名に変換します。 詳細説明: モード文字列から前後の空白を除去し、小文字に変換します。 スペルミス "seaerch" は "search" に修正され、古い "index" モードは "all" (search + guess) に変換されます。 サポートされていないモードが指定された場合はエラーを発生させます。 引数: :param mode: 入力されたモード文字列。 :type mode: str 戻り値: :returns: 標準化されたモード名。 :rtype: str 例外: :raises ValueError: サポートされていないモードが指定された場合。 """ m = mode.strip().lower() if m == "seaerch": return "search" if m == "index": # 後方互換性: 古い "index" モードは search + indexing を実行していました。 return "all" if m not in {"search", "guess", "all", "compare", "calc"}: raise ValueError(f"サポートされていないモード: {mode}") return m def parse_crystal_systems(text: str) -> list[str]: """ 概要: 結晶系指定の文字列をパースし、サポートされている結晶系のリストを返します。 詳細説明: カンマまたはセミコロンで区切られた結晶系名の文字列を処理します。 "auto" や "all" が指定された場合、全てのサポートされている結晶系 ("hexagonal", "orthorhombic", "tetragonal", "cubic") のリストを返します。エイリアスもサポートされます (例: "hex" -> "hexagonal")。 サポートされていない結晶系が指定された場合はエラーを発生させます。 引数: :param text: 結晶系指定の文字列。 :type text: str 戻り値: :returns: サポートされている結晶系のリスト。 :rtype: list[str] 例外: :raises ValueError: サポートされていない結晶系が指定された場合。 """ systems: list[str] = [] for raw in str(text).replace(";", ",").split(","): key = raw.strip().lower() if not key: continue if key not in CRYSTAL_SYSTEM_ALIASES: allowed = ", ".join(["auto"] + SUPPORTED_CRYSTAL_SYSTEMS) raise ValueError(f"サポートされていない結晶系: {raw}。許可されているもの: {allowed}") value = CRYSTAL_SYSTEM_ALIASES[key] if value == "auto": return list(SUPPORTED_CRYSTAL_SYSTEMS) if value not in systems: systems.append(value) return systems or list(SUPPORTED_CRYSTAL_SYSTEMS) def default_peak_file_for_input(infile: Path) -> Path: """ 概要: 入力ファイル名に基づいてデフォルトのピークファイル名を生成します。 詳細説明: 入力ファイルのステムをサニタイズし、"-peaks.xlsx" を追加した Pathオブジェクトを返します。 引数: :param infile: 入力データファイルへのパス。 :type infile: Path 戻り値: :returns: デフォルトのピークファイルパス。 :rtype: Path """ stem = sanitize_stem(infile.stem) return infile.with_name(f"{stem}-peaks.xlsx") def strip_peak_suffix(stem: str) -> str: """ 概要: ファイル名のステムからピーク関連のサフィックスを除去します。 詳細説明: ファイル名のステム (例: "sample-peaks") から、"-peaks", "_peaks", ".peaks", "-peak", "_peak" などのサフィックスを除去した文字列を返します。 引数: :param stem: ファイル名のステム文字列。 :type stem: str 戻り値: :returns: サフィックスが除去されたステム文字列。 :rtype: str """ for suffix in ("-peaks", "_peaks", ".peaks", "-peak", "_peak"): if stem.lower().endswith(suffix): return stem[: -len(suffix)] return stem def default_guess_file_for_peak_file(peak_file: Path, base_input: Path | None = None) -> Path: """ 概要: ピークファイル名に基づいてデフォルトの推定結果ファイル名を生成します。 詳細説明: 基本入力ファイル (base_input) が指定されている場合はそれを使用し、 そうでない場合はピークファイル名からピーク関連のサフィックスを除去し、 "-guess.xlsx" を追加したPathオブジェクトを返します。 引数: :param peak_file: ピーク情報ファイルへのパス。 :type peak_file: Path :param base_input: 元の入力データファイルへのパス（オプション）。 :type base_input: Path | None 戻り値: :returns: デフォルトの推定結果ファイルパス。 :rtype: Path """ if base_input is not None: stem = sanitize_stem(base_input.stem) return base_input.with_name(f"{stem}-guess.xlsx") stem = sanitize_stem(strip_peak_suffix(peak_file.stem)) return peak_file.with_name(f"{stem}-guess.xlsx") def peaks_to_frame(peaks: list[Peak]) -> pd.DataFrame: """ 概要: PeakオブジェクトのリストをPandasデータフレームに変換します。 詳細説明: Peakオブジェクトの各属性をデータフレームの列にマッピングし、 各ピークに一意のID ("peak_id") を割り当てます。 引数: :param peaks: Peakオブジェクトのリスト。 :type peaks: list[Peak] 戻り値: :returns: Peak情報を含むPandasデータフレーム。 :rtype: pandas.DataFrame """ rows = [] for i, p in enumerate(peaks, 1): rows.append({ "peak_id": i, "two_theta_deg": p.two_theta, "intensity": p.intensity, "intensity_raw": p.intensity_raw, "fwhm_deg": p.fwhm_deg, "d_A": p.d, "inv_d2_A-2": p.inv_d2, "q_A-1": p.q, "ka_role": p.ka_role, "ka_pair_peak_id": (p.ka_pair_index + 1) if p.ka_pair_index is not None else None, }) return pd.DataFrame(rows) def print_peak_table(peaks: list[Peak], title: str) -> None: """ 概要: ピーク情報を整形されたテーブル形式でコンソールに出力します。 詳細説明: 指定されたタイトルとともに、各ピークのID、2θ角度、強度、FWHM、Kα役割を 見やすい形式で標準出力に表示します。 引数: :param peaks: Peakオブジェクトのリスト。 :type peaks: list[Peak] :param title: テーブルのタイトル。 :type title: str """ print("") print(title) print(f"{'id':>3s} {'2theta':>9s} {'I':>11s} {'FWHM':>7s} {'role':>5s}") for i, p in enumerate(peaks, 1): print(f"{i:3d} {p.two_theta:9.3f} {p.intensity:11.2f} {p.fwhm_deg:7.3f} {p.ka_role or '-':>5s}") def save_workbook(path: Path, sheets: dict[str, pd.DataFrame]) -> None: """ 概要: 複数のPandasデータフレームをExcelワークブックとして保存します。 詳細説明: 辞書で指定されたシート名とデータフレームのペアを、 OpenPyXLライブラリを使用して単一のExcelファイルに書き込みます。 各シートの最初の行は列ヘッダーとして太字で表示され、行が固定されます。 欠損値はExcelの空セルとして扱われます。 引数: :param path: 保存先のExcelファイルパス。 :type path: Path :param sheets: シート名とPandasデータフレームの辞書。 :type sheets: dict[str, pandas.DataFrame] """ wb = Workbook() first = True for name, df in sheets.items(): ws = wb.active if first else wb.create_sheet(title=name[:31]) ws.title = name[:31] first = False for c, col in enumerate(df.columns, 1): cell = ws.cell(1, c, col) cell.font = Font(bold=True) for r, (_, row) in enumerate(df.iterrows(), 2): for c, val in enumerate(row, 1): ws.cell(r, c, None if (pd.isna(val) if not isinstance(val, str) else False) else val) ws.freeze_panes = "A2" wb.save(path) def _maybe_install_mplcursors(artists: list, hover_texts: list[str]) -> None: """ 概要: mplcursorsが利用可能な場合、Matplotlibのプロットにホバー注釈をアタッチします。 詳細説明: mplcursorsライブラリがインストールされていれば、指定されたアーティストに 対応するテキストがツールチップとして表示されるように設定します。 ライブラリがない場合は、メッセージを出力し機能は無効化されます。 引数: :param artists: Matplotlibアーティストオブジェクトのリスト。 :type artists: list :param hover_texts: 各アーティストに対応するホバーテキストのリスト。 :type hover_texts: list[str] """ if not artists: return try: import mplcursors # type: ignore except Exception: print("mplcursorsがインストールされていません。ホバー注釈は無効です。") print("以下でインストールできます: pip install mplcursors") return text_by_artist = {artist: text for artist, text in zip(artists, hover_texts)} cursor = mplcursors.cursor(artists, hover=True) @cursor.connect("add") def _on_add(sel): # pragma: no cover - interactive callback sel.annotation.set_text(text_by_artist.get(sel.artist, "")) sel.annotation.get_bbox_patch().set(alpha=0.92) def _peak_hover_text(p: Peak, hkl_text: str | None = None, calc_two_theta: float | None = None, resid: float | None = None) -> str: """ 概要: ピークのホバーテキストを生成します。 詳細説明: Peakオブジェクトの情報と、オプションでhkl指数、計算された2θ、残差を用いて、 Matplotlibのホバー注釈として表示する整形されたテキスト文字列を作成します。 引数: :param p: Peakオブジェクト。 :type p: Peak :param hkl_text: hkl指数を示す文字列（オプション）。 :type hkl_text: str | None :param calc_two_theta: 計算された2θ角度（オプション）。 :type calc_two_theta: float | None :param resid: 残差（オプション）。 :type resid: float | None 戻り値: :returns: ホバー注釈用の整形済みテキスト。 :rtype: str """ parts = [] if hkl_text: parts.append(f"hkl: {hkl_text}") parts.append(f"2theta(obs): {p.two_theta:.5f} deg") if calc_two_theta is not None and math.isfinite(calc_two_theta): parts.append(f"2theta(calc): {calc_two_theta:.5f} deg") if resid is not None and math.isfinite(resid): parts.append(f"resid: {resid:+.5f} deg") parts.append(f"intensity: {p.intensity:.6g}") parts.append(f"FWHM: {p.fwhm_deg:.5g} deg") if p.ka_role: parts.append(f"role: {p.ka_role}") return "\n".join(parts) def plot_results( x: np.ndarray, y: np.ndarray, info: dict[str, np.ndarray | float], peaks: list[Peak], outpath: Path | None, show: bool = False, save: bool = True, ) -> None: """ 概要: 入力データ、スムージングされたデータ、検出されたピークをプロットします。 詳細説明: 生のXRDデータ (x, y) とスムージングされたデータ (info["ysmooth"]) を線グラフで表示し、 検出されたピークを縦線で示します。Kα2ピークは異なる色で強調表示されます。 グラフは画像ファイルとして保存でき、Matplotlibウィンドウに表示することも可能です。 mplcursorsがインストールされていれば、ピークの縦線にホバー注釈が表示されます。 引数: :param x: 2θ角度のnumpy配列。 :type x: numpy.ndarray :param y: 強度のnumpy配列。 :type y: numpy.ndarray :param info: ピーク探索の診断情報を含む辞書（スムージングされたデータなど）。 :type info: dict[str, numpy.ndarray | float] :param peaks: 検出されたPeakオブジェクトのリスト。 :type peaks: list[Peak] :param outpath: プロットを保存するファイルパス（Noneの場合は保存しない）。 :type outpath: Path | None :param show: プロットウィンドウを表示するかどうか。 :type show: bool :param save: プロットをファイルに保存するかどうか。 :type save: bool """ if not show and not save: return ys = info["ysmooth"] fig, ax = plt.subplots(figsize=(12, 7)) ax.plot(x, y, lw=0.8, label="入力データ") ax.plot(x, ys, lw=1.0, label="スムージング後") hover_artists = [] hover_texts = [] ymin = float(np.nanmin(y)) if len(y) else 0.0 ymax = float(np.nanmax([np.nanmax(y), np.nanmax(ys)])) if len(y) else 1.0 for p in peaks: color = "tab:red" if p.ka_role == "ka2" else "tab:green" line = ax.axvline(p.two_theta, color=color, lw=0.9, alpha=0.85, picker=5) line.set_gid(f"peak_{p.two_theta:.5f}") hover_artists.append(line) hover_texts.append(_peak_hover_text(p)) ax.text(p.two_theta, p.intensity, f"{p.two_theta:.2f}", rotation=90, va="bottom", fontsize=7) ax.set_ylim(min(ymin, 0.0), ymax * 1.05 if ymax > 0 else ymax + 1.0) ax.set_xlabel("2Theta (度)") ax.set_ylabel("強度") ax.legend() fig.tight_layout() if save and outpath is not None: fig.savefig(outpath, dpi=160) if show: _maybe_install_mplcursors(hover_artists, hover_texts) plt.show() # plt.pause(0.1) # input("\nPress ENTER to terminate>>\n") if save and outpath is not None: plt.close(fig) def beta_from_params(system: str, params: dict[str, float]) -> np.ndarray: """ 概要: 格子定数の辞書から最小二乗法の係数 (beta) を計算します。 詳細説明: 与えられた結晶系と格子定数 (a, b, c) の辞書から、 inv_d2 (1 / d^2) の計算に必要な最小二乗係数 beta をnumpy配列として導出します。 引数: :param system: 結晶系名。 :type system: str :param params: 格子定数を含む辞書。 :type params: dict[str, float] 戻り値: :returns: beta係数のnumpy配列。 :rtype: numpy.ndarray 例外: :raises ValueError: サポートされていない結晶系が指定された場合。 """ system = system.lower() if system == "cubic": return np.array([1.0 / (params["a"] ** 2)], dtype=float) if system == "tetragonal": return np.array([1.0 / (params["a"] ** 2), 1.0 / (params["c"] ** 2)], dtype=float) if system == "hexagonal": return np.array([4.0 / (3.0 * params["a"] ** 2), 1.0 / (params["c"] ** 2)], dtype=float) if system == "orthorhombic": return np.array([1.0 / (params["a"] ** 2), 1.0 / (params["b"] ** 2), 1.0 / (params["c"] ** 2)], dtype=float) raise ValueError(system) def candidate_to_summary_row(rank: int, c: Candidate, method: str = "") -> dict[str, object]: """ 概要: Candidateオブジェクトをサマリー行（辞書）に変換します。 詳細説明: Candidateオブジェクトの情報を、Excelなどのサマリーシートに適した辞書形式に変換します。 ランキング、メソッド、結晶系、マッチ数、RMS誤差、格子定数を含みます。 引数: :param rank: 候補の順位。 :type rank: int :param c: Candidateオブジェクト。 :type c: Candidate :param method: 推定に使用されたメソッド（オプション）。 :type method: str 戻り値: :returns: サマリー行を表す辞書。 :rtype: dict[str, object] """ row: dict[str, object] = { "rank": rank, "method": method, "crystal_system": c.crystal_system, "score_matches": c.score_matches, "score_rms_rel": c.score_rms_rel, } row.update(c.params) return row def candidate_from_summary_row(row: pd.Series) -> Candidate: """ 概要: サマリー行（Pandas Series）からCandidateオブジェクトを生成します。 詳細説明: Excelファイルなどから読み込まれたサマリー行のデータを使って、 Candidateオブジェクトを再構築します。 結晶系と格子定数 (a, b, c, alpha, beta, gamma) を抽出し、 不足している角度はデフォルト値で補完します。 引数: :param row: サマリー行を表すPandas Series。 :type row: pandas.Series 戻り値: :returns: 生成されたCandidateオブジェクト。 :rtype: Candidate 例外: :raises ValueError: サポートされていない結晶系、または必要な格子定数が見つからない場合。 """ system = str(row.get("crystal_system", row.get("system", ""))).strip().lower() if system not in SUPPORTED_CRYSTAL_SYSTEMS: raise ValueError(f"候補サマリーにサポートされていない結晶系: {system}") params = { key: float(row[key]) for key in ["a", "b", "c", "alpha", "beta", "gamma"] if key in row.index and pd.notna(row[key]) } if "a" not in params: raise ValueError("候補サマリーには格子定数aが含まれていません。") if system in {"tetragonal", "hexagonal"} and "c" not in params: raise ValueError("候補サマリーには格子定数cが含まれていません。") if system == "orthorhombic" and not all(k in params for k in ["b", "c"]): raise ValueError("候補サマリーには格子定数bとcが含まれていません。") params.setdefault("b", params.get("a", float("nan"))) params.setdefault("c", params.get("a", float("nan"))) params.setdefault("alpha", 90.0) params.setdefault("beta", 90.0) params.setdefault("gamma", 120.0 if system == "hexagonal" else 90.0) ls_map = {"orthorhombic": 4, "tetragonal": 5, "cubic": 6, "hexagonal": 8} return Candidate( crystal_system=system, ls_code=ls_map[system], params=params, score_matches=int(row.get("score_matches", 0) or 0), score_rms_rel=float(row.get("score_rms_rel", float("nan"))), selected_matches=[], all_matches=[], ) def has_explicit_lattice_params(args: argparse.Namespace) -> bool: """ 概要: コマンドライン引数に明示的な格子定数が指定されているか確認します。 詳細説明: args.a, args.b, args.c のいずれかがNoneでない場合にTrueを返します。 引数: :param args: コマンドライン引数を格納したNamespaceオブジェクト。 :type args: argparse.Namespace 戻り値: :returns: 明示的な格子定数が指定されていればTrue、そうでなければFalse。 :rtype: bool """ return any(getattr(args, name, None) is not None for name in ["a", "b", "c"]) def candidate_from_lattice_args(args: argparse.Namespace) -> Candidate: """ 概要: コマンドライン引数で指定された格子定数からCandidateオブジェクトを生成します。 詳細説明: --a, --b, --c, --alpha, --beta, --gamma などの引数から 格子定数情報を抽出し、Candidateオブジェクトを作成します。 指定された結晶系に基づいて、必要な格子定数が全て提供されているかを検証します。 引数: :param args: コマンドライン引数を格納したNamespaceオブジェクト。 :type args: argparse.Namespace 戻り値: :returns: 生成されたCandidateオブジェクト。 :rtype: Candidate 例外: :raises ValueError: 結晶系の指定が適切でない、必要な格子定数が不足している、 または格子定数が正の値でない場合。 """ systems = parse_crystal_systems(args.crystal_system) if len(systems) != 1: raise ValueError("--crystal-system には auto またはカンマ区切りのリストではなく、単一の結晶系が必要です。") system = systems[0] a = getattr(args, "a", None) b = getattr(args, "b", None) c = getattr(args, "c", None) if a is None: raise ValueError("明示的な格子比較には --a が必要です。") params: dict[str, float] = { "a": float(a), "alpha": float(getattr(args, "alpha", 90.0) if getattr(args, "alpha", None) is not None else 90.0), "beta": float(getattr(args, "beta", 90.0) if getattr(args, "beta", None) is not None else 90.0), "gamma": float(getattr(args, "gamma", 120.0 if system == "hexagonal" else 90.0) if getattr(args, "gamma", None) is not None else (120.0 if system == "hexagonal" else 90.0)), } if system == "cubic": params["b"] = float(a) params["c"] = float(a if c is None else c) # 内部的には --b/--c が誤って供給されても立方晶として維持します。 params["b"] = params["a"] params["c"] = params["a"] elif system in {"tetragonal", "hexagonal"}: if c is None: raise ValueError(f"{system} の比較には --a に加えて --c が必要です。") params["b"] = float(a) params["c"] = float(c) if system == "hexagonal": params["gamma"] = 120.0 elif system == "orthorhombic": if b is None or c is None: raise ValueError("orthorhombic の比較には --a, --b, および --c が必要です。") params["b"] = float(b) params["c"] = float(c) else: raise ValueError(f"サポートされていない結晶系: {system}") for key in ["a", "b", "c"]: if not math.isfinite(params[key]) or params[key] <= 0: raise ValueError(f"格子定数 {key} は正の値でなければなりません: {params[key]}") ls_map = {"orthorhombic": 4, "tetragonal": 5, "cubic": 6, "hexagonal": 8} return Candidate( crystal_system=system, ls_code=ls_map[system], params=params, score_matches=0, score_rms_rel=float("nan"), selected_matches=[], all_matches=[], ) def filter_theoretical_lines_for_plot( theoretical_df: pd.DataFrame, mode: str = "near", near_deg: float = 0.5, ) -> pd.DataFrame: """ 概要: プロット表示のために理論的な回折線をフィルタリングします。 詳細説明: プロットモード (mode) に応じて、理論的な回折線データフレーム (theoretical_df) をフィルタリングします。 - "all": 全ての理論線を表示。 - "none": 全ての理論線を表示しない。 - "matched": 観測ピークにマッチした理論線のみ表示。 - "near": 観測ピークの近く（near_deg以内）にある理論線のみ表示。 引数: :param theoretical_df: 理論的な回折線情報を含むPandasデータフレーム。 :type theoretical_df: pandas.DataFrame :param mode: フィルタリングモード（"near", "matched", "all", "none"）。 :type mode: str :param near_deg: "near"モードで使用される2θ角度の許容範囲（度）。 :type near_deg: float 戻り値: :returns: フィルタリングされた理論回折線データフレーム。 :rtype: pandas.DataFrame 例外: :raises ValueError: サポートされていないプロットモードが指定された場合。 """ if theoretical_df.empty: return theoretical_df mode = str(mode).strip().lower() if mode == "all": return theoretical_df if mode == "none": return theoretical_df.iloc[0:0].copy() if mode == "matched": if "matched" not in theoretical_df.columns: return theoretical_df.iloc[0:0].copy() return theoretical_df[theoretical_df["matched"].astype(bool)].copy() if mode == "near": if "nearest_resid_two_theta_deg" not in theoretical_df.columns: if "matched" in theoretical_df.columns: return theoretical_df[theoretical_df["matched"].astype(bool)].copy() return theoretical_df.iloc[0:0].copy() resid = pd.to_numeric(theoretical_df["nearest_resid_two_theta_deg"], errors="coerce").abs() return theoretical_df[resid <= float(near_deg)].copy() raise ValueError(f"サポートされていない --plot-theory オプション: {mode}") def estimate_lattice_limit_from_peaks( peaks: list[Peak], wavelength: float = CU_KA1, factor: float = 1.2, ) -> tuple[float, float, float]: """ 概要: 観測ピークに基づいて格子定数の上限を推定します。 詳細説明: 観測されたピークの中で最も小さい2θ角度（Kα2ピークを除く）から、 対応する最大のd間隔を計算します。このd間隔に安全係数 (factor) を乗じて、 格子定数推定の妥当な上限値を決定します。この上限値は、低角度での回折線が 弱い、または欠落している可能性があるため、少し余裕を持たせるように設計されています。 引数: :param peaks: Peakオブジェクトのリスト。 :type peaks: list[Peak] :param wavelength: X線波長（オングストローム）。デフォルトはCU_KA1です。 :type wavelength: float :param factor: d間隔に乗じる安全係数。 :type factor: float 戻り値: :returns: (最小2θ角度, その角度でのd間隔, 推定された格子定数の上限) のタプル。 :rtype: tuple[float, float, float] 例外: :raises ValueError: 有効なピーク角度がない場合、またはd間隔の計算が不可能な場合。 """ usable = [ p for p in peaks if p.ka_role != "ka2" and math.isfinite(p.two_theta) and p.two_theta > 0.0 ] if not usable: usable = [p for p in peaks if math.isfinite(p.two_theta) and p.two_theta > 0.0] if not usable: raise ValueError("有限の正のピーク角度が利用できないため、格子定数の上限を推定できません。") min_two_theta = min(p.two_theta for p in usable) dmax = bragg_d(min_two_theta, wavelength) if not math.isfinite(dmax) or dmax <= 0.0: raise ValueError(f"最小の2θ角度={min_two_theta}から格子定数の上限を推定できません。") return float(min_two_theta), float(dmax), float(dmax * factor) def resolve_lattice_max_limit(args: argparse.Namespace, peaks: list[Peak]) -> float | None: """ 概要: 格子定数推定に使用される格子定数の上限値を解決します。 詳細説明: コマンドライン引数で --lattice-max が明示的に指定されている場合、その値を使用します。 指定がない場合、最も低角度の非Kα2ピークのd間隔に --lattice-max-factor を乗じて上限を推定します。 --lattice-max が0以下に設定されている場合、制限は無効化されます。 引数: :param args: コマンドライン引数を格納したNamespaceオブジェクト。 :type args: argparse.Namespace :param peaks: 観測されたPeakオブジェクトのリスト。 :type peaks: list[Peak] 戻り値: :returns: 解決された格子定数上限値（オングストローム）、または制限が無効化されている場合はNone。 :rtype: float | None 例外: :raises ValueError: --lattice-max-factor が正の値でない場合。 """ explicit = getattr(args, "lattice_max", None) if explicit is not None: explicit = float(explicit) if explicit <= 0.0: print("格子定数上限 : 無効 (--lattice-max <= 0)") return None print(f"格子定数上限 : {explicit:.6g} Å (ユーザー指定)") return explicit factor = float(getattr(args, "lattice_max_factor", 1.2)) if factor <= 0.0: raise ValueError("--lattice-max-factor は正の値でなければなりません。") min_tt, dmax, limit = estimate_lattice_limit_from_peaks(peaks, wavelength=args.wavelength, factor=factor) print(f"最小ピーク角度: {min_tt:.6g} 度 2θ") print(f"d(最小2θ) : {dmax:.6g} Å") print(f"格子定数上限 : {limit:.6g} Å (= d(最小2θ) x {factor:.6g})") return limit def params_within_lattice_max(params: dict[str, float], system: str, lattice_max: float | None) -> bool: """ 概要: 格子定数セットが指定された上限値内に収まっているかを確認します。 詳細説明: lattice_max がNoneでない場合、与えられた結晶系の主要な格子定数 (a, b, c) が lattice_max の値を超えていないかを検証します。 引数: :param params: 格子定数を含む辞書。 :type params: dict[str, float] :param system: 結晶系名。 :type system: str :param lattice_max: 格子定数の上限値（オングストローム）、またはNone。 :type lattice_max: float | None 戻り値: :returns: 全ての格子定数が上限値内であればTrue、そうでなければFalse。 :rtype: bool """ if lattice_max is None: return True keys = ["a"] if system in {"tetragonal", "hexagonal"}: keys = ["a", "c"] elif system == "orthorhombic": keys = ["a", "b", "c"] for key in keys: val = params.get(key, None) if val is None: continue try: v = float(val) except Exception: continue if math.isfinite(v) and v > lattice_max * (1.0 + 1.0e-12): return False return True def filter_candidates_by_lattice_max( candidates: list[Candidate], lattice_max: float | None, label: str = "candidates", is_print: bool = True, ) -> list[Candidate]: """ 概要: 格子定数候補を格子定数上限値でフィルタリングします。 詳細説明: lattice_max がNoneでない場合、その値を超える格子定数を持つCandidateオブジェクトを 候補リストから除外します。 除外された候補の数とラベルは、is_printがTrueの場合にコンソールに出力されます。 引数: :param candidates: Candidateオブジェクトのリスト。 :type candidates: list[Candidate] :param lattice_max: 格子定数の上限値（オングストローム）、またはNone。 :type lattice_max: float | None :param label: 出力メッセージに使用する候補のラベル。 :type label: str :param is_print: フィルタリング情報をコンソールに出力するかどうか。 :type is_print: bool 戻り値: :returns: フィルタリングされたCandidateオブジェクトのリスト。 :rtype: list[Candidate] """ if lattice_max is None: return candidates kept = [c for c in candidates if params_within_lattice_max(c.params, c.crystal_system, lattice_max)] n_removed = len(candidates) - len(kept) if is_print and n_removed > 0: print(f"格子定数上限 ({lattice_max:.6g} Å) で除外された候補: {n_removed} {label}") return kept def theoretical_lines_for_candidate( candidate: Candidate, hmax: int = 12, wavelength: float = CU_KA1, xmin: float | None = None, xmax: float | None = None, ) -> pd.DataFrame: """ 概要: 指定された格子定数候補に基づいて理論的な回折線を計算します。 詳細説明: Candidateオブジェクトの格子定数とhmaxで定義されるhkl指数範囲を用いて、 それぞれのhklに対する理論的な2θ角度、d間隔、inv_d2値を計算します。 結果はPandasデータフレームとして返され、オプションで2θ範囲 (xmin, xmax) でフィルタリングされます。 引数: :param candidate: 格子定数候補。 :type candidate: Candidate :param hmax: h, k, l指数の最大値。 :type hmax: int :param wavelength: X線波長（オングストローム）。デフォルトはCU_KA1です。 :type wavelength: float :param xmin: 最小2θ角度（度）のフィルタリング閾値（オプション）。 :type xmin: float | None :param xmax: 最大2θ角度（度）のフィルタリング閾値（オプション）。 :type xmax: float | None 戻り値: :returns: 理論的な回折線情報を含むPandasデータフレーム。 :rtype: pandas.DataFrame """ beta = beta_from_params(candidate.crystal_system, candidate.params) rows = [] for feat, (h, k, l) in system_rows(candidate.crystal_system, hmax=hmax): inv_d2 = predicted_inv_d2(candidate.crystal_system, feat, beta) tt = two_theta_from_inv_d2(inv_d2, wavelength) if tt is None or not math.isfinite(tt): continue if xmin is not None and tt < xmin: continue if xmax is not None and tt > xmax: continue d = 1.0 / math.sqrt(inv_d2) rows.append({ "h": h, "k": k, "l": l, "hkl": f"{h}{k}{l}", "two_theta_calc_deg": float(tt), "d_calc_A": float(d), "inv_d2_calc": float(inv_d2), }) return pd.DataFrame(rows).sort_values("two_theta_calc_deg").reset_index(drop=True) def compare_candidate_with_peaks( candidate: Candidate, peaks: list[Peak], hmax: int = 12, wavelength: float = CU_KA1, xmin: float | None = None, xmax: float | None = None, tolerance_deg: float = 0.25, ) -> tuple[pd.DataFrame, pd.DataFrame]: """ 概要: 格子定数候補と観測ピークを比較し、割り当てを行います。 詳細説明: 与えられた格子定数候補 (candidate) と観測ピーク (peaks) を比較し、 理論的な回折線と観測ピークを割り当てます。 各観測ピークに対して最も近い理論的な回折線を見つけ、 指定された許容誤差 (tolerance_deg) 内であればマッチとみなします。 理論的回折線のデータフレームと、観測ピークと理論線の割り当て結果のデータフレームを返します。 引数: :param candidate: 格子定数候補。 :type candidate: Candidate :param peaks: 観測されたPeakオブジェクトのリスト。 :type peaks: list[Peak] :param hmax: h, k, l指数の最大値。 :type hmax: int :param wavelength: X線波長（オングストローム）。デフォルトはCU_KA1です。 :type wavelength: float :param xmin: 最小2θ角度（度）のフィルタリング閾値（オプション）。 :type xmin: float | None :param xmax: 最大2θ角度（度）のフィルタリング閾値（オプション）。 :type xmax: float | None :param tolerance_deg: 2θ角度の割り当て許容誤差（度）。 :type tolerance_deg: float 戻り値: :returns: (理論的な回折線データフレーム, 観測ピークと理論線の割り当て結果データフレーム) のタプル。 :rtype: tuple[pandas.DataFrame, pandas.DataFrame] """ if xmin is None and peaks: xmin = max(0.0, min(p.two_theta for p in peaks) - 2.0) if xmax is None and peaks: xmax = max(p.two_theta for p in peaks) + 2.0 calc_df = theoretical_lines_for_candidate(candidate, hmax=hmax, wavelength=wavelength, xmin=xmin, xmax=xmax) if calc_df.empty: return calc_df, pd.DataFrame() calc_tt = calc_df["two_theta_calc_deg"].to_numpy(dtype=float) rows = [] used_calc: set[int] = set() for peak_id, p in enumerate(sorted(peaks, key=lambda pk: pk.two_theta), 1): diffs = calc_tt - p.two_theta order = np.argsort(np.abs(diffs)) best_idx = None for idx in order: if int(idx) not in used_calc: best_idx = int(idx) break if best_idx is None: continue used_calc.add(best_idx) calc = calc_df.iloc[best_idx] resid = p.two_theta - float(calc["two_theta_calc_deg"]) matched = abs(resid) <= tolerance_deg rows.append({ "peak_id": peak_id, "two_theta_obs_deg": p.two_theta, "intensity": p.intensity, "fwhm_deg": p.fwhm_deg, "ka_role": p.ka_role, "h": int(calc["h"]), "k": int(calc["k"]), "l": int(calc["l"]), "hkl": str(calc["hkl"]), "two_theta_calc_deg": float(calc["two_theta_calc_deg"]), "resid_two_theta_deg": float(resid), "abs_resid_two_theta_deg": float(abs(resid)), "matched": bool(matched), }) assign_df = pd.DataFrame(rows) calc_df = calc_df.copy() if not assign_df.empty: obs_by_hkl = assign_df.set_index("hkl") calc_df["nearest_peak_id"] = calc_df["hkl"].map(obs_by_hkl["peak_id"].to_dict()) calc_df["nearest_two_theta_obs_deg"] = calc_df["hkl"].map(obs_by_hkl["two_theta_obs_deg"].to_dict()) calc_df["nearest_intensity"] = calc_df["hkl"].map(obs_by_hkl["intensity"].to_dict()) calc_df["nearest_fwhm_deg"] = calc_df["hkl"].map(obs_by_hkl["fwhm_deg"].to_dict()) calc_df["nearest_resid_two_theta_deg"] = calc_df["hkl"].map(obs_by_hkl["resid_two_theta_deg"].to_dict()) calc_df["matched"] = calc_df["hkl"].map(obs_by_hkl["matched"].to_dict()).fillna(False) else: calc_df["matched"] = False return calc_df, assign_df def plot_compare_results( x: np.ndarray | None, y: np.ndarray | None, peaks: list[Peak], theoretical_df: pd.DataFrame, assignment_df: pd.DataFrame, candidate: Candidate, xmin: float, xmax: float, outpath: Path | None, show: bool = False, save: bool = True, ) -> None: """ 概要: 観測ピークと理論的な回折線の比較結果をプロットします。 詳細説明: 入力XRDデータ (x, y) があればそれをプロットし、観測されたピークを縦線で示します。 指定された格子定数候補に基づく理論的な回折線を、観測ピークの下部に棒グラフとして重ねて表示します。 理論線にはhkl指数が注釈され、マッチングされたものは強調されます。 グラフは画像ファイルとして保存でき、Matplotlibウィンドウに表示することも可能です。 mplcursorsがインストールされていれば、線にホバー注釈が表示されます。 引数: :param x: 2θ角度のnumpy配列（生のXRDデータ、オプション）。 :type x: numpy.ndarray | None :param y: 強度のnumpy配列（生のXRDデータ、オプション）。 :type y: numpy.ndarray | None :param peaks: 観測されたPeakオブジェクトのリスト。 :type peaks: list[Peak] :param theoretical_df: 理論的な回折線情報を含むPandasデータフレーム。 :type theoretical_df: pandas.DataFrame :param assignment_df: 観測ピークと理論線の割り当て結果データフレーム。 :type assignment_df: pandas.DataFrame :param candidate: 格子定数候補。 :type candidate: Candidate :param xmin: プロットのX軸最小値（度）。 :type xmin: float :param xmax: プロットのX軸最大値（度）。 :type xmax: float :param outpath: プロットを保存するファイルパス（Noneの場合は保存しない）。 :type outpath: Path | None :param show: プロットウィンドウを表示するかどうか。 :type show: bool :param save: プロットをファイルに保存するかどうか。 :type save: bool """ if not show and not save: return fig, ax = plt.subplots(figsize=(12, 7)) if x is not None and y is not None and len(x) and len(y): ax.plot(x, y, lw=0.8, label="入力XRD") ymax = float(np.nanmax(y)) if np.isfinite(y).any() else 1.0 else: ymax = max([p.intensity for p in peaks] + [1.0]) markerline, stemlines, baseline = ax.stem( [p.two_theta for p in peaks], [p.intensity for p in peaks], basefmt=" ", linefmt="C0-", markerfmt="C0o", label="観測ピーク", ) try: plt.setp(stemlines, linewidth=1.0) plt.setp(markerline, markersize=4) except Exception: pass hover_artists = [] hover_texts = [] assign_by_obs = {} if not assignment_df.empty: for _, row in assignment_df.iterrows(): assign_by_obs[round(float(row["two_theta_obs_deg"]), 6)] = row for p in sorted(peaks, key=lambda pk: pk.two_theta): key = round(float(p.two_theta), 6) row = assign_by_obs.get(key) hkl_text = None calc_tt = None resid = None if row is not None: hkl_text = f"{int(row['h'])}{int(row['k'])}{int(row['l'])}" calc_tt = float(row["two_theta_calc_deg"]) resid = float(row["resid_two_theta_deg"]) line = ax.axvline(p.two_theta, color="tab:green", lw=0.3, linestyle="--", alpha=0.75, picker=5) hover_artists.append(line) hover_texts.append(_peak_hover_text(p, hkl_text=hkl_text, calc_two_theta=calc_tt, resid=resid)) theory_height = ymax * 0.16 if ymax > 0 else 1.0 for _, row in theoretical_df.iterrows(): tt = float(row["two_theta_calc_deg"]) matched = bool(row.get("matched", False)) line = ax.vlines(tt, ymin=0, ymax=theory_height, colors="tab:red", lw=1.0 if matched else 0.7, alpha=0.85 if matched else 0.45) # mplcursorsはvlinesからのLineCollectionを受け入れます。 hkl_text = f"{int(row['h'])}{int(row['k'])}{int(row['l'])}" hover = [ f"hkl: {hkl_text}", f"2theta(calc): {tt:.5f} deg", ] if pd.notna(row.get("nearest_two_theta_obs_deg", np.nan)): hover.append(f"2theta(obs): {float(row['nearest_two_theta_obs_deg']):.5f} deg") if pd.notna(row.get("nearest_resid_two_theta_deg", np.nan)): hover.append(f"resid: {float(row['nearest_resid_two_theta_deg']):+.5f} deg") if pd.notna(row.get("nearest_intensity", np.nan)): hover.append(f"intensity: {float(row['nearest_intensity']):.6g}") if pd.notna(row.get("nearest_fwhm_deg", np.nan)): hover.append(f"FWHM: {float(row['nearest_fwhm_deg']):.5g} deg") hover_artists.append(line) hover_texts.append("\n".join(hover)) if matched: ax.text(tt, theory_height, hkl_text, rotation=90, va="bottom", fontsize=7) title_params = ", ".join(f"{k}={v:.5g}" for k, v in candidate.params.items() if k in {"a", "b", "c"}) ax.set_title(f"{candidate.crystal_system}: {title_params}") ax.set_xlabel("2Theta (度)") ax.set_ylabel("強度") ax.set_xlim([xmin, xmax]) ax.legend(loc="best") fig.tight_layout() if save and outpath is not None: fig.savefig(outpath, dpi=160) if show: _maybe_install_mplcursors(hover_artists, hover_texts) plt.pause(0.1) input("\n続行するにはEnterを押してください>>\n") if save and outpath is not None: plt.close(fig) def frange_values(vmin: float, vmax: float, step: float) -> np.ndarray: """ 概要: 指定された範囲とステップで浮動小数点数の配列を生成します。 詳細説明: numpy.arange のように動作しますが、浮動小数点数の丸め誤差を考慮して vmin から vmax までの値を step 間隔で生成します。 引数: :param vmin: 最小値。 :type vmin: float :param vmax: 最大値。 :type vmax: float :param step: ステップサイズ。 :type step: float 戻り値: :returns: 生成された浮動小数点数のnumpy配列。 :rtype: numpy.ndarray 例外: :raises ValueError: ステップが正でない、または最大値が最小値より小さい場合。 """ if step <= 0: raise ValueError("グリッドステップは正でなければなりません。") if vmax < vmin: raise ValueError("グリッドの最大値は最小値以上でなければなりません。") n = int(math.floor((vmax - vmin) / step + 0.5)) + 1 vals = vmin + step * np.arange(max(n, 1)) return vals[vals <= vmax + 1e-12] def _param_grid_for_system(system: str, ranges: dict[str, tuple[float, float, float]]) -> Iterable[dict[str, float]]: """ 概要: 特定の結晶系に対して格子定数のグリッドを生成します。 詳細説明: 与えられた結晶系 (system) と格子定数の範囲 (ranges) に基づいて、 可能な格子定数の組み合わせを辞書のイテラブルとして生成します。 各結晶系に応じて必要な格子定数 (a, b, c) と角度が設定されます。 引数: :param system: 結晶系名。 :type system: str :param ranges: 各格子定数（a, b, c）の (最小値, 最大値, ステップ) を含む辞書。 :type ranges: dict[str, tuple[float, float, float]] 戻り値: :returns: 格子定数辞書のイテラブル。 :rtype: Iterable[dict[str, float]] 例外: :raises ValueError: サポートされていない結晶系が指定された場合。 """ avec = frange_values(*ranges["a"]) if system == "cubic": for a in avec: yield {"a": float(a), "b": float(a), "c": float(a), "alpha": 90.0, "beta": 90.0, "gamma": 90.0} return cvec = frange_values(*ranges["c"]) if system == "tetragonal": for a in avec: for c in cvec: yield {"a": float(a), "b": float(a), "c": float(c), "alpha": 90.0, "beta": 90.0, "gamma": 90.0} return if system == "hexagonal": for a in avec: for c in cvec: yield {"a": float(a), "b": float(a), "c": float(c), "alpha": 90.0, "beta": 90.0, "gamma": 120.0} return if system == "orthorhombic": bvec = frange_values(*ranges["b"]) for a in avec: for b in bvec: for c in cvec: yield {"a": float(a), "b": float(b), "c": float(c), "alpha": 90.0, "beta": 90.0, "gamma": 90.0} return raise ValueError(system) def _ranges_from_args_for_system(args: argparse.Namespace, system: str) -> dict[str, tuple[float, float, float]] | None: """ 概要: コマンドライン引数から特定の結晶系に対する格子定数範囲を構築します。 詳細説明: argsオブジェクトに含まれる amin, amax, astep, bmin, bmax, bstep, cmin, cmax, cstep の値を解析し、指定された結晶系に必要な格子定数範囲の辞書を返します。 必要な範囲が完全に指定されていない場合はNoneを返します。 引数: :param args: コマンドライン引数を格納したNamespaceオブジェクト。 :type args: argparse.Namespace :param system: 結晶系名。 :type system: str 戻り値: :returns: 格子定数範囲の辞書、または不足している場合はNone。 :rtype: dict[str, tuple[float, float, float]] | None """ if args.amin is None or args.amax is None: return None ranges = {"a": (float(args.amin), float(args.amax), float(args.astep))} if system in {"tetragonal", "hexagonal", "orthorhombic"}: if args.cmin is None or args.cmax is None: return None ranges["c"] = (float(args.cmin), float(args.cmax), float(args.cstep)) if system == "orthorhombic": if args.bmin is None or args.bmax is None: return None ranges["b"] = (float(args.bmin), float(args.bmax), float(args.bstep)) return ranges def _ranges_around_candidate(candidate: Candidate, margin: float, args: argparse.Namespace) -> dict[str, tuple[float, float, float]]: """ 概要: 既存の格子定数候補の周囲に検索範囲を生成します。 詳細説明: 与えられたCandidateオブジェクトの格子定数に基づいて、 指定された相対マージン (margin) とステップサイズ (astep, bstep, cstep) を使用して、 新しい格子定数検索範囲の辞書を構築します。これは、グリッド検索のシードとして使用されます。 引数: :param candidate: 基準となるCandidateオブジェクト。 :type candidate: Candidate :param margin: 格子定数に適用する相対マージン。 :type margin: float :param args: コマンドライン引数を格納したNamespaceオブジェクト（ステップサイズ用）。 :type args: argparse.Namespace 戻り値: :returns: 格子定数検索範囲の辞書。 :rtype: dict[str, tuple[float, float, float]] """ params = candidate.params ranges: dict[str, tuple[float, float, float]] = {} for key, step_name in [("a", "astep"), ("b", "bstep"), ("c", "cstep")]: if key not in params or not math.isfinite(float(params[key])): continue value = float(params[key]) step = float(getattr(args, step_name)) ranges[key] = (max(0.1, value * (1.0 - margin)), value * (1.0 + margin), step) if candidate.crystal_system in {"cubic"}: ranges = {"a": ranges["a"]} elif candidate.crystal_system in {"tetragonal", "hexagonal"}: ranges = {"a": ranges["a"], "c": ranges["c"]} elif candidate.crystal_system == "orthorhombic": ranges = {"a": ranges["a"], "b": ranges["b"], "c": ranges["c"]} return ranges def _assign_peaks_for_params( peaks: list[Peak], system: str, params: dict[str, float], hmax: int, wavelength: float, tolerance_deg: float, ) -> tuple[list[dict], float, float, int]: """ 概要: 与えられた格子定数パラメータに対して観測ピークを割り当て、評価します。 詳細説明: 特定の結晶系 (system) と格子定数 (params) に基づいて理論的な回折線を生成し、 観測ピーク (peaks) をそれらの理論線に割り当てます。 割り当ては2θ角度の許容誤差 (tolerance_deg) 内で行われます。 割り当て結果、RMS残差 (2θと相対inv_d2)、およびマッチングされたピーク数を返します。 引数: :param peaks: 観測されたPeakオブジェクトのリスト。 :type peaks: list[Peak] :param system: 結晶系名。 :type system: str :param params: 格子定数を含む辞書。 :type params: dict[str, float] :param hmax: h, k, l指数の最大値。 :type hmax: int :param wavelength: X線波長（オングストローム）。 :type wavelength: float :param tolerance_deg: 2θ角度の割り当て許容誤差（度）。 :type tolerance_deg: float 戻り値: :returns: (割り当て結果のリスト, 2θのRMS残差, 相対inv_d2のRMS残差, マッチングされたピーク数) のタプル。 :rtype: tuple[list[dict], float, float, int] """ cand = Candidate(system, {"orthorhombic": 4, "tetragonal": 5, "cubic": 6, "hexagonal": 8}[system], params, 0, 0.0, [], []) calc_df = theoretical_lines_for_candidate(cand, hmax=hmax, wavelength=wavelength) if calc_df.empty: return [], float("inf"), float("inf"), 0 calc_tt = calc_df["two_theta_calc_deg"].to_numpy(dtype=float) used: set[int] = set() rows = [] residuals = [] rels = [] for pk in sorted(peaks, key=lambda p: p.two_theta): order = np.argsort(np.abs(calc_tt - pk.two_theta)) chosen = None for idx in order: if int(idx) not in used: chosen = int(idx) break if chosen is None: continue used.add(chosen) calc = calc_df.iloc[chosen] resid = pk.two_theta - float(calc["two_theta_calc_deg"]) matched = abs(resid) <= tolerance_deg if matched: residuals.append(float(resid)) rel = abs(float(calc["inv_d2_calc"]) - pk.inv_d2) / max(pk.inv_d2, 1e-12) rels.append(float(rel)) rows.append({ "two_theta_obs_deg": pk.two_theta, "intensity": pk.intensity, "ka_role": pk.ka_role, "fwhm_deg": pk.fwhm_deg, "inv_d2_obs": pk.inv_d2, "h": int(calc["h"]) if matched else None, "k": int(calc["k"]) if matched else None, "l": int(calc["l"]) if matched else None, "inv_d2_calc": float(calc["inv_d2_calc"]) if matched else None, "rel_error_inv_d2": rels[-1] if matched else None, "two_theta_calc_deg": float(calc["two_theta_calc_deg"]) if matched else None, "rel_error_2theta": abs(resid) / max(pk.two_theta, 1e-12) if matched else None, "resid_two_theta_deg": float(resid) if matched else None, "matched": bool(matched), }) nmatched = len(residuals) rms_tt = float(np.sqrt(np.mean(np.square(residuals)))) if residuals else float("inf") rms_rel = float(np.sqrt(np.mean(np.square(rels)))) if rels else float("inf") return rows, rms_tt, rms_rel, nmatched def _refit_params_from_assignment(system: str, rows: list[dict]) -> dict[str, float] | None: """ 概要: 割り当てられたピークから格子定数を再フィットします。 詳細説明: 割り当て結果 (rows) の中からマッチしたピークの inv_d2 値とhkl指数を用いて、 最小二乗法により格子定数を再計算します。 再計算された格子定数を含む辞書を返します。 フィットが不可能な場合（例: マッチングされたピークが少ない、係数が負になる）はNoneを返します。 引数: :param system: 結晶系名。 :type system: str :param rows: 割り当て結果のリスト。 :type rows: list[dict] 戻り値: :returns: 再フィットされた格子定数を含む辞書、またはNone。 :rtype: dict[str, float] | None """ matched = [r for r in rows if r.get("matched")] if not matched: return None y = np.array([float(r["inv_d2_obs"]) for r in matched], dtype=float) feats = [] for r in matched: h, k, l = int(r["h"]), int(r["k"]), int(r["l"]) if system == "cubic": feats.append((h*h + k*k + l*l,)) elif system == "tetragonal": feats.append((h*h + k*k, l*l)) elif system == "hexagonal": feats.append((h*h + h*k + k*k, l*l)) elif system == "orthorhombic": feats.append((h*h, k*k, l*l)) else: raise ValueError(system) X = np.array(feats, dtype=float) if len(y) < X.shape[1]: return None coef, _, rank, _ = np.linalg.lstsq(X, y, rcond=None) if rank < X.shape[1] or np.any(coef <= 0): return None return params_from_beta(system, coef) def score_system_grid( system: str, selected: list[Peak], all_peaks: list[Peak], ranges: dict[str, tuple[float, float, float]], hmax: int = 12, wavelength: float = CU_KA1, tolerance_deg: float = 0.25, refine: int = 2, keep: int = 50, penalty_unassigned: float = 0.5, lattice_max: float | None = None, ) -> list[Candidate]: """ 概要: グリッド検索法を用いて特定の結晶系の格子定数候補を評価します。 詳細説明: 指定された結晶系 (system) の格子定数範囲 (ranges) 内でグリッド検索を行い、 各格子定数セットに対して観測ピーク (selected) を割り当てて評価します。 割り当てられたピークから格子定数を繰り返し再フィット (refine回) し、 最もフィットの良い候補を絞り込みます。 未割り当てピークに対するペナルティ (penalty_unassigned) を適用し、 結果はスコアに基づいてソートされ、指定された数 (keep) の上位候補が返されます。 格子定数上限 (lattice_max) も考慮されます。 引数: :param system: 結晶系名。 :type system: str :param selected: 格子定数推定に使用する、フィルタリングされたPeakオブジェクトのリスト。 :type selected: list[Peak] :param all_peaks: 全ての観測されたPeakオブジェクトのリスト。 :type all_peaks: list[Peak] :param ranges: 各格子定数（a, b, c）の (最小値, 最大値, ステップ) を含む辞書。 :type ranges: dict[str, tuple[float, float, float]] :param hmax: h, k, l指数の最大値。 :type hmax: int :param wavelength: X線波長（オングストローム）。デフォルトはCU_KA1です。 :type wavelength: float :param tolerance_deg: 2θ角度の割り当て許容誤差（度）。 :type tolerance_deg: float :param refine: グリッド検索での割り当てと再フィットの反復回数。 :type refine: int :param keep: 内部的に保持する候補の最大数。 :type keep: int :param penalty_unassigned: 未割り当てピークに対するペナルティスコア。 :type penalty_unassigned: float :param lattice_max: 格子定数の上限値（オングストローム）、またはNone。 :type lattice_max: float | None 戻り値: :returns: 評価され、ソートされたCandidateオブジェクトのリスト。 :rtype: list[Candidate] """ candidates: list[Candidate] = [] ls_map = {"orthorhombic": 4, "tetragonal": 5, "cubic": 6, "hexagonal": 8} for params0 in _param_grid_for_system(system, ranges): if not params_within_lattice_max(params0, system, lattice_max): continue params = dict(params0) rows = [] for _ in range(max(refine, 1)): rows, _rms_tt, _rms_rel, nmatched = _assign_peaks_for_params( selected, system, params, hmax, wavelength, tolerance_deg ) if nmatched < len(beta_from_params(system, params)): rows = [] break refit = _refit_params_from_assignment(system, rows) if refit is None: rows = [] break params = refit if not rows: continue if not params_within_lattice_max(params, system, lattice_max): continue selected_rows, rms_tt, rms_rel, nmatched = _assign_peaks_for_params( selected, system, params, hmax, wavelength, tolerance_deg ) if nmatched == 0: continue all_rows, _all_rms_tt, all_rms_rel, all_nmatched = _assign_peaks_for_params( all_peaks, system, params, hmax, wavelength, tolerance_deg ) # ソート専用の private キーにペナルティ調整値を格納します。 penalty_score = rms_tt + penalty_unassigned * (len(selected) - nmatched) cand = Candidate( crystal_system=system, ls_code=ls_map[system], params=params, score_matches=int(all_nmatched), score_rms_rel=float(all_rms_rel if math.isfinite(all_rms_rel) else rms_rel), selected_matches=[r for r in selected_rows if r.get("matched")], all_matches=all_rows, ) cand.params = dict(cand.params) cand.params["_grid_penalty_score"] = float(penalty_score) candidates.append(cand) candidates.sort(key=lambda c: (c.params.get("_grid_penalty_score", float("inf")), -c.score_matches, c.score_rms_rel)) if len(candidates) > keep: candidates = candidates[:keep] for c in candidates: c.params.pop("_grid_penalty_score", None) candidates.sort(key=lambda c: (-c.score_matches, c.score_rms_rel)) return candidates def deduplicate_candidates(candidates: list[Candidate], ndigits: int = 4) -> list[Candidate]: """ 概要: 重複する格子定数候補をリストから除去します。 詳細説明: Candidateオブジェクトの結晶系と、a, b, cの格子定数（ndigits桁で丸められたもの） に基づいて重複を判断します。 リストの順序を維持しつつ、最初に出現したユニークな候補のみを返します。 引数: :param candidates: Candidateオブジェクトのリスト。 :type candidates: list[Candidate] :param ndigits: 格子定数を比較する際の丸め桁数。 :type ndigits: int 戻り値: :returns: 重複が除去されたCandidateオブジェクトのリスト。 :rtype: list[Candidate] """ seen = set() out = [] for c in candidates: key = (c.crystal_system,) + tuple(round(float(c.params.get(k, 0.0)), ndigits) for k in ["a", "b", "c"]) if key in seen: continue seen.add(key) out.append(c) return out def guess_candidates( args: argparse.Namespace, peaks: list[Peak], crystal_systems: list[str], ) -> tuple[list[Candidate], list[Peak]]: """ 概要: 与えられたピークと結晶系に基づいて格子定数候補を推定します。 詳細説明: まず、最も強度の高い非Kα2ピークを args.npeak の数だけ選択します。 その後、args.method に応じて「組み合わせ探索」または「グリッド検索」を実行します。 グリッド検索は、明示的な範囲が指定されていない場合、組み合わせ探索で得られた 上位候補をシードとして使用できます。 格子定数上限 (args.resolved_lattice_max) も適用されます。 引数: :param args: コマンドライン引数を格納したNamespaceオブジェクト。 :type args: argparse.Namespace :param peaks: 観測されたPeakオブジェクトのリスト。 :type peaks: list[Peak] :param crystal_systems: 探索する結晶系のリスト。 :type crystal_systems: list[str] 戻り値: :returns: (評価されソートされたCandidateオブジェクトのリスト, 推定に使用されたピークのリスト) のタプル。 :rtype: tuple[list[Candidate], list[Peak]] 例外: :raises ValueError: サポートされていない推定メソッドが指定された場合。 """ selected = [p for p in sorted(peaks, key=lambda p: p.intensity, reverse=True) if p.ka_role != "ka2"][:args.npeak] method = args.method.lower() if method not in SUPPORTED_METHODS: raise ValueError(f"サポートされていないメソッド: {args.method}") candidates: list[Candidate] = [] if method in {"combinatorial", "both", "grid"}: # --method grid で明示的な範囲がない場合、これらの候補がシードになります。 seed_candidates: list[Candidate] = [] for system in crystal_systems: cand = score_system(system, sorted(selected, key=lambda p: p.two_theta), sorted(peaks, key=lambda p: p.two_theta), hmax=args.hmax) if cand is not None: seed_candidates.append(cand) seed_candidates.sort(key=lambda c: (-c.score_matches, c.score_rms_rel)) seed_candidates = filter_candidates_by_lattice_max( seed_candidates, getattr(args, "resolved_lattice_max", None), label="組み合わせ探索シード候補" ) else: seed_candidates = [] if method == "combinatorial": candidates = seed_candidates elif method in {"grid", "both"}: if method == "both": candidates.extend(seed_candidates) for system in crystal_systems: ranges = _ranges_from_args_for_system(args, system) if ranges is not None: grid_cands = score_system_grid( system, selected, peaks, ranges, hmax=args.hmax, wavelength=args.wavelength, tolerance_deg=args.tolerance_deg, refine=args.grid_refine, keep=args.keep, penalty_unassigned=args.penalty_unassigned, lattice_max=getattr(args, "resolved_lattice_max", None), ) candidates.extend(grid_cands) continue # 明示的な範囲がなく、グリッド検索では同じシステムの上位組み合わせ候補を使用します。 seeds = [c for c in seed_candidates if c.crystal_system == system][:args.grid_seed_top] if not seeds and method == "grid": print(f"'{system}' の明示的なグリッド範囲もシード候補も見つかりませんでした。スキップしました。") for seed in seeds: ranges = _ranges_around_candidate(seed, args.grid_margin, args) grid_cands = score_system_grid( system, selected, peaks, ranges, hmax=args.hmax, wavelength=args.wavelength, tolerance_deg=args.tolerance_deg, refine=args.grid_refine, keep=args.keep, penalty_unassigned=args.penalty_unassigned, lattice_max=getattr(args, "resolved_lattice_max", None), ) candidates.extend(grid_cands) candidates = filter_candidates_by_lattice_max( candidates, getattr(args, "resolved_lattice_max", None), label="最終候補", is_print=False ) candidates.sort(key=lambda c: (-c.score_matches, c.score_rms_rel)) return deduplicate_candidates(candidates), selected def read_candidate_from_guess_file(path: Path, rank: int) -> Candidate: """ 概要: 推定結果ファイルから特定の順位の格子定数候補を読み込みます。 詳細説明: Excel形式の推定結果ファイルから "candidate_summary" シートを読み込み、 指定されたランク (rank) に対応する格子定数候補をCandidateオブジェクトとして返します。 ランクは1から始まります。 引数: :param path: 推定結果Excelファイルへのパス。 :type path: Path :param rank: 読み込む候補の順位（1から始まる）。 :type rank: int 戻り値: :returns: 指定された順位のCandidateオブジェクト。 :rtype: Candidate 例外: :raises ValueError: Excelファイルでない場合、"candidate_summary" シートがない場合、 または指定されたランクの候補が見つからない場合。 """ suffix = path.suffix.lower() if suffix not in {".xlsx", ".xls"}: raise ValueError("候補ランク比較にはExcel形式の推定ファイルか、ピークからの再計算が必要です。") xls = pd.ExcelFile(path) lower_map = {str(s).lower(): s for s in xls.sheet_names} if "candidate_summary" not in lower_map: raise ValueError(f"'{path}' に candidate_summary シートが見つかりませんでした。") df = pd.read_excel(path, sheet_name=lower_map["candidate_summary"]) if df.empty: raise ValueError(f"'{path}' の candidate_summary が空です。") if "rank" in df.columns: row_df = df[df["rank"] == rank] if row_df.empty: raise ValueError(f"'{path}' に候補ランク {rank} が見つかりませんでした。") row = row_df.iloc[0] else: if rank < 1 or rank > len(df): raise ValueError(f"候補ランク {rank} は 1..{len(df)} の範囲外です。") row = df.iloc[rank - 1] return candidate_from_summary_row(row) def try_read_raw_xy_for_plot(path: Path) -> tuple[np.ndarray | None, np.ndarray | None]: """ 概要: プロットのために生のX-Yデータを読み込みを試みます。 詳細説明: 指定されたファイルがピーク情報や推定結果のExcelワークブックでないことを確認し、 そうであれば read_xy_data を使って生のX-Yデータを読み込みます。 読み込みに失敗した場合や、ファイルが特定のワークブックである場合はNoneを返します。 引数: :param path: データファイルへのパス。 :type path: Path 戻り値: :returns: (Xデータnumpy配列, Yデータnumpy配列) のタプル、 または読み込みに失敗した場合は (None, None)。 :rtype: tuple[numpy.ndarray | None, numpy.ndarray | None] """ try: # ピーク/推定ワークブックを生のXRDデータとして扱わないようにします。 if path.suffix.lower() in {".xlsx", ".xls"}: xls = pd.ExcelFile(path) names = {str(s).lower() for s in xls.sheet_names} if names & {"peaks", "all_detected_peaks", "candidate_summary", "theoretical_lines"}: return None, None return read_xy_data(path) except Exception: return None, None def looks_like_raw_xrd_file(path: Path, min_rows: int = 50) -> bool: """ 概要: ファイルが生のXRDデータファイルのように見えるかどうかをヒューリスティックに判断します。 詳細説明: ファイルがピークリストではなく、高密度な生のXRD X-Yデータファイルであるかどうかを ヒューリスティックに判定します。 Excelファイルの場合、既知のピーク/推定シート名が含まれていないかを確認します。 ファイルの行数、列数、最初の数列に十分な数値データが含まれているかを確認します。 "two_theta_deg" のようなピークリスト特有の列名がないことも確認します。 引数: :param path: ファイルパス。 :type path: Path :param min_rows: 生のXRDファイルと見なすための最小行数。 :type min_rows: int 戻り値: :returns: ファイルが生のXRDデータファイルのように見える場合はTrue、そうでなければFalse。 :rtype: bool """ try: suffix = path.suffix.lower() if suffix in {".xlsx", ".xls"}: xls = pd.ExcelFile(path) names = {str(s).lower() for s in xls.sheet_names} if names & {"peaks", "all_detected_peaks", "candidate_summary", "theoretical_lines"}: return False df = pd.read_excel(path) else: df = pd.read_csv(path, sep=None, engine="python", comment="#") if df.empty or df.shape[1] < 2: return False peak_cols = { "two_theta_deg", "two_theta", "twotheta", "2theta", "2theta_deg", "theta2", "angle", "angle_deg", "position", "position_deg", "two_theta_obs_deg", "inv_d2", "inv_d2_obs", "inv_d2_a_2", } if any(_canonical_column_name(c) in peak_cols for c in df.columns): return False numeric_cols = 0 for col in df.columns[:4]: if pd.to_numeric(df[col], errors="coerce").notna().sum() >= min_rows: numeric_cols += 1 return len(df) >= min_rows and numeric_cols >= 2 except Exception: return False def build_argument_parser() -> argparse.ArgumentParser: """ 概要: コマンドライン引数を解析するためのArgumentParserオブジェクトを構築します。 詳細説明: ピーク探索、格子定数推定、比較などの機能に必要な全てのコマンドライン引数を定義します。 各引数には説明、型、デフォルト値、選択肢などが設定されています。 引数は、入力ファイル、モード、X線波長、ピーク探索オプション、インデクシング/比較オプション、 グリッド検索範囲、プロットオプションに分類されます。 戻り値: :returns: 構築されたArgumentParserオブジェクト。 :rtype: argparse.ArgumentParser """ p = argparse.ArgumentParser(description="粉末X線回折のためのピーク探索と格子定数推定") p.add_argument( "arg1", help=( "入力データファイル、ピーク情報ファイル、または推定ワークブック。 " "レガシーな使用法 'search infile' も受け入れられます。" ), ) p.add_argument( "arg2", nargs="?", help="レガシーな位置指定モードのための入力ファイル、例: 'search sample.TXT'", ) p.add_argument( "--mode", type=str, default="all", choices=["search", "guess", "all", "compare", "calc"], help=( "search: ピーク探索のみ; guess: ピークファイルから格子定数を推定; " "all: 探索 + 推定 (デフォルト); compare: ランク付けされた候補と観測ピークを比較; calc: ユーザー指定の格子定数を比較" ), ) p.add_argument( "--method", type=str, default="combinatorial", choices=SUPPORTED_METHODS, help=( "格子定数推定方法。combinatorial: 高速hkl組み合わせ探索; " "grid: 明示的な範囲または組み合わせシードを使用するグリッド検索; both: 両方の結果をマージ。" ), ) p.add_argument( "--crystal-system", "--system", type=str, default="auto", help=( "テストする結晶系。auto/all, cubic, tetragonal, hexagonal, " "orthorhombic、またはカンマ区切りの値を指定。デフォルト: auto" ), ) p.add_argument("--peak-file", type=Path, default=None, help="ピーク情報Excel/CSVファイルパス") p.add_argument("--guess-file", type=Path, default=None, help="--mode compare で使用される推定ワークブック") p.add_argument("--output-file", type=Path, default=None, help="guess/all/compare 結果の出力ワークブック") p.add_argument("--plot-file", type=Path, default=None, help="出力グラフパス") p.add_argument("--candidate-rank", type=int, default=1, help="--mode compare で使用される候補ランク") p.add_argument("--wavelength", type=float, default=CU_KA1, help=f"X線波長（オングストローム）。デフォルト: Cu Kα1 = {CU_KA1}") # --mode calc または推定ワークブックなしの --mode compare 用の明示的な格子定数。 p.add_argument("--a", "--lattice-a", dest="a", type=float, default=None, help="明示的な格子定数a（オングストローム）") p.add_argument("--b", "--lattice-b", dest="b", type=float, default=None, help="明示的な格子定数b（オングストローム）") p.add_argument("--c", "--lattice-c", dest="c", type=float, default=None, help="明示的な格子定数c（オングストローム）") p.add_argument("--alpha", type=float, default=None, help="明示的なα角（度）; 現在は出力に保存されます") p.add_argument("--beta", type=float, default=None, help="明示的なβ角（度）; 現在は出力に保存されます") p.add_argument("--gamma", type=float, default=None, help="明示的なγ角（度）; 現在は出力に保存されます") # 生XRD / ピーク探索オプション。 p.add_argument("--xmin", type=float, default=None) p.add_argument("--xmax", type=float, default=None) p.add_argument("--threshold", type=float, default=1000.0) p.add_argument("--nsmooth", type=int, default=11) p.add_argument("--norder", type=int, default=4) p.add_argument("--ydiff1-threshold", type=float, default=1.0e-2) p.add_argument("--fwhm-min-deg", type=float, default=0.06) p.add_argument("--fwhm-max-deg", type=float, default=1.0) # インデクシング / 比較オプション。 p.add_argument("--npeak", type=int, default=10) p.add_argument("--hmax", type=int, default=12) p.add_argument( "--lattice-max", "--max-lattice", dest="lattice_max", type=float, default=None, help=( "推定される格子定数の上限値（オングストローム）。 " "デフォルト: d(観測された非Kα2最小2θ) x --lattice-max-factor。 " "制限を無効にするには0以下に設定。" ), ) p.add_argument( "--lattice-max-factor", type=float, default=1.2, help="--lattice-max が省略された場合、d(観測された最小2θ) に乗じられる係数。デフォルト: 1.2", ) p.add_argument("--tolerance-deg", type=float, default=0.25, help="グリッド/比較での2θ割り当て許容誤差") p.add_argument("--penalty-unassigned", type=float, default=0.5) p.add_argument("--keep", type=int, default=200, help="内部的に保持されるグリッド候補の数") p.add_argument("--grid-refine", type=int, default=2, help="グリッド検索での割り当て/再フィット反復回数") p.add_argument("--grid-margin", type=float, default=0.05, help="グリッド検索のための組み合わせシード周辺の相対マージン") p.add_argument("--grid-seed-top", type=int, default=3, help="ローカルグリッド検索に使用されるシステムごとの組み合わせシードの数") # 明示的なグリッド範囲。省略した場合、グリッド検索は組み合わせシードを使用できます。 p.add_argument("--amin", type=float, default=None) p.add_argument("--amax", type=float, default=None) p.add_argument("--astep", type=float, default=0.02) p.add_argument("--bmin", type=float, default=None) p.add_argument("--bmax", type=float, default=None) p.add_argument("--bstep", type=float, default=0.02) p.add_argument("--cmin", type=float, default=None) p.add_argument("--cmax", type=float, default=None) p.add_argument("--cstep", type=float, default=0.02) p.add_argument("--lsq-script", type=Path, default=Path("lsq_latt2.py")) p.add_argument("--show", type=int, default=0, choices=[0, 1], help="Matplotlibグラフウィンドウを表示") p.add_argument("--save", "--save-plot", dest="save_plot", type=int, default=1, choices=[0, 1], help="グラフ画像を保存") p.add_argument( "--plot-theory", type=str, default="near", choices=["near", "matched", "all", "none"], help=( "比較/計算プロットで描画する計算された理論線の種類。 " "near は --plot-near-deg 内で観測ピークに割り当てられた線のみを描画; " "all は大きな単位セルでは非常に高密度になる可能性があります。デフォルト: near" ), ) p.add_argument("--plot-near-deg", type=float, default=0.5, help="--plot-theory near で使用される2θウィンドウ") p.add_argument("--no-show", action="store_true", help="レガシーオプション; --show 0 と同等") return p def resolve_mode_and_input(args: argparse.Namespace) -> tuple[str, Path]: """ 概要: コマンドライン引数から実行モードと入力ファイルを解決します。 詳細説明: 位置指定引数 arg1 と arg2、および --mode オプションを解析し、 スクリプトの実行モードと主要な入力ファイルパスを決定します。 レガシーな位置指定モード ("search sample.TXT") にも対応しています。 引数: :param args: コマンドライン引数を格納したNamespaceオブジェクト。 :type args: argparse.Namespace 戻り値: :returns: (解決されたモード名, 入力ファイルパス) のタプル。 :rtype: tuple[str, Path] 例外: :raises SystemExit: 予期しない位置指定引数が指定された場合。 """ first = str(args.arg1) if args.arg2 is not None and first.lower() in KNOWN_MODES: mode = normalize_mode(first) infile = Path(args.arg2) else: if args.arg2 is not None: raise SystemExit( "予期しない2番目の位置指定引数です。次のいずれかを使用してください:\n" " python guess_lattice_parameters.py sample.TXT --mode all\n" "またはレガシースタイル:\n" " python guess_lattice_parameters.py search sample.TXT" ) mode = normalize_mode(args.mode) infile = Path(args.arg1) return mode, infile def apply_x_range(x: np.ndarray, y: np.ndarray, xmin: float | None, xmax: float | None) -> tuple[np.ndarray, np.ndarray]: """ 概要: X-YデータにX軸範囲のフィルタリングを適用します。 詳細説明: X座標が xmin 以上 xmax 以下のデータポイントのみを保持するように numpy配列 x と y をフィルタリングします。 xmin または xmax がNoneの場合、その側のフィルタリングは適用されません。 引数: :param x: X座標のnumpy配列。 :type x: numpy.ndarray :param y: Y座標のnumpy配列。 :type y: numpy.ndarray :param xmin: 最小X値の閾値（オプション）。 :type xmin: float | None :param xmax: 最大X値の閾値（オプション）。 :type xmax: float | None 戻り値: :returns: フィルタリングされたX値とY値のnumpy配列のタプル。 :rtype: tuple[numpy.ndarray, numpy.ndarray] 例外: :raises ValueError: xmin/xmax適用後にデータポイントが残らない場合。 """ if xmin is not None: mask = x >= xmin x, y = x[mask], y[mask] if xmax is not None: mask = x <= xmax x, y = x[mask], y[mask] if len(x) == 0: raise ValueError("xmin/xmax 適用後にデータポイントが残りません。") return x, y def run_search(args: argparse.Namespace, infile: Path) -> tuple[list[Peak], Path, Path | None]: """ 概要: ピーク探索モードを実行します。 詳細説明: 入力XRDデータファイル (infile) を読み込み、ピーク探索アルゴリズム (peak_search_deriv3) を実行します。検出されたピークにKα2の役割を割り当て、結果をコンソールに表示します。 また、ピーク情報とプロットをファイルに保存します。 引数: :param args: コマンドライン引数を格納したNamespaceオブジェクト。 :type args: argparse.Namespace :param infile: 入力データファイルへのパス。 :type infile: Path 戻り値: :returns: (検出されたPeakオブジェクトのリスト, ピーク情報ファイルパス, プロットファイルパス) のタプル。 :rtype: tuple[list[Peak], Path, Path | None] """ x, y = read_xy_data(infile) x, y = apply_x_range(x, y, args.xmin, args.xmax) print(f"入力ファイル : {infile}") print("モード : search") print(f"2theta範囲 : {x.min():.3f} - {x.max():.3f}") print(f"閾値 : {args.threshold}") print(f"nsmooth : {args.nsmooth}") print(f"norder : {args.norder}") print(f"FWHM範囲 : {args.fwhm_min_deg} - {args.fwhm_max_deg} deg") peaks, info = peak_search_deriv3( x, y, nsmooth=args.nsmooth, norder=args.norder, threshold=args.threshold, ydiff1_threshold=args.ydiff1_threshold, fwhm_min_deg=args.fwhm_min_deg, fwhm_max_deg=args.fwhm_max_deg, is_print=True, ) assign_ka2(peaks) print_peak_table(peaks, "検出されたピーク") n_ka2 = sum(1 for p in peaks if p.ka_role == "ka2") print("") print(f"検出されたCu Kα2ピーク : {n_ka2}") print(f"非Kα2ピーク : {sum(1 for p in peaks if p.ka_role != 'ka2')}") stem = sanitize_stem(infile.stem) plot_path = args.plot_file if args.plot_file is not None else infile.with_name(f"{stem}-peaksearch.png") peak_file = args.peak_file if args.peak_file is not None else default_peak_file_for_input(infile) save_plot = bool(args.save_plot) show_plot = bool(args.show) and not bool(args.no_show) plot_results(x, y, info, peaks, plot_path, show=show_plot, save=save_plot) save_workbook(peak_file, {"peaks": peaks_to_frame(peaks)}) print("") print(f"保存されたピークファイル: {peak_file}") if save_plot: print(f"保存されたプロット : {plot_path}") return peaks, peak_file, plot_path if save_plot else None def run_guess(args: argparse.Namespace, peak_file: Path, base_input: Path | None = None) -> Path: """ 概要: 格子定数推定モードを実行します。 詳細説明: ピーク情報ファイル (peak_file) からピークを読み込み、指定された結晶系と 推定メソッド (args.method) に基づいて格子定数候補を推定します。 推定された候補はコンソールに表示され、最も良い候補に対しては lsq_latt2.py スクリプトを用いた精密化が試みられます。 結果はExcelワークブックとして保存されます。 引数: :param args: コマンドライン引数を格納したNamespaceオブジェクト。 :type args: argparse.Namespace :param peak_file: ピーク情報ファイルへのパス。 :type peak_file: Path :param base_input: 元の入力データファイルへのパス（オプション）。 :type base_input: Path | None 戻り値: :returns: 推定結果Excelファイルへのパス。 :rtype: Path """ crystal_systems = parse_crystal_systems(args.crystal_system) peaks = read_peaks_file(peak_file, wavelength=args.wavelength) print("") print(f"ピークファイル : {peak_file}") print("モード : guess") print(f"メソッド : {args.method}") print(f"結晶系 : {', '.join(crystal_systems)}") print(f"npeak : {args.npeak}") print(f"hmax : {args.hmax}") print_peak_table(peaks, "ファイルから読み込まれたピーク") n_ka2 = sum(1 for p in peaks if p.ka_role == "ka2") print("") print(f"ファイル/推定内のCu Kα2ピーク : {n_ka2}") print(f"非Kα2ピーク : {sum(1 for p in peaks if p.ka_role != 'ka2')}") print("") args.resolved_lattice_max = resolve_lattice_max_limit(args, peaks) candidates, selected = guess_candidates(args, peaks, crystal_systems) print("") print(f"推定に使用される最も強い非Kα2ピーク: {len(selected)}/{args.npeak}") print_peak_table(selected, "格子定数推定のために選択されたピーク") print("") print("格子定数候補") if not candidates: print("候補が見つかりませんでした。--npeak, --hmax を増やすか、--crystal-system auto を試してください。") for i, cand in enumerate(candidates[:10], 1): print(f"{i:2d}. {cand.crystal_system:12s} マッチ数={cand.score_matches:2d} rms_rel={cand.score_rms_rel:.6e} パラメータ={cand.params}") refined_summary = None refined_rows: list[dict] = [] if candidates and args.lsq_script.exists(): refined = refine_with_lsq(args.lsq_script, candidates[0], [r for r in candidates[0].all_matches if r["matched"]], wavelength=args.wavelength) if refined is not None: refined_summary, refined_rows = refined print("") print("精密化された格子定数") print(f" crystal_system : {refined_summary['crystal_system']}") print(f" a = {refined_summary['a']:.6f} ± {refined_summary['sigma_a']:.6f} Å") print(f" b = {refined_summary['b']:.6f} ± {refined_summary['sigma_b']:.6f} Å") print(f" c = {refined_summary['c']:.6f} ± {refined_summary['sigma_c']:.6f} Å") print(f" alpha = {refined_summary['alpha']:.6f} ± {refined_summary['sigma_alpha']:.6f} deg") print(f" beta = {refined_summary['beta']:.6f} ± {refined_summary['sigma_beta']:.6f} deg") print(f" gamma = {refined_summary['gamma']:.6f} ± {refined_summary['sigma_gamma']:.6f} deg") print(f" rms_rel_error_inv_d2 = {refined_summary['rms_rel_error_inv_d2']:.6e}") print(f" rms_rel_error_2theta = {refined_summary['rms_rel_error_2theta']:.6e}") elif candidates: print("") print(f"lsqスクリプトが見つからなかったため、精密化はスキップされました: {args.lsq_script}") sheets: dict[str, pd.DataFrame] = { "all_detected_peaks": peaks_to_frame(peaks), "selected_for_guess": peaks_to_frame(selected), } if candidates: cand_rows = [candidate_to_summary_row(i, c, method=args.method) for i, c in enumerate(candidates, 1)] sheets["candidate_summary"] = pd.DataFrame(cand_rows) sheets["all_peaks_assignment"] = pd.DataFrame(candidates[0].all_matches) sheets["selected_assignment"] = pd.DataFrame(candidates[0].selected_matches) if refined_summary: sheets["refined_summary"] = pd.DataFrame([refined_summary]) sheets["refined_matched"] = pd.DataFrame(refined_rows) out_xlsx = args.output_file if args.output_file is not None else default_guess_file_for_peak_file(peak_file, base_input) save_workbook(out_xlsx, sheets) print("") print(f"保存された推定ファイル: {out_xlsx}") return out_xlsx def _load_peaks_for_compare(args: argparse.Namespace, infile: Path) -> tuple[list[Peak], np.ndarray | None, np.ndarray | None, Path | None]: """ 概要: 比較モードのためにピークとオプションの生のXRDデータを読み込みます。 詳細説明: --peak-file が指定されていれば、そのファイルからピークを読み込みます。 入力ファイルが通常のXRDデータのように見える場合は、ピーク探索を実行してピークを取得します。 入力ファイルが推定ワークブックの場合、そこに含まれるピークシートを読み込みます。 いずれでもない場合、ファイルがピークリストであると仮定して読み込みを試み、 失敗した場合は生XRDデータと見なしてピーク探索を行います。 生のXRDデータがあれば、それも一緒に返します。 引数: :param args: コマンドライン引数を格納したNamespaceオブジェクト。 :type args: argparse.Namespace :param infile: 入力ファイルへのパス。 :type infile: Path 戻り値: :returns: (Peakオブジェクトのリスト, 生のXデータ, 生のYデータ, 使用されたピークファイルパス) のタプル。 :rtype: tuple[list[Peak], numpy.ndarray | None, numpy.ndarray | None, Path | None] """ if args.peak_file is not None: peaks = read_peaks_file(args.peak_file, wavelength=args.wavelength) x, y = try_read_raw_xy_for_plot(infile) return peaks, x, y, args.peak_file # 高密度な2列ファイルは通常、生のXRDプロファイルであり、ピークリストではありません。 if looks_like_raw_xrd_file(infile): peaks, peak_file, _plot = run_search(args, infile) x, y = try_read_raw_xy_for_plot(infile) return peaks, x, y, peak_file # 入力が推定ワークブックの場合、その all_detected_peaks シートを使用します。 if infile.suffix.lower() in {".xlsx", ".xls"}: try: xls = pd.ExcelFile(infile) names = {str(s).lower() for s in xls.sheet_names} if "all_detected_peaks" in names or "peaks" in names: peaks = read_peaks_file(infile, wavelength=args.wavelength) return peaks, None, None, infile except Exception: pass # それ以外の場合、まずピークファイルとして解釈を試みます。失敗した場合は生XRDとして扱い、探索を実行します。 try: peaks = read_peaks_file(infile, wavelength=args.wavelength) return peaks, None, None, infile except Exception: peaks, peak_file, _plot = run_search(args, infile) x, y = try_read_raw_xy_for_plot(infile) return peaks, x, y, peak_file def run_compare(args: argparse.Namespace, infile: Path) -> Path: """ 概要: 比較モードまたは計算モードを実行します。 詳細説明: 観測されたピークと、指定された格子定数候補（推定ファイルから読み込むか、 コマンドライン引数で明示的に指定されたもの）を比較します。 観測ピークと理論的回折線をマッチングし、結果をコンソールに表示します。 また、結果のデータフレームとプロットをファイルに保存します。 引数: :param args: コマンドライン引数を格納したNamespaceオブジェクト。 :type args: argparse.Namespace :param infile: 入力ファイルへのパス。 :type infile: Path 戻り値: :returns: 比較結果Excelファイルへのパス。 :rtype: Path 例外: :raises ValueError: 比較に利用可能なピークがない場合、 または候補が見つからない、またはランクが無効な場合。 """ explicit_lattice = has_explicit_lattice_params(args) mode_label = "calc" if normalize_mode(args.mode) == "calc" or explicit_lattice else "compare" print("") print(f"入力ファイル : {infile}") print(f"モード : {mode_label}") if explicit_lattice: print("候補元 : 明示的な格子定数") else: print(f"候補ランク : {args.candidate_rank}") peaks, x, y, peak_file_used = _load_peaks_for_compare(args, infile) if args.xmin is not None or args.xmax is not None: peaks = [p for p in peaks if (args.xmin is None or p.two_theta >= args.xmin) and (args.xmax is None or p.two_theta <= args.xmax)] if not peaks: raise ValueError("比較に利用可能なピークがありません。") guess_file = args.guess_file if guess_file is None and infile.suffix.lower() in {".xlsx", ".xls"}: try: xls = pd.ExcelFile(infile) if "candidate_summary" in {str(s).lower() for s in xls.sheet_names}: guess_file = infile except Exception: guess_file = None if explicit_lattice: candidate = candidate_from_lattice_args(args) elif guess_file is not None: candidate = read_candidate_from_guess_file(guess_file, args.candidate_rank) else: crystal_systems = parse_crystal_systems(args.crystal_system) if not hasattr(args, "resolved_lattice_max"): print("") args.resolved_lattice_max = resolve_lattice_max_limit(args, peaks) candidates, _selected = guess_candidates(args, peaks, crystal_systems) if not candidates: raise ValueError("比較のための候補が見つかりませんでした。--guess-file、明示的な--a/--b/--c を指定するか、推定オプションを調整してください。") if args.candidate_rank < 1 or args.candidate_rank > len(candidates): raise ValueError(f"候補ランク {args.candidate_rank} は 1..{len(candidates)} の範囲外です。") candidate = candidates[args.candidate_rank - 1] xmin = args.xmin xmax = args.xmax if xmin is None and x is not None and len(x): xmin = float(np.nanmin(x)) if xmax is None and x is not None and len(x): xmax = float(np.nanmax(x)) theoretical_df, assignment_df = compare_candidate_with_peaks( candidate, peaks, hmax=args.hmax, wavelength=args.wavelength, xmin=xmin, xmax=xmax, tolerance_deg=args.tolerance_deg, ) print("") print("使用された候補") print(f" crystal_system : {candidate.crystal_system}") for k in ["a", "b", "c", "alpha", "beta", "gamma"]: if k in candidate.params: print(f" {k:5s} = {candidate.params[k]:.8g}") if not assignment_df.empty: matched = assignment_df[assignment_df["matched"] == True] print("") print(f"比較された観測ピーク : {len(assignment_df)}") print(f"許容誤差内でマッチした数: {len(matched)}") if len(matched): rms = float(np.sqrt(np.mean(np.square(matched["resid_two_theta_deg"].to_numpy(dtype=float))))) print(f"RMS残差2θ : {rms:.6g} deg") cols = ["peak_id", "two_theta_obs_deg", "hkl", "two_theta_calc_deg", "resid_two_theta_deg", "intensity", "fwhm_deg", "matched"] with pd.option_context("display.max_rows", 300, "display.width", 180): print(assignment_df[cols].to_string(index=False)) stem = sanitize_stem(strip_peak_suffix((peak_file_used or infile).stem)) suffix = "calc" if explicit_lattice else f"compare-rank{args.candidate_rank}" out_xlsx = args.output_file if args.output_file is not None else (peak_file_used or infile).with_name(f"{stem}-{suffix}.xlsx") plot_path = args.plot_file if args.plot_file is not None else out_xlsx.with_suffix(".png") sheets = { "candidate_used": pd.DataFrame([candidate_to_summary_row(args.candidate_rank if not explicit_lattice else 0, candidate, method="explicit" if explicit_lattice else "")]), "observed_peaks": peaks_to_frame(peaks), "theoretical_lines": theoretical_df, "observed_vs_calc": assignment_df, } save_workbook(out_xlsx, sheets) save_plot = bool(args.save_plot) show_plot = bool(args.show) and not bool(args.no_show) plot_theoretical_df = filter_theoretical_lines_for_plot( theoretical_df, mode=args.plot_theory, near_deg=args.plot_near_deg if args.plot_near_deg is not None else args.tolerance_deg, ) plot_compare_results(x, y, peaks, plot_theoretical_df, assignment_df, candidate, xmin=xmin, xmax=xmax, outpath=plot_path, show=show_plot, save=save_plot) print("") print(f"理論線 : 計算された線 {len(theoretical_df)}本、描画された線 {len(plot_theoretical_df)}本 (--plot-theory {args.plot_theory})") print(f"保存された比較ファイル: {out_xlsx}") if save_plot: print(f"保存されたプロット : {plot_path}") if show_plot: print("ホバー注釈 : mplcursorsがインストールされている場合有効") return out_xlsx def main() -> int: """ 概要: スクリプトのメインエントリポイントです。 詳細説明: コマンドライン引数を解析し、指定されたモード（search, guess, all, compare, calc）に基づいて 適切な処理関数を呼び出します。エラーが発生した場合は、エラーメッセージを標準エラー出力に表示します。 戻り値: :returns: 成功した場合は0、エラーの場合は非0の終了コード。 :rtype: int 例外: :raises SystemExit: コマンドライン引数エラーや処理中の例外が発生した場合。 """ args = build_argument_parser().parse_args() try: if args.no_show: args.show = 0 mode, infile = resolve_mode_and_input(args) print(f"要求されたモード : {mode}") if mode == "search": run_search(args, infile) return 0 if mode == "guess": peak_file = args.peak_file if args.peak_file is not None else infile run_guess(args, peak_file) return 0 if mode in {"compare", "calc"}: if mode == "calc" and not has_explicit_lattice_params(args): raise ValueError("--mode calc には明示的な格子定数が必要です (例: --crystal-system cubic --a 3.905)") run_compare(args, infile) return 0 if mode == "all": _peaks, peak_file, _plot_path = run_search(args, infile) guess_file = run_guess(args, peak_file, base_input=infile) # --mode all は探索 + 推定のみです。選択されたランクを比較するには --mode compare と --guess-file を使用してください。 return 0 raise SystemExit(f"サポートされていないモード: {mode}") except Exception as exc: print(f"エラー: {exc}", file=sys.stderr) raise input("\n続行するにはEnterを押してください>>\n") if __name__ == "__main__": raise SystemExit(main())