Научный анализ данных с LabPlot в Python: обработка сигналов, аппроксимация спектральных пиков, визуализация и пакетная автоматизация

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__:<20} {self._name}"
       if isinstance(self, Column):
           s += f"  [{self.columnMode.name}, {self.rowCount()} rows, {self.plotDesignation.name}]"
       return "\n".join([s] + [c.tree(d+1) for c in self.children])
class Column(AbstractAspect):
   """LabPlot's fundamental data source: a typed vector + a plot designation."""
   def __init__(self, name, values=None, mode=ColumnMode.Double,
                designation=PlotDesignation.NoDesignation):
       super().__init__(name)
       self.columnMode, self.plotDesignation = mode, designation
       self._d = np.asarray([] if values is None else values, float)
   def values(self): return self._d
   def rowCount(self): return len(self._d)
   def clean(self): return self._d[np.isfinite(self._d)]
   def statistics(self):
       """The 20 quantities in LabPlot's Column Statistics dialog."""
       x = self.clean(); n = x.size
       if not n: return {}
       q1, med, q3 = np.percentile(x, [25, 50, 75]); iqr, pos = q3 - q1, x[x > 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": (x**2).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 < 2 or np.ptp(s) == 0 else "".join(
           b[i] for i in ((s - s.min()) / np.ptp(s) * 4).round().astype(int))
class Spreadsheet(AbstractAspect):
   def columns(self): return [c for c in self.children if isinstance(c, Column)]
   def column(self, k):
       cs = self.columns()
       return cs[k] if isinstance(k, int) else next(c for c in cs if c.name() == k)
   def columnCount(self): return len(self.columns())
   def rowCount(self): return max([c.rowCount() for c in self.columns()], default=0)
   def appendColumn(self, n, v, d=PlotDesignation.Y):
       return self.addChild(Column(n, v, designation=d))
   def toDataFrame(self):
       return pd.DataFrame({c.name(): c.values() for c in self.columns()})
   def info(self):
       print(f"  Spreadsheet '{self._name}': {self.rowCount()} rows x {self.columnCount()} cols")
       for c in self.columns():
           s = c.statistics()
           print(f"  {c.name():<12}{c.plotDesignation.name:<6}min {s['Minimum']:>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 < 0 else self.endRow]
       if not lines: raise ValueError("AsciiFilter: nothing to import")
       sep = (next((s for s in (",", ";", "\t", "|") if s in lines[0]), None)
              if self.separator == "auto" else self.separator)
       split = lambda l: [p.strip() for p in (l.split(sep) if sep else l.split()) if p.strip()]
       header = split(lines[0]) if self.headerEnabled else None
       rows = [split(l) for l in (lines[1:] if self.headerEnabled else lines)]
       ncol = max(map(len, rows))
       c0, c1 = self.startColumn - 1, ncol if self.endColumn < 0 else self.endColumn
       for j in range(c0, min(c1, ncol)):
           vals = []
           for r in rows:
               try: vals.append(float(r[j]))
               except (IndexError, ValueError): vals.append(np.nan)
           dataSource.addChild(Column(
               header[j] if header and j < len(header) else f"Column {j+1}", vals,
               designation=PlotDesignation.X if j == c0 else PlotDesignation.Y))
       return dataSource

class nsl_smooth:
   """Analysis -> 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 < 2: return float(np.trapezoid(y, x) if hasattr(np, "trapezoid") else np.trapz(y, x))
   tot, i = 0.0, 0
   while i + 2 <= n:
       h0, h1 = x[i+1] - x[i], x[i+2] - x[i+1]; hp, hd, hm = h1 + h0, h1 / h0, h1 * h0
       tot += hp / 6 * ((2-hd) * y[i] + hp**2 / hm * y[i+1] + (2 - 1/hd) * y[i+2]); i += 2
   return tot + ((x[n]-x[n-1]) * (y[n]+y[n-1]) / 2 if i < n else 0)
class nsl_int:
   """Analysis -> 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: m**2 / (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((f**2 - f0**2) / np.maximum(f * bw, eps))
           if type == "bandreject": r = 1 / np.maximum(r, eps)
       H = (r <= 1).astype(float) if form == "ideal" else 1 / np.sqrt(1 + r ** (2 * order))
       return np.fft.irfft(np.fft.rfft(y) * H, n=n)
class nsl_hilbert:
   """Analysis -> 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 <= i + 1: continue
           dx, dy = x[j]-x[i], y[j]-y[i]; den = math.hypot(dx, dy); sl = slice(i+1, j)
           d = (np.hypot(x[sl]-x[i], y[sl]-y[i]) if den == 0
                else np.abs(dy*(x[sl]-x[i]) - dx*(y[sl]-y[i])) / den)
           if d.size and d.max() > 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':<8}{'value':>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:<8}{v:>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]:<22}{it[i+1][1]:>14.6g}" if i + 1 < len(it) else ""
           print(f"  {it[i][0]:<22}{it[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((r**2).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 < 1 else np.inf
       logL = -.5*n * (math.log(2*math.pi) + math.log(sse/n) + 1); aic = 2*k - 2*logL
       gof = {"sum sq. residuals": sse, "mean squared error": sse/n, "root MSE": math.sqrt(sse/n),
              "mean abs. error": float(np.abs(r).mean()), "residual std dev": math.sqrt(sse/dof),
              "R^2": r2, "adjusted R^2": 1 - (1-r2)*(n-1)/dof, "chi^2": chisq,
              "reduced chi^2": red, "P > 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: d.pdf(t, *pr)}