Coverage for graphica/core/fit_models.py: 99.3%

1"""曲線フィットのモデルの表。種類の判定・パラメータ名・関数・初期値の推定をここ 1 か所に持つ。 

2 

3種類は表示名(例「ガウシアン (y = …)」)で渡され、保存ファイルと方法の文章にもそのまま入る。判定の順は、 

4プラグインの関数(名前の完全一致)→ カスタム数式 → 組み込み(語の部分一致、表の順:「2成分指数」を「指数関数」より、 

5「擬似フォークト」を「フォークト」より、「ボルツマン」を「シグモイド」より先に)。プラグインを先に見ないと、 

6「二重ガウシアン」のような名前が黙って組み込みのガウシアンになる。 

7初期値の式は結果の popt が 1 ビットも変わらないよう、演算の順まで変えないこと。 

8""" 

9import re 

10from dataclasses import dataclass 

11from typing import Any, Callable, Mapping 

12 

13import numpy as np 

14from scipy.special import wofz 

15 

16from graphica.core.safe_eval import DEFAULT_FUNCTIONS, safe_eval_formula 

17 

18CUSTOM_FORMULA_KEYWORD = "カスタム数式" 

19CUSTOM_FORMULA_LABEL = "カスタム数式..." 

20CUSTOM_FORMULA_MODEL_ID = "custom_formula" 

21PLUGIN_MODEL_ID_PREFIX = "plugin:" 

22 

23 

24def is_custom_formula_type(fit_type: str) -> bool: 

25 return CUSTOM_FORMULA_KEYWORD in fit_type 

26 

27 

28def _linear(x: Any, a: float, b: float) -> Any: 

29 return a * x + b 

30 

31 

32def _poly2(x: Any, a: float, b: float, c: float) -> Any: 

33 return a * x**2 + b * x + c 

34 

35 

36def _poly3(x: Any, a: float, b: float, c: float, d: float) -> Any: 

37 return a * x**3 + b * x**2 + c * x + d 

38 

39 

40def _exponential(x: Any, a: float, b: float) -> Any: 

41 return a * np.exp(b * x) 

42 

43 

44def _logarithmic(x: Any, a: float, b: float) -> Any: 

45 return a * np.log(x) + b 

46 

47 

48def _power(x: Any, a: float, b: float) -> Any: 

49 return a * np.power(x, b) 

50 

51 

52def _gaussian(x: Any, a: float, b: float, c: float, d: float) -> Any: 

53 return a * np.exp(-((x - b) ** 2) / (2 * c ** 2)) + d 

54 

55 

56def _sigmoid(x: Any, a: float, b: float, c: float) -> Any: 

57 return a / (1 + np.exp(-b * (x - c))) 

58 

59 

60def _two_exponentials(x: Any, a1: float, b1: float, a2: float, b2: float, c: float) -> Any: 

61 return a1 * np.exp(b1 * x) + a2 * np.exp(b2 * x) + c 

62 

63 

64def _lorentzian(x: Any, a: float, b: float, c: float, d: float) -> Any: 

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

66 

67 

68def _pseudo_voigt(x: Any, a: float, b: float, c: float, eta: float, d: float) -> Any: 

69 # 共通の中心 b と FWHM c を持つローレンツ型とガウス型を eta で混ぜる 

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

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

72 return a * (eta * lorentzian_shape + (1 - eta) * gaussian_shape) + d 

73 

74 

75def _voigt(x: Any, a: float, b: float, sigma: float, gamma: float, d: float) -> Any: 

76 # Faddeeva 関数による Voigt。この正規化で gamma→0 のとき a がピーク高さになる 

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

78 return a * np.real(wofz(z)) / (sigma * np.sqrt(2 * np.pi)) + d 

79 

80 

81def _boltzmann_sigmoid(x: Any, a1: float, a2: float, x0: float, dx: float) -> Any: 

82 return a2 + (a1 - a2) / (1 + np.exp((x - x0) / dx)) 

83 

84 

85def _hill(x: Any, vmax: float, k: float, n: float) -> Any: 

86 return (vmax * np.power(x, n)) / (np.power(k, n) + np.power(x, n)) 

87 

88 

89def _estimate_fwhm(x_arr: Any, y_arr: Any, amplitude: float) -> float: 

90 """半値を横切る X の幅を FWHM の粗い推定にする(裾の広い形で誤った局所解に落ちないため)。""" 

91 half_level = np.nanmin(y_arr) + amplitude / 2 

92 above_half = x_arr[y_arr >= half_level] 

93 if len(above_half) == 0: 93 ↛ 94line 93 didn't jump to line 94 because the condition on line 93 was never true

94 return (np.nanmax(x_arr) - np.nanmin(x_arr)) / 4 or 1.0 

95 return (above_half.max() - above_half.min()) or 1.0 

96 

97 

98def _p0_linear(x: Any, y: Any) -> list[Any]: 

99 return [1.0, 0.0] 

100 

101 

102def _p0_poly2(x: Any, y: Any) -> list[Any]: 

103 return [1.0, 1.0, 0.0] 

104 

105 

106def _p0_poly3(x: Any, y: Any) -> list[Any]: 

107 return [1.0, 1.0, 1.0, 0.0] 

108 

109 

110def _p0_two_exponentials(x: Any, y: Any) -> list[Any]: 

111 amplitude = (np.nanmax(y) - np.nanmin(y)) / 2 or 1.0 

112 x_span = (np.nanmax(x) - np.nanmin(x)) or 1.0 

113 # 2成分が同じ初期値だと縮退して収束しないので、符号の違う率から始める 

114 rate0 = 2.0 / x_span 

115 return [amplitude, rate0, amplitude, -rate0, np.nanmin(y)] 

116 

117 

118def _p0_exponential(x: Any, y: Any) -> list[Any]: 

119 amplitude = np.nanmean(np.abs(y)) or 1.0 

120 return [amplitude, 0.01] 

121 

122 

123def _p0_logarithmic(x: Any, y: Any) -> list[Any]: 

124 return [1.0, 0.0] 

125 

126 

127def _p0_power(x: Any, y: Any) -> list[Any]: 

128 return [1.0, 1.0] 

129 

130 

131def _peak_amplitude_and_center(x: Any, y: Any) -> tuple[Any, Any]: 

132 amplitude = (np.nanmax(y) - np.nanmin(y)) or 1.0 

133 center = x[np.nanargmax(y)] if len(x) else 0.0 

134 return amplitude, center 

135 

136 

137def _p0_gaussian(x: Any, y: Any) -> list[Any]: 

138 amplitude, center = _peak_amplitude_and_center(x, y) 

139 width = (np.nanmax(x) - np.nanmin(x)) / 4 or 1.0 

140 return [amplitude, center, width, np.nanmin(y)] 

141 

142 

143def _p0_lorentzian(x: Any, y: Any) -> list[Any]: 

144 amplitude, center = _peak_amplitude_and_center(x, y) 

145 # c は HWHM なので FWHM の半分 

146 fwhm0 = _estimate_fwhm(x, y, amplitude) 

147 return [amplitude, center, fwhm0 / 2, np.nanmin(y)] 

148 

149 

150def _p0_pseudo_voigt(x: Any, y: Any) -> list[Any]: 

151 amplitude, center = _peak_amplitude_and_center(x, y) 

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

153 fwhm0 = _estimate_fwhm(x, y, amplitude) 

154 return [amplitude, center, fwhm0, 0.5, np.nanmin(y)] 

155 

156 

157def _p0_voigt(x: Any, y: Any) -> list[Any]: 

158 amplitude, center = _peak_amplitude_and_center(x, y) 

159 fwhm0 = _estimate_fwhm(x, y, amplitude) 

160 # FWHM をガウス成分とローレンツ成分に大まかに割り振る経験的な初期値 

161 sigma0 = (fwhm0 / 2.355) or 1.0 

162 gamma0 = (fwhm0 / 4) or 1.0 

163 # _voigt の正規化に合わせ、ピーク高さが振幅に近くなるよう a をスケールする 

164 return [amplitude * sigma0 * np.sqrt(2 * np.pi), center, sigma0, gamma0, np.nanmin(y)] 

165 

166 

167def _p0_boltzmann_sigmoid(x: Any, y: Any) -> list[Any]: 

168 order = np.argsort(x) 

169 x_sorted, y_sorted = x[order], y[order] 

170 y_start = y_sorted[0] 

171 y_end = y_sorted[-1] 

172 mid_level = (y_start + y_end) / 2 

173 # 遷移の中心は Y の中点を最初に横切る X から推定する(X の平均だと局所解に落ちやすい) 

174 crossing_mask = y_sorted < mid_level if y_start >= y_end else y_sorted > mid_level 

175 crossing_indices = np.flatnonzero(crossing_mask) 

176 x0 = x_sorted[crossing_indices[0]] if len(crossing_indices) else np.nanmean(x) 

177 dx0 = (np.nanmax(x) - np.nanmin(x)) / 10 or 1.0 

178 return [y_start, y_end, x0, dx0] 

179 

180 

181def _p0_sigmoid(x: Any, y: Any) -> list[Any]: 

182 amplitude = np.nanmax(y) or 1.0 

183 return [amplitude, 1.0, np.nanmean(x)] 

184 

185 

186def _p0_hill(x: Any, y: Any) -> list[Any]: 

187 vmax0 = np.nanmax(y) or 1.0 

188 positive_x = x[x > 0] 

189 k0 = np.nanmedian(positive_x) if len(positive_x) else 1.0 

190 return [vmax0, k0, 1.0] 

191 

192 

193def _require_positive_x(message: str) -> Callable[[Any], None]: 

194 def check(x: Any) -> None: 

195 if np.any(x <= 0): 

196 raise ValueError(message) 

197 return check 

198 

199 

200def _require_non_negative_x(x: Any) -> None: 

201 if np.any(x < 0): 

202 raise ValueError("ヒル式は X >= 0 のデータにのみ使用できます。") 

203 

204 

205@dataclass(frozen=True) 

206class FitModel: 

207 model_id: str # 保存用の ID(fit_result['fit_model_id'])。一度決めたら変えない 

208 keyword: str # 種類の名前にこの語が含まれればこのモデル 

209 menu_label: str # ダイアログの選択肢(= 保存される fit_type) 

210 param_names: tuple[str, ...] 

211 func: Callable[..., Any] 

212 initial_guess: Callable[[Any, Any], list[Any]] 

213 check_data: Callable[[Any], None] | None = None # 使えないデータなら ValueError 

214 

215 

216# 判定の順。ダイアログの並びは MENU_ORDER。 

217BUILTIN_FIT_MODELS: tuple[FitModel, ...] = ( 

218 FitModel("linear", "線形", "線形 (y = ax + b)", ("a", "b"), _linear, _p0_linear), 

219 FitModel("poly2", "2次多項式", "2次多項式 (y = ax^2 + bx + c)", ("a", "b", "c"), _poly2, _p0_poly2), 

220 FitModel("poly3", "3次多項式", "3次多項式 (y = ax^3 + bx^2 + cx + d)", ("a", "b", "c", "d"), _poly3, _p0_poly3), 

221 FitModel("two_exponentials", "2成分指数", "2成分指数関数 (y = a1*exp(b1*x) + a2*exp(b2*x) + c)", ("a1", "b1", "a2", "b2", "c"), 

222 _two_exponentials, _p0_two_exponentials), 

223 FitModel("exponential", "指数関数", "指数関数 (y = a * exp(bx))", ("a", "b"), _exponential, _p0_exponential), 

224 FitModel("logarithmic", "対数", "対数 (y = a * ln(x) + b)", ("a", "b"), _logarithmic, _p0_logarithmic, 

225 _require_positive_x("対数フィットは X > 0 のデータにのみ使用できます。")), 

226 FitModel("power", "べき乗", "べき乗 (y = a * x^b)", ("a", "b"), _power, _p0_power, 

227 _require_positive_x("べき乗フィットは X > 0 のデータにのみ使用できます。")), 

228 FitModel("gaussian", "ガウシアン", "ガウシアン (y = a * exp(-(x-b)^2 / (2c^2)) + d)", ("a", "b", "c", "d"), 

229 _gaussian, _p0_gaussian), 

230 FitModel("lorentzian", "ローレンツ", "ローレンツ関数 (y = a / (1 + ((x-b)/c)^2) + d)", ("a", "b", "c", "d"), 

231 _lorentzian, _p0_lorentzian), 

232 FitModel("pseudo_voigt", "擬似フォークト", "擬似フォークト関数 (y = a*(η/(1+((x-b)/c)^2) + (1-η)*exp(-4ln2*((x-b)/c)^2)) + d)", 

233 ("a", "b", "c", "eta", "d"), _pseudo_voigt, _p0_pseudo_voigt), 

234 FitModel("voigt", "フォークト", "フォークト関数 (y = a*Re[wofz((x-b+iγ)/(σ√2))] / (σ√(2π)) + d)", 

235 ("a", "b", "sigma", "gamma", "d"), _voigt, _p0_voigt), 

236 FitModel("boltzmann_sigmoid", "ボルツマン", "ボルツマンシグモイド (y = a2 + (a1-a2) / (1 + exp((x-x0)/dx)))", 

237 ("a1", "a2", "x0", "dx"), _boltzmann_sigmoid, _p0_boltzmann_sigmoid), 

238 FitModel("sigmoid", "シグモイド", "シグモイド (y = a / (1 + exp(-b(x-c))))", ("a", "b", "c"), _sigmoid, _p0_sigmoid), 

239 FitModel("hill", "ヒル", "ヒルの式 (y = vmax*x^n / (k^n + x^n))", ("vmax", "k", "n"), _hill, _p0_hill, 

240 _require_non_negative_x), 

241) 

242 

243MENU_ORDER: tuple[str, ...] = ( 

244 "線形", "2次多項式", "3次多項式", "指数関数", "対数", "べき乗", "ガウシアン", "ローレンツ", 

245 "擬似フォークト", "フォークト", "2成分指数", "ボルツマン", "シグモイド", "ヒル", 

246) 

247 

248# プラグインの関数名がこれらと同じだと、保存済みの組み込みの種類を取り違えるので、登録を断る 

249RESERVED_FIT_TYPE_NAMES: tuple[str, ...] = (CUSTOM_FORMULA_KEYWORD,) + tuple(m.keyword for m in BUILTIN_FIT_MODELS) 

250 

251 

252def builtin_menu_labels() -> list[str]: 

253 by_keyword = {m.keyword: m for m in BUILTIN_FIT_MODELS} 

254 return [by_keyword[keyword].menu_label for keyword in MENU_ORDER] 

255 

256 

257_RESERVED_FORMULA_NAMES = set(DEFAULT_FUNCTIONS.keys()) | {'x'} 

258 

259 

260def extract_formula_params(formula: str) -> list[str]: 

261 """x でも既知の関数名でもない識別子を、出現順にパラメータとして取り出す。""" 

262 params = [] 

263 for name in re.findall(r'[a-zA-Z_][a-zA-Z_0-9]*', formula): 

264 if name in _RESERVED_FORMULA_NAMES or name in params: 

265 continue 

266 params.append(name) 

267 if not params: 

268 raise ValueError("数式にフィットパラメータ(x以外の文字)が見つかりません。") 

269 return params 

270 

271 

272def build_custom_fit_func(formula: str, param_names: list[str]) -> Callable[..., Any]: 

273 def custom_func(x: Any, *params: float) -> Any: 

274 variables = {'x': x} 

275 variables.update(zip(param_names, params)) 

276 try: 

277 return safe_eval_formula(formula, variables) 

278 except Exception as e: 

279 raise ValueError(f"数式の評価に失敗しました: {e}") from e 

280 return custom_func 

281 

282 

283@dataclass 

284class ResolvedFitModel: 

285 """種類の名前から決まったモデル。param_names はプラグインのときは登録されたリストそのもの。""" 

286 param_names: list[str] 

287 func: Callable[..., Any] 

288 initial_guess: Callable[[Any, Any], list[Any]] 

289 check_data: Callable[[Any], None] | None = None 

290 

291 

292def _custom_formula_or_raise(custom_formula: str | None) -> str: 

293 if not custom_formula or not custom_formula.strip(): 

294 raise ValueError("カスタム数式が入力されていません。") 

295 return custom_formula 

296 

297 

298_BUILTIN_BY_ID = {m.model_id: m for m in BUILTIN_FIT_MODELS} 

299 

300 

301def _identify(fit_type: str, model_id: str | None, 

302 plugin_functions: Mapping[str, dict[str, Any]]) -> tuple[str, Any] | None: 

303 """("custom", None) / ("builtin", FitModel) / ("plugin", 名前)。ID が分かればそれを先に使う(表示名が変わっても同じモデル)。""" 

304 if model_id: 

305 if model_id == CUSTOM_FORMULA_MODEL_ID: 

306 return ("custom", None) 

307 if model_id in _BUILTIN_BY_ID: 

308 return ("builtin", _BUILTIN_BY_ID[model_id]) 

309 plugin_name = model_id[len(PLUGIN_MODEL_ID_PREFIX):] 

310 if model_id.startswith(PLUGIN_MODEL_ID_PREFIX) and plugin_name in plugin_functions: 

311 return ("plugin", plugin_name) 

312 if fit_type in plugin_functions: 

313 return ("plugin", fit_type) 

314 if is_custom_formula_type(fit_type): 

315 return ("custom", None) 

316 for model in BUILTIN_FIT_MODELS: 

317 if model.keyword in fit_type: 

318 return ("builtin", model) 

319 return None 

320 

321 

322def fit_model_id(fit_type: str, plugin_functions: Mapping[str, dict[str, Any]]) -> str | None: 

323 found = _identify(fit_type, None, plugin_functions) 

324 if found is None: 

325 return None 

326 kind, target = found 

327 if kind == "custom": 

328 return CUSTOM_FORMULA_MODEL_ID 

329 if kind == "builtin": 

330 return str(target.model_id) 

331 return PLUGIN_MODEL_ID_PREFIX + target 

332 

333 

334def fit_type_label(fit_result: Mapping[str, Any], default: Any = None) -> Any: 

335 """保存したフィット結果の種類の表示名。組み込みは ID から今の表示名を引く(表示名を変えても古いファイルとそろう)。""" 

336 model = _BUILTIN_BY_ID.get(fit_result.get('fit_model_id') or "") 

337 return model.menu_label if model is not None else fit_result.get('fit_type', default) 

338 

339 

340def resolve_fit_param_names(fit_type: str, custom_formula: str | None, 

341 plugin_functions: Mapping[str, dict[str, Any]], model_id: str | None = None) -> list[str]: 

342 """パラメータ名だけを決める(カスタム数式の関数は作らない)。""" 

343 found = _identify(fit_type, model_id, plugin_functions) 

344 if found is None: 

345 raise ValueError(f"不明なフィットタイプ: {fit_type}") 

346 kind, target = found 

347 if kind == "custom": 

348 return extract_formula_params(_custom_formula_or_raise(custom_formula)) 

349 if kind == "builtin": 

350 return list(target.param_names) 

351 return list(plugin_functions[target]["params"]) 

352 

353 

354def resolve_fit_model(fit_type: str, custom_formula: str | None, 

355 plugin_functions: Mapping[str, dict[str, Any]], model_id: str | None = None) -> ResolvedFitModel: 

356 found = _identify(fit_type, model_id, plugin_functions) 

357 if found is None: 

358 raise ValueError(f"不明なフィットタイプ: {fit_type}") 

359 kind, target = found 

360 if kind == "custom": 

361 params = extract_formula_params(_custom_formula_or_raise(custom_formula)) 

362 # *params 形式ではパラメータ数を推定できないので p0 で数を伝える 

363 return ResolvedFitModel(params, build_custom_fit_func(custom_formula or "", params), 

364 lambda x, y: [1.0] * len(params)) 

365 if kind == "builtin": 

366 return ResolvedFitModel(list(target.param_names), target.func, target.initial_guess, target.check_data) 

367 entry = plugin_functions[target] 

368 return ResolvedFitModel(entry["params"], entry["func"], _plugin_initial_guess(entry)) 

369 

370 

371def _plugin_initial_guess(entry: dict[str, Any]) -> Callable[[Any, Any], list[Any]]: 

372 def initial_guess(x: Any, y: Any) -> list[Any]: 

373 plugin_p0 = entry["p0"] 

374 if callable(plugin_p0): 

375 return list(plugin_p0(x, y)) 

376 if plugin_p0 is not None: 

377 return list(plugin_p0) 

378 return [1.0] * len(entry["params"]) 

379 return initial_guess