"""
最小二乗法を用いた格子定数計算スクリプト。
このスクリプトは、X線回折データなどの回折データから結晶の格子定数を導出します。
様々な格子系に対応し、観測データの前処理、重み付け、最小二乗フィッティング、
そして最終的な格子定数（逆格子および実格子）の導出を行います。
外れ値の自動除去機能も備えています。
詳細説明:
入力ファイルは特定のフォーマットに従い、複数の計算ジョブを連続して処理できます。
各ジョブでは、格子系、補正フラグ、重み付けモード、回折位置の変換モードなどが
設定され、それに従って観測データが処理されます。
重み付き最小二乗法によりパラメータが決定され、その結果から逆格子定数と実格子定数、
およびそれらの標準偏差が計算されます。計算結果は整形されたレポートファイルとして出力されます。
関連リンク:
:doc:`lsq_latt2_usage`
"""
from __future__ import annotations
import argparse
import math
from dataclasses import dataclass, field
from pathlib import Path
from typing import Callable, Iterable, TextIO
import numpy as np
DEG_PER_RAD = 180.0 / math.pi
RAD_PER_DEG = math.pi / 180.0
ANGLE_EPS = 1e-12
PERTURB_SCALE = 1e-4
OUTLIER_SIGMA = 3.0
[ドキュメント]
class LatticeCalculationError(RuntimeError):
"""
格子定数計算におけるエラーを示す例外。
この例外は、入力データが無効な場合、または最小二乗問題の解決中に
計算上の問題（例: 特異行列）が発生した場合に発生します。
"""
[ドキュメント]
@dataclass(slots=True)
class Observation:
"""
一つの回折観測データを保持するデータクラス。
詳細説明:
このクラスは、特定の回折面(hkl)に対応する生の回折位置、変換された位置、
計算で使用される重み、デザイン行列の行、および観測時の波長などの
詳細な情報を含みます。
:param serial_index: int: 観測データの通し番号。
:param h: int: ミラー指数 h。
:param k: int: ミラー指数 k。
:param l: int: ミラー指数 l。
:param raw_position: float: 生の回折位置データ。
:param transformed_position: float: 変換された回折位置データ（最小二乗法に使用）。
:param weight: float: この観測データに対する重み。
:param design_row: np.ndarray: この観測データに対応するデザイン行列の行。
:param wavelength: float: この観測で使用されたX線の波長。
"""
serial_index: int
h: int
k: int
l: int
raw_position: float
transformed_position: float
weight: float
design_row: np.ndarray
wavelength: float
[ドキュメント]
@dataclass(slots=True)
class LeastSquaresResult:
"""
最小二乗法の結果を保持するデータクラス。
詳細説明:
最小二乗法によって計算された格子定数パラメータとその標準偏差、
観測データへのフィット値、残差、観測の標準偏差、および
計算に残された観測データのリストと共分散行列を格納します。
:param parameters: np.ndarray: 最小二乗法によって決定されたパラメータの配列。
:param parameter_sigma: np.ndarray: 各パラメータの標準偏差の配列。
:param fitted_y: np.ndarray: 各観測データに対してフィットされた予測値。
:param residuals: np.ndarray: 各観測データに対する残差 (観測値 - フィット値)。
:param observation_sigma: np.ndarray: 各観測データの標準偏差。
:param kept_observations: list[Observation]: 外れ値除去後に残った観測データのリスト。
:param covariance: np.ndarray: パラメータ間の共分散行列。
"""
parameters: np.ndarray
parameter_sigma: np.ndarray
fitted_y: np.ndarray
residuals: np.ndarray
observation_sigma: np.ndarray
kept_observations: list[Observation]
covariance: np.ndarray
[ドキュメント]
@dataclass(slots=True)
class CellConstants:
"""
格子定数（逆格子および実格子）を保持するデータクラス。
詳細説明:
逆格子と実格子の両方について、長さ、角度（度数）、コサイン、
体積およびそれぞれの標準偏差を格納します。
デフォルト値としてゼロまたは標準的な値（角度90度）が設定されています。
:param reciprocal_lengths: np.ndarray: 逆格子の長さ (a*, b*, c*) の配列。デフォルトは全て0.0。
:param reciprocal_length_sigma: np.ndarray: 逆格子の長さの標準偏差の配列。デフォルトは全て0.0。
:param reciprocal_angles_deg: np.ndarray: 逆格子の角度 (alpha*, beta*, gamma*) の配列（度数）。デフォルトは全て90.0。
:param reciprocal_angle_sigma_deg: np.ndarray: 逆格子の角度の標準偏差の配列（度数）。デフォルトは全て0.0。
:param reciprocal_cosines: np.ndarray: 逆格子の角度のコサイン (cos(alpha*), cos(beta*), cos(gamma*)) の配列。デフォルトは全て0.0。
:param reciprocal_cosine_sigma: np.ndarray: 逆格子の角度のコサインの標準偏差の配列。デフォルトは全て0.0。
:param reciprocal_volume: float: 逆格子の体積。デフォルトは0.0。
:param reciprocal_volume_sigma: float: 逆格子の体積の標準偏差。デフォルトは0.0。
:param direct_lengths: np.ndarray: 実格子の長さ (a, b, c) の配列。デフォルトは全て0.0。
:param direct_length_sigma: np.ndarray: 実格子の長さの標準偏差の配列。デフォルトは全て0.0。
:param direct_angles_deg: np.ndarray: 実格子の角度 (alpha, beta, gamma) の配列（度数）。デフォルトは全て90.0。
:param direct_angle_sigma_deg: np.ndarray: 実格子の角度の標準偏差の配列（度数）。デフォルトは全て0.0。
:param direct_cosines: np.ndarray: 実格子の角度のコサイン (cos(alpha), cos(beta), cos(gamma)) の配列。デフォルトは全て0.0。
:param direct_cosine_sigma: np.ndarray: 実格子の角度のコサインの標準偏差の配列。デフォルトは全て0.0。
:param direct_volume: float: 実格子の体積。デフォルトは0.0。
:param direct_volume_sigma: float: 実格子の体積の標準偏差。デフォルトは0.0。
"""
reciprocal_lengths: np.ndarray = field(default_factory=lambda: np.zeros(3, dtype=float))
reciprocal_length_sigma: np.ndarray = field(default_factory=lambda: np.zeros(3, dtype=float))
reciprocal_angles_deg: np.ndarray = field(default_factory=lambda: np.full(3, 90.0, dtype=float))
reciprocal_angle_sigma_deg: np.ndarray = field(default_factory=lambda: np.zeros(3, dtype=float))
reciprocal_cosines: np.ndarray = field(default_factory=lambda: np.zeros(3, dtype=float))
reciprocal_cosine_sigma: np.ndarray = field(default_factory=lambda: np.zeros(3, dtype=float))
reciprocal_volume: float = 0.0
reciprocal_volume_sigma: float = 0.0
direct_lengths: np.ndarray = field(default_factory=lambda: np.zeros(3, dtype=float))
direct_length_sigma: np.ndarray = field(default_factory=lambda: np.zeros(3, dtype=float))
direct_angles_deg: np.ndarray = field(default_factory=lambda: np.full(3, 90.0, dtype=float))
direct_angle_sigma_deg: np.ndarray = field(default_factory=lambda: np.zeros(3, dtype=float))
direct_cosines: np.ndarray = field(default_factory=lambda: np.zeros(3, dtype=float))
direct_cosine_sigma: np.ndarray = field(default_factory=lambda: np.zeros(3, dtype=float))
direct_volume: float = 0.0
direct_volume_sigma: float = 0.0
[ドキュメント]
@dataclass(slots=True)
class JobConfig:
"""
格子定数計算ジョブの設定を保持するデータクラス。
詳細説明:
このクラスは、計算に使用する格子系の種類、適用する補正の種類、
観測データの重み付け方法、回折位置の変換方法、およびプロファイルモードの
設定情報をカプセル化します。
:param lattice_system: int: 格子系の種類を示す整数コード（例: 1=三斜晶, 4=直方晶など）。
:param correction_flags: list[int]: 5つの補正フラグのリスト。各フラグは特定の補正の適用有無を示す。
:param weight_mode: int: 観測データに適用する重み付けモード。
:param position_mode: int: 生の回折位置データを変換するモード。
:param profile_mode: int: プロファイル分析モード（このコードでは未使用だが、入力形式の一部）。
"""
lattice_system: int
correction_flags: list[int]
weight_mode: int
position_mode: int
profile_mode: int
[ドキュメント]
@dataclass(slots=True)
class JobData:
"""
一つの格子定数計算ジョブに関する全てのデータを保持するデータクラス。
詳細説明:
ジョブのタイトル、設定、処理する観測データのリスト、および
参照として使用される波長を含みます。
:param title: str: ジョブのタイトル。
:param config: JobConfig: このジョブの設定情報。
:param observations: list[Observation]: このジョブで処理される観測データのリスト。
:param reference_wavelength: float: このジョブの参照波長。
"""
title: str
config: JobConfig
observations: list[Observation]
reference_wavelength: float
[ドキュメント]
class LineReader:
"""
テキストストリームから行を読み込むためのヘルパークラス。
詳細説明:
このクラスは、テキストファイルを一行ずつ読み込み、現在の行番号を追跡し、
ファイルの予期せぬ終端を検出する機能を提供します。
空行でない行を確実に読み込むためのメソッドも備えています。
:param stream: TextIO: 読み込むテキストストリーム。
"""
def __init__(self, stream: TextIO):
"""
LineReaderの新しいインスタンスを初期化する。
:param stream: TextIO: 読み込むテキストストリーム。
"""
self.stream = stream
self.line_number = 0
[ドキュメント]
def read_line(self, *, allow_eof: bool = False) -> str | None:
"""
ストリームから一行を読み込む。
:param allow_eof: bool: Trueの場合、ファイルの終端に達してもEOFErrorを発生させず、Noneを返す。デフォルトはFalse。
:returns: str | None: 読み込まれた行（末尾の改行文字は除去）。EOFでallow_eofがTrueの場合はNone。
:raises EOFError: `allow_eof`がFalseのときに、ファイルの予期せぬ終端に達した場合。
"""
line = self.stream.readline()
if line == "":
if allow_eof:
return None
raise EOFError("unexpected end of file")
self.line_number += 1
return line.rstrip("\n")
[ドキュメント]
def read_nonempty_line(self) -> str:
"""
ストリームから空でない行を読み込む。
:returns: str: 読み込まれた空でない行（末尾の改行文字は除去）。
:raises EOFError: ファイルの終端に達した場合。
:raises AssertionError: 読み込まれた行がNoneであった場合 (内部的なガード)。
"""
line = self.read_line()
assert line is not None
return line
[ドキュメント]
class SigmaPropagator:
"""
誤差伝播の計算を行うクラス。
詳細説明:
このクラスは、複数の入力値とそれぞれの標準偏差が与えられたときに、
それらの入力値に依存する関数が返す値とその標準偏差を、
数値微分（摂動法）を用いて計算します。
:param values: Iterable[float]: 関数に渡す入力値のシーケンス。
:param sigmas: Iterable[float]: 各入力値の標準偏差のシーケンス。
"""
def __init__(self, values: Iterable[float], sigmas: Iterable[float]):
"""
SigmaPropagatorの新しいインスタンスを初期化する。
:param values: Iterable[float]: 関数に渡す入力値のシーケンス。
:param sigmas: Iterable[float]: 各入力値の標準偏差のシーケンス。
"""
self.values = np.array(list(values), dtype=float)
self.sigmas = np.array(list(sigmas), dtype=float)
[ドキュメント]
def evaluate(self, func: Callable[[np.ndarray], float]) -> tuple[float, float]:
"""
与えられた関数を評価し、その値と標準偏差を計算する。
詳細説明:
入力値の各要素を微小に摂動させ、その変化から数値微分を計算します。
これらの数値微分と各入力値の標準偏差を用いて、誤差伝播の法則により
関数の結果の標準偏差を導出します。
:param func: Callable[[np.ndarray], float]: 誤差伝播を計算する関数。
NumPy配列を受け取り、単一のfloatを返す必要があります。
:returns: tuple[float, float]: 関数の評価値と、その標準偏差のタプル。
"""
base = func(self.values.copy())
variance = 0.0
for idx, sigma in enumerate(self.sigmas):
if abs(sigma) < ANGLE_EPS:
continue
shifted = self.values.copy()
shifted[idx] += PERTURB_SCALE * sigma
deriv_sigma = (func(shifted) - base) / PERTURB_SCALE
variance += deriv_sigma * deriv_sigma
return base, math.sqrt(max(variance, 0.0))
# ---------- elementary math helpers ----------
[ドキュメント]
def clamp_unit(x: float) -> float:
"""
値を -1.0 から 1.0 の範囲にクランプする。
詳細説明:
主に `math.acos` の引数が有効な範囲内にあることを保証するために使用されます。
入力値が -1.0 未満であれば -1.0 に、1.0 を超えれば 1.0 に丸めます。
:param x: float: クランプする値。
:returns: float: -1.0 から 1.0 の範囲にクランプされた値。
"""
return max(-1.0, min(1.0, x))
[ドキュメント]
def safe_sqrt(x: float) -> float:
"""
非負の引数に対して安全に平方根を計算する。
詳細説明:
`math.sqrt` は負の引数に対して `ValueError` を発生させますが、
この関数は引数が負の場合でも0.0として処理することでエラーを回避し、
常に非負の値を返します。
:param x: float: 平方根を計算する値。
:returns: float: xの平方根。xが負の場合は0.0。
"""
return math.sqrt(max(0.0, x))
[ドキュメント]
def acos_deg(cosine: float) -> float:
"""
コサイン値から角度を度数で計算する。
詳細説明:
入力されたコサイン値が有効な範囲 [-1, 1] 内にあることを保証するために、
`clamp_unit` 関数で値をクランプしてから `math.acos` を適用します。
結果はラジアンから度数に変換されます。
:param cosine: float: コサイン値。
:returns: float: 度数で表された角度。
"""
return math.degrees(math.acos(clamp_unit(cosine)))
[ドキュメント]
def reciprocal_volume_from_cosines(cosines: np.ndarray) -> float:
"""
逆格子のコサイン値から逆格子単位胞体積の計算に用いる係数を計算する。
詳細説明:
逆格子単位胞体積 V* は、V* = a*b*c* * sqrt(1 - cos^2(alpha*) - cos^2(beta*) - cos^2(gamma*) + 2cos(alpha*)cos(beta*)cos(gamma*))
で与えられます。この関数は、上記式の `sqrt(...)` の部分を計算します。
:param cosines: np.ndarray: 逆格子の角度のコサイン値 (cos_alpha_star, cos_beta_star, cos_gamma_star) を含むNumPy配列。
:returns: float: 逆格子単位胞体積計算に使用される平方根の値。
"""
ca, cb, cg = cosines
term = 1.0 - ca * ca - cb * cb - cg * cg + 2.0 * ca * cb * cg
return safe_sqrt(term)
[ドキュメント]
def direct_cosine_from_reciprocal(values: np.ndarray) -> float:
"""
逆格子のコサイン値から実格子のコサイン値を計算する。
詳細説明:
逆格子の角度 (alpha*, beta*, gamma*) のコサイン値から、
対応する実格子の角度 (alpha, beta, gamma) のコサイン値を導出します。
例えば、実格子の cos(alpha) は逆格子の cos(alpha*), cos(beta*), cos(gamma*) から計算されます。
:param values: np.ndarray: (cos_alpha_star, cos_beta_star, cos_gamma_star) を含むNumPy配列。
:returns: float: 対応する実格子のコサイン値。
"""
cos_alpha_star, cos_beta_star, cos_gamma_star = values
sin_beta_star = safe_sqrt(1.0 - cos_beta_star * cos_beta_star)
sin_gamma_star = safe_sqrt(1.0 - cos_gamma_star * cos_gamma_star)
denom = sin_beta_star * sin_gamma_star
if denom <= ANGLE_EPS:
return 0.0
return (cos_beta_star * cos_gamma_star - cos_alpha_star) / denom
[ドキュメント]
def direct_length_from_reciprocal(values: np.ndarray) -> float:
"""
逆格子の長さとコサイン値から実格子の長さを計算する。
詳細説明:
逆格子定数 (a*, cos_alpha*, cos_beta*, cos_gamma*) から、
対応する実格子の長さ (a) を導出します。
例: 実格子の長さaは、逆格子のa*, cos(alpha*), cos(beta*), cos(gamma*) から計算されます。
:param values: np.ndarray: (a_star, cos_alpha_star, cos_beta_star, cos_gamma_star) を含むNumPy配列。
:returns: float: 対応する実格子の長さ。
"""
a_star, cos_alpha_star, cos_beta_star, cos_gamma_star = values
if abs(a_star) <= ANGLE_EPS:
return 0.0
sin_alpha_star = safe_sqrt(1.0 - cos_alpha_star * cos_alpha_star)
v_star = reciprocal_volume_from_cosines(np.array([cos_alpha_star, cos_beta_star, cos_gamma_star]))
if v_star <= ANGLE_EPS:
return 0.0
return sin_alpha_star / (a_star * v_star)
[ドキュメント]
def direct_volume_from_reciprocal(values: np.ndarray) -> float:
"""
逆格子の長さとコサイン値から実格子の体積を計算する。
詳細説明:
逆格子定数 (a*, b*, c*, cos_alpha*, cos_beta*, cos_gamma*) から
実格子の体積を導出します。実格子の体積は逆格子の体積の逆数です。
:param values: np.ndarray: (a_star, b_star, c_star, cos_alpha_star, cos_beta_star, cos_gamma_star) を含むNumPy配列。
:returns: float: 実格子の体積。
"""
a_star, b_star, c_star, cos_alpha_star, cos_beta_star, cos_gamma_star = values
v_star = a_star * b_star * c_star * reciprocal_volume_from_cosines(
np.array([cos_alpha_star, cos_beta_star, cos_gamma_star])
)
if v_star <= ANGLE_EPS:
return 0.0
return 1.0 / v_star
[ドキュメント]
def trigonal_direct_length(values: np.ndarray) -> float:
"""
三方晶系の逆格子定数から実格子の長さを計算する。
詳細説明:
三方晶系において、逆格子のa*とcos(alpha*)から実格子の長さaを導出します。
等方的な格子定数を持つ三方晶系（菱面体晶）の特殊なケースに適用されます。
:param values: np.ndarray: (a_star, cos_alpha_star) を含むNumPy配列。
:returns: float: 三方晶系における実格子の長さa。
"""
a_star, cos_alpha_star = values
if abs(a_star) <= ANGLE_EPS:
return 0.0
numerator = 1.0 - cos_alpha_star * cos_alpha_star
denominator = 1.0 - 3.0 * cos_alpha_star * cos_alpha_star + 2.0 * cos_alpha_star**3
if denominator <= ANGLE_EPS or numerator < 0.0:
return 0.0
return safe_sqrt(numerator / denominator) / a_star
[ドキュメント]
def trigonal_direct_cosine(values: np.ndarray) -> float:
"""
三方晶系の逆格子のコサイン値から実格子のコサイン値を計算する。
詳細説明:
三方晶系において、逆格子のcos(alpha*)から実格子のcos(alpha)を導出します。
:param values: np.ndarray: (cos_alpha_star,) を含むNumPy配列。
:returns: float: 三方晶系における実格子のコサイン値cos(alpha)。
"""
(cos_alpha_star,) = values
denominator = 1.0 - cos_alpha_star * cos_alpha_star
if abs(denominator) <= ANGLE_EPS:
return 0.0
return cos_alpha_star * (cos_alpha_star - 1.0) / denominator
[ドキュメント]
def trigonal_direct_volume(values: np.ndarray) -> float:
"""
三方晶系の逆格子定数から実格子の体積を計算する。
詳細説明:
三方晶系において、逆格子のa*とcos(alpha*)から実格子の体積Vを導出します。
:param values: np.ndarray: (a_star, cos_alpha_star) を含むNumPy配列。
:returns: float: 三方晶系における実格子の体積V。
"""
a_star, cos_alpha_star = values
term = 1.0 - 3.0 * cos_alpha_star * cos_alpha_star + 2.0 * cos_alpha_star**3
if term <= 0.0:
return 0.0
v_star = a_star**3 * safe_sqrt(term)
if v_star <= ANGLE_EPS:
return 0.0
return 1.0 / v_star
# ---------- model construction ----------
[ドキュメント]
def base_design_terms(h: int, k: int, l: int) -> np.ndarray:
"""
回折指数 (h, k, l) からデザイン行列の基本的な項を計算する。
詳細説明:
格子定数計算におけるデザイン行列を構築するための基本要素として、
h^2, k^2, l^2, 2kl, 2lh, 2hk の6つの項を計算します。
これらの項は、一般的な三斜晶系を含む様々な格子系のデザイン行列の基礎となります。
:param h: int: ミラー指数h。
:param k: int: ミラー指数k。
:param l: int: ミラー指数l。
:returns: np.ndarray: 6つのデザイン行列項を含むNumPy配列。
"""
return np.array(
[
float(h * h),
float(k * k),
float(l * l),
2.0 * float(k * l),
2.0 * float(l * h),
2.0 * float(h * k),
],
dtype=float,
)
[ドキュメント]
def lattice_design_row(lattice_system: int, h: int, k: int, l: int) -> np.ndarray:
"""
指定された格子系とミラー指数に基づいて、デザイン行列の行を構築する。
詳細説明:
各格子系（三斜晶、単斜晶、直方晶、正方晶、六方晶、立方晶、菱面体晶）は、
結晶学的に異なる格子定数の制約を持ちます。この関数は、その制約に基づいて、
`base_design_terms` から必要な項を選択・結合し、デザイン行列の1行を生成します。
格子系のコードと対応するパラメータ:
- 1: 三斜晶 (triclinic): a*, b*, c*, cos(alpha*), cos(beta*), cos(gamma*)
- 2: 単斜晶 (monoclinic, unique b): a*, b*, c*, cos(beta*)
- 3: 単斜晶 (monoclinic, unique c): a*, b*, c*, cos(gamma*) (※コードではgamma*として実装されているが、入力形式によってはalpha*を意味する場合もある)
- 4: 直方晶 (orthorhombic): a*, b*, c*
- 5: 正方晶 (tetragonal): a*, c*
- 6: 立方晶 (cubic): a*
- 7: 菱面体晶 (rhombohedral, trigonal): a*, cos(alpha*) (等軸菱面体)
- 8: 六方晶 (hexagonal): a*, c*
:param lattice_system: int: 格子系の種類を示す整数コード（1-8）。
:param h: int: ミラー指数h。
:param k: int: ミラー指数k。
:param l: int: ミラー指数l。
:returns: np.ndarray: 指定された格子系のデザイン行列の行。
:raises LatticeCalculationError: 未対応の格子系が指定された場合。
"""
x = base_design_terms(h, k, l)
if lattice_system == 1: # Triclinic
return x.copy()
if lattice_system == 2: # Monoclinic (unique b axis, beta*)
return np.array([x[0], x[1], x[2], x[4]], dtype=float) # h^2, k^2, l^2, 2lh
if lattice_system == 3: # Monoclinic (unique c axis, gamma*) or (unique a axis, alpha*)
return np.array([x[0], x[1], x[2], x[5]], dtype=float) # h^2, k^2, l^2, 2hk
if lattice_system == 4: # Orthorhombic
return np.array([x[0], x[1], x[2]], dtype=float) # h^2, k^2, l^2
if lattice_system == 5: # Tetragonal
return np.array([x[0] + x[1], x[2]], dtype=float) # h^2+k^2, l^2
if lattice_system == 6: # Cubic
return np.array([x[0] + x[1] + x[2]], dtype=float) # h^2+k^2+l^2
if lattice_system == 7: # Rhombohedral (trigonal)
return np.array([x[0] + x[1] + x[2], x[3] + x[4] + x[5]], dtype=float) # h^2+k^2+l^2, 2(kl+lh+hk)
if lattice_system == 8: # Hexagonal
return np.array([x[0] + x[1] + 0.5 * x[5], x[2]], dtype=float) # h^2+k^2+hk, l^2
raise LatticeCalculationError(f"unsupported lattice system: {lattice_system}")
[ドキュメント]
def calc_weight(weight_mode: int, transformed: float, sigma_or_weight: float | None) -> float:
"""
指定された重みモードと変換された位置データに基づいて観測の重みを計算する。
詳細説明:
最小二乗法では、各観測データの信頼度に応じて重みを設定することで、
より信頼性の高いデータに大きな影響を与えることができます。
この関数は、異なる重み付けスキーム（一定、測定誤差に基づくなど）を実装しています。
重みモードの種類:
- 0: シグマ・シータに基づく重み。
- 1: シグマ・ツーシータに基づく重み。
- 2: 一定の重み (1.0)。
- 3: 直接指定された重み。
- 4: シグマ値の逆二乗に基づく重み (1/sigma^2)。
:param weight_mode: int: 重み計算モード（0-4）。
:param transformed: float: 変換された回折位置データ。
:param sigma_or_weight: float | None: モードによって、回折角度の標準偏差または直接指定される重み値。Noneの場合もある。
:returns: float: 計算された観測の重み。
:raises LatticeCalculationError: 重みモードに必要な追加データが提供されない場合、または未対応のモードが指定された場合。
"""
bas = transformed * transformed # transformed position squared
if weight_mode == 0:
if sigma_or_weight is None:
raise LatticeCalculationError("weight mode 0 requires an additional sigma line")
sigma_theta = 0.25 * sigma_or_weight
sigma_theta = max(sigma_theta, 0.25) # Ensure minimum sigma_theta
denom = bas * (1.0 - bas)
# Original Fortran: 1/(sigma_theta^2 * 4 * transformed^2 * (1-transformed^2))
# This implementation: sigma_theta / denom, which might be a simplified constant, or a typo from 1/(sigma_theta^2 * X)
# However, following the existing code logic:
# The weight should be proportional to 1/variance, so usually 1/(sigma^2).
# It's possible that `sigma_theta` here represents the numerator factor for 1/variance.
# Given the existing code is `sigma_theta / denom`
return 0.0 if abs(denom) <= ANGLE_EPS else sigma_theta / denom
if weight_mode == 1:
if sigma_or_weight is None:
raise LatticeCalculationError("weight mode 1 requires an additional sigma line")
sigma_two_theta = 2.0 * sigma_or_weight # sigma_two_theta is related to two_theta
denom = bas * (1.0 - bas) * sigma_two_theta * sigma_two_theta
# Original Fortran: 1/(sigma_2theta^2 * 4 * transformed^2 * (1-transformed^2))
# This implementation: 0.25 / denom, which is consistent with the Fortran `1 / (4 * denom)` if denom already includes sigma^2
return 0.0 if abs(denom) <= ANGLE_EPS else 0.25 / denom
if weight_mode == 2:
return 1.0
if weight_mode == 3:
if sigma_or_weight is None:
raise LatticeCalculationError("weight mode 3 requires an additional weight line")
return sigma_or_weight if sigma_or_weight > 0.0 else 1.0
if weight_mode == 4:
if sigma_or_weight is None:
raise LatticeCalculationError("weight mode 4 requires an additional sigma line")
# Weight = 1/sigma^2
return 0.0 if sigma_or_weight == 0.0 else 1.0 / (sigma_or_weight * sigma_or_weight)
raise LatticeCalculationError(f"unsupported weight mode: {weight_mode}")
[ドキュメント]
def correction_terms(flags: list[int], transformed: float) -> list[float]:
"""
指定された補正フラグと変換された位置データに基づいて補正項を計算する。
詳細説明:
回折データの測定には、様々な系統誤差（例:吸収、分極、ローレンツ因子など）
が影響を与える可能性があります。この関数は、これらの誤差を補正するための
特定の補正項を計算し、最小二乗法のデザイン行列に追加できるようにします。
:param flags: list[int]: 5つの補正フラグのリスト。各フラグは特定の補正の有無を示す（0:なし、1:あり）。
:param transformed: float: 変換された回折位置データ（通常はsin(θ)）。
:returns: list[float]: 計算された補正項のリスト。
"""
bas = transformed * transformed # sin(theta)^2
theta = math.asin(clamp_unit(transformed)) # theta in radians
terms: list[float] = []
if flags[0]:
# Correction 1: absorption? or 1/sin(2θ)? The original Fortran code suggests a zero shift correction
# A.k.a. "dθ_c" in some contexts
terms.append((math.pi / 2.0 - theta) * math.tan(math.pi / 2.0 - theta))
if flags[1]:
# Correction 2: polarization factor, or L.P. factor component
# sin^2(2θ) term, or 4 * sin^2(θ) * cos^2(θ)
terms.append(4.0 * bas * (1.0 - bas))
if flags[2]:
# Correction 3: usually related to sample transparency or beam divergence
# The term 1/transformed + 1/theta is unusual.
# If transformed = sin(theta), then 1/sin(theta) + 1/theta.
# This term is only added if transformed and theta are not extremely small.
if abs(transformed) <= ANGLE_EPS or abs(theta) <= ANGLE_EPS:
terms.append(0.0)
else:
terms.append(4.0 * bas * (1.0 - bas) * (1.0 / transformed + 1.0 / theta))
if flags[3]:
# Correction 4: Often related to preferred orientation or sample height error
# sin(2θ) term
terms.append(math.sin(2.0 * theta))
if flags[4]:
# Correction 5: Another component, possibly related to preferred orientation or asymmetric peak broadening
# sin(2θ) * cos(θ) term
terms.append(math.sin(2.0 * theta) * math.cos(theta))
return terms
[ドキュメント]
def read_numeric_line(reader: LineReader, expected: int) -> list[float]:
"""
リーダーから一行を読み込み、指定された数の浮動小数点数に変換する。
詳細説明:
入力ファイルから設定値やパラメータを読み込む際に使用されるヘルパー関数。
読み込んだ行をスペースで分割し、最初の`expected`個の要素を浮動小数点数に変換します。
:param reader: LineReader: 入力ファイルを読み込むためのLineReaderインスタンス。
:param expected: int: 期待される数値の最小数。
:returns: list[float]: 読み込まれた数値のリスト。
:raises LatticeCalculationError: 期待されるよりも少ない数値が読み込まれた場合。
"""
raw = reader.read_nonempty_line().split()
if len(raw) < expected:
raise LatticeCalculationError(
f"line {reader.line_number}: expected at least {expected} numeric values, got {len(raw)}"
)
return [float(x) for x in raw[:expected]]
[ドキュメント]
def parse_jobs(stream: TextIO) -> list[JobData]:
"""
入力ストリームから複数の格子定数計算ジョブを解析する。
詳細説明:
入力ファイルは、一連の計算ジョブで構成されます。各ジョブはタイトル、
計算設定、および複数の回折観測データを含みます。この関数は、
ストリームの終端に達するか、特定の終了コード（`h=1000`で次のジョブ、`job_code=0`で全体終了）
が読み込まれるまで、ジョブを順次解析して `JobData` オブジェクトのリストを生成します。
入力ファイルの主要なセクション:
1. ジョブタイトル
2. ジョブ設定パラメータ（格子系、補正フラグ、モードなど）
3. 波長と装置半径
4. 回折観測データのブロック
- h k l raw_position
- (position_mode=3の場合) zero_shift
- (weight_mode in {0,1,3,4}の場合) sigma_or_weight
- 特殊行: `2000 wavelength radius` で波長と半径を更新
- 特殊行: `1000 ...` で現在のジョブの終了を示す
5. ジョブ終了コード (0: 全体の終了, 1: 次のジョブへ)
:param stream: TextIO: 入力データを含むテキストストリーム。
:returns: list[JobData]: 解析されたJobDataオブジェクトのリスト。
:raises EOFError: ファイルの予期せぬ終端に達した場合。
:raises LatticeCalculationError: 入力ファイルの形式が無効な場合（例: ヘッダー情報不足、反射行のデータ不足）。
"""
reader = LineReader(stream)
jobs: list[JobData] = []
while True:
title = reader.read_line(allow_eof=True)
if title is None:
break # EOF encountered
title = title.strip()
if title == "":
continue # Skip empty title lines
# Read job configuration parameters
params = reader.read_nonempty_line().split()
if len(params) < 9:
raise LatticeCalculationError(f"line {reader.line_number}: invalid job header")
lattice_system, *raw_flags, weight_mode, position_mode, profile_mode = map(int, params[:9])
correction_flags = raw_flags[:5] # Take first 5 flags
config = JobConfig(
lattice_system=lattice_system,
correction_flags=correction_flags,
weight_mode=weight_mode,
position_mode=position_mode,
profile_mode=profile_mode,
)
# Read initial wavelength and radius
wavelength, radius = read_numeric_line(reader, 2)
reference_wavelength = wavelength
observations: list[Observation] = []
serial_index = 0
# Read reflection data
line = reader.read_nonempty_line()
while True:
parts = line.split()
if len(parts) < 4:
raise LatticeCalculationError(f"line {reader.line_number}: invalid reflection line")
h, k, l = map(int, parts[:3])
raw_position = float(parts[3])
if not (h < 1000 or h == 2000): # Check for job termination (h >= 1000 and h != 2000)
break
if h == 2000: # Special line: update wavelength and radius
wavelength, radius = read_numeric_line(reader, 2)
if len(observations) == 0: # If this is the first wavelength for the job
reference_wavelength = wavelength
line = reader.read_nonempty_line()
continue
serial_index += 1
zero_shift = None
if position_mode == 3:
zero_shift = float(reader.read_nonempty_line().strip())
transformed = transformed_position(raw_position, position_mode, wavelength, radius, zero_shift)
sigma_or_weight = None
if weight_mode in {0, 1, 3, 4}:
sigma_or_weight = float(reader.read_nonempty_line().strip())
weight = calc_weight(weight_mode, transformed, sigma_or_weight)
# Construct design row
row = lattice_design_row(lattice_system, h, k, l)
row = np.concatenate([row, np.array(correction_terms(correction_flags, transformed), dtype=float)])
# The y-value for least squares, typically proportional to 1/d^2
# Here it's 4.0 / (lambda^2) * sin^2(theta) = 1/d^2, if transformed is sin(theta)
# The original code comment says "transformed_position (for least square)" and is (4.0/lambda^2)*sin^2(theta)
# If transformed is already sin(theta), then this is 1/d^2.
y_value = (4.0 / (wavelength * wavelength)) * transformed * transformed
observations.append(
Observation(
serial_index=serial_index,
h=h,
k=k,
l=l,
raw_position=raw_position,
transformed_position=y_value,
weight=weight,
design_row=row,
wavelength=wavelength,
)
)
line = reader.read_nonempty_line()
jobs.append(JobData(title=title, config=config, observations=observations, reference_wavelength=reference_wavelength))
# Read job termination code
job_code = int(reader.read_nonempty_line().strip())
if job_code == 0:
break # All jobs finished
return jobs
# ---------- least squares ----------
[ドキュメント]
def weighted_least_squares(observations: list[Observation], output: TextIO) -> LeastSquaresResult:
"""
観測データに対して重み付き最小二乗法を実行する。
詳細説明:
この関数は、与えられた観測データに対し、反復的に重み付き最小二乗法を適用します。
各反復では、現在のパラメータで残差を計算し、`OUTLIER_SIGMA` を超える残差を持つ観測データを
外れ値として除去します。外れ値がなくなるか、自由度が0になるまでこのプロセスを繰り返します。
最終的に、最適なパラメータ、その標準偏差、適合値、残差、および共分散行列が返されます。
:param observations: list[Observation]: 観測データのリスト。
:param output: TextIO: 外れ値の除去情報を書き込むための出力ストリーム。
:returns: LeastSquaresResult: 最小二乗計算の結果を含むオブジェクト。
:raises LatticeCalculationError: 観測データが提供されない場合、観測数よりもパラメータ数が多い場合、または正規方程式が特異な場合。
"""
if not observations:
raise LatticeCalculationError("no observations were provided")
current = observations[:]
n_params = current[0].design_row.size
while True:
y = np.array([obs.transformed_position for obs in current], dtype=float)
w = np.array([obs.weight for obs in current], dtype=float)
x = np.vstack([obs.design_row for obs in current])
if len(current) < n_params:
raise LatticeCalculationError("number of observations is smaller than number of parameters")
# Construct normal equations: (X^T W X) beta = X^T W y
# X^T W X is 'normal'
# X^T W y is 'rhs'
normal = x.T @ (w[:, None] * x)
rhs = x.T @ (w * y)
try:
covariance_unscaled = np.linalg.inv(normal)
except np.linalg.LinAlgError as exc:
raise LatticeCalculationError("normal matrix is singular") from exc
params = covariance_unscaled @ rhs
fitted = x @ params
residuals = y - fitted
dof = len(current) - n_params # Degrees of freedom
if dof == 0:
# If no degrees of freedom, cannot calculate scale factor for sigma
sigma_obs = np.zeros(len(current), dtype=float)
sigma_params = np.zeros(n_params, dtype=float)
covariance = covariance_unscaled.copy() # No scaling
return LeastSquaresResult(params, sigma_params, fitted, residuals, sigma_obs, current, covariance)
weighted_ss = float(np.sum(w * residuals * residuals)) # Weighted sum of squares of residuals
scale = math.sqrt(weighted_ss / dof) # Scale factor for standard deviation
# Calculate observation sigma: sigma_obs = scale / sqrt(weight)
sigma_obs = np.where(w > 0.0, scale / np.sqrt(w), 0.0)
# Identify outliers
outlier_mask = np.abs(residuals) > OUTLIER_SIGMA * sigma_obs
if np.any(outlier_mask):
# Outliers detected, remove them and re-run least squares
removed = [obs for obs, flag in zip(current, outlier_mask) if flag]
for obs, pred in zip(current, fitted):
# The original code has an empty loop here, suggesting debug output
# or a placeholder for detailed logging of all observations.
pass
for obs in removed:
# Log removed observations
pred_d = 1.0 / math.sqrt(max(obs.transformed_position, ANGLE_EPS)) # d_obs from original transformed position
output.write(
"This observation is eliminated.\n"
f"No. = {obs.serial_index:2d} hkl=({obs.h:d} {obs.k:d} {obs.l:d}) raw={obs.raw_position:10.3f}\n"
)
current = [obs for obs, flag in zip(current, outlier_mask) if not flag]
continue # Repeat the loop with cleaned observations
# No outliers, final result
covariance = covariance_unscaled * (scale * scale) # Scale covariance matrix
sigma_params = np.sqrt(np.maximum(np.diag(covariance), 0.0)) # Parameter standard deviations
return LeastSquaresResult(params, sigma_params, fitted, residuals, sigma_obs, current, covariance)
# ---------- post-processing ----------
[ドキュメント]
def derive_cell_constants(lattice_system: int, params: np.ndarray, sigma: np.ndarray) -> CellConstants:
"""
最小二乗法で得られたパラメータから逆格子および実格子の定数を導出する。
詳細説明:
最小二乗法によって決定された格子パラメータとそれらの標準偏差を受け取り、
指定された格子系の結晶学的な関係に基づいて、最終的な逆格子定数
(a*, b*, c*, alpha*, beta*, gamma*, V*) および実格子定数
(a, b, c, alpha, beta, gamma, V) を計算します。
各定数に対して誤差伝播法を用いて標準偏差も同時に計算されます。
:param lattice_system: int: 格子系の種類を示す整数コード（1-8）。
:param params: np.ndarray: 最小二乗法で得られたパラメータの配列。
:param sigma: np.ndarray: パラメータの標準偏差の配列。
:returns: CellConstants: 導出された格子定数を含むオブジェクト。
:raises LatticeCalculationError: 未対応の格子系が指定された場合。
"""
cell = CellConstants()
def reciprocal_to_direct_general(rec_lengths: np.ndarray, rec_length_sigma: np.ndarray, rec_cos: np.ndarray, rec_cos_sigma: np.ndarray) -> None:
"""
一般的な逆格子定数から実格子定数を導出し、CellConstantsオブジェクトに格納するヘルパー関数。
"""
cell.reciprocal_lengths[:] = rec_lengths
cell.reciprocal_length_sigma[:] = rec_length_sigma
cell.reciprocal_cosines[:] = rec_cos
cell.reciprocal_cosine_sigma[:] = rec_cos_sigma
for i in range(3):
# Calculate reciprocal angles and their sigmas
cell.reciprocal_angles_deg[i] = acos_deg(rec_cos[i])
cell.reciprocal_angle_sigma_deg[i] = SigmaPropagator([rec_cos[i]], [rec_cos_sigma[i]]).evaluate(lambda v: acos_deg(v[0]))[1]
# Propagate error for direct volume
prop = SigmaPropagator(
[*rec_lengths, *rec_cos],
[*rec_length_sigma, *rec_cos_sigma],
)
cell.direct_volume, cell.direct_volume_sigma = prop.evaluate(direct_volume_from_reciprocal)
# Calculate reciprocal volume (inverse of direct volume)
if cell.direct_volume > 0.0:
cell.reciprocal_volume = 1.0 / cell.direct_volume
# Error propagation for 1/x: d(1/x) = -1/x^2 dx => sigma(1/x) = dx/x^2
cell.reciprocal_volume_sigma = cell.direct_volume_sigma / (cell.direct_volume * cell.direct_volume)
for i in range(3):
j = (i + 1) % 3
k = (i + 2) % 3
# Calculate direct cosines and their sigmas
cell.direct_cosines[i], cell.direct_cosine_sigma[i] = SigmaPropagator(
[rec_cos[i], rec_cos[j], rec_cos[k]],
[rec_cos_sigma[i], rec_cos_sigma[j], rec_cos_sigma[k]],
).evaluate(direct_cosine_from_reciprocal)
# Calculate direct angles and their sigmas
cell.direct_angles_deg[i] = acos_deg(cell.direct_cosines[i])
cell.direct_angle_sigma_deg[i] = SigmaPropagator(
[cell.direct_cosines[i]], [cell.direct_cosine_sigma[i]]
).evaluate(lambda v: acos_deg(v[0]))[1]
# Calculate direct lengths and their sigmas
cell.direct_lengths[i], cell.direct_length_sigma[i] = SigmaPropagator(
[rec_lengths[i], rec_cos[i], rec_cos[j], rec_cos[k]],
[rec_length_sigma[i], rec_cos_sigma[i], rec_cos_sigma[j], rec_cos_sigma[k]],
).evaluate(direct_length_from_reciprocal)
if lattice_system == 1: # Triclinic
# Params: a*^2, b*^2, c*^2, 2b*c*cos(alpha*), 2c*a*cos(beta*), 2a*b*cos(gamma*)
rec_lengths = np.sqrt(np.maximum(params[:3], 0.0))
# Error propagation for sqrt(x): d(sqrt(x)) = 0.5 * dx / sqrt(x)
rec_length_sigma = np.divide(sigma[:3], 2.0 * np.maximum(rec_lengths, ANGLE_EPS), out=np.zeros(3), where=rec_lengths > ANGLE_EPS)
rec_cos = np.array([
params[3] / (rec_lengths[1] * rec_lengths[2]) if rec_lengths[1] * rec_lengths[2] > ANGLE_EPS else 0.0,
params[4] / (rec_lengths[2] * rec_lengths[0]) if rec_lengths[2] * rec_lengths[0] > ANGLE_EPS else 0.0,
params[5] / (rec_lengths[0] * rec_lengths[1]) if rec_lengths[0] * rec_lengths[1] > ANGLE_EPS else 0.0,
])
rec_cos_sigma = np.zeros(3, dtype=float)
for i in range(3):
j = (i + 1) % 3
k = (i + 2) % 3
# Error propagation for x/(y*z)
rec_cos_sigma[i] = SigmaPropagator(
[params[3 + i], rec_lengths[j], rec_lengths[k]],
[sigma[3 + i], rec_length_sigma[j], rec_length_sigma[k]],
).evaluate(lambda v: 0.0 if abs(v[1] * v[2]) <= ANGLE_EPS else v[0] / (v[1] * v[2]))[1]
reciprocal_to_direct_general(rec_lengths, rec_length_sigma, rec_cos, rec_cos_sigma)
return cell
if lattice_system == 2: # Monoclinic (unique b axis, beta*)
# Parameters: a*^2, b*^2, c*^2, 2c*a*cos(beta*)
# Missing 2b*c*cos(alpha*) (0), 2a*b*cos(gamma*) (0)
expanded_p = np.array([params[0], params[1], params[2], 0.0, params[3], 0.0], dtype=float)
expanded_s = np.array([sigma[0], sigma[1], sigma[2], 0.0, sigma[3], 0.0], dtype=float)
return derive_cell_constants(1, expanded_p, expanded_s)
if lattice_system == 3: # Monoclinic (unique c axis, gamma*) (as implemented)
# Parameters: a*^2, b*^2, c*^2, 2a*b*cos(gamma*)
# Missing 2b*c*cos(alpha*) (0), 2c*a*cos(beta*) (0)
expanded_p = np.array([params[0], params[1], params[2], 0.0, 0.0, params[3]], dtype=float)
expanded_s = np.array([sigma[0], sigma[1], sigma[2], 0.0, 0.0, sigma[3]], dtype=float)
return derive_cell_constants(1, expanded_p, expanded_s)
if lattice_system == 4: # Orthorhombic
# Parameters: a*^2, b*^2, c*^2 (all angles are 90 deg, cosines are 0)
rec_lengths = np.sqrt(np.maximum(params[:3], 0.0))
rec_length_sigma = np.divide(sigma[:3], 2.0 * np.maximum(rec_lengths, ANGLE_EPS), out=np.zeros(3), where=rec_lengths > ANGLE_EPS)
reciprocal_to_direct_general(rec_lengths, rec_length_sigma, np.zeros(3), np.zeros(3)) # Reciprocal cosines and their sigmas are 0
return cell
if lattice_system == 5: # Tetragonal
# Parameters: a*^2, c*^2 (a*=b*, all angles 90 deg)
a_star = safe_sqrt(params[0])
c_star = safe_sqrt(params[1])
da_star = 0.0 if a_star <= ANGLE_EPS else 0.5 * sigma[0] / a_star
dc_star = 0.0 if c_star <= ANGLE_EPS else 0.5 * sigma[1] / c_star
cell.reciprocal_lengths[:] = [a_star, a_star, c_star]
cell.reciprocal_length_sigma[:] = [da_star, da_star, dc_star]
# Direct lengths a = 1/a*, c = 1/c*
a = 0.0 if a_star <= ANGLE_EPS else 1.0 / a_star
c = 0.0 if c_star <= ANGLE_EPS else 1.0 / c_star
cell.direct_lengths[:] = [a, a, c]
cell.direct_length_sigma[:] = [da_star / max(a_star * a_star, ANGLE_EPS), da_star / max(a_star * a_star, ANGLE_EPS), dc_star / max(c_star * c_star, ANGLE_EPS)]
cell.direct_volume = a * a * c # For tetragonal, direct_volume = a^2 * c
cell.direct_volume_sigma = SigmaPropagator([a, a, c], [da_star / max(a_star * a_star, ANGLE_EPS), da_star / max(a_star * a_star, ANGLE_EPS), dc_star / max(c_star * c_star, ANGLE_EPS)]).evaluate(lambda v: v[0] * v[1] * v[2])[1] # Propagate error for a*a*c
if cell.direct_volume > 0.0:
cell.reciprocal_volume = 1.0 / cell.direct_volume
cell.reciprocal_volume_sigma = cell.direct_volume_sigma / (cell.direct_volume * cell.direct_volume)
return cell
if lattice_system == 6: # Cubic
# Parameters: a*^2 (a*=b*=c*, all angles 90 deg)
a_star = safe_sqrt(params[0])
da_star = 0.0 if a_star <= ANGLE_EPS else 0.5 * sigma[0] / a_star
cell.reciprocal_lengths[:] = [a_star, a_star, a_star]
cell.reciprocal_length_sigma[:] = [da_star, da_star, da_star]
# Direct length a = 1/a*
a = 0.0 if a_star <= ANGLE_EPS else 1.0 / a_star
da = da_star / max(a_star * a_star, ANGLE_EPS)
cell.direct_lengths[:] = [a, a, a]
cell.direct_length_sigma[:] = [da, da, da]
cell.direct_volume = a**3
cell.direct_volume_sigma = 3.0 * a * a * da # For a^3, sigma is 3*a^2*da
if cell.direct_volume > 0.0:
cell.reciprocal_volume = 1.0 / cell.direct_volume
cell.reciprocal_volume_sigma = cell.direct_volume_sigma / (cell.direct_volume * cell.direct_volume)
return cell
if lattice_system == 7: # Rhombohedral (trigonal)
# Parameters: (a*^2 + b*^2 + c*^2), 2(b*c*cos(alpha*) + c*a*cos(beta*) + a*b*cos(gamma*))
# With a*=b*=c* and alpha*=beta*=gamma*, this simplifies:
# P0 = 3*a*^2
# P1 = 6*a*^2*cos(alpha*)
# So a*^2 = P0/3 => a* = sqrt(P0/3)
# cos(alpha*) = P1/P0
# a_star = safe_sqrt(params[0] / 3.0) # params[0] is 3*a*^2 for rhombohedral
a_star = safe_sqrt(params[0])
da_star = SigmaPropagator([params[0]], [sigma[0]]).evaluate(lambda v: safe_sqrt(v[0] / 3.0))[1]
cos_alpha_star = 0.0 if abs(params[0]) <= ANGLE_EPS else params[1] / params[0]
dcos_alpha_star = SigmaPropagator([params[0], params[1]], [sigma[0], sigma[1]]).evaluate(
lambda v: 0.0 if abs(v[0]) <= ANGLE_EPS else v[1] / v[0]
)[1]
cell.reciprocal_lengths[:] = [a_star, a_star, a_star]
cell.reciprocal_length_sigma[:] = [da_star, da_star, da_star]
cell.reciprocal_cosines[:] = [cos_alpha_star, cos_alpha_star, cos_alpha_star]
cell.reciprocal_cosine_sigma[:] = [dcos_alpha_star, dcos_alpha_star, dcos_alpha_star]
angle_star = acos_deg(cos_alpha_star)
dangle_star = SigmaPropagator([cos_alpha_star], [dcos_alpha_star]).evaluate(lambda v: acos_deg(v[0]))[1]
cell.reciprocal_angles_deg[:] = [angle_star, angle_star, angle_star]
cell.reciprocal_angle_sigma_deg[:] = [dangle_star, dangle_star, dangle_star]
# Direct lengths, cosines, angles, and volume for trigonal
cell.direct_lengths[0], cell.direct_length_sigma[0] = SigmaPropagator([a_star, cos_alpha_star], [da_star, dcos_alpha_star]).evaluate(trigonal_direct_length)
cell.direct_lengths[:] = [cell.direct_lengths[0]] * 3
cell.direct_length_sigma[:] = [cell.direct_length_sigma[0]] * 3
cell.direct_cosines[0], cell.direct_cosine_sigma[0] = SigmaPropagator([cos_alpha_star], [dcos_alpha_star]).evaluate(trigonal_direct_cosine)
cell.direct_cosines[:] = [cell.direct_cosines[0]] * 3
cell.direct_cosine_sigma[:] = [cell.direct_cosine_sigma[0]] * 3
angle = acos_deg(cell.direct_cosines[0])
dangle = SigmaPropagator([cell.direct_cosines[0]], [cell.direct_cosine_sigma[0]]).evaluate(lambda v: acos_deg(v[0]))[1]
cell.direct_angles_deg[:] = [angle, angle, angle]
cell.direct_angle_sigma_deg[:] = [dangle, dangle, dangle]
cell.direct_volume, cell.direct_volume_sigma = SigmaPropagator([a_star, cos_alpha_star], [da_star, dcos_alpha_star]).evaluate(trigonal_direct_volume)
if cell.direct_volume > 0.0:
cell.reciprocal_volume = 1.0 / cell.direct_volume
cell.reciprocal_volume_sigma = cell.direct_volume_sigma / (cell.direct_volume * cell.direct_volume)
return cell
if lattice_system == 8: # Hexagonal
# Parameters: (a*^2 + b*^2 + a*b*), c*^2
# Which is 3*a*^2 for hexagonal a* = b* and gamma*=120deg.
# So P0 = a*^2 + a*^2 + a*^2*(-0.5)*2 = 3*a*^2
# No, a*^2+k*^2+hk gives a_star^2 (from base_design_terms for h^2+k^2+2hk/2)
# So params[0] is a*^2, params[1] is c*^2
a_star = safe_sqrt(params[0])
c_star = safe_sqrt(params[1])
da_star = 0.0 if a_star <= ANGLE_EPS else 0.5 * sigma[0] / a_star
dc_star = 0.0 if c_star <= ANGLE_EPS else 0.5 * sigma[1] / c_star
cell.reciprocal_lengths[:] = [a_star, a_star, c_star]
cell.reciprocal_length_sigma[:] = [da_star, da_star, dc_star]
# Direct lengths for hexagonal: a = 2/(sqrt(3)a*), c=1/c*
a = 0.0 if a_star <= ANGLE_EPS else 2.0 / (math.sqrt(3.0) * a_star)
c = 0.0 if c_star <= ANGLE_EPS else 1.0 / c_star
# Error propagation for 1/x: d(1/x) = dx/x^2
da = 0.0 if a_star <= ANGLE_EPS else da_star * a / a_star # (2/sqrt(3)) * (-1/a*^2) * da* = (2/sqrt(3)a*) * (-1/a*) * da* = a * (-1/a*) * da*
dc = 0.0 if c_star <= ANGLE_EPS else dc_star * c / c_star # c * (-1/c*) * dc*
cell.direct_lengths[:] = [a, a, c]
cell.direct_length_sigma[:] = [da, da, dc]
# Angles for hexagonal
cell.direct_cosines[2] = -0.5 # cos(120 deg)
cell.reciprocal_cosines[2] = 0.5 # cos(60 deg)
cell.direct_angles_deg[2] = 120.0
cell.reciprocal_angles_deg[2] = 60.0
# Direct volume for hexagonal V = sqrt(3)/2 * a^2 * c
cell.direct_volume = math.sqrt(3.0) * 0.5 * a * a * c
cell.direct_volume_sigma = SigmaPropagator([a, a, c], [da, da, dc]).evaluate(lambda v: math.sqrt(3.0) * 0.5 * v[0] * v[1] * v[2])[1]
if cell.direct_volume > 0.0:
cell.reciprocal_volume = 1.0 / cell.direct_volume
cell.reciprocal_volume_sigma = cell.direct_volume_sigma / (cell.direct_volume * cell.direct_volume)
return cell
raise LatticeCalculationError(f"unsupported lattice system: {lattice_system}")
# ---------- reporting ----------
[ドキュメント]
def write_job_report(output: TextIO, job: JobData, fit: LeastSquaresResult, cell: CellConstants) -> None:
"""
計算結果をフォーマットされたレポートとして出力ストリームに書き込む。
詳細説明:
この関数は、単一の計算ジョブに対する詳細なレポートを生成します。
レポートには、ジョブのタイトル、設定、最小二乗法で得られたパラメータと誤差、
各観測データに対する計算値と残差、そして最終的に導出された逆格子および実格子の
定数とその標準偏差が含まれます。
:param output: TextIO: レポートを書き込むための出力ストリーム。
:param job: JobData: 計算に使用されたジョブデータ。
:param fit: LeastSquaresResult: 最小二乗法の結果。
:param cell: CellConstants: 導出された格子定数。
:returns: None
"""
output.write(f"\n {job.title}\n")
cfg = job.config
output.write(" HLS LE(1) LE(2) LE(3) LE(4) LE(5) IW IO IP\n")
output.write(
f"{cfg.lattice_system:7d}{cfg.correction_flags[0]:7d}{cfg.correction_flags[1]:7d}{cfg.correction_flags[2]:7d}"
f"{cfg.correction_flags[3]:7d}{cfg.correction_flags[4]:7d}{cfg.weight_mode:7d}{cfg.position_mode:7d}{cfg.profile_mode:7d}\n"
)
output.write(f"Wave Length ={job.reference_wavelength:f}\n")
output.write(" P(I) SIGP(I)\n")
for value, sigma in zip(fit.parameters, fit.parameter_sigma):
output.write(f" {value:10.7f} {sigma:10.7f}\n")
output.write(" H K L Dobs. Dcal. BO BC Resid. d(2Q)\n")
for obs, fitted, resid in zip(fit.kept_observations, fit.fitted_y, fit.residuals):
# Calculate d_obs and d_cal from transformed positions
# transformed_position = 4 * sin^2(theta) / wavelength^2 = 1/d^2
# So d = 1 / sqrt(transformed_position)
d_obs = 1.0 / math.sqrt(max(obs.transformed_position, ANGLE_EPS))
d_cal = 1.0 / math.sqrt(max(fitted, ANGLE_EPS))
# Calculate calculated 2theta from d_cal
# 2d sin(theta) = wavelength => sin(theta) = wavelength / (2d)
theta_calc = clamp_unit(obs.wavelength / (2.0 * d_cal)) if d_cal > 0 else 1.0
two_theta_calc = 2.0 * math.degrees(math.asin(theta_calc))
# Print observation details, fitted values, residuals, and 2theta difference
output.write(
f" ({obs.serial_index:2d}){obs.h:5d}{obs.k:5d}{obs.l:5d}{d_obs:8.4f}{d_cal:8.4f}"
f"{obs.raw_position:8.3f}{two_theta_calc:8.3f}{resid:10.5f}{(obs.raw_position - two_theta_calc):10.5f}\n"
)
output.write("Reciprocal cell constant\n")
output.write(" a* b* c*\n")
for j in range(3):
output.write(f"{cell.reciprocal_lengths[j]:10.7f}({cell.reciprocal_length_sigma[j]:10.7f}) ")
output.write("\n alpha* beta* gamma*\n")
for j in range(3):
output.write(f"{cell.reciprocal_angles_deg[j]:10.6f}({cell.reciprocal_angle_sigma_deg[j]:10.6f}) ")
output.write("\n cos a* cos b* cos c*\n")
for j in range(3):
output.write(f"{cell.reciprocal_cosines[j]:10.6f}({cell.reciprocal_cosine_sigma[j]:10.6f}) ")
output.write(f"\n V*={cell.reciprocal_volume:.10g}({cell.reciprocal_volume_sigma:.10g})\n")
output.write("\nDirect cell constant\n")
output.write(" a b c\n")
for j in range(3):
output.write(f"{cell.direct_lengths[j]:8.5f}({cell.direct_length_sigma[j]:8.5f}) ")
output.write("\n alpha beta gamma\n")
for j in range(3):
output.write(f"{cell.direct_angles_deg[j]:8.4f}({cell.direct_angle_sigma_deg[j]:8.4f}) ")
output.write("\n cos a cos b cos c\n")
for j in range(3):
output.write(f"{cell.direct_cosines[j]:8.4f}({cell.direct_cosine_sigma[j]:8.4f}) ")
output.write(f"\n V={cell.direct_volume:.10g}({cell.direct_volume_sigma:.10g})\n")
# ---------- CLI ----------
[ドキュメント]
def build_argument_parser() -> argparse.ArgumentParser:
"""
コマンドライン引数を解析するためのパーサーを構築する。
詳細説明:
この関数は、スクリプトが受け入れるコマンドライン引数を定義します。
具体的には、入力ファイルパス、出力ファイルパス、およびバナー出力を抑制するための
`--silent` フラグを設定します。
:returns: argparse.ArgumentParser: 構築されたArgumentParserオブジェクト。
"""
parser = argparse.ArgumentParser(description="Least-square calculation of lattice constants from reflection data")
parser.add_argument("input_file", type=Path, help="input file")
parser.add_argument("output_file", type=Path, help="output report file")
parser.add_argument("--silent", action="store_true", help="suppress banner output")
return parser
[ドキュメント]
def main() -> int:
"""
スクリプトのエントリポイント。
詳細説明:
コマンドライン引数を解析し、入力ファイルから格子定数計算ジョブを読み込みます。
各ジョブに対して重み付き最小二乗計算を実行し、得られた結果から格子定数を導出します。
最終的なレポートは指定された出力ファイルに書き込まれます。
:returns: int: 終了コード（成功の場合は0）。
"""
parser = build_argument_parser()
args = parser.parse_args()
if not args.silent:
print("Least Square Calculation for Lattice Constant")
print("Refactored Python version")
print(f"Input : {args.input_file}")
print(f"Output: {args.output_file}")
with args.input_file.open("r", encoding="utf-8") as fin:
jobs = parse_jobs(fin)
with args.output_file.open("w", encoding="utf-8") as fout:
fout.write("Least Square Calculation for Lattice Constant\n")
fout.write("Refactored Python version\n\n")
for job in jobs:
fit = weighted_least_squares(job.observations, fout)
cell = derive_cell_constants(job.config.lattice_system, fit.parameters, fit.parameter_sigma)
write_job_report(fout, job, fit, cell)
fout.write("\nLattice Constant Calculation was finished!\n")
if not args.silent:
print("Lattice Constant Calculation was finished!")
return 0
if __name__ == "__main__":
main()