{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# HSI preprocessing in Python: bad bands, noise and scaling\n", "\n", "A self-contained, CPU-only teaching notebook. All spectra and targets below are synthetic, generated with seed 20261002. No real-world prediction gain is claimed. The spectrum axis is in nm and original signal is fractional reflectance.\n", "\n", "Tested with Python 3.12.14, NumPy 2.3.5, SciPy 1.17.0 and scikit-learn 1.8.0. Install the numerical requirements in your own environment if needed; this notebook does not install anything or access the network.\n", "\n", "## 1. Explicit contracts and reusable functions\n", "\n", "Rows are observations; columns are bands. Metadata decide bad bands. Retain the original uniform coordinate grid, detrend valid runs when requested, apply the selected row-wise or training-fitted normalisation, then filter contiguous runs separately. The functions reject missing values and irregular grids instead of silently imputing or bridging gaps. SNV uses sample standard deviation (ddof=1)." ], "id": "hsi-preprocessing-00" }, { "cell_type": "code", "metadata": {}, "source": [ "\"\"\"Auditable HSI preprocessing on synthetic spectra; CPU only.\n", "\n", "The example is not a chickpea-study reproduction or a performance benchmark.\n", "Axis contract: rows are observations, columns are wavelength samples in nm.\n", "\"\"\"\n", "from __future__ import annotations\n", "import argparse\n", "import json\n", "from pathlib import Path\n", "import numpy as np\n", "from scipy.signal import savgol_filter, savgol_coeffs\n", "from sklearn.base import BaseEstimator, TransformerMixin\n", "from sklearn.decomposition import PCA\n", "from sklearn.linear_model import Ridge\n", "from sklearn.model_selection import GroupShuffleSplit\n", "from sklearn.pipeline import Pipeline\n", "from sklearn.preprocessing import StandardScaler\n", "from sklearn.utils.validation import check_is_fitted\n", "\n", "\n", "def canonicalise(X, wavelength, *, unit=\"nm\", sort=False):\n", " \"\"\"Return float data, increasing nm wavelengths, original-column permutation.\n", "\n", " Sorting is an explicit operation; apply the returned permutation to every\n", " band-level annotation as well. No units are guessed from numerical values.\n", " Missing data are rejected rather than silently filled.\n", " \"\"\"\n", " X, w = np.asarray(X, dtype=float), np.asarray(wavelength, dtype=float)\n", " if X.ndim != 2 or w.ndim != 1 or X.shape[1] != w.size or w.size < 2:\n", " raise ValueError(\"Expected X[n_observations, n_bands] and matching wavelengths.\")\n", " if not np.isfinite(X).all() or not np.isfinite(w).all():\n", " raise ValueError(\"Non-finite data: apply an explicit missing-data policy first.\")\n", " factors = {\"nm\": 1.0, \"um\": 1000.0, \"\u00b5m\": 1000.0}\n", " if unit not in factors:\n", " raise ValueError(\"Declare wavelength unit as nm, um or \u00b5m.\")\n", " w = w * factors[unit]\n", " if np.any(w <= 0):\n", " raise ValueError(\"Wavelengths must be positive.\")\n", " order = np.argsort(w, kind=\"stable\") if sort else np.arange(w.size)\n", " w, X = w[order], X[:, order]\n", " if np.any(np.diff(w) <= 0):\n", " raise ValueError(\"Wavelengths must be strictly increasing and unique; reorder explicitly.\")\n", " return X, w, order\n", "\n", "\n", "def uniform_step(wavelength_nm):\n", " \"\"\"Reject deleted-band gaps and irregular sampling; never use median delta.\"\"\"\n", " w = np.asarray(wavelength_nm, dtype=float)\n", " if w.ndim != 1 or len(w) < 2 or not np.isfinite(w).all():\n", " raise ValueError(\"At least two finite wavelengths are required.\")\n", " d = np.diff(w)\n", " if d[0] <= 0 or not np.allclose(d, d[0], rtol=1e-6, atol=1e-9):\n", " raise ValueError(\"Non-uniform wavelength spacing or deleted-band gap: use separate contiguous segments; do not smooth across the gap.\")\n", " return float(d[0])\n", "\n", "\n", "def snv(X, *, minimum_sd=1e-12):\n", " \"\"\"Per-row SNV using sample SD (ddof=1), across retained bands only.\"\"\"\n", " X = np.asarray(X, dtype=float)\n", " if X.ndim != 2 or X.shape[1] < 2 or not np.isfinite(X).all():\n", " raise ValueError(\"SNV needs at least two finite retained bands per row.\")\n", " scale = X.std(axis=1, ddof=1, keepdims=True)\n", " if np.any(scale <= minimum_sd):\n", " raise ValueError(\"SNV is undefined or unstable for a flat/near-flat spectrum.\")\n", " return (X - X.mean(axis=1, keepdims=True)) / scale\n", "\n", "\n", "def retained_segments(keep):\n", " keep = np.asarray(keep)\n", " if keep.ndim != 1 or keep.dtype.kind != \"b\" or not keep.any():\n", " raise ValueError(\"keep must be a non-empty boolean band mask.\")\n", " ids = np.flatnonzero(keep)\n", " return np.split(ids, np.flatnonzero(np.diff(ids) != 1) + 1)\n", "\n", "\n", "NORMALISATIONS = (\"none\", \"minmax\", \"zscore\", \"snv\", \"l2\", \"msc\", \"feature_zscore\")\n", "OPERATIONS = (\"none\", \"smooth\", \"derivative\", \"derivative2\", \"moving_average\")\n", "\n", "\n", "def baseline_correct(X, wavelength_nm, keep, degree=0):\n", " \"\"\"Subtract a polynomial fitted per spectrum and contiguous segment.\n", "\n", " degree=0 means no detrending; degree=1 or 2 subtracts the corresponding\n", " least-squares polynomial. Coordinates are centred/scaled for conditioning.\n", " This can remove real broad spectral structure as well as a nuisance baseline.\n", " \"\"\"\n", " X, w, _ = canonicalise(X, wavelength_nm)\n", " uniform_step(w) # deliberately strict contract for this combined workflow\n", " keep = np.asarray(keep)\n", " if keep.shape != w.shape:\n", " raise ValueError(\"Band mask and wavelengths must have the same shape.\")\n", " runs = retained_segments(keep)\n", " if degree not in (0, 1, 2):\n", " raise ValueError(\"Detrending degree must be 0 (off), 1 or 2.\")\n", " out = np.full_like(X, np.nan)\n", " out[:, keep] = X[:, keep]\n", " if degree:\n", " for ids in runs:\n", " if len(ids) <= degree:\n", " raise ValueError(\"Too few bands for this detrending degree.\")\n", " t = (w[ids] - (w[ids[0]] + w[ids[-1]]) / 2) / ((w[ids[-1]] - w[ids[0]]) / 2)\n", " design = np.vander(t, N=degree+1, increasing=True)\n", " coefficients = np.linalg.lstsq(design, X[:, ids].T, rcond=None)[0]\n", " out[:, ids] -= (design @ coefficients).T\n", " return out\n", "\n", "\n", "def fit_training_statistics(X, wavelength_nm, keep, *, detrend=0):\n", " \"\"\"Learn only from training rows, after the same fixed mask/detrending.\n", "\n", " Returned metadata guards against applying the reference to a different band\n", " schema or detrending recipe. Sliders/test rows must never update this object.\n", " \"\"\"\n", " keep = np.asarray(keep)\n", " Z = baseline_correct(X, wavelength_nm, keep, detrend)[:, keep]\n", " if Z.shape[0] < 2:\n", " raise ValueError(\"Use at least two training observations for this example.\")\n", " scale = Z.std(axis=0, ddof=0)\n", " scale[scale <= 1e-12] = 1.0 # constant features remain centred, not divided by 0\n", " mean = Z.mean(axis=0)\n", " return {\"reference\": mean.tolist(), \"feature_mean\": mean.tolist(),\n", " \"feature_scale\": scale.tolist(), \"n_training_rows\": len(Z),\n", " \"wavelength_nm\": np.asarray(wavelength_nm)[keep].tolist(), \"detrend\": detrend}\n", "\n", "\n", "def msc(X, reference, *, minimum_slope=1e-12):\n", " \"\"\"Regress each row x = a + b * reference, return (x-a)/b.\n", "\n", " This teaching implementation rejects non-positive or near-zero slopes and\n", " nearly flat references instead of silently generating an unstable result.\n", " It does not guarantee that physical scattering has been identified.\n", " \"\"\"\n", " X, reference = np.asarray(X, float), np.asarray(reference, float)\n", " if X.ndim != 2 or reference.shape != (X.shape[1],) or not np.isfinite(X).all() or not np.isfinite(reference).all():\n", " raise ValueError(\"MSC needs finite rows and a matching finite training reference.\")\n", " centred_reference = reference - reference.mean()\n", " denominator = centred_reference @ centred_reference\n", " if denominator <= 1e-24:\n", " raise ValueError(\"MSC reference is flat or near-flat.\")\n", " slope = ((X - X.mean(axis=1, keepdims=True)) @ centred_reference) / denominator\n", " if np.any(slope <= minimum_slope):\n", " raise ValueError(\"MSC slope is non-positive or near zero; inspect the spectrum/reference.\")\n", " intercept = X.mean(axis=1) - slope * reference.mean()\n", " return (X - intercept[:, None]) / slope[:, None]\n", "\n", "\n", "def normalise_rows(X, method=\"none\", *, statistics=None):\n", " \"\"\"Explicit axis semantics: six row/reference transforms plus train z-score.\"\"\"\n", " X = np.asarray(X, float)\n", " if X.ndim != 2 or X.shape[1] < 2 or not np.isfinite(X).all():\n", " raise ValueError(\"Normalisation needs finite retained spectra.\")\n", " if method not in NORMALISATIONS:\n", " raise ValueError(\"Unknown normalisation method.\")\n", " if method == \"none\":\n", " return X.copy()\n", " if method == \"snv\":\n", " return snv(X)\n", " if method in (\"zscore\", \"minmax\", \"l2\"):\n", " if method == \"zscore\":\n", " centre, scale = X.mean(1, keepdims=True), X.std(1, ddof=0, keepdims=True)\n", " elif method == \"minmax\":\n", " centre = X.min(1, keepdims=True)\n", " scale = X.max(1, keepdims=True) - centre\n", " else:\n", " centre, scale = 0.0, np.linalg.norm(X, axis=1, keepdims=True)\n", " if np.any(scale <= 1e-12):\n", " raise ValueError(f\"{method} has a zero or near-zero denominator.\")\n", " return (X-centre)/scale\n", " if statistics is None:\n", " raise ValueError(f\"{method} requires explicit training-fitted statistics; never fit on test rows.\")\n", " if method == \"msc\":\n", " return msc(X, statistics[\"reference\"])\n", " mean, scale = np.asarray(statistics[\"feature_mean\"], float), np.asarray(statistics[\"feature_scale\"], float)\n", " if mean.shape != (X.shape[1],) or scale.shape != mean.shape or not np.isfinite(mean).all() or not np.isfinite(scale).all() or np.any(scale <= 0):\n", " raise ValueError(\"Invalid training feature statistics.\")\n", " return (X-mean)/scale\n", "\n", "\n", "def preprocess(X, wavelength_nm, keep, *, normalise=False, normalisation=None,\n", " detrend=0, statistics=None, operation=\"smooth\", window=11, polyorder=2):\n", " \"\"\"Mask \u2192 detrend \u2192 normalise/scatter-correct \u2192 segment-wise filter.\n", "\n", " Input retains the original complete uniform grid; output preserves rejected\n", " bands as NaN. Individual normalisation formulae do not require a uniform grid,\n", " but this combined teaching workflow intentionally enforces one to protect its\n", " filters. Detrending and local filters never cross a missing gap. Row normalisation\n", " uses all retained bands together; detrending/filtering use each run.\n", " \"\"\"\n", " X, w, _ = canonicalise(X, wavelength_nm)\n", " delta = uniform_step(w)\n", " keep = np.asarray(keep)\n", " if keep.shape != w.shape:\n", " raise ValueError(\"Band mask and wavelengths must have the same shape.\")\n", " runs = retained_segments(keep)\n", " method = normalisation if normalisation is not None else (\"snv\" if normalise else \"none\")\n", " if normalise and method != \"snv\":\n", " raise ValueError(\"Use either the old normalise=True SNV flag or a named method, not conflicting values.\")\n", " if operation not in OPERATIONS:\n", " raise ValueError(\"Unknown filtering operation.\")\n", " if not isinstance(window, (int, np.integer)) or window % 2 != 1 or window < 3:\n", " raise ValueError(\"This teaching implementation requires an odd window >= 3.\")\n", " if not isinstance(polyorder, (int, np.integer)) or not 0 <= polyorder < window:\n", " raise ValueError(\"Polynomial order must be an integer in [0, window).\")\n", " derivative = {\"derivative\":1,\"derivative2\":2}.get(operation,0)\n", " if derivative > polyorder:\n", " raise ValueError(\"The derivative order cannot exceed the polynomial degree.\")\n", " if method in (\"msc\", \"feature_zscore\"):\n", " if statistics is None:\n", " raise ValueError(f\"{method} requires explicit training-fitted statistics.\")\n", " schema = np.asarray(statistics.get(\"wavelength_nm\", []), float)\n", " if schema.shape != w[keep].shape or not np.array_equal(schema, w[keep]) or statistics.get(\"detrend\") != detrend:\n", " raise ValueError(\"Training statistics do not match this wavelength/mask/detrending schema.\")\n", " out = baseline_correct(X, w, keep, detrend)\n", " out[:, keep] = normalise_rows(out[:, keep], method, statistics=statistics)\n", " if operation != \"none\":\n", " for ids in runs:\n", " if len(ids) < window:\n", " raise ValueError(f\"Segment of {len(ids)} bands is shorter than window={window}; reduce the window or exclude the segment explicitly.\")\n", " if operation == \"moving_average\":\n", " before = out[:, ids].copy()\n", " half = window//2\n", " for i, band in enumerate(ids):\n", " # At edges shrink the centred window; no padding or wraparound.\n", " out[:, band] = before[:,max(0,i-half):min(len(ids),i+half+1)].mean(axis=1)\n", " else:\n", " out[:, ids] = savgol_filter(out[:, ids], window, polyorder,\n", " deriv=derivative, delta=delta, axis=1, mode=\"interp\")\n", " return out\n", "\n", "\n", "class SpectralPreprocessor(TransformerMixin, BaseEstimator):\n", " \"\"\"Row-wise transforms; optional train-reference MSC or feature scaling.\"\"\"\n", " def __init__(self, wavelength_nm, keep, normalise=False, normalisation=None,\n", " detrend=0, operation=\"smooth\", window=11):\n", " self.wavelength_nm = wavelength_nm\n", " self.keep = keep\n", " self.normalise = normalise\n", " self.normalisation = normalisation\n", " self.detrend = detrend\n", " self.operation = operation\n", " self.window = window\n", "\n", " def fit(self, X, y=None):\n", " X = np.asarray(X)\n", " method = self.normalisation if self.normalisation is not None else (\"snv\" if self.normalise else \"none\")\n", " self.statistics_ = fit_training_statistics(X, self.wavelength_nm, self.keep,\n", " detrend=self.detrend) if method in (\"msc\", \"feature_zscore\") else None\n", " self._apply(X)\n", " self.n_features_in_ = X.shape[1]\n", " return self\n", "\n", " def _apply(self, X):\n", " return preprocess(X, self.wavelength_nm, self.keep, normalise=self.normalise,\n", " normalisation=self.normalisation, detrend=self.detrend, statistics=self.statistics_,\n", " operation=self.operation, window=self.window)[:, np.asarray(self.keep)]\n", "\n", " def transform(self, X):\n", " check_is_fitted(self, \"n_features_in_\")\n", " if np.asarray(X).ndim != 2 or np.asarray(X).shape[1] != self.n_features_in_:\n", " raise ValueError(\"Prediction spectra must use the training wavelength schema.\")\n", " return self._apply(X)\n", "\n", "\n", "def synthetic_spectra(n_groups=30, repeats=2, seed=20261002):\n", " \"\"\"Known toy curves plus controlled offset/gain/noise; no field measurements.\"\"\"\n", " rng = np.random.default_rng(seed)\n", " w = np.arange(900., 1701., 5.)\n", " groups = np.repeat(np.arange(n_groups), repeats)\n", " latent = rng.uniform(0.35, 1.0, n_groups)[groups]\n", " baseline = .58 + .00008 * (w - 1300)\n", " clean = baseline[None, :] - .12 * latent[:, None] * np.exp(-.5*((w-1190)/28)**2)\n", " clean -= .16 * np.exp(-.5*((w-1510)/48)**2)[None, :]\n", " gains = rng.uniform(.90, 1.10, n_groups)[groups, None]\n", " offsets = rng.uniform(-.035, .035, n_groups)[groups, None]\n", " clean = gains * clean + offsets\n", " X = clean + rng.normal(0, .008, clean.shape)\n", " keep = ~((w >= 1340) & (w <= 1390))\n", " # Deliberately injected corruption, not an atmospheric-band prescription.\n", " X[:, ~keep] += .09 * np.sin(np.arange((~keep).sum()) * 2.2)[None, :]\n", " y = latent # dimensionless teaching target, not composition or abundance\n", " return X, w, keep, clean, y, groups\n", "\n", "\n", "def training_example():\n", " X, w, keep, clean, y, groups = synthetic_spectra()\n", " train, test = next(GroupShuffleSplit(n_splits=1, test_size=.25, random_state=7).split(X, y, groups))\n", " model = Pipeline([\n", " (\"spectral\", SpectralPreprocessor(w, keep, normalise=False)),\n", " (\"scale\", StandardScaler()),\n", " (\"pca\", PCA(n_components=5, svd_solver=\"full\")),\n", " (\"regressor\", Ridge(alpha=1.0)),\n", " ])\n", " model.fit(X[train], y[train])\n", " prediction = model.predict(X[test])\n", " return model, (X, w, keep, y, groups, train, test, prediction)\n", "\n", "\n", "PRESETS = {\n", " \"clean\": {\"gain\":1.0,\"offset\":0.0,\"noise\":0.0,\"curve\":0.0},\n", " \"noise\": {\"gain\":1.0,\"offset\":0.0,\"noise\":0.025,\"curve\":0.0},\n", " \"scatter\": {\"gain\":1.3,\"offset\":0.08,\"noise\":0.003,\"curve\":0.0},\n", " \"baseline\": {\"gain\":1.0,\"offset\":0.0,\"noise\":0.003,\"curve\":0.12},\n", " \"mixed\": {\"gain\":1.15,\"offset\":0.04,\"noise\":0.012,\"curve\":0.06},\n", "}\n", "\n", "\n", "def playground_data():\n", " X,w,keep,clean,y,groups = synthetic_spectra()\n", " train,test = next(GroupShuffleSplit(n_splits=1,test_size=.25,random_state=7).split(X,y,groups))\n", " rng = np.random.default_rng(711)\n", " samples = [{\"label\":f\"Synthetic sample {i+1}\",\"clean\":clean[test[i*2]].tolist(),\n", " \"noise_vector\":rng.normal(size=len(w)).tolist()} for i in range(3)]\n", " corruption = np.zeros_like(w)\n", " corruption[~keep] = .09*np.sin(np.arange((~keep).sum())*2.2)\n", " learned = {str(int(mask)): {str(degree): fit_training_statistics(X[train],w,\n", " keep if mask else np.ones_like(keep),detrend=degree) for degree in (0,1,2)} for mask in (False,True)}\n", " kernels = {str(window): {str(deriv): [savgol_coeffs(window,2,deriv=deriv,\n", " delta=5.,pos=pos,use=\"dot\").tolist() for pos in range(window)]\n", " for deriv in (0,1,2)} for window in (5,11,21)}\n", " return {\"version\":2,\"status\":\"Synthetic teaching data only; no sensor or model benchmark\",\n", " \"seed\":20261002,\"wavelength_nm\":w.tolist(),\"keep\":keep.tolist(),\"samples\":samples,\n", " \"bad_interval_nm\":[1340,1390],\"corruption\":corruption.tolist(),\"delta_nm\":5,\n", " \"polyorder\":2,\"kernels\":kernels,\"presets\":PRESETS,\"training_statistics\":learned,\n", " \"training_groups\":len(set(groups[train])),\"training_rows\":len(train)}\n", "\n", "\n", "def synthetic_observation(data, sample=0, preset=\"mixed\", controls=None):\n", " p = dict(PRESETS[preset] if controls is None else controls)\n", " if not all(np.isfinite(p[k]) for k in (\"gain\",\"offset\",\"noise\",\"curve\")) or p[\"gain\"]<=0 or p[\"noise\"]<0:\n", " raise ValueError(\"Synthetic controls must be finite, with positive gain and non-negative noise.\")\n", " s=data[\"samples\"][sample];w=np.asarray(data[\"wavelength_nm\"])\n", " return p[\"gain\"]*np.asarray(s[\"clean\"])+p[\"offset\"]+p[\"noise\"]*np.asarray(s[\"noise_vector\"])+p[\"curve\"]*((w-1300)/400)**2+np.asarray(data[\"corruption\"])\n", "\n", "\n", "def export_assets(destination):\n", " dest=Path(destination);dest.mkdir(parents=True,exist_ok=True)\n", " data=playground_data()\n", " (dest/\"hsi-preprocessing-data.json\").write_text(json.dumps(data,separators=(\",\",\":\")))\n", " golden=[]\n", " # Every method/order/filter combination, with representative distinct inputs.\n", " for mask in (False,True):\n", " for degree in (0,1,2):\n", " for method in NORMALISATIONS:\n", " for op in OPERATIONS:\n", " for window in (5,11,21):\n", " sample=(degree+window)%3\n", " preset=list(PRESETS)[(window+degree+NORMALISATIONS.index(method)+OPERATIONS.index(op))%len(PRESETS)]\n", " X=synthetic_observation(data,sample,preset)[None,:]\n", " keep=np.asarray(data[\"keep\"]) if mask else np.ones(len(data[\"keep\"]),bool)\n", " out=preprocess(X,data[\"wavelength_nm\"],keep,normalisation=method,\n", " detrend=degree,statistics=data[\"training_statistics\"][str(int(mask))][str(degree)],\n", " operation=op,window=window)[0]\n", " golden.append({\"sample\":sample,\"preset\":preset,\"mask\":mask,\"normalisation\":method,\n", " \"detrend\":degree,\"operation\":op,\"window\":window,\n", " \"values\":[None if not np.isfinite(v) else float(v) for v in out]})\n", " (dest/\"hsi-preprocessing-golden.json\").write_text(json.dumps(golden,separators=(\",\",\":\")))\n", " return len(golden)\n", "\n", "\n", "def main(argv=None):\n", " parser=argparse.ArgumentParser(description=__doc__)\n", " parser.add_argument(\"--export\",type=Path,help=\"Write browser fixtures to this directory\")\n", " parser.add_argument(\"--normalisation\",choices=NORMALISATIONS,default=\"none\")\n", " parser.add_argument(\"--detrend\",type=int,choices=(0,1,2),default=0)\n", " parser.add_argument(\"--operation\",choices=OPERATIONS,default=\"smooth\")\n", " parser.add_argument(\"--window\",type=int,choices=(5,11,21),default=11)\n", " parser.add_argument(\"--preset\",choices=tuple(PRESETS),default=\"mixed\")\n", " parser.add_argument(\"--sample\",type=int,choices=(1,2,3),default=1)\n", " parser.add_argument(\"--gain\",type=float,help=\"Override preset positive multiplicative gain\")\n", " parser.add_argument(\"--offset\",type=float,help=\"Override preset additive reflectance offset\")\n", " parser.add_argument(\"--noise\",type=float,help=\"Override preset noise standard deviation\")\n", " parser.add_argument(\"--curve\",type=float,help=\"Override preset quadratic curvature amplitude\")\n", " parser.add_argument(\"--include-flagged\",action=\"store_true\")\n", " parser.add_argument(\"--output\",type=Path,help=\"Write displayed synthetic spectrum as CSV\")\n", " parser.add_argument(\"--self-test\",action=\"store_true\",help=\"Run the companion unittest suite\")\n", " args=parser.parse_args(argv)\n", " if args.self_test:\n", " import unittest\n", " suite=unittest.defaultTestLoader.discover(str(Path(__file__).parent),pattern=\"test_preprocessing.py\")\n", " if not unittest.TextTestRunner(verbosity=2).run(suite).wasSuccessful():\n", " return 1\n", " data=playground_data();w=np.asarray(data[\"wavelength_nm\"])\n", " keep=np.asarray(data[\"keep\"]) if not args.include_flagged else np.ones(len(w),bool)\n", " synthetic_controls=dict(PRESETS[args.preset])\n", " for key in (\"gain\",\"offset\",\"noise\",\"curve\"):\n", " if getattr(args,key) is not None: synthetic_controls[key]=getattr(args,key)\n", " X=synthetic_observation(data,args.sample-1,args.preset,synthetic_controls)\n", " result=preprocess(X[None,:],w,keep,normalisation=args.normalisation,detrend=args.detrend,\n", " statistics=data[\"training_statistics\"][str(int(not args.include_flagged))][str(args.detrend)],\n", " operation=args.operation,window=args.window)[0]\n", " model,(_,_,_,_,groups,train,test,prediction)=training_example()\n", " print(f\"Synthetic CPU example: {len(w)} wavelengths, {keep.sum()} retained\")\n", " print(f\"Recipe: mask -> detrend({args.detrend}) -> {args.normalisation} -> {args.operation}({args.window})\")\n", " print(f\"Training reference: {data['training_rows']} rows / {data['training_groups']} independent synthetic groups; never refitted by presets\")\n", " print(f\"Group-disjoint train/test observations: {len(train)}/{len(test)}; finite predictions: {np.isfinite(prediction).all()}\")\n", " print(\"No real-world performance claim. Flagged interval is injected corruption, not a universal mask.\")\n", " if args.output:\n", " args.output.parent.mkdir(parents=True,exist_ok=True)\n", " np.savetxt(args.output,np.column_stack([w,X,~np.asarray(data['keep']),keep,result]),delimiter=\",\",\n", " header=\"wavelength_nm,original_reflectance_fraction,flagged_band,used_band,processed_value\",comments=\"\")\n", " unit = \"reflectance fraction\" if args.normalisation in (\"none\",\"msc\") else \"dimensionless\"\n", " if args.detrend and args.normalisation in (\"none\",\"msc\"):\n", " unit = \"reflectance residual fraction\"\n", " if args.operation == \"derivative\": unit += \" per nm\"\n", " if args.operation == \"derivative2\": unit += \" per nm^2\"\n", " metadata = {\"data_status\":\"synthetic teaching example\",\"preset\":args.preset,\n", " \"sample\":args.sample,\"synthetic_controls\":synthetic_controls,\"normalisation\":args.normalisation,\"detrend\":args.detrend,\n", " \"operation\":args.operation,\"window_bands\":args.window,\"delta_nm\":5,\n", " \"processed_unit\":unit,\"include_flagged\":args.include_flagged,\n", " \"order\":[\"mask\",\"detrend\",\"normalise/scatter-correct\",\"filter\"],\n", " \"training_rows\":data[\"training_rows\"],\"training_groups\":data[\"training_groups\"],\n", " \"seed\":data[\"seed\"],\"missing_output\":\"NaN marks excluded bands\"}\n", " args.output.with_suffix(\".metadata.json\").write_text(json.dumps(metadata,indent=2)+\"\\n\")\n", " print(f\"Wrote {args.output} and recipe metadata; output units: {unit}\")\n", " if args.export:\n", " print(f\"Exported {export_assets(args.export)} Python/JavaScript reference cases\")\n", " return 0\n", "\n" ], "execution_count": 1, "outputs": [], "id": "hsi-preprocessing-01" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Generate the controlled teaching data\n", "\n", "Two repeat observations belong to each synthetic sample group. The 1340\u20131390 nm corruption is injected deliberately; it is not a recommended exclusion interval for a real sensor." ], "id": "hsi-preprocessing-02" }, { "cell_type": "code", "metadata": {}, "source": [ "X, wavelength_nm, keep, clean, target, groups = synthetic_spectra()\n", "print(\"Shape:\", X.shape)\n", "print(\"Wavelength range and step (nm):\", wavelength_nm[0], wavelength_nm[-1], uniform_step(wavelength_nm))\n", "print(\"Retained / excluded bands:\", int(keep.sum()), int((~keep).sum()))\n", "print(\"Valid runs (nm):\", [(wavelength_nm[r[0]], wavelength_nm[r[-1]]) for r in retained_segments(keep)])" ], "execution_count": 2, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Shape: (60, 161)\n", "Wavelength range and step (nm): 900.0 1700.0 5.0\n", "Retained / excluded bands: 150 11\n", "Valid runs (nm): [(np.float64(900.0), np.float64(1335.0)), (np.float64(1395.0), np.float64(1700.0))]\n" ] } ], "id": "hsi-preprocessing-03" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. Verify that compressing a gap is rejected\n", "\n", "The step from 1335 to 1395 nm is 60 nm, not 5 nm. Catch the expected error; do not substitute a median spacing." ], "id": "hsi-preprocessing-04" }, { "cell_type": "code", "metadata": {}, "source": [ "try:\n", " uniform_step(wavelength_nm[keep])\n", "except ValueError as error:\n", " print(\"Expected validation error:\", error)\n", "else:\n", " raise AssertionError(\"Deleted-band gap was not detected\")" ], "execution_count": 3, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Expected validation error: Non-uniform wavelength spacing or deleted-band gap: use separate contiguous segments; do not smooth across the gap.\n" ] } ], "id": "hsi-preprocessing-05" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4. Compare SNV, smoothing and a first derivative\n", "\n", "Mask \u2192 SNV \u2192 smoothing is a different recipe from mask \u2192 smoothing \u2192 SNV. First derivatives use delta=5 nm. Excluded output bands remain NaN. No-data NaNs in original input are deliberately rejected by the minimal reader contract; use a documented real-data policy first." ], "id": "hsi-preprocessing-06" }, { "cell_type": "code", "metadata": {}, "source": [ "smooth = preprocess(X, wavelength_nm, keep, operation=\"smooth\")\n", "normalised = preprocess(X, wavelength_nm, keep, normalise=True, operation=\"none\")\n", "derivative = preprocess(X, wavelength_nm, keep, operation=\"derivative\")\n", "np.testing.assert_allclose(normalised[:, keep].mean(axis=1), 0, atol=1e-13)\n", "np.testing.assert_allclose(normalised[:, keep].std(axis=1, ddof=1), 1, atol=1e-13)\n", "assert np.isnan(smooth[:, ~keep]).all()\n", "print(\"First retained value: raw / smooth / SNV / dR per nm\")\n", "print(X[0,0], smooth[0,0], normalised[0,0], derivative[0,0])\n", "print(\"SNV row invariants and explicit masked gaps verified\")" ], "execution_count": 4, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "First retained value: raw / smooth / SNV / dR per nm\n", "0.571614933186864 0.5631569858442029 0.08745685448918786 0.0005199068021641358\n", "SNV row invariants and explicit masked gaps verified\n" ] } ], "id": "hsi-preprocessing-07" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 5. Check the derivative against a known function\n", "\n", "This verifies calculus, units and edge handling. It is not a sensor-performance experiment." ], "id": "hsi-preprocessing-08" }, { "cell_type": "code", "metadata": {}, "source": [ "quadratic = (.000003 * (wavelength_nm - 1200)**2 + .0007 * wavelength_nm)[None, :]\n", "expected = (.000006 * (wavelength_nm - 1200) + .0007)[None, :]\n", "actual = preprocess(quadratic, wavelength_nm, keep, operation=\"derivative\")\n", "np.testing.assert_allclose(actual[:, keep], expected[:, keep], atol=2e-14)\n", "print(\"Maximum derivative error:\", np.max(np.abs(actual[:, keep] - expected[:, keep])))" ], "execution_count": 5, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Maximum derivative error: 1.0581813203458523e-16\n" ] } ], "id": "hsi-preprocessing-09" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 6. Fit scaler, PCA and model only on training groups\n", "\n", "The target is an arbitrary parameter controlling a synthetic band depth, not a chemical concentration or a mixing fraction. No accuracy score is presented. The purpose is a correctly scoped pipeline." ], "id": "hsi-preprocessing-10" }, { "cell_type": "code", "metadata": {}, "source": [ "model, example = training_example()\n", "X, wavelength_nm, keep, target, groups, train, test, prediction = example\n", "assert not (set(groups[train]) & set(groups[test]))\n", "training_features = model.named_steps[\"spectral\"].transform(X[train])\n", "np.testing.assert_allclose(model.named_steps[\"scale\"].mean_, training_features.mean(axis=0))\n", "scaled_training = model.named_steps[\"scale\"].transform(training_features)\n", "np.testing.assert_allclose(model.named_steps[\"pca\"].mean_, scaled_training.mean(axis=0), atol=1e-13)\n", "print(\"Group-disjoint train/test observations:\", len(train), len(test))\n", "print(\"Finite held-out predictions:\", bool(np.isfinite(prediction).all()))\n", "print(\"Training-only scaler and PCA statistics verified\")" ], "execution_count": 6, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Group-disjoint train/test observations: 44 16\n", "Finite held-out predictions: True\n", "Training-only scaler and PCA statistics verified\n" ] } ], "id": "hsi-preprocessing-11" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 7. Refit within each cross-validation fold\n", "\n", "Keep all related observations together. Tuning preprocessing belongs inside the appropriate validation loop, never on a final test set." ], "id": "hsi-preprocessing-12" }, { "cell_type": "code", "metadata": {}, "source": [ "from sklearn.model_selection import GroupKFold, cross_validate\n", "cv = GroupKFold(n_splits=3)\n", "result = cross_validate(model, X, target, groups=groups, cv=cv, return_estimator=True)\n", "for fitted, (training_ids, validation_ids) in zip(result[\"estimator\"], cv.split(X, target, groups)):\n", " Z = fitted.named_steps[\"spectral\"].transform(X[training_ids])\n", " np.testing.assert_allclose(fitted.named_steps[\"scale\"].mean_, Z.mean(axis=0))\n", "print(\"All 3 grouped folds independently fitted their preprocessing\")" ], "execution_count": 7, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "All 3 grouped folds independently fitted their preprocessing\n" ] } ], "id": "hsi-preprocessing-13" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 8. Compare the complete normalisation family\n", "\n", "The formulas act across different axes. These checks inspect output invariants immediately after normalisation, with no later filter." ], "id": "hsi-preprocessing-14" }, { "cell_type": "code", "metadata": {}, "source": [ "training_statistics = fit_training_statistics(X[train], wavelength_nm, keep)\n", "for method in NORMALISATIONS:\n", " transformed = preprocess(X[test], wavelength_nm, keep,\n", " normalisation=method, statistics=training_statistics, operation=\"none\")[:, keep]\n", " print(method, \"first-row mean:\", round(float(transformed[0].mean()), 6),\n", " \"range:\", round(float(transformed[0].min()), 6), round(float(transformed[0].max()), 6))\n", "print(\"The feature-wise z-score does not standardise each row to zero mean.\")" ], "execution_count": 8, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "none first-row mean: 0.538651 range: 0.416989 0.61549\n", "minmax first-row mean: 0.612904 range: 0.0 1.0\n", "zscore first-row mean: -0.0 range: -2.924322 1.846933\n", "snv first-row mean: -0.0 range: -2.914558 1.840767\n", "l2 first-row mean: 0.081407 range: 0.06302 0.09302\n", "msc first-row mean: 0.544481 range: 0.410629 0.62902\n", "feature_zscore first-row mean: -0.171715 range: -0.963317 1.024716\n", "The feature-wise z-score does not standardise each row to zero mean.\n" ] } ], "id": "hsi-preprocessing-15" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 9. Keep the MSC reference frozen at prediction time\n", "\n", "MSC regresses each spectrum against a reference learned only from training rows. The constant/near-zero denominator and non-positive slope guards are part of this teaching implementation." ], "id": "hsi-preprocessing-16" }, { "cell_type": "code", "metadata": {}, "source": [ "msc_transformer = SpectralPreprocessor(wavelength_nm, keep,\n", " normalisation=\"msc\", detrend=1, operation=\"none\").fit(X[train])\n", "reference_before = np.array(msc_transformer.statistics_[\"reference\"])\n", "msc_result = msc_transformer.transform(X[test])\n", "np.testing.assert_array_equal(reference_before, msc_transformer.statistics_[\"reference\"])\n", "reference = reference_before\n", "perfect_affine = 1.3 * reference[None, :] + .07\n", "np.testing.assert_allclose(msc(perfect_affine, reference), reference[None, :], atol=1e-13)\n", "print(\"MSC reference fixed during prediction; exact positive affine example recovered\")" ], "execution_count": 9, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "MSC reference fixed during prediction; exact positive affine example recovered\n" ] } ], "id": "hsi-preprocessing-17" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 10. Test detrending, moving average and the second derivative\n", "\n", "Broad fitted baselines can also remove genuine signal. A second derivative has units per wavelength squared. The moving-average edge window shrinks; it does not cross a gap." ], "id": "hsi-preprocessing-18" }, { "cell_type": "code", "metadata": {}, "source": [ "polynomial = (.5 + .0002 * wavelength_nm + 1e-7 * wavelength_nm**2)[None, :]\n", "residual = baseline_correct(polynomial, wavelength_nm, keep, degree=2)\n", "np.testing.assert_allclose(residual[:, keep], 0, atol=1e-13)\n", "second = preprocess(polynomial, wavelength_nm, keep, operation=\"derivative2\")\n", "np.testing.assert_allclose(second[:, keep], 2e-7, atol=3e-14)\n", "averaged = preprocess(X[:1], wavelength_nm, keep, operation=\"moving_average\", window=5)\n", "np.testing.assert_allclose(averaged[0, 0], X[0, :3].mean())\n", "print(\"Quadratic detrend, second derivative per nm\u00b2 and truncated moving-average edge verified\")" ], "execution_count": 10, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Quadratic detrend, second derivative per nm\u00b2 and truncated moving-average edge verified\n" ] } ], "id": "hsi-preprocessing-19" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 11. Reproduce a playground preset or custom controls\n", "\n", "The browser and Python use the same seeded noise vector, independent training reference and recipe. No hardware data or account is required." ], "id": "hsi-preprocessing-20" }, { "cell_type": "code", "metadata": {}, "source": [ "playground = playground_data()\n", "current = synthetic_observation(playground, sample=0, preset=\"scatter\")\n", "reference_stats = playground[\"training_statistics\"][\"1\"][\"0\"]\n", "processed = preprocess(current[None, :], playground[\"wavelength_nm\"],\n", " np.asarray(playground[\"keep\"]), normalisation=\"msc\", statistics=reference_stats,\n", " operation=\"derivative2\", window=11)\n", "print(\"Preset:\", \"scatter\", \"; recipe: mask \u2192 MSC \u2192 second derivative\")\n", "print(\"Fixed training rows:\", playground[\"training_rows\"], \"; retained output values:\", int(np.isfinite(processed).sum()))\n", "print(\"CLI: python hsi_preprocessing.py --preset scatter --normalisation msc --operation derivative2 --output example.csv\")" ], "execution_count": 11, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Preset: scatter ; recipe: mask \u2192 MSC \u2192 second derivative\n", "Fixed training rows: 44 ; retained output values: 150\n", "CLI: python hsi_preprocessing.py --preset scatter --normalisation msc --operation derivative2 --output example.csv\n" ] } ], "id": "hsi-preprocessing-21" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Further checks and reading\n", "\n", "The companion script bundle includes 25 unit tests, covering reordered metadata, missing data, deleted-band gaps, SNV invariants, derivative units, corrupt-band isolation, pipeline cloning and train-only fitting. Run `python -m unittest -v test_preprocessing.py` after extracting that bundle.\n", "\n", "- [SciPy savgol_filter](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.savgol_filter.html)\n", "- [Barnes, Dhanoa and Lister (1989), SNV](https://doi.org/10.1366/0003702894202201)\n", "- [scikit-learn preprocessing pitfalls](https://scikit-learn.org/stable/common_pitfalls.html)\n", "- [scikit-learn Pipeline](https://scikit-learn.org/stable/modules/generated/sklearn.pipeline.Pipeline.html)\n", "\n", "Record the original measurement representation, units, masks, transform order, filter settings, split identifiers and dependency versions when adapting the workflow. Never copy the synthetic exclusion interval as a universal sensor rule." ], "id": "hsi-preprocessing-22" } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.12.14" } }, "nbformat": 4, "nbformat_minor": 5 }