Coverage for graphica/core/analysis.py: 95.8%

1from typing import Any, Callable

2import numpy as np

3from scipy import sparse

4from scipy.sparse.linalg import spsolve

5from scipy.integrate import simpson, cumulative_trapezoid, cumulative_simpson

6from scipy.interpolate import CubicSpline

7from scipy.ndimage import uniform_filter1d, median_filter, gaussian_filter1d

8from scipy.optimize import curve_fit

9from scipy.signal import find_peaks, peak_widths, savgol_filter, correlate, correlation_lags

10from scipy.special import wofz

11from scipy.stats import gaussian_kde

13from graphica.core.fit_models import RESERVED_FIT_TYPE_NAMES, fit_model_id, resolve_fit_model, resolve_fit_param_names

15CURVE_FIT_MAX_ITERATIONS = 5000

17# プラグインのフィット関数。{name: {"func": f(x, *params), "params": [...], "p0": list | callable | None}}

18_PLUGIN_FIT_FUNCTIONS: dict[str, dict[str, Any]] = {}

21def register_fit_function(name: str, func: Callable[..., Any], param_names: list[str],

22 p0: list[float] | Callable[..., Any] | None = None) -> None:

23 """p0 は初期値のリストか (x_data, y_data) -> list を返す関数。省略時は全て 1.0。"""

24 if not name or not name.strip(): 24 ↛ 25line 24 didn't jump to line 25 because the condition on line 24 was never true

25 raise ValueError("フィット関数名が空です。")

26 if name in RESERVED_FIT_TYPE_NAMES:

27 raise ValueError(f"'{name}' は組み込みのフィットタイプ名と衝突します。")

28 if name in _PLUGIN_FIT_FUNCTIONS:

29 raise ValueError(f"フィット関数 '{name}' は既に登録されています。")

30 if not param_names:

31 raise ValueError("param_names が空です。")

32 _PLUGIN_FIT_FUNCTIONS[name] = {"func": func, "params": list(param_names), "p0": p0}

35def get_plugin_fit_type_names() -> list[str]:

36 return list(_PLUGIN_FIT_FUNCTIONS.keys())

39def get_fit_param_names(fit_type: str, custom_formula: str | None = None, model_id: str | None = None) -> list[str]:

40 """フィットせずにパラメータ名を返す(フィットの前に入力欄を組み立てるため)。model_id が分かればそれを優先する。"""

41 return resolve_fit_param_names(fit_type, custom_formula, _PLUGIN_FIT_FUNCTIONS, model_id)

44def get_fit_model_id(fit_type: str) -> str | None:

45 """保存用の安定した ID(組み込みは 'gaussian' など、カスタム数式、'plugin:<名前>')。分からなければ None。"""

46 return fit_model_id(fit_type, _PLUGIN_FIT_FUNCTIONS)

49_ROBUST_LOSS_FUNCTIONS = ('linear', 'soft_l1', 'huber')

52def _run_curve_fit_with_overrides(fit_func: Callable[..., Any], params_info: list[str], p0: Any, x_data: Any, y_data: Any,

53 sigma: Any, p0_overrides: dict[str, float] | None, fixed_params: dict[str, float] | None,

54 bounds: dict[str, tuple[float, float]] | None, fit_type_label: str,

55 loss: str = 'linear') -> tuple[np.ndarray, np.ndarray]:

56 """p0 の上書き・固定・範囲拘束を適用して curve_fit を実行する。

58 popt と pcov は固定パラメータも含めたフルサイズで返す(固定分の pcov は 0)。

59 """

60 if loss not in _ROBUST_LOSS_FUNCTIONS:

61 raise ValueError(f"未知の損失関数です(loss): '{loss}' (使用可能: {_ROBUST_LOSS_FUNCTIONS})")

62 if len(x_data) < len(params_info):

63 raise ValueError(

64 f"データ点数 ({len(x_data)}) がフィットに必要なパラメータ数 "

65 f"({len(params_info)}) より少ないため、フィッティングできません。"

66 )

68 p0_overrides = p0_overrides or {}

69 fixed_params = fixed_params or {}

70 bounds = bounds or {}

72 for name_dict, label in (

73 (p0_overrides, "p0_overrides"), (fixed_params, "fixed_params"), (bounds, "bounds"),

74 ):

75 for pname in name_dict:

76 if pname not in params_info:

77 raise ValueError(

78 f"未知のパラメータ名です({label}): '{pname}' "

79 f"(このフィットタイプのパラメータ: {params_info})"

80 )

82 if fixed_params and len(fixed_params) >= len(params_info):

83 raise ValueError(

84 "すべてのパラメータを固定することはできません"

85 "(最適化する自由パラメータが1つも残りません)。"

86 )

88 p0 = list(p0)

89 for i, name in enumerate(params_info):

90 if name in p0_overrides:

91 p0[i] = float(p0_overrides[name])

93 free_indices = [i for i, name in enumerate(params_info) if name not in fixed_params]

94 fixed_indices = [i for i, name in enumerate(params_info) if name in fixed_params]

95 fixed_values = {i: float(fixed_params[params_info[i]]) for i in fixed_indices}

97 if fixed_indices:

98 # curve_fit には自由パラメータだけを渡し、固定値は関数の中で元の位置に挿し込む

99 original_fit_func = fit_func

100

101 def fit_func_for_curve_fit(x: Any, *free_args: float) -> Any:

102 full_params: list[float | None] = [None] * len(params_info)

103 for i in fixed_indices:

104 full_params[i] = fixed_values[i]

105 for idx, i in enumerate(free_indices):

106 full_params[i] = free_args[idx]

107 return original_fit_func(x, *full_params)

108

109 p0_for_curve_fit = [p0[i] for i in free_indices]

110 else:

111 fit_func_for_curve_fit = fit_func

112 p0_for_curve_fit = p0

113

114 curve_fit_kwargs = {

115 "p0": p0_for_curve_fit,

116 "sigma": sigma,

117 "absolute_sigma": sigma is not None,

118 }

119 if bounds:

120 lower, upper = [], []

121 for i in free_indices:

122 lo, hi = bounds.get(params_info[i], (-np.inf, np.inf))

123 lower.append(lo)

124 upper.append(hi)

125 # bounds は p0 が境界の厳密に内側にあることを要求する(等しいだけで例外)。

126 # 「初期値=下限」のような自然な入力で落ちないよう、わずかに内側へずらす。

127 for idx in range(len(p0_for_curve_fit)):

128 lo, hi = lower[idx], upper[idx]

129 val = p0_for_curve_fit[idx]

130 if val <= lo or val >= hi:

131 span = hi - lo

132 nudge = span * 1e-6 if np.isfinite(span) and span > 0 else max(abs(val), 1.0) * 1e-6 or 1e-9

133 p0_for_curve_fit[idx] = min(max(val, lo + nudge), hi - nudge)

134 curve_fit_kwargs["bounds"] = (lower, upper)

135 # bounds 付きは trf 法になり、maxfev ではなく max_nfev しか受け付けない

136 curve_fit_kwargs["max_nfev"] = CURVE_FIT_MAX_ITERATIONS

137 else:

138 curve_fit_kwargs["maxfev"] = CURVE_FIT_MAX_ITERATIONS

139

140 if loss != 'linear':

141 # 既定の 'lm' は loss を受け付けないので trf にし、maxfev も max_nfev に付け替える

142 curve_fit_kwargs["loss"] = loss

143 curve_fit_kwargs.setdefault("method", "trf")

144 if "maxfev" in curve_fit_kwargs:

145 curve_fit_kwargs["max_nfev"] = curve_fit_kwargs.pop("maxfev")

146

147 try:

148 popt_free, pcov_free = curve_fit(fit_func_for_curve_fit, x_data, y_data, **curve_fit_kwargs)

149 except RuntimeError as e:

150 raise RuntimeError(

151 f"フィッティングが収束しませんでした（{fit_type_label}）。データの分布が"

152 f"このモデルに適していない可能性があります。詳細: {e}"

153 ) from e

154

155 if fixed_indices:

156 popt = np.empty(len(params_info))

157 for i in fixed_indices:

158 popt[i] = fixed_values[i]

159 for idx, i in enumerate(free_indices):

160 popt[i] = popt_free[idx]

161

162 pcov = np.zeros((len(params_info), len(params_info)))

163 for row_idx, i in enumerate(free_indices):

164 for col_idx, j in enumerate(free_indices):

165 pcov[i, j] = pcov_free[row_idx, col_idx]

166 else:

167 popt, pcov = popt_free, pcov_free

168

169 return popt, pcov

170

171

172def evaluate_saved_fit(fit_result: dict[str, Any], x: Any) -> Any:

173 """保存したフィット結果(Dataset.fit_result)の曲線の値を x で計算する。式が分からない・計算できないときは None。

174

175 Args:

176 fit_result: calculate_curve_fit の結果から作った、データセットに保存する形の辞書。

177 x: 値を求める X(数値の並び)。

178

179 Returns:

180 x と同じ長さの float の ndarray、または None。

181 """

182 try:

183 model = resolve_fit_model(fit_result.get('fit_type', ''), fit_result.get('custom_formula'),

184 _PLUGIN_FIT_FUNCTIONS, fit_result.get('fit_model_id'))

185 x_values = np.asarray(x, dtype=float)

186 values = np.asarray(model.func(x_values, *fit_result.get('params', [])), dtype=float)

187 return np.broadcast_to(values, x_values.shape).copy()

188 except (ValueError, TypeError, ArithmeticError, KeyError, IndexError):

189 return None

190

191

192def calculate_curve_fit(x_data: Any, y_data: Any, fit_type: str, custom_formula: str | None = None, sigma: Any = None,

193 x_range: tuple[float, float] | None = None, p0_overrides: dict[str, float] | None = None,

194 fixed_params: dict[str, float] | None = None,

195 bounds: dict[str, tuple[float, float]] | None = None, loss: str = 'linear',

196 model_id: str | None = None) -> dict[str, Any]:

197 """曲線フィット。popt / pcov / perr / x_fit・y_fit / r_squared / residuals などの dict を返す。

198

199 x_range は両端を含み、範囲外の点は p0 の推定にも使わない。sigma は absolute_sigma=True で渡す。

200 """

201 x_data = np.asarray(x_data)

202 y_data = np.asarray(y_data)

203 if sigma is not None:

204 sigma = np.asarray(sigma)

205 if x_range is not None:

206 x_min, x_max = x_range

207 range_mask = (x_data >= x_min) & (x_data <= x_max)

208 x_data, y_data = x_data[range_mask], y_data[range_mask]

209 if sigma is not None:

210 sigma = sigma[range_mask]

211

212 # NaN が残ると curve_fit は収束失敗ではない ValueError を出し、p0 の推定も落ちる

213 nan_mask = np.isnan(x_data) | np.isnan(y_data)

214 if sigma is not None:

215 nan_mask |= np.isnan(sigma)

216 if nan_mask.any():

217 x_data, y_data = x_data[~nan_mask], y_data[~nan_mask]

218 if sigma is not None: 218 ↛ 219line 218 didn't jump to line 219 because the condition on line 218 was never true

219 sigma = sigma[~nan_mask]

220

221 if len(x_data) == 0:

222 raise ValueError("有効なデータ点がありません(すべて欠損値です)。フィッティングできません。")

223

224 model = resolve_fit_model(fit_type, custom_formula, _PLUGIN_FIT_FUNCTIONS, model_id)

225 if model.check_data is not None:

226 model.check_data(x_data)

227 fit_func, params_info = model.func, model.param_names

228 p0 = model.initial_guess(x_data, y_data)

229

230 popt, pcov = _run_curve_fit_with_overrides(

231 fit_func, params_info, p0, x_data, y_data, sigma,

232 p0_overrides, fixed_params, bounds, fit_type, loss=loss,

233 )

234

235 # 退化したフィットでは pcov の対角が負や inf になる。例外にせず nan / inf のまま返す

236 perr = np.sqrt(np.diag(pcov))

237

238 x_fit = np.linspace(x_data.min(), x_data.max(), 200)

239 y_fit = fit_func(x_fit, *popt)

240

241 residuals = y_data - fit_func(x_data, *popt)

242 ss_res = np.sum(residuals ** 2)

243 ss_tot = np.sum((y_data - np.mean(y_data)) ** 2)

244 r_squared = 1.0 if ss_tot == 0 else 1.0 - (ss_res / ss_tot)

245

246 return {

247 'popt': popt,

248 'pcov': pcov,

249 'perr': perr,

250 'param_names': params_info,

251 # 信頼帯の計算用。保存はしない(Dataset.fit_result には入れない)

252 'fit_func': fit_func,

253 'x_fit': x_fit,

254 'y_fit': y_fit,

255 'r_squared': r_squared,

256 'residuals': residuals,

257 'x_data_used': x_data,

258 'y_data_used': y_data,

259 'loss': loss,

260 }

261

262

263# 多峰分離の成分関数。ベースラインを全成分で共有するので、単峰版と違いオフセットを持たない。

264

265def _gaussian_component(x: Any, a: float, b: float, c: float) -> Any:

266 return a * np.exp(-((x - b) ** 2) / (2 * c ** 2))

267

268

269def _lorentzian_component(x: Any, a: float, b: float, c: float) -> Any:

270 return a / (1 + ((x - b) / c) ** 2)

271

272

273def _pseudo_voigt_component(x: Any, a: float, b: float, c: float, eta: float) -> Any:

274 lorentzian_shape = 1 / (1 + ((x - b) / c) ** 2)

275 gaussian_shape = np.exp(-4 * np.log(2) * ((x - b) / c) ** 2)

276 return a * (eta * lorentzian_shape + (1 - eta) * gaussian_shape)

277

278

279def _voigt_component(x: Any, a: float, b: float, sigma: float, gamma: float) -> Any:

280 z = ((x - b) + 1j * gamma) / (sigma * np.sqrt(2))

281 return a * np.real(wofz(z)) / (sigma * np.sqrt(2 * np.pi))

282

283

284def _constant_baseline(x: Any, c: float) -> Any:

285 return np.full_like(np.asarray(x, dtype=float), c)

286

287

288def _linear_baseline(x: Any, m: float, b: float) -> Any:

289 return m * np.asarray(x, dtype=float) + b

290

291

292_MULTI_PEAK_COMPONENT_TYPES: dict[str, dict[str, Any]] = {

293 'gaussian': {'label': 'ガウシアン', 'func': _gaussian_component, 'param_names': ['a', 'b', 'c']},

294 'lorentzian': {'label': 'ローレンツ', 'func': _lorentzian_component, 'param_names': ['a', 'b', 'c']},

295 'pseudo_voigt': {'label': '擬似フォークト', 'func': _pseudo_voigt_component, 'param_names': ['a', 'b', 'c', 'eta']},

296 'voigt': {'label': 'フォークト', 'func': _voigt_component, 'param_names': ['a', 'b', 'sigma', 'gamma']},

297}

298

299

300def calculate_multi_peak_fit(x_data: Any, y_data: Any, component_type: str, initial_guesses: list[dict[str, float]],

301 baseline_type: str = 'constant', x_range: tuple[float, float] | None = None,

302 sigma: Any = None, p0_overrides: dict[str, float] | None = None,

303 fixed_params: dict[str, float] | None = None,

304 bounds: dict[str, tuple[float, float]] | None = None) -> dict[str, Any]:

305 """同じ種類の N 成分とベースラインを同時にフィットする。

306

307 initial_guesses は [{'center', 'height', 'width'(FWHM)}, ...]。戻り値は calculate_curve_fit() の

308 dict に component_type / n_components / baseline_type / components を足したもの。

309 """

310 if component_type not in _MULTI_PEAK_COMPONENT_TYPES:

311 raise ValueError(f"不明な成分タイプ: {component_type}")

312 if not initial_guesses:

313 raise ValueError("少なくとも1つのピークの初期値が必要です。")

314 if baseline_type not in ('none', 'constant', 'linear'):

315 raise ValueError(f"不明なベースラインタイプ: {baseline_type}")

316

317 x_data = np.asarray(x_data, dtype=float)

318 y_data = np.asarray(y_data, dtype=float)

319 if sigma is not None: 319 ↛ 320line 319 didn't jump to line 320 because the condition on line 319 was never true

320 sigma = np.asarray(sigma, dtype=float)

321 if x_range is not None:

322 x_min, x_max = x_range

323 range_mask = (x_data >= x_min) & (x_data <= x_max)

324 x_data, y_data = x_data[range_mask], y_data[range_mask]

325 if sigma is not None: 325 ↛ 326line 325 didn't jump to line 326 because the condition on line 325 was never true

326 sigma = sigma[range_mask]

327

328 # calculate_curve_fit() と同じ理由で NaN を除く

329 nan_mask = np.isnan(x_data) | np.isnan(y_data)

330 if sigma is not None: 330 ↛ 331line 330 didn't jump to line 331 because the condition on line 330 was never true

331 nan_mask |= np.isnan(sigma)

332 if nan_mask.any():

333 x_data, y_data = x_data[~nan_mask], y_data[~nan_mask]

334 if sigma is not None: 334 ↛ 335line 334 didn't jump to line 335 because the condition on line 334 was never true

335 sigma = sigma[~nan_mask]

336

337 if len(x_data) == 0: 337 ↛ 338line 337 didn't jump to line 338 because the condition on line 337 was never true

338 raise ValueError("有効なデータ点がありません(すべて欠損値です)。フィッティングできません。")

339

340 component_info = _MULTI_PEAK_COMPONENT_TYPES[component_type]

341 component_func = component_info['func']

342 component_param_names = component_info['param_names']

343 n_component_params = len(component_param_names)

344 n_components = len(initial_guesses)

345

346 params_info: list[str] = []

347 p0: list[float] = []

348 for i, guess in enumerate(initial_guesses, start=1):

349 height = float(guess['height'])

350 center = float(guess['center'])

351 width = float(guess.get('width') or 1.0) or 1.0

352 params_info.extend(f"{p}{i}" for p in component_param_names)

353 if component_type == 'gaussian':

354 p0.extend([height, center, width])

355 elif component_type == 'lorentzian':

356 # width は FWHM、c は HWHM

357 p0.extend([height, center, (width / 2) or 1.0])

358 elif component_type == 'pseudo_voigt':

359 # この定義の c は FWHM そのもの

360 p0.extend([height, center, width, 0.5])

361 else: # voigt

362 # 単峰のフォークトと同じ初期値の割り振り

363 sigma0 = (width / 2.355) or 1.0

364 gamma0 = (width / 4) or 1.0

365 amplitude0 = height * sigma0 * np.sqrt(2 * np.pi)

366 p0.extend([amplitude0, center, sigma0, gamma0])

367

368 baseline_func: Callable[..., Any] | None

369 if baseline_type == 'constant':

370 params_info.append('baseline_c')

371 p0.append(float(np.nanmin(y_data)))

372 baseline_func = _constant_baseline

373 elif baseline_type == 'linear':

374 params_info.extend(['baseline_m', 'baseline_b'])

375 p0.extend([0.0, float(np.nanmin(y_data))])

376 baseline_func = _linear_baseline

377 else:

378 baseline_func = None

379

380 def fit_func(x: Any, *params: float) -> Any:

381 total = np.zeros_like(np.asarray(x, dtype=float))

382 for i in range(n_components):

383 comp_params = params[i * n_component_params:(i + 1) * n_component_params]

384 total = total + component_func(x, *comp_params)

385 if baseline_func is not None:

386 baseline_params = params[n_components * n_component_params:]

387 total = total + baseline_func(x, *baseline_params)

388 return total

389

390 fit_type_label = f"多峰分離({component_info['label']} x{n_components})"

391 popt, pcov = _run_curve_fit_with_overrides(

392 fit_func, params_info, p0, x_data, y_data, sigma,

393 p0_overrides, fixed_params, bounds, fit_type_label,

394 )

395

396 perr = np.sqrt(np.diag(pcov))

397 x_fit = np.linspace(x_data.min(), x_data.max(), 200)

398 y_fit = fit_func(x_fit, *popt)

399 residuals = y_data - fit_func(x_data, *popt)

400 ss_res = np.sum(residuals ** 2)

401 ss_tot = np.sum((y_data - np.mean(y_data)) ** 2)

402 r_squared = 1.0 if ss_tot == 0 else 1.0 - (ss_res / ss_tot)

403

404 components = []

405 for i in range(n_components):

406 start = i * n_component_params

407 end = start + n_component_params

408 components.append({

409 'type': component_type,

410 'param_names': list(params_info[start:end]),

411 'params': [float(v) for v in popt[start:end]],

412 })

413

414 return {

415 'popt': popt,

416 'pcov': pcov,

417 'perr': perr,

418 'param_names': params_info,

419 'fit_func': fit_func,

420 'x_fit': x_fit,

421 'y_fit': y_fit,

422 'r_squared': r_squared,

423 'residuals': residuals,

424 'x_data_used': x_data,

425 'y_data_used': y_data,

426 'component_type': component_type,

427 'n_components': n_components,

428 'baseline_type': baseline_type,

429 'components': components,

430 }

431

432

433def fit_curve_task(x_data: Any, y_data: Any, fit_type: str, custom_formula: str | None = None, sigma: Any = None,

434 x_range: tuple[float, float] | None = None, p0_overrides: dict[str, float] | None = None,

435 fixed_params: dict[str, float] | None = None,

436 bounds: dict[str, tuple[float, float]] | None = None, loss: str = 'linear',

437 report_progress: Callable[..., Any] | None = None,

438 is_cancelled: Callable[[], bool] | None = None) -> dict[str, Any]:

439 """TaskRunner 用。curve_fit は進捗も中断もできないので、report_progress と is_cancelled は受け取るだけ。"""

440 return calculate_curve_fit(

441 x_data, y_data, fit_type, custom_formula=custom_formula, sigma=sigma, x_range=x_range,

442 p0_overrides=p0_overrides, fixed_params=fixed_params, bounds=bounds, loss=loss,

443 )

444

445

446def multi_peak_fit_task(x_data: Any, y_data: Any, component_type: str, initial_guesses: list[dict[str, float]],

447 baseline_type: str = 'constant', x_range: tuple[float, float] | None = None,

448 sigma: Any = None, p0_overrides: dict[str, float] | None = None,

449 fixed_params: dict[str, float] | None = None,

450 bounds: dict[str, tuple[float, float]] | None = None,

451 report_progress: Callable[..., Any] | None = None,

452 is_cancelled: Callable[[], bool] | None = None) -> dict[str, Any]:

453 """TaskRunner 用。fit_curve_task() と同じく進捗と中断は受け取るだけ。"""

454 return calculate_multi_peak_fit(

455 x_data, y_data, component_type, initial_guesses, baseline_type=baseline_type,

456 x_range=x_range, sigma=sigma, p0_overrides=p0_overrides, fixed_params=fixed_params, bounds=bounds,

457 )

458

459

460_INTEGRAL_METHODS = ("trapezoid", "simpson")

461

462

463def _trapezoid_integrate(y: Any, x: Any) -> float:

464 """NumPy 2.0 で trapz が trapezoid になった。古い版でも動くようにする。"""

465 trapezoid_func = getattr(np, "trapezoid", None)

466 if trapezoid_func is None: 466 ↛ 467line 466 didn't jump to line 467 because the condition on line 466 was never true

467 trapezoid_func = np.trapz

468 return float(trapezoid_func(y, x))

469

470

471def calculate_interval_integral(x_data: Any, y_data: Any, x_range: tuple[float, float], method: str = "trapezoid",

472 subtract_baseline: bool = False) -> dict[str, Any]:

473 """x_range(両端を含む)で y を x について積分する。

474

475 subtract_baseline は範囲の両端を結ぶ直線を引いてから積分する。両端の Y は、その X 位置で補間した値。

476 """

477 if method not in _INTEGRAL_METHODS:

478 raise ValueError(f"未知の積分方法です: {method}")

479 if x_range is None or len(x_range) != 2: 479 ↛ 480line 479 didn't jump to line 480 because the condition on line 479 was never true

480 raise ValueError("積分範囲(x_range)を指定してください。")

481

482 x_min, x_max = float(x_range[0]), float(x_range[1])

483 if x_min >= x_max:

484 raise ValueError("積分範囲の最小値は最大値より小さい値である必要があります。")

485

486 x_data = np.asarray(x_data, dtype=float)

487 y_data = np.asarray(y_data, dtype=float)

488

489 # NaN があると積分値が nan になる

490 nan_mask = np.isnan(x_data) | np.isnan(y_data)

491 if nan_mask.any():

492 x_data, y_data = x_data[~nan_mask], y_data[~nan_mask]

493

494 if len(x_data) == 0:

495 raise ValueError("有効なデータ点がありません(すべて欠損値です)。積分できません。")

496

497 # trapezoid / simpson は X が単調増加であることを前提にする

498 order = np.argsort(x_data)

499 x_sorted, y_sorted = x_data[order], y_data[order]

500

501 x_min_data, x_max_data = float(x_sorted[0]), float(x_sorted[-1])

502 if x_min < x_min_data or x_max > x_max_data:

503 raise ValueError(

504 f"積分範囲はデータのX範囲({x_min_data:.6g} 〜 {x_max_data:.6g})内で指定してください。"

505 )

506

507 range_mask = (x_sorted >= x_min) & (x_sorted <= x_max)

508 x_in, y_in = x_sorted[range_mask], y_sorted[range_mask]

509

510 if len(x_in) < 2:

511 raise ValueError(f"積分範囲内に十分なデータ点がありません(最低2点必要、現在{len(x_in)}点)。")

512

513 baseline_used = None

514 y_for_integration = y_in

515 if subtract_baseline:

516 y_at_min = float(np.interp(x_min, x_sorted, y_sorted))

517 y_at_max = float(np.interp(x_max, x_sorted, y_sorted))

518 baseline_used = y_at_min + (y_at_max - y_at_min) * (x_in - x_min) / (x_max - x_min)

519 y_for_integration = y_in - baseline_used

520

521 if method == "trapezoid":

522 integral = _trapezoid_integrate(y_for_integration, x_in)

523 else:

524 integral = float(simpson(y_for_integration, x=x_in))

525

526 return {

527 'integral': integral,

528 'method': method,

529 'x_range': (x_min, x_max),

530 'subtract_baseline': subtract_baseline,

531 'x_used': x_in,

532 'y_used': y_for_integration,

533 'y_raw_used': y_in,

534 'baseline_used': baseline_used,

535 'n_points': len(x_in),

536 }

537

538

539def calculate_cumulative_integral(x_data: Any, y_data: Any, method: str = "trapezoid") -> dict[str, Any]:

540 """各点までの累積積分 ∫[x_min, x_i] y dx を返す(先頭は 0)。"""

541 if method not in _INTEGRAL_METHODS:

542 raise ValueError(f"未知の積分方法です: {method}")

543

544 x_data = np.asarray(x_data, dtype=float)

545 y_data = np.asarray(y_data, dtype=float)

546

547 nan_mask = np.isnan(x_data) | np.isnan(y_data)

548 if nan_mask.any():

549 x_data, y_data = x_data[~nan_mask], y_data[~nan_mask]

550

551 if len(x_data) < 2:

552 raise ValueError("有効なデータ点が不足しています(最低2点必要です)。")

553

554 order = np.argsort(x_data)

555 x_sorted, y_sorted = x_data[order], y_data[order]

556

557 if method == "trapezoid":

558 y_cumulative = cumulative_trapezoid(y_sorted, x_sorted, initial=0.0)

559 else:

560 y_cumulative = cumulative_simpson(y_sorted, x=x_sorted, initial=0.0)

561

562 return {

563 'x_used': x_sorted,

564 'y_cumulative': y_cumulative,

565 'method': method,

566 'n_points': len(x_sorted),

567 }

568

569

570def _peak_detection_signal_and_kwargs(x_data: Any, y_data: Any, peak_type: str,

571 settings: dict[str, Any]) -> tuple[np.ndarray, dict[str, Any]]:

572 """calculate_peaks と calculate_peak_quantification の共通の前処理。

573

574 谷は -y_data で探すので height も符号を反転する(利用者は谷でも「Y < -10」を -10 と入力する)。

575 height=None は高さで絞らない。

576 """

577 y_data_to_find = y_data

578 height = settings["height"]

579

580 if "下に凸" in peak_type:

581 y_data_to_find = -y_data

582 if height is not None:

583 height = -height

584

585 kwargs = {"height": height}

586

587 if settings["prominence"] is not None:

588 kwargs["prominence"] = settings["prominence"]

589

590 if len(x_data) > 1 and settings["distance_x"] > 0:

591 sorted_x = np.sort(x_data)

592 avg_x_diff = np.mean(np.diff(sorted_x))

593 if avg_x_diff > 0: 593 ↛ 596line 593 didn't jump to line 596 because the condition on line 593 was always true

594 kwargs["distance"] = max(1, int(np.ceil(settings["distance_x"] / avg_x_diff)))

595 else:

596 kwargs["distance"] = 1

597 else:

598 kwargs["distance"] = 1

599

600 return y_data_to_find, kwargs

601

602

603def calculate_peaks(x_data: Any, y_data: Any, peak_type: str, settings: dict[str, Any]) -> tuple[np.ndarray, np.ndarray]:

604 y_data_to_find, kwargs = _peak_detection_signal_and_kwargs(x_data, y_data, peak_type, settings)

605

606 peak_indices, _ = find_peaks(y_data_to_find, **kwargs)

607

608 if len(peak_indices) == 0:

609 return np.array([]), np.array([])

610

611 return x_data[peak_indices], y_data[peak_indices]

612

613

614def calculate_peak_quantification(x_data: Any, y_data: Any, peak_type: str, settings: dict[str, Any]) -> dict[str, Any]:

615 """ピーク/谷の位置に加え、FWHM・面積・重心を返す。

616

617 谷は反転した信号の上で計算するので、面積は谷でも正の値(検出方向への突出量)になる。

618 ローカル基線は peak_widths(rel_height=1.0) の左右の裾を結ぶ直線(裾の高さが違う

619 非対称なピークでも面積と重心が定まる)。重心の重みは生の Y ではなく基線からの高さ

620 (大きなオフセットの上のピークで重心が引っ張られない)。x_data の並びがそのまま

621 ピークの前後関係になるので、ソートは呼び出し側で行う。

622 """

623 x_data = np.asarray(x_data, dtype=float)

624 y_data = np.asarray(y_data, dtype=float)

625

626 y_signal, kwargs = _peak_detection_signal_and_kwargs(x_data, y_data, peak_type, settings)

627 peak_indices, _ = find_peaks(y_signal, **kwargs)

628

629 empty = np.array([])

630 if len(peak_indices) == 0:

631 return {'peak_x': empty, 'peak_y': empty, 'fwhm': empty, 'area': empty, 'centroid': empty}

632

633 idx_axis = np.arange(len(x_data))

634

635 # 幅はインデックス単位で出るので、x_data に補間して X の単位にする(等間隔でなくてよい)

636 _, _, left_ips_half, right_ips_half = peak_widths(y_signal, peak_indices, rel_height=0.5)

637 x_left_half = np.interp(left_ips_half, idx_axis, x_data)

638 x_right_half = np.interp(right_ips_half, idx_axis, x_data)

639 fwhm = x_right_half - x_left_half

640

641 # rel_height=1.0 は prominence を測った高さでの幅、つまりピークの裾

642 _, _, left_ips_full, right_ips_full = peak_widths(y_signal, peak_indices, rel_height=1.0)

643

644 area = np.empty(len(peak_indices))

645 centroid = np.empty(len(peak_indices))

646

647 for i, pk in enumerate(peak_indices):

648 li, ri = left_ips_full[i], right_ips_full[i]

649

650 x_left = np.interp(li, idx_axis, x_data)

651 x_right = np.interp(ri, idx_axis, x_data)

652 y_left = np.interp(li, idx_axis, y_signal)

653 y_right = np.interp(ri, idx_axis, y_signal)

654

655 li_idx, ri_idx = int(np.floor(li)), int(np.ceil(ri))

656 inner_idx = np.arange(max(li_idx, 0), min(ri_idx, len(x_data) - 1) + 1)

657 xs_inner, ys_inner = x_data[inner_idx], y_signal[inner_idx]

658 inner_mask = (xs_inner >= x_left) & (xs_inner <= x_right)

659 xs_inner, ys_inner = xs_inner[inner_mask], ys_inner[inner_mask]

660

661 xs = np.concatenate(([x_left], xs_inner, [x_right]))

662 ys = np.concatenate(([y_left], ys_inner, [y_right]))

663 xs, unique_order = np.unique(xs, return_index=True)

664 ys = ys[unique_order]

665

666 baseline = np.interp(xs, [x_left, x_right], [y_left, y_right])

667 y_above = ys - baseline

668

669 area[i] = np.trapezoid(y_above, xs)

670

671 weights = np.clip(y_above, 0, None)

672 total_weight = np.sum(weights)

673 if total_weight > 0: 673 ↛ 676line 673 didn't jump to line 676 because the condition on line 673 was always true

674 centroid[i] = np.sum(xs * weights) / total_weight

675 else:

676 centroid[i] = x_data[pk]

677

678 return {

679 'peak_x': x_data[peak_indices],

680 'peak_y': y_data[peak_indices],

681 'fwhm': fwhm,

682 'area': area,

683 'centroid': centroid,

684 }

685

686

687def calculate_savgol(x_data: Any, y_data: Any, window_length: int, polyorder: int,

688 deriv: int = 0) -> tuple[np.ndarray, np.ndarray]:

689 """Savitzky-Golay による平滑化(deriv=0)または微分(deriv=1, 2)。

690

691 delta に X の間隔の中央値を渡すので、微分は dy/dx の大きさになる(X がほぼ等間隔である前提)。

692 """

693 if window_length % 2 == 0:

694 raise ValueError("窓幅(window_length)は奇数である必要があります。")

695 if polyorder >= window_length:

696 raise ValueError("多項式の次数は窓幅より小さくする必要があります。")

697 if window_length > len(y_data):

698 raise ValueError(f"窓幅({window_length})がデータ点数({len(y_data)})を超えています。")

699

700 order = np.argsort(x_data)

701 x_sorted, y_sorted = x_data[order], y_data[order]

702 diffs = np.diff(x_sorted)

703 dx = float(np.median(diffs)) if len(diffs) > 0 and np.median(diffs) > 0 else 1.0

704

705 y_result = savgol_filter(y_sorted, window_length, polyorder, deriv=deriv, delta=dx)

706 return x_sorted, y_result

707

708

709# 線の平滑化(Dataset.smoothing_method)。描画用で元データは変えない。

710# どれも X でソートして実データ点のまま返す(CubicSpline と違い点を補間で増やさない)。

711

712def calculate_moving_average_smooth(x_data: Any, y_data: Any, window: int = 5) -> tuple[np.ndarray, np.ndarray]:

713 """window は点数。データ点数を超えたら切り詰める。端は端の値を延長する。"""

714 x_data = np.asarray(x_data, dtype=float)

715 y_data = np.asarray(y_data, dtype=float)

716 order = np.argsort(x_data)

717 x_sorted, y_sorted = x_data[order], y_data[order]

718

719 n = len(y_sorted)

720 w = max(1, min(int(window), n)) if n > 0 else 1

721 if w <= 1:

722 return x_sorted, y_sorted.copy()

723 return x_sorted, uniform_filter1d(y_sorted, size=w, mode='nearest')

724

725

726def calculate_median_smooth(x_data: Any, y_data: Any, window: int = 5) -> tuple[np.ndarray, np.ndarray]:

727 x_data = np.asarray(x_data, dtype=float)

728 y_data = np.asarray(y_data, dtype=float)

729 order = np.argsort(x_data)

730 x_sorted, y_sorted = x_data[order], y_data[order]

731

732 n = len(y_sorted)

733 w = max(1, min(int(window), n)) if n > 0 else 1

734 if w <= 1: 734 ↛ 735line 734 didn't jump to line 735 because the condition on line 734 was never true

735 return x_sorted, y_sorted.copy()

736 return x_sorted, median_filter(y_sorted, size=w, mode='nearest')

737

738

739def calculate_gaussian_smooth(x_data: Any, y_data: Any, sigma: float = 2.0) -> tuple[np.ndarray, np.ndarray]:

740 """sigma はデータ点のインデックス単位(X の間隔ではない)。"""

741 x_data = np.asarray(x_data, dtype=float)

742 y_data = np.asarray(y_data, dtype=float)

743 order = np.argsort(x_data)

744 x_sorted, y_sorted = x_data[order], y_data[order]

745

746 if len(y_sorted) < 2 or sigma <= 0:

747 return x_sorted, y_sorted.copy()

748 return x_sorted, gaussian_filter1d(y_sorted, sigma=float(sigma), mode='nearest')

749

750

751# ベースライン補正。どの手法も (ソート済みの x, ベースライン, 差し引いた y) を返す。

752

753def _sort_xy_for_baseline(x_data: Any, y_data: Any) -> tuple[np.ndarray, np.ndarray]:

754 x_data = np.asarray(x_data, dtype=float)

755 y_data = np.asarray(y_data, dtype=float)

756 order = np.argsort(x_data)

757 return x_data[order], y_data[order]

758

759

760def calculate_baseline_als(x_data: Any, y_data: Any, lam: float = 1e5, p: float = 0.01,

761 niter: int = 10) -> tuple[np.ndarray, np.ndarray, np.ndarray]:

762 """Asymmetric Least Squares (Eilers & Boelens 2005)。

763

764 lam が大きいほど滑らか、p が小さいほどピークを避けて下側を通る(目安 0.001〜0.1)。

765 """

766 if lam <= 0:

767 raise ValueError("lam(平滑化パラメータ)は正の値である必要があります。")

768 if not (0 < p < 1):

769 raise ValueError("p(非対称重み)は0より大きく1より小さい値である必要があります。")

770 if niter < 1:

771 raise ValueError("反復回数(niter)は1以上である必要があります。")

772

773 x_sorted, y_sorted = _sort_xy_for_baseline(x_data, y_data)

774 n_points = len(y_sorted)

775 if n_points < 3:

776 raise ValueError(f"ALSベースライン補正には少なくとも3点のデータが必要です(現在{n_points}点)。")

777

778 # 2階差分の二乗和へのペナルティ

779 diff_matrix = sparse.diags([1, -2, 1], [0, -1, -2], shape=(n_points, n_points - 2))

780 penalty = lam * (diff_matrix @ diff_matrix.transpose())

781

782 weights = np.ones(n_points)

783 baseline = y_sorted.copy()

784 for _ in range(niter):

785 weight_matrix = sparse.diags(weights, 0, shape=(n_points, n_points))

786 baseline = spsolve((weight_matrix + penalty).tocsc(), weights * y_sorted)

787 weights = p * (y_sorted > baseline) + (1 - p) * (y_sorted <= baseline)

788

789 return x_sorted, baseline, y_sorted - baseline

790

791

792def calculate_baseline_polynomial(x_data: Any, y_data: Any, degree: int = 3,

793 iterations: int = 10) -> tuple[np.ndarray, np.ndarray, np.ndarray]:

794 """反復多項式フィット(Lieber & Mahadevan-Jansen 2003 の ModPoly)。

795

796 フィット曲線を上回る点をフィット値で置き換えて再フィットし、ピークを少しずつ削る。

797 """

798 if degree < 0:

799 raise ValueError("多項式の次数は0以上である必要があります。")

800 if iterations < 1:

801 raise ValueError("反復回数(iterations)は1以上である必要があります。")

802

803 x_sorted, y_sorted = _sort_xy_for_baseline(x_data, y_data)

804 if degree >= len(y_sorted):

805 raise ValueError(f"多項式の次数({degree})がデータ点数({len(y_sorted)})以上です。")

806

807 work_y = y_sorted.copy()

808 baseline = work_y

809 for _ in range(iterations):

810 coeffs = np.polyfit(x_sorted, work_y, degree)

811 baseline = np.polyval(coeffs, x_sorted)

812 work_y = np.minimum(work_y, baseline)

813

814 return x_sorted, baseline, y_sorted - baseline

815

816

817def calculate_baseline_rubberband(x_data: Any, y_data: Any) -> tuple[np.ndarray, np.ndarray, np.ndarray]:

818 """ラバーバンド法。下側凸包の頂点を区分線形でつないだものをベースラインにする。"""

819 x_sorted, y_sorted = _sort_xy_for_baseline(x_data, y_data)

820 n_points = len(x_sorted)

821 if n_points < 3:

822 raise ValueError(f"ラバーバンド法には少なくとも3点のデータが必要です(現在{n_points}点)。")

823

824 # Andrew の monotone chain の下側だけ

825 hull_indices: list[int] = []

826 for i in range(n_points):

827 while len(hull_indices) >= 2:

828 o, a = hull_indices[-2], hull_indices[-1]

829 cross = ((x_sorted[a] - x_sorted[o]) * (y_sorted[i] - y_sorted[o])

830 - (y_sorted[a] - y_sorted[o]) * (x_sorted[i] - x_sorted[o]))

831 if cross <= 0:

832 hull_indices.pop()

833 else:

834 break

835 hull_indices.append(i)

836

837 x_hull = x_sorted[hull_indices]

838 y_hull = y_sorted[hull_indices]

839 baseline = np.interp(x_sorted, x_hull, y_hull)

840

841 return x_sorted, baseline, y_sorted - baseline

842

843

844def calculate_baseline_manual(x_data: Any, y_data: Any, anchor_x: Any,

845 method: str = "linear") -> tuple[np.ndarray, np.ndarray, np.ndarray]:

846 """指定した X 位置でデータを補間した点を、線形("linear")か 3 次スプライン("spline")で結ぶ。"""

847 if method not in ("linear", "spline"):

848 raise ValueError(f"未知の補間方法です: {method}")

849

850 x_sorted, y_sorted = _sort_xy_for_baseline(x_data, y_data)

851

852 anchor_x = np.unique(np.asarray(anchor_x, dtype=float))

853 if len(anchor_x) < 2:

854 raise ValueError("アンカー点は重複しない値で2点以上指定してください。")

855 if method == "spline" and len(anchor_x) < 3:

856 raise ValueError("スプライン補間には3点以上のアンカー点が必要です。")

857

858 x_min, x_max = x_sorted[0], x_sorted[-1]

859 if anchor_x[0] < x_min or anchor_x[-1] > x_max:

860 raise ValueError(

861 f"アンカー点はデータのX範囲({x_min:.6g} 〜 {x_max:.6g})内で指定してください。"

862 )

863

864 anchor_y = np.interp(anchor_x, x_sorted, y_sorted)

865

866 if method == "linear":

867 baseline = np.interp(x_sorted, anchor_x, anchor_y)

868 else:

869 baseline = CubicSpline(anchor_x, anchor_y)(x_sorted)

870

871 return x_sorted, baseline, y_sorted - baseline

872

873

874def calculate_confidence_band(x_eval: Any, fit_func: Callable[..., Any], popt: Any, pcov: Any, residuals: Any,

875 confidence: float = 0.95, band_type: str = "confidence") -> dict[str, Any]:

876 """非線形回帰の信頼帯・予測帯をデルタ法(線形化)で近似する。

877

878 各点で J(x) @ pcov @ J(x).T を分散とする。予測帯はさらに残差の平均二乗誤差を足す。

879 固定パラメータは pcov の行と列が 0 なので不確かさに寄与しない。

880 """

881 if band_type not in ("confidence", "prediction"):

882 raise ValueError(f"未知のband_typeです: {band_type}")

883 if not (0.0 < confidence < 1.0):

884 raise ValueError("confidence(信頼水準)は0より大きく1より小さい値である必要があります。")

885

886 x_eval = np.asarray(x_eval, dtype=float)

887 popt = np.asarray(popt, dtype=float)

888 pcov = np.asarray(pcov, dtype=float)

889 residuals = np.asarray(residuals, dtype=float)

890

891 n_params = len(popt)

892 n_data = len(residuals)

893 dof = n_data - n_params

894 if dof < 1:

895 raise ValueError(

896 f"自由度(データ点数{n_data} - パラメータ数{n_params})が1未満のため、"

897 "信頼帯・予測帯を計算できません。"

898 )

899

900 from scipy.stats import t as _t_dist

901 t_value = _t_dist.ppf(1.0 - (1.0 - confidence) / 2.0, dof)

902

903 y_center = fit_func(x_eval, *popt)

904 jacobian = np.empty((len(x_eval), n_params))

905 for i in range(n_params):

906 # 刻みはパラメータの大きさに比例させる(極端な値でも数値誤差が出にくい)

907 step = max(abs(popt[i]), 1.0) * 1e-6

908 popt_plus, popt_minus = popt.copy(), popt.copy()

909 popt_plus[i] += step

910 popt_minus[i] -= step

911 jacobian[:, i] = (fit_func(x_eval, *popt_plus) - fit_func(x_eval, *popt_minus)) / (2 * step)

912

913 # diag(J @ pcov @ J.T) だけを求める

914 var_yhat = np.einsum('ij,jk,ik->i', jacobian, pcov, jacobian)

915 # 数値誤差でわずかに負になることがある

916 var_yhat = np.clip(var_yhat, 0.0, None)

917

918 if band_type == "prediction":

919 mse = np.sum(residuals ** 2) / dof

920 variance = var_yhat + mse

921 else:

922 variance = var_yhat

923

924 margin = t_value * np.sqrt(variance)

925

926 return {

927 'y_center': y_center,

928 'y_lower': y_center - margin,

929 'y_upper': y_center + margin,

930 'confidence': confidence,

931 'band_type': band_type,

932 }

933

934

935def calculate_resample_to_grid(x_data: Any, y_data: Any, target_x: Any, method: str = "linear",

936 extrapolate: bool = False) -> np.ndarray:

937 """(x_data, y_data) を target_x の格子に補間する。

938

939 extrapolate=False(既定)は元の X の範囲外を NaN にする。True のとき "linear" は両端の傾きで

940 外挿する(np.interp は端の値で止めるだけで外挿しない)。

941 """

942 if method not in ("linear", "cubic"):

943 raise ValueError(f"未知の補間方法です: {method}")

944

945 x_data = np.asarray(x_data, dtype=float)

946 y_data = np.asarray(y_data, dtype=float)

947 target_x = np.asarray(target_x, dtype=float)

948

949 nan_mask = np.isnan(x_data) | np.isnan(y_data)

950 if nan_mask.any():

951 x_data, y_data = x_data[~nan_mask], y_data[~nan_mask]

952

953 if len(x_data) == 0:

954 raise ValueError("有効なデータ点がありません(すべて欠損値です)。リサンプリングできません。")

955

956 order = np.argsort(x_data)

957 x_sorted, y_sorted = x_data[order], y_data[order]

958 # 補間は X が厳密に単調増加であることを要求するので、同じ X は最後の値を残してまとめる。

959 # np.unique は最初の出現を返すので、反転してから除く。

960 x_rev, y_rev = x_sorted[::-1], y_sorted[::-1]

961 x_sorted, rev_unique_indices = np.unique(x_rev, return_index=True)

962 y_sorted = y_rev[rev_unique_indices]

963

964 n_points = len(x_sorted)

965 if method == "linear" and n_points < 2:

966 raise ValueError(f"線形補間には少なくとも2点のデータが必要です(現在{n_points}点)。")

967 if method == "cubic" and n_points < 4:

968 raise ValueError(f"3次スプライン補間には少なくとも4点のデータが必要です(現在{n_points}点)。")

969

970 x_min, x_max = x_sorted[0], x_sorted[-1]

971

972 if method == "linear":

973 if extrapolate:

974 # np.interp は範囲外を端の値で止めるだけなので、両端の傾きで外挿する

975 result = np.interp(target_x, x_sorted, y_sorted)

976 below = target_x < x_min

977 if below.any(): 977 ↛ 980line 977 didn't jump to line 980 because the condition on line 977 was always true

978 slope = (y_sorted[1] - y_sorted[0]) / (x_sorted[1] - x_sorted[0])

979 result[below] = y_sorted[0] + slope * (target_x[below] - x_min)

980 above = target_x > x_max

981 if above.any(): 981 ↛ 997line 981 didn't jump to line 997 because the condition on line 981 was always true

982 slope = (y_sorted[-1] - y_sorted[-2]) / (x_sorted[-1] - x_sorted[-2])

983 result[above] = y_sorted[-1] + slope * (target_x[above] - x_max)

984 else:

985 result = np.interp(target_x, x_sorted, y_sorted)

986 out_of_range = (target_x < x_min) | (target_x > x_max)

987 result = result.astype(float)

988 result[out_of_range] = np.nan

989 else:

990 spline = CubicSpline(x_sorted, y_sorted, extrapolate=extrapolate)

991 result = spline(target_x)

992 if not extrapolate:

993 out_of_range = (target_x < x_min) | (target_x > x_max)

994 result = np.asarray(result, dtype=float)

995 result[out_of_range] = np.nan

996

997 return result

998

999

1000def calculate_average_duplicate_x(x_data: Any, y_data: Any) -> dict[str, Any]:

1001 """同じ X の行の Y を平均する。重複の「除去」は index ラベルが要るので呼び出し側で行う。"""

1002 x_data = np.asarray(x_data, dtype=float)

1003 y_data = np.asarray(y_data, dtype=float)

1004

1005 nan_mask = np.isnan(x_data) | np.isnan(y_data)

1006 if nan_mask.any():

1007 x_data, y_data = x_data[~nan_mask], y_data[~nan_mask]

1008

1009 if len(x_data) == 0:

1010 raise ValueError("有効なデータ点がありません(すべて欠損値です)。")

1011

1012 order = np.argsort(x_data, kind='stable')

1013 x_sorted, y_sorted = x_data[order], y_data[order]

1014

1015 unique_x, inverse, counts = np.unique(x_sorted, return_inverse=True, return_counts=True)

1016 # NumPy の版によって inverse の形が (N,) か (N, 1) になる

1017 inverse = np.asarray(inverse).reshape(-1)

1018 sums = np.zeros(len(unique_x))

1019 np.add.at(sums, inverse, y_sorted)

1020 y_averaged = sums / counts

1021

1022 return {

1023 'x_used': unique_x,

1024 'y_averaged': y_averaged,

1025 'group_sizes': counts,

1026 'n_duplicate_groups': int(np.sum(counts > 1)),

1027 'n_points_in': len(x_sorted),

1028 'n_points_out': len(unique_x),

1029 }

1030

1031

1032# 外れ値の検出。マスクへの適用は利用者が選ぶので、ここでは検出だけ。

1033# is_outlier は入力と同じ長さと並び(NaN は False)。visible_df.index[is_outlier] でそのまま引けるように。

1034

1035def sample_standard_deviation(values: Any) -> float:

1036 """標本標準偏差(n−1 で割る)。有効な値が 2 未満なら NaN。

1037

1038 アプリ内の「標準偏差」はすべてこれで計算する。場所によって n と n−1 が混ざると同じデータで値が違って見える。

1039 """

1040 arr = np.asarray(values, dtype=float)

1041 arr = arr[~np.isnan(arr)]

1042 if arr.size < 2:

1043 return float('nan')

1044 return float(np.std(arr, ddof=1))

1045

1046

1047def calculate_zscore_outliers(y_data: Any, threshold: float = 3.0) -> dict[str, Any]:

1048 if threshold <= 0:

1049 raise ValueError("しきい値は正の値である必要があります。")

1050

1051 y = np.asarray(y_data, dtype=float)

1052 is_outlier = np.zeros(len(y), dtype=bool)

1053 z_scores = np.full(len(y), np.nan)

1054 valid = ~np.isnan(y)

1055

1056 if valid.sum() >= 2: 1056 ↛ 1063line 1056 didn't jump to line 1063 because the condition on line 1056 was always true

1057 mean = np.mean(y[valid])

1058 std = sample_standard_deviation(y[valid])

1059 if std > 0:

1060 z_scores[valid] = (y[valid] - mean) / std

1061 is_outlier[valid] = np.abs(z_scores[valid]) > threshold

1062

1063 return {

1064 'is_outlier': is_outlier,

1065 'z_scores': z_scores,

1066 'threshold': threshold,

1067 'n_outliers': int(np.sum(is_outlier)),

1068 }

1069

1070

1071def calculate_iqr_outliers(y_data: Any, multiplier: float = 1.5) -> dict[str, Any]:

1072 """[Q1 - multiplier*IQR, Q3 + multiplier*IQR] の外を外れ値とする。"""

1073 if multiplier <= 0:

1074 raise ValueError("係数は正の値である必要があります。")

1075

1076 y = np.asarray(y_data, dtype=float)

1077 is_outlier = np.zeros(len(y), dtype=bool)

1078 valid = ~np.isnan(y)

1079 lower_bound = upper_bound = None

1080

1081 if valid.sum() >= 4:

1082 q1, q3 = np.percentile(y[valid], [25, 75])

1083 iqr = q3 - q1

1084 lower_bound = float(q1 - multiplier * iqr)

1085 upper_bound = float(q3 + multiplier * iqr)

1086 is_outlier[valid] = (y[valid] < lower_bound) | (y[valid] > upper_bound)

1087

1088 return {

1089 'is_outlier': is_outlier,

1090 'lower_bound': lower_bound,

1091 'upper_bound': upper_bound,

1092 'multiplier': multiplier,

1093 'n_outliers': int(np.sum(is_outlier)),

1094 }

1095

1096

1097def calculate_lttb_downsample(x_data: Any, y_data: Any, n_out: int) -> np.ndarray:

1098 """Largest-Triangle-Three-Buckets(Steinarsson 2013)による表示用の間引き。

1099

1100 x_data はソート済みであること。選んだ点のインデックス(先頭と末尾を含む昇順)を返すので、

1101 同じインデックスを誤差などほかの配列にも使える。n_out 以下なら間引かない。

1102 """

1103 n = len(x_data)

1104 if n_out < 3 or n <= n_out:

1105 return np.arange(n)

1106

1107 x_data = np.asarray(x_data, dtype=float)

1108 y_data = np.asarray(y_data, dtype=float)

1109

1110 selected = np.empty(n_out, dtype=np.int64)

1111 selected[0] = 0

1112 selected[-1] = n - 1

1113

1114 bucket_edges = np.linspace(1, n - 1, n_out - 1).astype(np.int64)

1115 a = 0

1116 for i in range(n_out - 2):

1117 bucket_start, bucket_end = bucket_edges[i], bucket_edges[i + 1]

1118 if bucket_end <= bucket_start: 1118 ↛ 1119line 1118 didn't jump to line 1119 because the condition on line 1118 was never true

1119 bucket_end = bucket_start + 1

1120 next_start = bucket_end

1121 next_end = bucket_edges[i + 2] if i + 2 < len(bucket_edges) else n

1122 if next_end <= next_start: 1122 ↛ 1123line 1122 didn't jump to line 1123 because the condition on line 1122 was never true

1123 next_end = min(next_start + 1, n)

1124 avg_x = np.mean(x_data[next_start:next_end])

1125 avg_y = np.mean(y_data[next_start:next_end])

1126

1127 cand_x = x_data[bucket_start:bucket_end]

1128 cand_y = y_data[bucket_start:bucket_end]

1129 ax_, ay_ = x_data[a], y_data[a]

1130 areas = np.abs(

1131 (ax_ - avg_x) * (cand_y - ay_) - (ax_ - cand_x) * (avg_y - ay_)

1132 )

1133 best_local = np.argmax(areas)

1134 chosen = bucket_start + best_local

1135 selected[i + 1] = chosen

1136 a = chosen

1137

1138 return selected

1139

1140

1141def calculate_histogram(data: Any, bins: Any = 'auto', density: bool = False) -> dict[str, Any]:

1142 """ビン中心と度数(density=True なら確率密度)を返す。"""

1143 data = np.asarray(data, dtype=float)

1144 data = data[~np.isnan(data)]

1145 if len(data) == 0:

1146 raise ValueError("有効なデータ点がありません(すべて欠損値です)。")

1147 counts, bin_edges = np.histogram(data, bins=bins, density=density)

1148 bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2

1149 return {

1150 'bin_centers': bin_centers,

1151 'bin_edges': bin_edges,

1152 'counts': counts,

1153 'n_points_used': len(data),

1154 }

1155

1156

1157def calculate_kde(data: Any, n_points: int = 200, bw_method: Any = None) -> dict[str, Any]:

1158 data = np.asarray(data, dtype=float)

1159 data = data[~np.isnan(data)]

1160 if len(data) < 2:

1161 raise ValueError("カーネル密度推定には少なくとも2点の有効なデータが必要です。")

1162 try:

1163 kde = gaussian_kde(data, bw_method=bw_method)

1164 x_grid = np.linspace(data.min(), data.max(), n_points)

1165 density = kde(x_grid)

1166 except np.linalg.LinAlgError:

1167 raise ValueError("データにばらつきが無いため、カーネル密度推定を計算できません。") from None

1168 return {'x_grid': x_grid, 'density': density, 'n_points_used': len(data)}

1169

1170

1171def calculate_error_propagation(operation: str, value_a: Any, error_a: Any, value_b: Any, error_b: Any) -> np.ndarray:

1172 """独立な誤差を仮定した誤差伝播。

1173

1174 除算は、分子がゼロのときに NaN にならないよう、比ではなく分母で割る形の式にしている。

1175 """

1176 a = np.asarray(value_a, dtype=float)

1177 ea = np.asarray(error_a, dtype=float)

1178 b = np.asarray(value_b, dtype=float)

1179 eb = np.asarray(error_b, dtype=float)

1180

1181 if operation in ("A - B", "B - A", "A + B"):

1182 return np.sqrt(ea ** 2 + eb ** 2)

1183 if operation == "A × B":

1184 return np.sqrt((ea * b) ** 2 + (eb * a) ** 2)

1185 with np.errstate(divide='ignore', invalid='ignore'):

1186 if operation == "A ÷ B":

1187 return np.sqrt((ea / b) ** 2 + (a * eb / b ** 2) ** 2)

1188 else:

1189 return np.sqrt((eb / a) ** 2 + (b * ea / a ** 2) ** 2)

1190

1191

1192# 共通グリッドの点数の上限(大きいデータでフリーズしないため)

1193_XCORR_MAX_GRID_POINTS = 20000

1194

1195

1196def calculate_cross_correlation_alignment(x_a: Any, y_a: Any, x_b: Any, y_b: Any) -> dict[str, float]:

1197 """B を A に重ねるための X のシフト量(x_b + shift)を相互相関のピークから求める。

1198

1199 両者の範囲全体を覆う等間隔グリッドに補間し(範囲外は 0)、平均を引いてから相関を取る。

1200 """

1201 x_a = np.asarray(x_a, dtype=float)

1202 y_a = np.asarray(y_a, dtype=float)

1203 x_b = np.asarray(x_b, dtype=float)

1204 y_b = np.asarray(y_b, dtype=float)

1205

1206 valid_a = ~(np.isnan(x_a) | np.isnan(y_a))

1207 valid_b = ~(np.isnan(x_b) | np.isnan(y_b))

1208 x_a, y_a = x_a[valid_a], y_a[valid_a]

1209 x_b, y_b = x_b[valid_b], y_b[valid_b]

1210 if len(x_a) < 2 or len(x_b) < 2:

1211 raise ValueError("有効なデータ点が不足しています(それぞれ最低2点必要です)。")

1212

1213 order_a = np.argsort(x_a)

1214 x_a, y_a = x_a[order_a], y_a[order_a]

1215 order_b = np.argsort(x_b)

1216 x_b, y_b = x_b[order_b], y_b[order_b]

1217

1218 step_a = np.median(np.diff(x_a))

1219 step_b = np.median(np.diff(x_b))

1220 candidates = [s for s in (step_a, step_b) if s and s > 0]

1221 if not candidates: 1221 ↛ 1222line 1221 didn't jump to line 1222 because the condition on line 1221 was never true

1222 raise ValueError("X間隔を推定できません(Xの値が一定、または重複しています)。")

1223 step = min(candidates)

1224

1225 lo = min(x_a.min(), x_b.min())

1226 hi = max(x_a.max(), x_b.max())

1227 n_points = int(np.ceil((hi - lo) / step)) + 1

1228 if n_points > _XCORR_MAX_GRID_POINTS: 1228 ↛ 1229line 1228 didn't jump to line 1229 because the condition on line 1228 was never true

1229 step = (hi - lo) / _XCORR_MAX_GRID_POINTS

1230 n_points = _XCORR_MAX_GRID_POINTS + 1

1231 grid = np.linspace(lo, hi, n_points)

1232

1233 # 範囲外を含む全体で平均を引くので中心化は厳密でないが、ピーク位置を探すには足りる

1234 sig_a = np.interp(grid, x_a, y_a, left=0.0, right=0.0)

1235 sig_b = np.interp(grid, x_b, y_b, left=0.0, right=0.0)

1236 sig_a = sig_a - sig_a.mean()

1237 sig_b = sig_b - sig_b.mean()

1238

1239 corr = correlate(sig_a, sig_b, mode='full')

1240 lags = correlation_lags(len(sig_a), len(sig_b), mode='full')

1241 best_idx = int(np.argmax(corr))

1242 shift = float(lags[best_idx] * step)

1243

1244 return {'shift': shift, 'grid_step': float(step), 'correlation_peak': float(corr[best_idx])}

1245

1246

1247def assign_peak_label_levels(peak_x_values: Any, x_axis_span: float, proximity_ratio: float = 0.05,

1248 max_levels: int = 4) -> list[int]:

1249 """X が近いピークのラベルに互い違いの段番号を割り当てる(ずらす量は呼び出し側で決める)。"""

1250 peak_x_values = list(peak_x_values)

1251 if x_axis_span <= 0 or not peak_x_values:

1252 return [0] * len(peak_x_values)

1253

1254 threshold = x_axis_span * proximity_ratio

1255 levels = []

1256 prev_x = None

1257 current_level = 0

1258 for x in peak_x_values:

1259 if prev_x is not None and abs(x - prev_x) < threshold:

1260 current_level = (current_level + 1) % max_levels

1261 else:

1262 current_level = 0

1263 levels.append(current_level)

1264 prev_x = x

1265 return levels

1266

1267

1268def split_dataframe_by_column(df: Any, split_col: str) -> dict[str, Any]:

1269 """区分列の値ごとに DataFrame を分ける(long 形式のデータ)。

1270

1271 グループはファイル中の初出順(測定順などの並びに意味があることが多い)。

1272 区分列が NaN の行は除き、その数を n_dropped で返す。各グループの index は 0 から振り直す。

1273 """

1274 if split_col not in df.columns:

1275 raise ValueError(f"列 '{split_col}' がデータに存在しません。")

1276

1277 values = df[split_col]

1278 keep = values.notna()

1279 n_dropped = int((~keep).sum())

1280 filtered = df[keep]

1281

1282 groups = []

1283 for label, sub_df in filtered.groupby(split_col, sort=False, observed=True):

1284 groups.append((_format_group_label(label), sub_df.reset_index(drop=True)))

1285

1286 return {'groups': groups, 'n_dropped': n_dropped, 'n_groups': len(groups)}

1287

1288

1289def _format_group_label(value: Any) -> str:

1290 """numpy のスカラは str() すると 'np.float64(1.5)' になるので、Python の型にしてから文字列にする。"""

1291 if hasattr(value, 'item'): 1291 ↛ 1292line 1291 didn't jump to line 1292 because the condition on line 1291 was never true

1292 try:

1293 value = value.item()

1294 except (ValueError, TypeError):

1295 pass

1296 if isinstance(value, float):

1297 return f"{value:g}"

1298 return str(value)