#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ adaptive_gaussian_ridge_practical.py Adaptive variable-width Gaussian basis regression with Chebyshev background, weighted Ridge regularization, and optional pruning. Main features ------------- - Read x-y data from Excel with --infile. - Use built-in demo data when --infile demo. - Fit y = background + sum_i w_i Gaussian_i. - Background and Gaussian coefficients can have different Ridge penalties. - Wide Gaussian bases can be penalized more strongly. - Optional pruning removes weak Gaussian bases step by step. - Export calculation results to Excel for plotting: * input data * fitted total spectrum * background spectrum * Gaussian sum spectrum * basis parameters * pruning history * settings/summary Examples -------- python adaptive_gaussian_ridge.py --infile demo --mode prune --n-centers 60 --bg-order 3 python adaptive_gaussian_ridge.py --infile data.xlsx --xcol Energy --ycol Intensity --mode prune python adaptive_gaussian_ridge.py --infile data.xlsx --sheet Sheet1 --xcol 0 --ycol 1 --outfile result.xlsx """ import argparse import os import sys import traceback from dataclasses import dataclass import numpy as np import pandas as pd import matplotlib.pyplot as plt from scipy.signal import savgol_filter from numpy.polynomial.chebyshev import chebvander sigma_slope_eps = 1.0e-6 # ============================================================================= # Demo data # ============================================================================= def true_background(x): return 0.10 + 0.08 * x - 0.025 * x**2 def true_peaks(x): return ( 1.2 * np.exp(-0.5 * ((x - 1.0) / 0.25) ** 2) - 0.8 * np.exp(-0.5 * ((x + 1.5) / 0.15) ** 2) + 0.25 * np.exp(-0.5 * ((x - 2.0) / 0.45) ** 2) ) def true_function(x): return true_background(x) + true_peaks(x) def generate_demo_data(n=180, noise=0.06, seed=1): rng = np.random.default_rng(seed) x = np.linspace(-3.0, 3.0, n) y_true = true_function(x) y = y_true + rng.normal(0.0, noise, size=n) return x, y, y_true # ============================================================================= # I/O helpers # ============================================================================= def safe_import_message(exc, package_hint): print("ERROR: required package import/use failed.") traceback.print_exc() print(f"\nPlease install required package, e.g.\n pip install {package_hint}") input("\nPress ENTER to terminate>>\n") sys.exit(1) def parse_column_selector(selector, columns): """Return column name from a selector that may be a name or zero-based index.""" if selector is None: return None if isinstance(selector, int): return columns[selector] s = str(selector).strip() if s in columns: return s # zero-based integer column index is allowed try: idx = int(s) return columns[idx] except Exception: pass raise ValueError( f"Column '{selector}' was not found. Available columns: {list(columns)}" ) def read_input_data(args): """Read x,y from Excel or generate demo data.""" infile = str(args.infile).strip() if infile.lower() == "demo": x, y, y_true = generate_demo_data( n=args.n_data, noise=args.noise, seed=args.seed, ) src = "demo" df_input = pd.DataFrame({ "x": x, "y": y, "y_true_demo": y_true, "true_background_demo": true_background(x), "true_gaussian_sum_demo": true_peaks(x), }) return x, y, df_input, src if not os.path.exists(infile): raise FileNotFoundError(f"Input Excel file was not found: {infile}") try: if args.sheet: df = pd.read_excel(infile, sheet_name=args.sheet) else: df = pd.read_excel(infile) except Exception as exc: safe_import_message(exc, "pandas openpyxl") if df.empty: raise ValueError("Input sheet is empty.") columns = list(df.columns) if args.xcol is None: xcol = columns[0] else: xcol = parse_column_selector(args.xcol, columns) if args.ycol is None: ycol = columns[1] if len(columns) >= 2 else None else: ycol = parse_column_selector(args.ycol, columns) if ycol is None: raise ValueError("Could not determine y column. Please specify --ycol.") work = df[[xcol, ycol]].copy() work.columns = ["x", "y"] work["x"] = pd.to_numeric(work["x"], errors="coerce") work["y"] = pd.to_numeric(work["y"], errors="coerce") work = work.dropna(subset=["x", "y"]).copy() if len(work) < 5: raise ValueError("Too few valid numeric x-y data points after removing NaN rows.") work = work.sort_values("x").reset_index(drop=True) # Merge duplicated x by average y. if work["x"].duplicated().any(): work = work.groupby("x", as_index=False)["y"].mean() x = work["x"].to_numpy(dtype=float) y = work["y"].to_numpy(dtype=float) df_input = work.copy() df_input["source_xcol"] = xcol df_input["source_ycol"] = ycol return x, y, df_input, infile def make_default_outfile(infile, suffix="_adaptive_gaussian_ridge.xlsx"): if str(infile).strip().lower() == "demo": return "demo" + suffix stem, _ = os.path.splitext(os.path.basename(infile)) return stem + suffix # ============================================================================= # Basis functions # ============================================================================= def scale_x_to_chebyshev_domain(x, xmin=None, xmax=None): x = np.asarray(x, dtype=float) if xmin is None: xmin = np.min(x) if xmax is None: xmax = np.max(x) if xmax <= xmin: raise ValueError("xmax must be larger than xmin.") return 2.0 * (x - xmin) / (xmax - xmin) - 1.0 def background_design_matrix(x, order, xmin, xmax): if order < 0: return np.empty((len(x), 0)) xr = scale_x_to_chebyshev_domain(x, xmin, xmax) return chebvander(xr, order) def gaussian_design_matrix(x, centers, sigmas): x = np.asarray(x, dtype=float)[:, None] centers = np.asarray(centers, dtype=float)[None, :] sigmas = np.asarray(sigmas, dtype=float)[None, :] return np.exp(-0.5 * ((x - centers) / sigmas) ** 2) def design_matrix(x, centers, sigmas, bg_order, xmin, xmax): phi_bg = background_design_matrix(x, bg_order, xmin, xmax) phi_g = gaussian_design_matrix(x, centers, sigmas) return np.column_stack([phi_bg, phi_g]) def adaptive_centers_from_derivatives( x, y, n_centers=45, eps=0.15, w_grad=0.3, w_curv=1.0, smooth_window=21, smooth_order=3, ): x = np.asarray(x, dtype=float) y = np.asarray(y, dtype=float) if len(x) < 5: raise ValueError("At least 5 data points are required.") # Robust Savitzky-Golay window adjustment. smooth_window = int(smooth_window) smooth_order = int(smooth_order) if smooth_window < smooth_order + 2: smooth_window = smooth_order + 2 if smooth_window % 2 == 0: smooth_window += 1 if smooth_window >= len(x): smooth_window = len(x) if len(x) % 2 == 1 else len(x) - 1 if smooth_window <= smooth_order: smooth_order = max(1, smooth_window - 2) y_smooth = savgol_filter(y, smooth_window, smooth_order) dy = np.gradient(y_smooth, x) d2y = np.gradient(dy, x) grad_score = np.abs(dy) curv_score = np.abs(d2y) grad_score /= np.max(grad_score) + 1.0e-30 curv_score /= np.max(curv_score) + 1.0e-30 density = eps + w_grad * grad_score + w_curv * curv_score dx = np.diff(x) u = np.zeros_like(x) u[1:] = np.cumsum(0.5 * (density[1:] + density[:-1]) * dx) u_centers = np.linspace(u[0], u[-1], n_centers) centers = np.interp(u_centers, u, x) return centers, density, y_smooth, dy, d2y def estimate_local_amplitude( x, y_smooth, centers, window_points=21, p_low=10.0, p_high=90.0, floor_frac=1.0e-4, ): """ Estimate local amplitude around each center without using background. The amplitude is estimated from local percentile width of smoothed y. This avoids using absolute |y|, which becomes problematic when the spectrum has a large offset or large background intensity. """ x = np.asarray(x, dtype=float) y_smooth = np.asarray(y_smooth, dtype=float) centers = np.asarray(centers, dtype=float) n = len(x) window_points = int(window_points) if window_points < 3: window_points = 3 if window_points % 2 == 0: window_points += 1 if window_points > n: window_points = n if n % 2 == 1 else n - 1 half = window_points // 2 global_scale = np.percentile(y_smooth, 95.0) - np.percentile(y_smooth, 5.0) amp_floor = floor_frac * max(abs(global_scale), 1.0e-30) amp = np.empty(len(centers), dtype=float) for i, c in enumerate(centers): idx = int(np.argmin(np.abs(x - c))) i0 = max(0, idx - half) i1 = min(n, idx + half + 1) yw = y_smooth[i0:i1] local_low = np.percentile(yw, p_low) local_high = np.percentile(yw, p_high) local_median = np.median(yw) # local percentile half-width amp_width = 0.5 * abs(local_high - local_low) # deviation from local median y_c = np.interp(c, x, y_smooth) amp_dev = abs(y_c - local_median) # use the larger local scale, but keep a small floor amp[i] = max(amp_width, amp_dev, amp_floor) return amp def estimate_variable_sigma( centers, alpha=1.6, m=2, sigma_min=None, sigma_max=None, method="hybrid_slope", x=None, y_smooth=None, dy=None, alpha_slope=1.0, sigma_amp_window=21, slope_eps=1.0e-6, return_details=False, ): """ Estimate Gaussian widths. method: "distance": sigma = alpha * local center spacing. The local spacing is estimated from +/- m neighboring centers. "nearest": sigma = alpha * nearest left/right center distance. "slope": sigma = alpha_slope * |y_smooth| / (|dy| + slope_eps). "hybrid_slope": sigma = min(sigma_nearest, sigma_slope), then clipped. Notes: The slope-based estimate deliberately does not use background subtraction. This avoids feedback from an unstable background fit when many Gaussian bases are initially present. """ centers = np.asarray(centers, dtype=float) n = len(centers) if n < 2: raise ValueError("At least two centers are required to estimate sigma.") method = str(method).lower() # Conventional distance estimate using +/- m neighbors. sigma_dist = np.empty(n) for j in range(n): i0 = max(0, j - m) i1 = min(n - 1, j + m) if i1 == i0: local_spacing = centers[1] - centers[0] else: local_spacing = (centers[i1] - centers[i0]) / (i1 - i0) sigma_dist[j] = alpha * local_spacing # Nearest-neighbor distance estimate. This is stricter for isolated bases. sigma_nearest = np.empty(n) for j in range(n): if j == 0: d = centers[1] - centers[0] elif j == n - 1: d = centers[-1] - centers[-2] else: d = min(centers[j] - centers[j - 1], centers[j + 1] - centers[j]) sigma_nearest[j] = alpha * d # Slope estimate. Large values are allowed at flat regions and later rejected # by min(distance, slope), sigma_max, and sigma-dependent Ridge penalty. sigma_slope = np.full(n, np.nan, dtype=float) if method in ("slope", "hybrid_slope"): if x is None or y_smooth is None or dy is None: raise ValueError("x, y_smooth, and dy are required for slope-based sigma.") # y_c = np.interp(centers, x, y_smooth) # dy_c = np.interp(centers, x, dy) # sigma_slope = alpha_slope * np.abs(y_c) / (np.abs(dy_c) + slope_eps) amp_c = estimate_local_amplitude( x, y_smooth, centers, window_points=sigma_amp_window, ) dy_c = np.interp(centers, x, dy) sigma_slope = alpha_slope * amp_c / (np.abs(dy_c) + sigma_slope_eps) if method == "distance": sigmas = sigma_dist.copy() elif method == "nearest": sigmas = sigma_nearest.copy() elif method == "slope": sigmas = sigma_slope.copy() elif method == "hybrid_slope": sigmas = np.minimum(sigma_nearest, sigma_slope) else: raise ValueError( f"Unknown sigma method: {method}. " "Use distance, nearest, slope, or hybrid_slope." ) if sigma_min is not None: sigmas = np.maximum(sigmas, sigma_min) if sigma_max is not None: sigmas = np.minimum(sigmas, sigma_max) sigmas = np.maximum(sigmas, 1.0e-30) details = { "sigma_dist": sigma_dist, "sigma_nearest": sigma_nearest, "sigma_slope": sigma_slope, "sigma_raw": sigmas.copy(), } if return_details: return sigmas, details return sigmas # ============================================================================= # Weighted Ridge # ============================================================================= @dataclass class WeightedRidgeResult: coef_: np.ndarray penalty_diag_: np.ndarray def predict(self, xmat): return xmat @ self.coef_ def make_penalty_diag( n_bg, sigmas, lambda_bg=1.0e-8, lambda_gauss=1.0e-3, sigma_penalty_power=2.0, ): sigmas = np.asarray(sigmas, dtype=float) n_g = len(sigmas) p = n_bg + n_g penalty = np.zeros(p, dtype=float) if n_bg > 0: penalty[:n_bg] = lambda_bg if n_g > 0: sigma_ref = np.median(sigmas) sigma_ref = max(sigma_ref, 1.0e-30) scale = (sigmas / sigma_ref) ** sigma_penalty_power penalty[n_bg:] = lambda_gauss * scale return penalty def fit_weighted_ridge( x, y, centers, sigmas, bg_order, xmin, xmax, lambda_bg=1.0e-8, lambda_gauss=1.0e-3, sigma_penalty_power=2.0, ): xmat = design_matrix(x, centers, sigmas, bg_order, xmin, xmax) n_bg = bg_order + 1 if bg_order >= 0 else 0 penalty_diag = make_penalty_diag( n_bg, sigmas, lambda_bg=lambda_bg, lambda_gauss=lambda_gauss, sigma_penalty_power=sigma_penalty_power, ) amat = xmat.T @ xmat + np.diag(penalty_diag) bvec = xmat.T @ y try: coef = np.linalg.solve(amat, bvec) except np.linalg.LinAlgError: print("WARNING: solve failed. Falling back to np.linalg.lstsq.") coef = np.linalg.lstsq(amat, bvec, rcond=None)[0] return WeightedRidgeResult(coef, penalty_diag) def predict(model, x, centers, sigmas, bg_order, xmin, xmax): xmat = design_matrix(x, centers, sigmas, bg_order, xmin, xmax) return model.predict(xmat) def predict_background(model, x, bg_order, xmin, xmax): phi_bg = background_design_matrix(x, bg_order, xmin, xmax) n_bg = bg_order + 1 if bg_order >= 0 else 0 if n_bg <= 0: return np.zeros(len(x), dtype=float) coef_bg = model.coef_[:n_bg] return phi_bg @ coef_bg def calc_metrics(y_ref, y_pred): mae = np.mean(np.abs(y_pred - y_ref)) rmse = np.sqrt(np.mean((y_pred - y_ref) ** 2)) rss = np.sum((y_pred - y_ref) ** 2) return mae, rmse, rss def ridge_effective_df( x, centers, sigmas, bg_order, xmin, xmax, lambda_bg, lambda_gauss, sigma_penalty_power, ): xmat = design_matrix(x, centers, sigmas, bg_order, xmin, xmax) n_bg = bg_order + 1 if bg_order >= 0 else 0 penalty_diag = make_penalty_diag( n_bg, sigmas, lambda_bg=lambda_bg, lambda_gauss=lambda_gauss, sigma_penalty_power=sigma_penalty_power, ) amat = xmat.T @ xmat + np.diag(penalty_diag) try: return np.trace(xmat @ np.linalg.solve(amat, xmat.T)) except np.linalg.LinAlgError: return np.nan def calc_aic_bic(y_ref, y_pred, df): n = len(y_ref) rss = np.sum((y_pred - y_ref) ** 2) rss = max(rss, 1.0e-300) aic = n * np.log(rss / n) + 2.0 * df bic = n * np.log(rss / n) + np.log(n) * df return aic, bic # ============================================================================= # Pruning # ============================================================================= def prune_basis_by_ridge_weight( x, y, centers0, sigmas0, bg_order, xmin, xmax, lambda_bg=1.0e-8, lambda_gauss=1.0e-3, sigma_penalty_power=2.0, prune_frac=0.10, min_basis=5, metric="mae", allowed_ratio=1.02, ): centers = np.asarray(centers0, dtype=float).copy() sigmas = np.asarray(sigmas0, dtype=float).copy() n_bg = bg_order + 1 if bg_order >= 0 else 0 results = [] step = 0 while len(centers) >= min_basis: model = fit_weighted_ridge( x, y, centers, sigmas, bg_order, xmin, xmax, lambda_bg=lambda_bg, lambda_gauss=lambda_gauss, sigma_penalty_power=sigma_penalty_power, ) y_pred = predict(model, x, centers, sigmas, bg_order, xmin, xmax) mae, rmse, rss = calc_metrics(y, y_pred) df = ridge_effective_df( x, centers, sigmas, bg_order, xmin, xmax, lambda_bg=lambda_bg, lambda_gauss=lambda_gauss, sigma_penalty_power=sigma_penalty_power, ) aic, bic = calc_aic_bic(y, y_pred, df) results.append( { "step": step, "n_basis": len(centers), "mae": mae, "rmse": rmse, "rss": rss, "df": df, "aic": aic, "bic": bic, "centers": centers.copy(), "sigmas": sigmas.copy(), "model": model, } ) if len(centers) == min_basis: break gaussian_weights = np.abs(model.coef_[n_bg:]) n_remove = max(1, int(np.ceil(prune_frac * len(centers)))) n_remove = min(n_remove, len(centers) - min_basis) remove_idx = np.argsort(gaussian_weights)[:n_remove] keep = np.ones(len(centers), dtype=bool) keep[remove_idx] = False centers = centers[keep] sigmas = sigmas[keep] step += 1 key = metric.lower() best_value = min(r[key] for r in results) candidates = [r for r in results if r[key] <= allowed_ratio * best_value] selected = min(candidates, key=lambda r: r["n_basis"]) return results, selected def print_prune_results(results, selected, metric): print() print("===== Pruning history =====") print(" step n_gauss MAE RMSE df AIC BIC") for r in results: mark = " <== selected" if r is selected else "" print( f"{r['step']:5d} {r['n_basis']:7d} " f"{r['mae']:10.6g} {r['rmse']:10.6g} " f"{r['df']:10.4f} {r['aic']:11.4f} {r['bic']:11.4f}" f"{mark}" ) print() print(f"Selection metric = {metric}") print(f"Selected n_gauss = {selected['n_basis']}") print(f"Selected MAE = {selected['mae']:.6g}") print(f"Selected RMSE = {selected['rmse']:.6g}") # ============================================================================= # Output # ============================================================================= def build_output_tables( args, source, x, y, df_input, centers0, sigmas0, density, y_smooth, dy, d2y, centers, sigmas, model, xmin, xmax, prune_results, ): x_fit = np.linspace(x.min(), x.max(), args.n_fit) y_fit = predict(model, x_fit, centers, sigmas, args.bg_order, xmin, xmax) y_bg_fit = predict_background(model, x_fit, args.bg_order, xmin, xmax) y_gauss_fit = y_fit - y_bg_fit # Components on input x as well. y_pred_input = predict(model, x, centers, sigmas, args.bg_order, xmin, xmax) y_bg_input = predict_background(model, x, args.bg_order, xmin, xmax) y_gauss_input = y_pred_input - y_bg_input mae_input, rmse_input, rss_input = calc_metrics(y, y_pred_input) df_data = pd.DataFrame({ "x": x, "y_input": y, "y_smooth_for_density": y_smooth, "dy_smooth": dy, "d2y_smooth": d2y, "basis_density": density, "y_fit_at_input_x": y_pred_input, "background_at_input_x": y_bg_input, "gaussian_sum_at_input_x": y_gauss_input, "residual": y - y_pred_input, }) df_fit = pd.DataFrame({ "x": x_fit, "y_fit": y_fit, "background": y_bg_fit, "gaussian_sum": y_gauss_fit, }) if str(args.infile).strip().lower() == "demo": df_fit["true_function_demo"] = true_function(x_fit) df_fit["true_background_demo"] = true_background(x_fit) df_fit["true_gaussian_sum_demo"] = true_peaks(x_fit) n_bg = args.bg_order + 1 if args.bg_order >= 0 else 0 coef_g_for_components = model.coef_[n_bg:] phi_g_fit = gaussian_design_matrix(x_fit, centers, sigmas) df_components = pd.DataFrame({"x": x_fit}) for i in range(len(centers)): df_components[f"G{i:03d}"] = phi_g_fit[:, i] * coef_g_for_components[i] coef_bg = model.coef_[:n_bg] coef_g = model.coef_[n_bg:] penalty_bg = model.penalty_diag_[:n_bg] penalty_g = model.penalty_diag_[n_bg:] df_bg_coef = pd.DataFrame({ "component": [f"Chebyshev_T{i}" for i in range(n_bg)], "coef": coef_bg, "penalty": penalty_bg, }) # Initial vs selected Gaussian bases. selected_set = set(np.round(centers, 14)) # Sigma diagnostic columns for the initial basis. These are useful for # checking whether the final sigma was limited by the nearest-center rule, # by the slope rule, or by sigma_min/sigma_max. try: _, sigma_details = estimate_variable_sigma( centers0, alpha=args.sigma_alpha, sigma_amp_window=args.sigma_amp_window, m=args.sigma_m, sigma_min=args.sigma_min, sigma_max=args.sigma_max, method=args.sigma_method, x=x, y_smooth=y_smooth, dy=dy, alpha_slope=args.sigma_alpha_slope, slope_eps=args.sigma_slope_eps, return_details=True, ) except Exception: sigma_details = { "sigma_dist": np.full(len(centers0), np.nan), "sigma_nearest": np.full(len(centers0), np.nan), "sigma_slope": np.full(len(centers0), np.nan), } df_basis_initial = pd.DataFrame({ "basis_index_initial": np.arange(len(centers0)), "center_initial": centers0, "sigma_initial": sigmas0, "sigma_dist_pm_m": sigma_details["sigma_dist"], "sigma_nearest": sigma_details["sigma_nearest"], "sigma_slope": sigma_details["sigma_slope"], "selected_by_center_match": [bool(np.round(c, 14) in selected_set) for c in centers0], }) df_basis_selected = pd.DataFrame({ "basis_index_selected": np.arange(len(centers)), "center": centers, "sigma": sigmas, "coef": coef_g, "abs_coef": np.abs(coef_g), "penalty": penalty_g, "area_like_coef_sigma_sqrt2pi": coef_g * sigmas * np.sqrt(2.0 * np.pi), }) df_prune = None if prune_results is not None: df_prune = pd.DataFrame([ { "step": r["step"], "n_gaussian_basis": r["n_basis"], "mae": r["mae"], "rmse": r["rmse"], "rss": r["rss"], "effective_df": r["df"], "aic": r["aic"], "bic": r["bic"], "selected": len(r["centers"]) == len(centers), } for r in prune_results ]) df_summary = pd.DataFrame({ "key": [ "source", "mode", "bg_order", "n_input_data", "n_centers_initial", "n_gaussian_selected", "lambda_bg", "lambda_gauss", "sigma_penalty_power", "eps", "w_grad", "w_curv", "sigma_method", "sigma_alpha", "sigma_m", "sigma_alpha_slope", "sigma_slope_eps", "sigma_min", "sigma_max", "smooth_window", "smooth_order", "prune_frac", "min_basis", "metric", "allowed_ratio", "mae_input", "rmse_input", "rss_input", ], "value": [ source, args.mode, args.bg_order, len(x), len(centers0), len(centers), args.lambda_bg, args.lambda_gauss, args.sigma_penalty_power, args.eps, args.w_grad, args.w_curv, args.sigma_method, args.sigma_alpha, args.sigma_m, args.sigma_alpha_slope, args.sigma_slope_eps, args.sigma_min, args.sigma_max, args.smooth_window, args.smooth_order, args.prune_frac, args.min_basis, args.metric, args.allowed_ratio, mae_input, rmse_input, rss_input, ], }) return { "InputData": df_input, "DataAndFit": df_data, "FitSpectrum": df_fit, "GaussianComponents": df_components, "BackgroundCoef": df_bg_coef, "InitialGaussianBasis": df_basis_initial, "SelectedGaussianBasis": df_basis_selected, "PruneHistory": df_prune, "Summary": df_summary, } def write_excel_output(outfile, tables): try: with pd.ExcelWriter(outfile, engine="openpyxl") as writer: for sheet_name, df in tables.items(): if df is None: continue df.to_excel(writer, sheet_name=sheet_name, index=False) # Simple readable widths. ws = writer.book[sheet_name] for col_cells in ws.columns: col_letter = col_cells[0].column_letter max_len = 0 for cell in col_cells[:200]: if cell.value is not None: max_len = max(max_len, len(str(cell.value))) ws.column_dimensions[col_letter].width = min(max(max_len + 2, 10), 32) ws.freeze_panes = "A2" except Exception as exc: safe_import_message(exc, "pandas openpyxl") def plot_results(args, x, y, centers0, sigmas0, density, y_smooth, centers, sigmas, model, xmin, xmax, prune_results): """ Split visualization into three windows: Window A: diagnostics 1. raw data + smoothed data 3. density + sigma + selected centers 5. error vs number of Gaussian bases Window B: main fit raw data + background + Gaussian sum + total fit Window C: decomposition raw data + total fit + Gaussian sum + individual Gaussian components Gi(x) The Gi(x) curves do not include background. """ x_fit = np.linspace(x.min(), x.max(), args.n_fit) y_fit = predict(model, x_fit, centers, sigmas, args.bg_order, xmin, xmax) y_bg_fit = predict_background(model, x_fit, args.bg_order, xmin, xmax) y_gauss_fit = y_fit - y_bg_fit n_bg = args.bg_order + 1 if args.bg_order >= 0 else 0 coef_g = model.coef_[n_bg:] phi_g = gaussian_design_matrix(x_fit, centers, sigmas) g_each = phi_g * coef_g[None, :] # Recompute sigma diagnostic candidates for plotting, when possible. try: _, sigma_details = estimate_variable_sigma( centers0, alpha=args.sigma_alpha, sigma_amp_window=args.sigma_amp_window, m=args.sigma_m, sigma_min=args.sigma_min, sigma_max=args.sigma_max, method=args.sigma_method, x=x, y_smooth=y_smooth, dy=np.gradient(y_smooth, x), alpha_slope=args.sigma_alpha_slope, slope_eps=args.sigma_slope_eps, return_details=True, ) except Exception: sigma_details = None # ------------------------------------------------------------------------- # Window A: diagnostics # ------------------------------------------------------------------------- nrows_a = 3 fig_a, axes = plt.subplots(nrows_a, 1, figsize=(8, 9), sharex=False) fig_a.suptitle("Diagnostics: data, basis density, sigma, pruning history") ax = axes[0] ax.scatter(x, y, s=16, label="input data") ax.plot(x, y_smooth, "-", label="smoothed for density/sigma") ax.set_xlabel("x") ax.set_ylabel("y") ax.legend() ax.grid(True) ax = axes[1] ax.plot(x, density, "-", label="basis density") ax.set_xlabel("x") ax.set_ylabel("density") ax.grid(True) ax2 = ax.twinx() if sigma_details is not None: ax2.plot(centers0, sigma_details["sigma_nearest"], ":", label="sigma nearest candidate") if np.isfinite(sigma_details["sigma_slope"]).any(): ax2.plot(centers0, sigma_details["sigma_slope"], "--", label="sigma slope candidate") ax2.plot(centers0, sigmas0, "o-", label="initial sigma") ax2.plot(centers, sigmas, "s", label="selected sigma") ax2.set_ylabel("sigma") for c in centers: ax.axvline(c, color="0.85", linewidth=0.6, zorder=0) lines1, labels1 = ax.get_legend_handles_labels() lines2, labels2 = ax2.get_legend_handles_labels() ax.legend(lines1 + lines2, labels1 + labels2, loc="best") ax = axes[2] if args.mode == "prune" and prune_results is not None: n_basis = np.array([r["n_basis"] for r in prune_results]) mae = np.array([r["mae"] for r in prune_results]) rmse = np.array([r["rmse"] for r in prune_results]) ax.plot(n_basis, mae, "o-", label="MAE") ax.plot(n_basis, rmse, "s-", label="RMSE") ax.axvline(len(centers), linestyle="--", label="selected") ax.invert_xaxis() ax.set_xlabel("number of Gaussian basis") ax.set_ylabel("error") ax.legend() else: y_pred_input = predict(model, x, centers, sigmas, args.bg_order, xmin, xmax) residual = y - y_pred_input ax.plot(x, residual, ".-", label="residual") ax.axhline(0.0, linestyle="--", linewidth=1) ax.set_xlabel("x") ax.set_ylabel("residual") ax.legend() ax.grid(True) fig_a.tight_layout() # ------------------------------------------------------------------------- # Window B: main fit # ------------------------------------------------------------------------- fig_b, ax = plt.subplots(figsize=(9, 5)) fig_b.suptitle("Main fit: data, background, Gaussian sum, total fit") ax.scatter(x, y, s=16, label="input data") ax.plot(x_fit, y_bg_fit, ":", label="background") ax.plot(x_fit, y_gauss_fit, "-", label="Gaussian sum") ax.plot(x_fit, y_fit, "-", linewidth=2, label="total fit = background + Gaussian sum") if str(args.infile).strip().lower() == "demo": ax.plot(x_fit, true_function(x_fit), "--", label="true demo function") ax.plot(x_fit, true_background(x_fit), "--", label="true demo background") ax.plot(x_fit, true_peaks(x_fit), "--", label="true demo Gaussian sum") ax.set_xlabel("x") ax.set_ylabel("y") ax.legend() ax.grid(True) fig_b.tight_layout() # ------------------------------------------------------------------------- # Window C: decomposition # ------------------------------------------------------------------------- fig_c, ax = plt.subplots(figsize=(9, 5)) fig_c.suptitle("Decomposition: individual Gaussian components Gi(x), no background") ax.scatter(x, y, s=14, label="input data") ax.plot(x_fit, y_fit, "-", linewidth=2, label="total fit") ax.plot(x_fit, y_gauss_fit, "-", linewidth=2, label="Gaussian sum") # Plot selected Gaussian components. Suppress individual labels to avoid # unreadable legends when many bases remain. for i in range(g_each.shape[1]): ax.plot(x_fit, g_each[:, i], "--", linewidth=1, alpha=0.55) ax.axhline(0.0, linestyle=":", linewidth=1) ax.set_xlabel("x") ax.set_ylabel("y or component amplitude") ax.legend() ax.grid(True) fig_c.tight_layout() if args.save_plot: base, ext = os.path.splitext(args.plotfile) if ext == "": ext = ".png" files = [ f"{base}_diagnostics{ext}", f"{base}_mainfit{ext}", f"{base}_decomposition{ext}", ] for fig, filename in zip([fig_a, fig_b, fig_c], files): fig.savefig(filename, dpi=200) print(f"Saved plot: {filename}") if args.show: plt.show(block=False) plt.pause(0.1) input("Press ENTER to close figures>>\n") plt.close("all") else: plt.close(fig_a) plt.close(fig_b) plt.close(fig_c) # ============================================================================= # Main # ============================================================================= def main(): parser = argparse.ArgumentParser( description="Adaptive Gaussian basis + background + weighted Ridge + pruning for spectra" ) parser.add_argument("--infile", type=str, default="demo", help="Input Excel file. Use 'demo' for generated data.") parser.add_argument("--sheet", type=str, default="", help="Excel sheet name or index. Default: ''") parser.add_argument("--xcol", type=str, default=None, help="x column name or zero-based column index. Default: first column") parser.add_argument("--ycol", type=str, default=None, help="y column name or zero-based column index. Default: second column") parser.add_argument("--outfile", type=str, default=None, help="Output Excel file. Default: created from infile name.") parser.add_argument("--mode", type=str, default="prune", choices=["fit", "prune"]) parser.add_argument("--n-data", type=int, default=180, help="Number of generated data points for --infile demo.") parser.add_argument("--n-fit", type=int, default=1000, help="Number of points in exported fitted spectra.") parser.add_argument("--n-centers", type=int, default=60) parser.add_argument("--noise", type=float, default=0.06) parser.add_argument("--seed", type=int, default=1) parser.add_argument("--bg-order", type=int, default=3, help="Chebyshev background order. Use -1 for no background.") parser.add_argument("--lambda-bg", type=float, default=1.0e-8) parser.add_argument("--lambda-gauss", type=float, default=1.0e-3) parser.add_argument("--sigma-penalty-power", type=float, default=2.0) parser.add_argument("--eps", type=float, default=0.15) parser.add_argument("--w-grad", type=float, default=0.3) parser.add_argument("--w-curv", type=float, default=1.0) parser.add_argument("--sigma-method", type=str, default="nearest", # parser.add_argument("--sigma-method", type=str, default="hybrid_slope", choices=["distance", "nearest", "slope", "hybrid_slope"], help="Sigma estimation method. hybrid_slope uses min(nearest-distance, slope estimate).") parser.add_argument("--sigma-alpha", type=float, default=1.6, help="Distance-based sigma multiplier.") parser.add_argument("--sigma-m", type=int, default=2, help="Neighbor count for conventional distance sigma.") parser.add_argument("--sigma-alpha-slope", type=float, default=1.0, help="Slope-based sigma multiplier.") parser.add_argument( "--sigma-amp-window", type=int, default=21, help="Window points for local amplitude estimation used in slope-based sigma.", ) parser.add_argument("--sigma-slope-eps", type=float, default=1.0e-6, help="Small denominator for |y|/(|dy/dx|+eps) sigma estimate.") parser.add_argument("--sigma-min", type=float, default=None) parser.add_argument("--sigma-max", type=float, default=None) parser.add_argument("--smooth-window", type=int, default=21) parser.add_argument("--smooth-order", type=int, default=3) parser.add_argument("--prune-frac", type=float, default=0.10) parser.add_argument("--min-basis", type=int, default=5) parser.add_argument("--metric", type=str, default="mae", choices=["mae", "rmse"]) parser.add_argument("--allowed-ratio", type=float, default=1.02) parser.add_argument("--show", type=int, default=1, choices=[0, 1]) parser.add_argument("--save-plot", type=int, default=0, choices=[0, 1]) parser.add_argument("--plotfile", type=str, default="adaptive_gaussian_ridge_plot.png") args = parser.parse_args() args.show = bool(args.show) args.save_plot = bool(args.save_plot) if args.outfile is None: args.outfile = make_default_outfile(args.infile) # Make output path relative to current directory unless user gives directory. if os.path.dirname(args.outfile): os.makedirs(os.path.dirname(args.outfile), exist_ok=True) x, y, df_input, source = read_input_data(args) xmin = np.min(x) xmax = np.max(x) centers0, density, y_smooth, dy, d2y = adaptive_centers_from_derivatives( x, y, n_centers=args.n_centers, eps=args.eps, w_grad=args.w_grad, w_curv=args.w_curv, smooth_window=args.smooth_window, smooth_order=args.smooth_order, ) sigmas0 = estimate_variable_sigma( centers0, alpha=args.sigma_alpha, sigma_amp_window=args.sigma_amp_window, m=args.sigma_m, sigma_min=args.sigma_min, sigma_max=args.sigma_max, method=args.sigma_method, x=x, y_smooth=y_smooth, dy=dy, alpha_slope=args.sigma_alpha_slope, slope_eps=args.sigma_slope_eps, ) if args.mode == "fit": model = fit_weighted_ridge( x, y, centers0, sigmas0, args.bg_order, xmin, xmax, lambda_bg=args.lambda_bg, lambda_gauss=args.lambda_gauss, sigma_penalty_power=args.sigma_penalty_power, ) centers = centers0 sigmas = sigmas0 prune_results = None else: prune_results, selected = prune_basis_by_ridge_weight( x, y, centers0, sigmas0, args.bg_order, xmin, xmax, lambda_bg=args.lambda_bg, lambda_gauss=args.lambda_gauss, sigma_penalty_power=args.sigma_penalty_power, prune_frac=args.prune_frac, min_basis=args.min_basis, metric=args.metric, allowed_ratio=args.allowed_ratio, ) print_prune_results(prune_results, selected, args.metric) centers = selected["centers"] sigmas = selected["sigmas"] model = selected["model"] y_pred_input = predict(model, x, centers, sigmas, args.bg_order, xmin, xmax) mae_input, rmse_input, _ = calc_metrics(y, y_pred_input) print() print("===== Adaptive Gaussian Weighted Ridge Regression =====") print(f"source = {source}") print(f"outfile = {args.outfile}") print(f"mode = {args.mode}") print(f"bg_order = {args.bg_order}") print(f"n_input_data = {len(x)}") print(f"n_centers_initial = {len(centers0)}") print(f"n_gauss_selected = {len(centers)}") print(f"lambda_bg = {args.lambda_bg}") print(f"lambda_gauss = {args.lambda_gauss}") print(f"sigma_penalty_power = {args.sigma_penalty_power}") print(f"sigma_method = {args.sigma_method}") print(f"sigma_alpha = {args.sigma_alpha}") print(f"sigma_alpha_slope = {args.sigma_alpha_slope}") print(f"sigma_min = {args.sigma_min}") print(f"sigma_max = {args.sigma_max}") print(f"MAE input = {mae_input:.6g}") print(f"RMSE input = {rmse_input:.6g}") tables = build_output_tables( args, source, x, y, df_input, centers0, sigmas0, density, y_smooth, dy, d2y, centers, sigmas, model, xmin, xmax, prune_results, ) write_excel_output(args.outfile, tables) print(f"Saved Excel output: {args.outfile}") plot_results( args, x, y, centers0, sigmas0, density, y_smooth, centers, sigmas, model, xmin, xmax, prune_results, ) if __name__ == "__main__": try: main() except Exception: print("ERROR: program terminated unexpectedly.") traceback.print_exc() input("\nPress ENTER to terminate>>\n") sys.exit(1)