#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ lsq_func_error_argparse.py Generic nonlinear least-squares fitting program. Features -------- - Model function is given from command line as a string, e.g. --func "p[0]*exp(-((x[0]-p[1])/p[2])**2)" - Multiple independent variables are supported through --xlabels. - Parameters can be optimized with scipy.optimize.minimize. - Parameter standard errors are estimated from numerical Jacobian. - Confidence / uncertainty bands are calculated by the delta method. - Results are saved to Excel files. Notes ----- This script is intended for local research use. The expression evaluator uses eval(), but with restricted globals. Do not pass untrusted expressions. """ import sys import os import argparse import types import traceback from typing import List, Tuple, Dict, Any, Optional import numpy as np import pandas as pd import matplotlib.pyplot as plt from scipy.optimize import minimize from scipy.stats import t as student_t # ---------------------------------------------------------------------- # Optional tkProg/tklib helpers # ---------------------------------------------------------------------- try: from tklib.tkapplication import tkApplication except Exception: tkApplication = None try: from tklib.tkvariousdata import tkVariousData except Exception: tkVariousData = None # ---------------------------------------------------------------------- # Small utilities # ---------------------------------------------------------------------- def pause_if_needed(pause: int) -> None: if int(pause): input("\nPress ENTER to terminate >> ") def terminate(message: str, pause: int = 1, code: int = 1) -> None: if message: print(message) pause_if_needed(pause) sys.exit(code) def replace_path(path: str, suffix: str) -> str: """Return dirname/filebody + suffix.""" d = os.path.dirname(path) b = os.path.splitext(os.path.basename(path))[0] return os.path.join(d, b + suffix) def parse_float_list(text: str, name: str) -> List[float]: try: return [float(s.strip()) for s in text.split(',') if s.strip() != ''] except Exception as exc: raise ValueError(f"Could not parse {name} as comma-separated floats: {text}") from exc def parse_str_list(text: str) -> List[str]: return [s.strip() for s in text.split(',') if s.strip() != ''] def parse_ranges(text: str, n: int, default_min: float = -1.0e100, default_max: float = 1.0e100) -> List[Tuple[float, float]]: """ Parse ranges such as '-1:1, 0:100'. If text is empty or '*', use full range for all variables. """ if text is None or str(text).strip() in ('', '*'): return [(default_min, default_max) for _ in range(n)] items = parse_str_list(text) if len(items) == 1 and n > 1: items = items * n if len(items) != n: raise ValueError(f"Number of ranges must be 1 or nvar={n}: {text}") out = [] for item in items: if ':' not in item: raise ValueError(f"Range must be xmin:xmax format: {item}") a, b = item.split(':', 1) xmin = default_min if a.strip() in ('', '*') else float(a) xmax = default_max if b.strip() in ('', '*') else float(b) out.append((xmin, xmax)) return out def parse_bounds(text: str, n: int) -> Optional[List[Tuple[Optional[float], Optional[float]]]]: """ Parse optimizer bounds such as '0:*, *:10, -5:5'. Use '*' or empty for None. """ if text is None or str(text).strip() in ('', '*', 'none', 'None'): return None items = parse_str_list(text) if len(items) == 1 and n > 1: items = items * n if len(items) != n: raise ValueError(f"Number of bounds must be 1 or nparam={n}: {text}") out = [] for item in items: if ':' not in item: raise ValueError(f"Bound must be lower:upper format: {item}") a, b = item.split(':', 1) lo = None if a.strip() in ('', '*') else float(a) hi = None if b.strip() in ('', '*') else float(b) out.append((lo, hi)) return out def parse_fixed(text: str, n: int) -> np.ndarray: """ Parse fixed parameter indices. Example: --fix 0,2 fixes p[0] and p[2]. """ fixed = np.zeros(n, dtype=bool) if text is None or str(text).strip() in ('', '*', 'none', 'None'): return fixed for s in parse_str_list(text): idx = int(s) if idx < 0 or idx >= n: raise ValueError(f"Fixed parameter index out of range: {idx}") fixed[idx] = True return fixed # ---------------------------------------------------------------------- # Argument parser # ---------------------------------------------------------------------- def initialize() -> types.SimpleNamespace: parser = argparse.ArgumentParser( description='Generic nonlinear least-squares fitting with parameter errors and confidence bands' ) # Basic IO and model parser.add_argument('--mode', type=str, default='fit', choices=['fit', 'sim'], help='fit or sim') parser.add_argument('--infile', type=str, default='peak.xlsx', help='Input Excel/CSV/text file') parser.add_argument('--func', type=str, default='p[0]*exp(-((x[0]-p[1])/p[2])**2)', help='Model function. Use x[0], x[1], ... and p[0], p[1], ...') parser.add_argument('--p0', type=str, default='1.3,0.6,0.1', help='Initial parameters, comma separated') parser.add_argument('--xlabels', type=str, default='0', help='Independent variable labels or column indices, comma separated') parser.add_argument('--ylabel', type=str, default='1', help='Dependent variable label or column index') parser.add_argument('--fit_range', type=str, default='*', help='Fit range for each x as xmin:xmax, comma separated. Example: -1:1,0:10') # Optimizer parser.add_argument('--method', type=str, default='BFGS', help='scipy.optimize.minimize method') parser.add_argument('--maxiter', type=int, default=1000, help='Maximum optimizer iterations') parser.add_argument('--tol', type=float, default=1.0e-8, help='Optimizer tolerance') parser.add_argument('--bounds', type=str, default='*', help='Bounds for parameters as lower:upper, comma separated. Use with L-BFGS-B/SLSQP/etc.') parser.add_argument('--fix', type=str, default='', help='Fixed parameter indices, comma separated. Example: 0,2') # Numerical derivatives parser.add_argument('--jac_absstep', type=float, default=1.0e-8, help='Absolute step for numerical derivative') parser.add_argument('--jac_relstep', type=float, default=1.0e-5, help='Relative step for numerical derivative') parser.add_argument('--use_jac', type=int, default=1, choices=[0, 1], help='Use numerical Jacobian for optimizer, 0/1') # Calculation grid and output parser.add_argument('--xcalmin', type=str, default='*', help='Calculation x min for 1D plot/grid') parser.add_argument('--xcalmax', type=str, default='*', help='Calculation x max for 1D plot/grid') parser.add_argument('--ncal', type=int, default=201, help='Number of calculation grid points for 1D plot') parser.add_argument('--out_prefix', type=str, default='', help='Output prefix. Default: input file body') parser.add_argument('--xlsm_template', type=str, default='', help='Reserved for compatibility') # Uncertainty and plot options, Launcher friendly 0/1 flags parser.add_argument('--confidence', type=float, default=0.682689492, help='Confidence level. 0.6827 is about 1 sigma') parser.add_argument('--plot', type=int, default=1, choices=[0, 1], help='Show plot, 0/1') parser.add_argument('--nplot_interval', type=int, default=10, help='Update live fitting plot every N iterations. 0 disables live update') parser.add_argument('--plot_pause', type=float, default=0.01, help='Pause time for live plot update') parser.add_argument('--plot_ci', type=int, default=1, choices=[0, 1], help='Plot uncertainty bands, 0/1') parser.add_argument('--plot_sigma_param', type=int, default=1, choices=[0, 1], help='Plot parameter uncertainty band, 0/1') parser.add_argument('--plot_sigma_pred', type=int, default=1, choices=[0, 1], help='Plot prediction uncertainty band, 0/1') parser.add_argument('--plot_sigma_combined', type=int, default=0, choices=[0, 1], help='Plot param + measured-y scatter band, 0/1') parser.add_argument('--figsize', type=float, nargs=2, default=[10.0, 5.0], help='Figure size') parser.add_argument('--fontsize', type=int, default=14, help='Font size') parser.add_argument('--fontsize_legend', type=int, default=10, help='Legend font size') parser.add_argument('--pause', type=int, default=1, choices=[0, 1], help='Pause before termination, 0/1') cfg = parser.parse_args() cfg.p0 = np.asarray(parse_float_list(cfg.p0, 'p0'), dtype=float) cfg.xlabels_list = parse_str_list(cfg.xlabels) cfg.nparam = len(cfg.p0) cfg.fixed_mask = parse_fixed(cfg.fix, cfg.nparam) cfg.active_mask = ~cfg.fixed_mask cfg.bounds_full = parse_bounds(cfg.bounds, cfg.nparam) if np.sum(cfg.active_mask) == 0: terminate('Error: all parameters are fixed.', cfg.pause) if cfg.out_prefix.strip() == '': cfg.out_prefix = os.path.splitext(cfg.infile)[0] cfg.output_fitting_path = cfg.out_prefix + '-fit.xlsx' cfg.output_parameter_path = cfg.out_prefix + '-parameters.xlsx' cfg.output_convergence_path = cfg.out_prefix + '-convergence.xlsx' return cfg # ---------------------------------------------------------------------- # Data loading # ---------------------------------------------------------------------- def _column_by_label_or_index(df: pd.DataFrame, key: str) -> Tuple[str, np.ndarray]: key = str(key).strip() if key in df.columns: label = key else: try: idx = int(key) except Exception as exc: raise KeyError(f"Column [{key}] was not found as label or integer index") from exc if idx < 0 or idx >= len(df.columns): raise IndexError(f"Column index out of range: {idx}") label = df.columns[idx] return str(label), pd.to_numeric(df[label], errors='coerce').to_numpy(dtype=float) def read_table(path: str) -> pd.DataFrame: if not os.path.exists(path): terminate(f"Error: input file does not exist: {path}", 1) ext = os.path.splitext(path)[1].lower() if ext in ['.xlsx', '.xlsm', '.xls']: return pd.read_excel(path, engine='openpyxl') if ext in ['.csv']: return pd.read_csv(path) # Whitespace or tab separated text fallback try: return pd.read_csv(path, sep=None, engine='python') except Exception: return pd.read_csv(path, delim_whitespace=True) def load_data(cfg: types.SimpleNamespace) -> Dict[str, Any]: df = read_table(cfg.infile) xlabels = [] xlist = [] for key in cfg.xlabels_list: label, arr = _column_by_label_or_index(df, key) xlabels.append(label) xlist.append(arr) ylabel, y = _column_by_label_or_index(df, cfg.ylabel) Xraw = np.vstack(xlist) ranges = parse_ranges(cfg.fit_range, len(xlabels)) mask = np.isfinite(y) for j in range(len(xlabels)): xmin, xmax = ranges[j] mask &= np.isfinite(Xraw[j]) mask &= (xmin <= Xraw[j]) & (Xraw[j] <= xmax) X = Xraw[:, mask] yfit = y[mask] if len(yfit) == 0: terminate('Error: no valid data points in the specified fit range.', cfg.pause) return { 'df': df, 'xlabels': xlabels, 'ylabel': ylabel, 'x': X, 'y': yfit, 'ranges': ranges, } # ---------------------------------------------------------------------- # Safe-ish expression evaluation # ---------------------------------------------------------------------- def build_eval_env() -> Dict[str, Any]: allowed = { '__builtins__': {}, 'np': np, 'numpy': np, 'pi': np.pi, 'e': np.e, 'exp': np.exp, 'log': np.log, 'log10': np.log10, 'sqrt': np.sqrt, 'abs': np.abs, 'sin': np.sin, 'cos': np.cos, 'tan': np.tan, 'arcsin': np.arcsin, 'arccos': np.arccos, 'arctan': np.arctan, 'sinh': np.sinh, 'cosh': np.cosh, 'tanh': np.tanh, 'minimum': np.minimum, 'maximum': np.maximum, 'where': np.where, 'clip': np.clip, } return allowed def cal_y_one(x_one: np.ndarray, p: np.ndarray, cfg: types.SimpleNamespace, env: Dict[str, Any]) -> float: val = eval(cfg.func, env, {'x': x_one, 'p': p}) return float(val) def cal_y_array(x: np.ndarray, p: np.ndarray, cfg: types.SimpleNamespace, env: Dict[str, Any]) -> np.ndarray: """Evaluate model for x with shape (nvar, ndata).""" y = np.empty(x.shape[1], dtype=float) for i in range(x.shape[1]): y[i] = cal_y_one(x[:, i], p, cfg, env) return y # ---------------------------------------------------------------------- # Parameter packing/unpacking for fixed parameters # ---------------------------------------------------------------------- def pack_active(p_full: np.ndarray, active_mask: np.ndarray) -> np.ndarray: return np.asarray(p_full[active_mask], dtype=float) def unpack_active(q: np.ndarray, p_template: np.ndarray, active_mask: np.ndarray) -> np.ndarray: p = np.asarray(p_template, dtype=float).copy() p[active_mask] = q return p def active_bounds(bounds_full, active_mask): if bounds_full is None: return None return [b for b, active in zip(bounds_full, active_mask) if active] # ---------------------------------------------------------------------- # Objective and derivatives # ---------------------------------------------------------------------- def residual_vector(p: np.ndarray, x: np.ndarray, y: np.ndarray, cfg: types.SimpleNamespace, env: Dict[str, Any]) -> np.ndarray: return y - cal_y_array(x, p, cfg, env) def objective_q(q: np.ndarray, p_template: np.ndarray, x: np.ndarray, y: np.ndarray, cfg: types.SimpleNamespace, env: Dict[str, Any], history: Optional[Dict[str, list]] = None) -> float: p = unpack_active(q, p_template, cfg.active_mask) r = residual_vector(p, x, y, cfg, env) rss = float(r @ r) mse = rss / len(y) if history is not None: history['last_p'] = p.copy() history['last_mse'] = mse if not np.isfinite(mse): return 1.0e300 return mse def step_for_param(value: float, cfg: types.SimpleNamespace) -> float: return cfg.jac_absstep + cfg.jac_relstep * max(abs(float(value)), 1.0) def gradient_q(q: np.ndarray, p_template: np.ndarray, x: np.ndarray, y: np.ndarray, cfg: types.SimpleNamespace, env: Dict[str, Any]) -> np.ndarray: g = np.empty_like(q, dtype=float) for k in range(len(q)): h = step_for_param(q[k], cfg) qm = q.copy() qp = q.copy() qm[k] -= h qp[k] += h fm = objective_q(qm, p_template, x, y, cfg, env) fp = objective_q(qp, p_template, x, y, cfg, env) g[k] = (fp - fm) / (2.0 * h) return g def numerical_model_jacobian(p: np.ndarray, x: np.ndarray, cfg: types.SimpleNamespace, env: Dict[str, Any]) -> np.ndarray: """ Jacobian of model y_model with respect to active parameters. Returns J with shape (ndata, n_active). """ active_indices = np.where(cfg.active_mask)[0] ndata = x.shape[1] J = np.empty((ndata, len(active_indices)), dtype=float) for col, ip in enumerate(active_indices): h = step_for_param(p[ip], cfg) pm = p.copy() pp = p.copy() pm[ip] -= h pp[ip] += h ym = cal_y_array(x, pm, cfg, env) yp = cal_y_array(x, pp, cfg, env) J[:, col] = (yp - ym) / (2.0 * h) return J def covariance_from_jacobian(J: np.ndarray, residuals: np.ndarray) -> Dict[str, Any]: """ Estimate parameter covariance from nonlinear least squares linearization. cov = sigma^2 * inv(J^T J), where residual = y - f(x,p). """ ndata, npar = J.shape rss = float(residuals @ residuals) dof = ndata - npar result = { 'rss': rss, 'ndata': ndata, 'npar_active': npar, 'dof': dof, 'sigma2_resid': np.nan, 'sigma_resid': np.nan, 'cov_active': None, 'std_active': None, 'corr_active': None, 'valid': False, 'message': '', } if dof <= 0: result['message'] = ( f'WARNING: ndata={ndata} <= n_active_parameters={npar}. ' 'Parameter errors and confidence bands are not estimated.' ) return result JTJ = J.T @ J try: JTJ_inv = np.linalg.pinv(JTJ) except Exception as exc: result['message'] = f'WARNING: could not invert J^T J: {exc}' return result sigma2 = rss / dof cov = sigma2 * JTJ_inv diag = np.diag(cov) if np.any(diag < 0): result['message'] = 'WARNING: covariance has negative diagonal elements; errors are not reliable.' return result std = np.sqrt(diag) denom = np.outer(std, std) with np.errstate(divide='ignore', invalid='ignore'): corr = cov / denom result.update({ 'sigma2_resid': sigma2, 'sigma_resid': np.sqrt(sigma2), 'cov_active': cov, 'std_active': std, 'corr_active': corr, 'valid': True, 'message': 'OK', }) return result def expand_active_vector(v_active: Optional[np.ndarray], cfg: types.SimpleNamespace, fill: float = np.nan) -> np.ndarray: v = np.full(cfg.nparam, fill, dtype=float) if v_active is not None: v[cfg.active_mask] = v_active v[cfg.fixed_mask] = 0.0 return v def expand_active_matrix(m_active: Optional[np.ndarray], cfg: types.SimpleNamespace, fill: float = np.nan) -> np.ndarray: m = np.full((cfg.nparam, cfg.nparam), fill, dtype=float) if m_active is not None: idx = np.where(cfg.active_mask)[0] for ia, i in enumerate(idx): for ja, j in enumerate(idx): m[i, j] = m_active[ia, ja] for i in np.where(cfg.fixed_mask)[0]: m[i, i] = 0.0 return m # ---------------------------------------------------------------------- # Bands # ---------------------------------------------------------------------- def compute_uncertainty_bands(xcal: np.ndarray, p: np.ndarray, cfg: types.SimpleNamespace, env: Dict[str, Any], cov_active: np.ndarray, sigma2_resid: float, confidence: float) -> Dict[str, np.ndarray]: y_mean = cal_y_array(xcal, p, cfg, env) Jcal = numerical_model_jacobian(p, xcal, cfg, env) var_param = np.sum((Jcal @ cov_active) * Jcal, axis=1) var_param = np.maximum(var_param, 0.0) sigma_param = np.sqrt(var_param) sigma_pred = np.sqrt(var_param + sigma2_resid) dof = max(1, xcal.shape[1] - Jcal.shape[1]) # central two-sided interval multiplier tval = student_t.ppf(0.5 + confidence / 2.0, dof) return { 'y_mean': y_mean, 'sigma_param': sigma_param, 'sigma_pred': sigma_pred, 'ci_param': tval * sigma_param, 'ci_pred': tval * sigma_pred, 'tval': np.full_like(y_mean, tval, dtype=float), } def make_calculation_grid(x: np.ndarray, cfg: types.SimpleNamespace) -> np.ndarray: nvar, ndata = x.shape if nvar == 1: if cfg.xcalmin in ('*', '', None): xmin = float(np.min(x[0])) else: xmin = float(cfg.xcalmin) if cfg.xcalmax in ('*', '', None): xmax = float(np.max(x[0])) else: xmax = float(cfg.xcalmax) return np.asarray([np.linspace(xmin, xmax, cfg.ncal)], dtype=float) # For multi-dimensional models, a 1D line is ambiguous. # Use the original x rows for uncertainty output. return x.copy() # ---------------------------------------------------------------------- # Saving and plotting # ---------------------------------------------------------------------- def save_outputs(cfg, data, p_opt, p_std, cov_full, corr_full, y_initial, y_fit, history, cov_info, bands, xcal): x = data['x'] y = data['y'] xlabels = data['xlabels'] ylabel = data['ylabel'] print() print(f"Save fitting results to [{cfg.output_fitting_path}]") fit_dict = {} for j, label in enumerate(xlabels): fit_dict[label] = x[j] fit_dict[ylabel] = y fit_dict[f'{ylabel}(initial)'] = y_initial fit_dict[f'{ylabel}(fit)'] = y_fit fit_dict[f'{ylabel}(residual)'] = y - y_fit df_fit = pd.DataFrame(fit_dict) cal_dict = {} for j, label in enumerate(xlabels): cal_dict[f'{label}(cal)'] = xcal[j] if bands is not None: cal_dict[f'{ylabel}(mean)'] = bands['y_mean'] cal_dict[f'{ylabel}(sigma_param)'] = bands['sigma_param'] cal_dict[f'{ylabel}(sigma_pred)'] = bands['sigma_pred'] cal_dict[f'{ylabel}(ci_param)'] = bands['ci_param'] cal_dict[f'{ylabel}(ci_pred)'] = bands['ci_pred'] cal_dict['t_multiplier'] = bands['tval'] df_cal = pd.DataFrame(cal_dict) df_conv = pd.DataFrame(history) param_rows = [] for i in range(cfg.nparam): param_rows.append({ 'index': i, 'initial': cfg.p0[i], 'final': p_opt[i], 'std': p_std[i], 'fixed': int(cfg.fixed_mask[i]), }) df_param = pd.DataFrame(param_rows) df_cov = pd.DataFrame(cov_full, columns=[f'p{i}' for i in range(cfg.nparam)], index=[f'p{i}' for i in range(cfg.nparam)]) df_corr = pd.DataFrame(corr_full, columns=[f'p{i}' for i in range(cfg.nparam)], index=[f'p{i}' for i in range(cfg.nparam)]) df_summary = pd.DataFrame([ {'key': 'func', 'value': cfg.func}, {'key': 'method', 'value': cfg.method}, {'key': 'ndata', 'value': cov_info['ndata']}, {'key': 'n_active_parameters', 'value': cov_info['npar_active']}, {'key': 'dof', 'value': cov_info['dof']}, {'key': 'RSS', 'value': cov_info['rss']}, {'key': 'MSE', 'value': cov_info['rss'] / max(cov_info['ndata'], 1)}, {'key': 'sigma2_resid', 'value': cov_info['sigma2_resid']}, {'key': 'sigma_resid', 'value': cov_info['sigma_resid']}, {'key': 'confidence', 'value': cfg.confidence}, {'key': 'covariance_status', 'value': cov_info['message']}, ]) with pd.ExcelWriter(cfg.output_fitting_path, engine='openpyxl') as writer: df_fit.to_excel(writer, sheet_name='fit_data', index=False) df_cal.to_excel(writer, sheet_name='calculation_grid', index=False) df_conv.to_excel(writer, sheet_name='convergence', index=False) df_summary.to_excel(writer, sheet_name='summary', index=False) print(f"Save parameter results to [{cfg.output_parameter_path}]") with pd.ExcelWriter(cfg.output_parameter_path, engine='openpyxl') as writer: df_param.to_excel(writer, sheet_name='parameters', index=False) df_cov.to_excel(writer, sheet_name='covariance') df_corr.to_excel(writer, sheet_name='correlation') df_summary.to_excel(writer, sheet_name='summary', index=False) print(f"Save convergence history to [{cfg.output_convergence_path}]") df_conv.to_excel(cfg.output_convergence_path, index=False) class LiveFitPlotter: """Live plotter for monitoring nonlinear fitting progress.""" def __init__(self, cfg, data, y_initial): self.cfg = cfg self.data = data self.y_initial = y_initial self.enabled = bool(cfg.plot) and int(cfg.nplot_interval) > 0 self.fig = None self.ax_fit = None self.ax_conv = None self.line_current = None self.line_conv = None self.nvar = data['x'].shape[0] self.xplot = None if not self.enabled: return # Interactive mode keeps the window responsive during scipy.optimize.minimize(). plt.ion() self.fig, (self.ax_fit, self.ax_conv) = plt.subplots(1, 2, figsize=cfg.figsize) x = data['x'] y = data['y'] ylabel = data['ylabel'] xlabels = data['xlabels'] if self.nvar == 1: self.xplot = x[0] self.ax_fit.plot(x[0], y, 'o', label='data', markersize=3) self.ax_fit.plot(x[0], y_initial, '--', label='initial', linewidth=0.8) (self.line_current,) = self.ax_fit.plot(x[0], y_initial, '-', label='current', linewidth=1.0) self.ax_fit.set_xlabel(xlabels[0], fontsize=cfg.fontsize) else: self.xplot = np.arange(len(y)) self.ax_fit.plot(self.xplot, y, 'o', label='data', markersize=3) self.ax_fit.plot(self.xplot, y_initial, '--', label='initial', linewidth=0.8) (self.line_current,) = self.ax_fit.plot(self.xplot, y_initial, '-', label='current', linewidth=1.0) self.ax_fit.set_xlabel('index', fontsize=cfg.fontsize) self.ax_fit.set_ylabel(ylabel, fontsize=cfg.fontsize) self.ax_fit.tick_params(labelsize=cfg.fontsize) self.ax_fit.legend(fontsize=cfg.fontsize_legend) (self.line_conv,) = self.ax_conv.plot([], [], 'o-', linewidth=0.8, markersize=3) self.ax_conv.set_xlabel('# of iteration', fontsize=cfg.fontsize) self.ax_conv.set_ylabel('MSE', fontsize=cfg.fontsize) self.ax_conv.set_yscale('log') self.ax_conv.tick_params(labelsize=cfg.fontsize) self.fig.suptitle('Fitting progress') self.fig.tight_layout() self.fig.canvas.draw_idle() plt.pause(cfg.plot_pause) def update(self, iteration, p_current, mse_history, env): if not self.enabled: return if iteration % int(self.cfg.nplot_interval) != 0: return if self.fig is None or not plt.fignum_exists(self.fig.number): self.enabled = False return x = self.data['x'] y_current = cal_y_array(x, p_current, self.cfg, env) self.line_current.set_ydata(y_current) self.ax_fit.relim() self.ax_fit.autoscale_view() iters = np.asarray(mse_history['iter'], dtype=float) mses = np.asarray(mse_history['MSE'], dtype=float) mask = np.isfinite(mses) & (mses > 0) self.line_conv.set_data(iters, mses) if len(iters) > 0: self.ax_conv.set_xlim(0, max(1.0, float(np.max(iters)) + 1.0)) if np.any(mask): ymin = float(np.min(mses[mask])) ymax = float(np.max(mses[mask])) if ymin == ymax: ymin *= 0.8 ymax *= 1.2 self.ax_conv.set_ylim(ymin * 0.8, ymax * 1.2) self.fig.canvas.draw_idle() plt.pause(self.cfg.plot_pause) def finalize(self, p_final, mse_history, env): if not self.enabled: return # Always show the final state, even if it is not exactly on the interval. final_iter = int(mse_history['iter'][-1]) if len(mse_history['iter']) else 0 old_interval = self.cfg.nplot_interval self.cfg.nplot_interval = 1 self.update(final_iter, p_final, mse_history, env) self.cfg.nplot_interval = old_interval def plot_results(cfg, data, p_opt, y_initial, y_fit, p_std, bands, xcal): if not cfg.plot: return x = data['x'] y = data['y'] xlabels = data['xlabels'] ylabel = data['ylabel'] nvar = x.shape[0] if cfg.plot_ci: fig, (ax1, ax2) = plt.subplots(1, 2, figsize=cfg.figsize) else: fig, ax1 = plt.subplots(1, 1, figsize=cfg.figsize) ax2 = None if nvar == 1: ax1.plot(x[0], y, 'o', label='data', markersize=3) ax1.plot(x[0], y_initial, '--', label='initial', linewidth=0.8) if bands is not None: ax1.plot(xcal[0], bands['y_mean'], '-', label='fit', linewidth=1.0) if cfg.plot_sigma_pred: ax1.fill_between(xcal[0], bands['y_mean'] - bands['ci_pred'], bands['y_mean'] + bands['ci_pred'], alpha=0.25, label=f'prediction CI ({cfg.confidence:g})') if cfg.plot_sigma_param: ax1.fill_between(xcal[0], bands['y_mean'] - bands['ci_param'], bands['y_mean'] + bands['ci_param'], alpha=0.35, label=f'parameter CI ({cfg.confidence:g})') else: ax1.plot(x[0], y_fit, '-', label='fit', linewidth=1.0) ax1.set_xlabel(xlabels[0], fontsize=cfg.fontsize) else: idx = np.arange(len(y)) ax1.plot(idx, y, 'o', label='data', markersize=3) ax1.plot(idx, y_initial, '--', label='initial', linewidth=0.8) ax1.plot(idx, y_fit, '-', label='fit', linewidth=1.0) if bands is not None and cfg.plot_sigma_param: ax1.fill_between(idx, bands['y_mean'] - bands['ci_param'], bands['y_mean'] + bands['ci_param'], alpha=0.35, label=f'parameter CI ({cfg.confidence:g})') if bands is not None and cfg.plot_sigma_pred: ax1.fill_between(idx, bands['y_mean'] - bands['ci_pred'], bands['y_mean'] + bands['ci_pred'], alpha=0.25, label=f'prediction CI ({cfg.confidence:g})') ax1.set_xlabel('index', fontsize=cfg.fontsize) ax1.set_ylabel(ylabel, fontsize=cfg.fontsize) ax1.tick_params(labelsize=cfg.fontsize) ax1.legend(fontsize=cfg.fontsize_legend) if ax2 is not None: idx = np.arange(cfg.nparam) ax2.errorbar(idx, p_opt, yerr=p_std, fmt='o', capsize=3, label='parameter ±1σ') ax2.set_xlabel('parameter index', fontsize=cfg.fontsize) ax2.set_ylabel('parameter value', fontsize=cfg.fontsize) ax2.tick_params(labelsize=cfg.fontsize) ax2.legend(fontsize=cfg.fontsize_legend) plt.tight_layout() plt.show() # ---------------------------------------------------------------------- # Main # ---------------------------------------------------------------------- def main() -> None: cfg = initialize() print('Generic nonlinear least-squares fitting') print(f' infile : {cfg.infile}') print(f' func : {cfg.func}') print(f' xlabels : {cfg.xlabels_list}') print(f' ylabel : {cfg.ylabel}') print(f' method : {cfg.method}') print(f' p0 : {cfg.p0}') print(f' fixed mask : {cfg.fixed_mask.astype(int)}') try: data = load_data(cfg) env = build_eval_env() x = data['x'] y = data['y'] print(f' ndata : {len(y)}') print(f' fit ranges : {data["ranges"]}') y_initial = cal_y_array(x, cfg.p0, cfg, env) if cfg.mode == 'sim': print('Simulation mode: no fitting is performed.') history = {'iter': [0], 'MSE': [float(np.mean((y - y_initial)**2))]} p_std = np.full(cfg.nparam, np.nan) cov_full = np.full((cfg.nparam, cfg.nparam), np.nan) corr_full = np.full((cfg.nparam, cfg.nparam), np.nan) cov_info = {'ndata': len(y), 'npar_active': int(np.sum(cfg.active_mask)), 'dof': np.nan, 'rss': float(np.sum((y - y_initial)**2)), 'sigma2_resid': np.nan, 'sigma_resid': np.nan, 'message': 'simulation mode'} save_outputs(cfg, data, cfg.p0, p_std, cov_full, corr_full, y_initial, y_initial, history, cov_info, None, x) plot_results(cfg, data, cfg.p0, y_initial, y_initial, p_std, None, x) pause_if_needed(cfg.pause) return q0 = pack_active(cfg.p0, cfg.active_mask) b_active = active_bounds(cfg.bounds_full, cfg.active_mask) history = {'iter': [], 'MSE': []} hist_state = {'iter': 0, 'last_p': cfg.p0.copy(), 'last_mse': np.nan} live_plotter = LiveFitPlotter(cfg, data, y_initial) grad_methods = {'CG', 'BFGS', 'L-BFGS-B', 'TNC', 'SLSQP', 'trust-constr'} no_grad_methods = {'Nelder-Mead', 'Powell', 'COBYLA'} method_upper_name = cfg.method.strip() use_jac = bool(cfg.use_jac) and (method_upper_name not in no_grad_methods) def callback(q): mse = objective_q(q, cfg.p0, x, y, cfg, env, hist_state) hist_state['iter'] += 1 history['iter'].append(hist_state['iter']) history['MSE'].append(mse) print(f"callback {hist_state['iter']:5d}: MSE={mse:.8e} p={hist_state['last_p']}") live_plotter.update(hist_state['iter'], hist_state['last_p'], history, env) print() print('Minimize:') kwargs = dict( fun=lambda q: objective_q(q, cfg.p0, x, y, cfg, env, hist_state), x0=q0, method=cfg.method, callback=callback, tol=cfg.tol, options={'maxiter': cfg.maxiter, 'disp': True}, ) if b_active is not None: kwargs['bounds'] = b_active if use_jac: kwargs['jac'] = lambda q: gradient_q(q, cfg.p0, x, y, cfg, env) res = minimize(**kwargs) p_opt = unpack_active(res.x, cfg.p0, cfg.active_mask) y_fit = cal_y_array(x, p_opt, cfg, env) residuals = y - y_fit mse_final = float(np.mean(residuals**2)) print() print('Optimization result:') print(f' success : {res.success}') print(f' message : {res.message}') print(f' nit : {getattr(res, "nit", "-")}') print(f' MSE : {mse_final:.12g}') print(f' p : {p_opt}') J = numerical_model_jacobian(p_opt, x, cfg, env) cov_info = covariance_from_jacobian(J, residuals) print() print('Uncertainty estimate:') print(f" status : {cov_info['message']}") print(f" ndata : {cov_info['ndata']}") print(f" n_active_param: {cov_info['npar_active']}") print(f" dof : {cov_info['dof']}") print(f" RSS : {cov_info['rss']:.12g}") print(f" sigma_resid : {cov_info['sigma_resid']}") p_std = expand_active_vector(cov_info['std_active'], cfg) cov_full = expand_active_matrix(cov_info['cov_active'], cfg) corr_full = expand_active_matrix(cov_info['corr_active'], cfg) print() print('Fitted parameters:') for i in range(cfg.nparam): if cfg.fixed_mask[i]: print(f' p[{i}] = {p_opt[i]:.12g} fixed') else: print(f' p[{i}] = {p_opt[i]:.12g} +- {p_std[i]:.6g}') xcal = make_calculation_grid(x, cfg) bands = None if cov_info['valid']: bands = compute_uncertainty_bands(xcal, p_opt, cfg, env, cov_info['cov_active'], cov_info['sigma2_resid'], cfg.confidence) save_outputs(cfg, data, p_opt, p_std, cov_full, corr_full, y_initial, y_fit, history, cov_info, bands, xcal) live_plotter.finalize(p_opt, history, env) plot_results(cfg, data, p_opt, y_initial, y_fit, p_std, bands, xcal) pause_if_needed(cfg.pause) except SystemExit: raise except Exception: print('\nERROR:') traceback.print_exc() pause_if_needed(cfg.pause) sys.exit(1) if __name__ == '__main__': main()