AI 服務暫時不可用，以下為來源正文，待恢復後補全翻譯。
In this tutorial, we explore a LabPlot-inspired scientific data analysis workflow in Python while preserving the structure and terminology of LabPlot’s aspect tree, analysis kernels, plotting system, and project model. We build reusable components to import tabular data, compute descriptive statistics, smooth and differentiate signals, perform Fourier analysis and filtering, detect peaks, integrate curves, reduce data, and fit nonlinear models with detailed statistical diagnostics. We then apply these tools to a realistic spectroscopy example: removing periodic interference, identifying overlapping peaks, fitting a multi-Gaussian model, inspecting residuals, visualizing results through themed worksheets, exporting figures, and saving project data in LabPlot-compatible .lml-style files. Finally, we extend the same workflow to batch processing so we can analyze multiple temperature-dependent spectra and fit secondary trends across the resulting measurements. Copy CodeCopiedUse a different Browser import os, sys, gzip, bz2, lzma, time, math, textwrap, warnings import xml.etree.ElementTree as ET from dataclasses import dataclass, field from enum import Enum import numpy as np, pandas as pd, matplotlib, matplotlib.pyplot as plt from matplotlib.ticker import AutoMinorLocator import scipy from scipy import signal, stats, optimize warnings.filterwarnings("ignore", category=RuntimeWarning) np.random.seed(20260815) IN_COLAB = "google.colab" in sys.modules OUT = "/content/labplot_out" if IN_COLAB else os.path.join(os.getcwd(), "labplot_out") os.makedirs(OUT, exist_ok=True) try: from pylabplot import *; HAVE_SDK = True except Exception: HAVE_SDK = False banner = lambda t: print("\n" + "=" * 76 + f"\n {t}\n" + "=" * 76) banner("environment") print(f" numpy {np.version} | scipy {scipy.version} | mpl {matplotlib.version} | " f"colab={IN_COLAB} | pylabplot={'yes' if HAVE_SDK else 'no -> emulation'}\n -> {OUT}") class PlotDesignation(Enum): NoDesignation = 0; X = 1; Y = 2; Z = 3 XError = 4; XErrorMinus = 5; XErrorPlus = 6 YError = 7; YErrorMinus = 8; YErrorPlus = 9 class ColumnMode(Enum): Double = 0; Text = 1; Integer = 2; BigInt = 3; DateTime = 4 class AbstractAspect: def init(self, name, comment=""): self._name, self.comment, self.parent, self.children = name, comment, None, [] def name(self): return self._name def addChild(self, a): a.parent = self; self.children.append(a); return a def tree(self, d=0): s = " " * d + f"{'|- ' if d else ''}{type(self).name: 0] h = 2 * iqr / n(1/3) if iqr > 0 else 0 c, _ = np.histogram(x, bins=int(np.clip(np.ptp(x)/h, 1, 1000)) if h else 10) p = c[c > 0] / c.sum(); v, k = np.unique(np.round(x, 12), return_counts=True) return {"Count": n, "Minimum": x.min(), "Maximum": x.max(), "Arithmetic mean": x.mean(), "Geometric mean": stats.gmean(pos) if pos.size else np.nan, "Harmonic mean": stats.hmean(pos) if pos.size else np.nan, "Contraharmonic mean": (x2).sum() / x.sum() if x.sum() else np.nan, "Mode": v[k.argmax()] if k.max() > 1 else np.nan, "First quartile": q1, "Median": med, "Third quartile": q3, "Interquartile range": iqr, "Trimean": (q1 + 2*med + q3) / 4, "Variance": x.var(ddof=1), "Standard deviation": x.std(ddof=1), "Skewness": stats.skew(x), "Mean absolute deviation": np.abs(x - x.mean()).mean(), "Median absolute deviation": np.median(np.abs(x - med)), "Kurtosis": stats.kurtosis(x, fisher=False), "Entropy": float(-(p * np.log2(p)).sum())} def sparkline(self, w=26): """LabPlot 2.11+ draws these in the column header; text version.""" b, x = "_.-~^", self.clean() s = x[np.linspace(0, x.size-1, min(w, x.size)).astype(int)] if x.size > 1 else x return "" if s.size 9.4g} max " f"{s['Maximum']:>9.4g} mean {s['Arithmetic mean']:>9.4g} {c.sparkline()}") class Project(AbstractAspect): XML_VERSION = 15 def init(self, name="project", author=""): super().init(name); self.author, self.version = author, "2.12.1" def spreadsheets(self): return [c for c in self.children if isinstance(c, Spreadsheet)] class AsciiFilter: """LabPlot's text import: separator auto-detect, comments, row/col limits.""" def init(self, separator="auto", commentCharacter="#", headerEnabled=True, startRow=1, endRow=-1, startColumn=1, endColumn=-1): self.separator, self.commentCharacter = separator, commentCharacter self.headerEnabled, self.startRow, self.endRow = headerEnabled, startRow, endRow self.startColumn, self.endColumn = startColumn, endColumn def readDataFromFile(self, path, dataSource): with open(path, encoding="utf-8", errors="replace") as fh: lines = [l.rstrip("\n") for l in fh if l.strip() and not l.lstrip().startswith(self.commentCharacter)] lines = lines[self.startRow - 1: None if self.endRow Smooth (Savitzky-Golay; LabPlot also offers moving average/percentile).""" @staticmethod def savitzky_golay(y, points=11, order=3, deriv=0): points += points % 2 == 0 return signal.savgol_filter(y, points, min(order, points-1), deriv=deriv, mode="interp") class nsl_diff: """Analysis -> Differentiate: order 1..6; SG differentiation for noisy data.""" @staticmethod def derive(x, y, order=1, smooth_points=0, sg_order=3): if smooth_points: return nsl_smooth.savitzky_golay(y, smooth_points, sg_order, deriv=order) \ / np.gradient(x) order out = np.asarray(y, float) for _ in range(order): out = np.gradient(out, x, edge_order=2) return out def simpson(x, y): """Composite Simpson on a non-uniform grid (Cartwright's formula).""" n = len(x) - 1 if n Integrate: rectangle / trapezoid / Simpson, cumulative.""" @staticmethod def integrate(x, y, method="trapezoid", absolute=False): yy = np.abs(y) if absolute else np.asarray(y, float) seg = np.diff(x) * (yy[:-1] if method == "rectangle" else (yy[:-1] + yy[1:]) / 2) cum = np.r_[0.0, np.cumsum(seg)] return cum * simpson(x, yy) / cum[-1] if method == "simpson" and cum[-1] else cum class nsl_dft: """Analysis -> Fourier Transform: amplitude/magnitude/power/dB, 5 windows.""" WIN = {"rectangular": np.ones, "hann": lambda n: signal.windows.hann(n, sym=False), "hamming": lambda n: signal.windows.hamming(n, sym=False), "blackman": lambda n: signal.windows.blackman(n, sym=False), "flattop": lambda n: signal.windows.flattop(n, sym=False)} @staticmethod def transform(x, y, output="amplitude", window="rectangular"): n, dt = len(y), float(np.mean(np.diff(x))) w = nsl_dft.WIN[window](n); cg = w.mean() Y = np.fft.rfft(y * w); f = np.fft.rfftfreq(n, dt); m = np.abs(Y) v = {"magnitude": lambda: m, "power": lambda: m2 / (n*cg)2, "phase": lambda: np.angle(Y), "amplitude": lambda: np.r_[m[0]/(n*cg), 2*m[1:]/(n*cg)], "dB": lambda: 20*np.log10(np.maximum(m / (m.max() or 1), 1e-16))}[output]() return f, v class nsl_filter: """Analysis -> Fourier Filter: low/high/band pass + band reject; ideal or Butterworth.""" @staticmethod def apply(x, y, type="lowpass", form="butterworth", cutoff=.1, cutoff2=.3, order=3): n = len(y); f = np.fft.rfftfreq(n, float(np.mean(np.diff(x)))); eps = 1e-30 if type == "lowpass": r = f / cutoff elif type == "highpass": r = cutoff / np.maximum(f, eps) else: f0, bw = math.sqrt(cutoff * cutoff2), cutoff2 - cutoff r = np.abs((f2 - f02) / np.maximum(f * bw, eps)) if type == "bandreject": r = 1 / np.maximum(r, eps) H = (r Hilbert Transform (LabPlot 2.9+).""" @staticmethod def transform(y, output="envelope"): a = signal.hilbert(y) return {"imag": a.imag, "real": a.real, "envelope": np.abs(a), "phase": np.unwrap(np.angle(a))}[output] class nsl_geom: """Analysis -> Data Reduction: Douglas-Peucker, iterative (no recursion limit).""" @staticmethod def douglas_peucker(x, y, tol): n = len(x); keep = np.zeros(n, bool); keep[[0, -1]] = True; stack = [(0, n-1)] while stack: i, j = stack.pop() if j tol: k = i + 1 + int(d.argmax()); keep[k] = True; stack += [(i, k), (k, j)] return np.flatnonzero(keep) class nsl_peak: """Analysis -> Peak Find (LabPlot 2.11+); seeds multi-peak fits.""" @staticmethod def find(x, y, prominence=None, distance=None): pk, pr = signal.find_peaks(y, prominence=prominence, distance=distance) w = signal.peak_widths(y, pk, rel_height=.5)[0] if pk.size else np.array([]) return pk, {"positions": x[pk], "heights": y[pk], "prominences": pr.get("prominences", np.array([])), "fwhm": w * float(np.mean(np.diff(x)))} class nsl_fit_model: """LabPlot's model catalogue (Basic / Peak / Growth / Distribution).""" @staticmethod def gaussian(x, a, mu, s): return a / (math.sqrt(2 * np.pi) * s) * np.exp(-(x - mu) 2 / (2 * s 2)) @staticmethod def lorentz(x, a, mu, g): return a / np.pi * (g / 2) / ((x - mu) 2 + (g / 2) 2) @dataclass class FitResult: names: list; values: np.ndarray; errors: np.ndarray; t: np.ndarray; p: np.ndarray margin: np.ndarray; gof: dict; dof: int; nfev: int; status: str elapsed: float; unweighted: bool residuals: np.ndarray = field(repr=False, default=None) cov: np.ndarray = field(repr=False, default=None) def report(self, title="Fit result"): print(f"\n{'-'*76}\n {title}\n{'-'*76}") print(f" {self.status} | nfev {self.nfev} | dof {self.dof} | {self.elapsed*1e3:.1f} ms") print(f"\n {'param':13}{'error':>11}{'err%':>8}{'t':>8}{'P>|t|':>10}{'95% CI':>26}") for i, n in enumerate(self.names): v, e, m = self.values[i], self.errors[i], self.margin[i] print(f" {n:13.6g}{e:>11.4g}{abs(100*e/v) if v else np.inf:>7.2f}%" f"{self.t[i]:>8.1f}{self.p[i]:>10.2g}{f'[{v-m:.5g},{v+m:.5g}]':>26}") print("\n goodness of fit") it = list(self.gof.items()) for i in range(0, len(it), 2): r = f"{it[i+1][0]:14.6g}" if i + 1 14.6g} {r}") if self.unweighted: print(" note: no y-errors given, so chi^2 == SSE and 'P > chi^2' is not a real\n" " test. Pass yerr= for a meaningful reduced chi^2.") print("-" * 76) class nsl_fit: @staticmethod def fit(model, x, y, p0, yerr=None, bounds=None, paramNames=None, conf=.95): """GSL's multifit_nlinear == scipy least_squares(method='lm').""" t0 = time.perf_counter() x, y, p0 = np.asarray(x, float), np.asarray(y, float), np.asarray(p0, float) sig = np.ones_like(y) if yerr is None else np.asarray(yerr, float) res = lambda p: (model(x, *p) - y) / sig kw = dict(max_nfev=500*len(p0), ({"method": "lm"} if bounds is None else {"bounds": bounds})) out = optimize.least_squares(res, p0, kw) p = out.x; n, k = len(y), len(p); dof = max(n - k, 1); r = y - model(x, *p) sse = float((r2).sum()); chisq = float(((r / sig)2).sum()); red = chisq / dof try: cov = np.linalg.inv(out.jac.T @ out.jac) except np.linalg.LinAlgError: cov = np.linalg.pinv(out.jac.T @ out.jac) cov = cov * (red if yerr is None else 1.0); err = np.sqrt(np.abs(np.diag(cov))) tv = np.divide(p, err, out=np.full_like(p, np.inf), where=err > 0) sst = float(((y - y.mean())2).sum()); r2 = 1 - sse/sst if sst else np.nan F = (r2 / max(k-1, 1)) / ((1-r2) / dof) if r2 chi^2": stats.chi2.sf(chisq, dof), "F statistic": F, "P > F": stats.f.sf(F, max(k-1, 1), dof), "log-likelihood": logL, "AIC": aic, "AICc": aic + 2*k*(k+1)/max(n-k-1, 1), "BIC": k*math.log(n) - 2*logL} return FitResult(paramNames or [f"p{i}" for i in range(k)], p, err, tv, 2 * stats.t.sf(np.abs(tv), dof), stats.t.ppf(.5 + conf/2, dof) * err, gof, dof, int(out.nfev), out.message, time.perf_counter() - t0, yerr is None, r, cov) @staticmethod def confidenceBand(model, x, res, level=.95, eps=1e-7): """Delta method sqrt(diag(J C J^T)) * t -- LabPlot's CI overlay.""" p = res.values; J = np.empty((len(x), len(p))) for i in range(len(p)): dp = np.zeros_like(p); dp[i] = eps * max(abs(p[i]), 1) J[:, i] = (model(x, *(p + dp)) - model(x, *(p - dp))) / (2 * dp[i]) v = np.einsum("ij,jk,ik->i", J, res.cov, J) return stats.t.ppf(.5 + level/2, res.dof) * np.sqrt(np.maximum(v, 0)) @staticmethod def distributionFitML(data, dist="norm"): """nsl_fit_algorithm_ml -- max-likelihood distribution fit (the SDK demo).""" d = getattr(stats, dist); pr = d.fit(data); ks = stats.kstest(data, dist, args=pr) ll = float(d.logpdf(data, *pr).sum()) return {"params": pr, "logLik": ll, "AIC": 2 * len(pr) - 2 * ll, "KS_stat": ks.statistic, "KS_p": ks.pvalue, "pdf": lambda t: [truncated for AI cost control]