{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Sensor-response integration before model adaptation\n",
        "\n",
        "A self-contained CPU notebook. All spectra are original analytic constructions; all source and target responses are idealized Gaussians. No DOFA/SpectralEarth model is loaded or evaluated. Python 3.10+ standard library only. JupyterLab is optional.\n",
        "\n",
        "The task: predict target band averages from already-integrated source measurements, then compare against synthetic truth available only because we constructed it."
      ],
      "id": "sensor-00"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 1. Define the numerical model\n",
        "\n",
        "The code below is embedded so this notebook works as a standalone download. FWHM and centres are in nm. Responses are normalized and truncated at \u00b14\u03c3; composite midpoint cells are at most 0.5 nm. The 99.9% coverage threshold is an explicit numerical choice, not confidence. No source extrapolation or interpolation across invalid neighbours is allowed."
      ],
      "id": "sensor-01"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "\"\"\"Synthetic spectral-response lab. Python >=3.10, standard library only; CPU.\"\"\"\n",
        "import csv\n",
        "import json\n",
        "import math\n",
        "from pathlib import Path\n",
        "\n",
        "DEFAULTS = dict(sourceFwhm=8, spacing=5, targetFwhm=80, featureSigma=8, shift=0, coverage=\"full\")\n",
        "PRESETS = {\"baseline\": {}, \"centre\": {\"targetFwhm\":120}, \"broad\": {\"sourceFwhm\":100,\"spacing\":20,\"targetFwhm\":12}, \"sparse\":{\"spacing\":60,\"targetFwhm\":20}, \"gap\":{\"coverage\":\"gap\"}, \"shift\":{\"shift\":20,\"targetFwhm\":20}}\n",
        "SQ = 2 * math.sqrt(2 * math.log(2))\n",
        "\n",
        "def truth(x, sigma=8):\n",
        "    return .48 + .00004*(x-2150) - .22*math.exp(-.5*((x-2200)/sigma)**2) - .055*math.exp(-.5*((x-2040)/45)**2)\n",
        "\n",
        "def integrate(fn, centre, fwhm, step=.5):\n",
        "    \"\"\"Normalized midpoint quadrature. Null for <99.9% valid response coverage.\n",
        "\n",
        "    Model assumption: nonnegative Gaussian SRF truncated at +/-4 sigma.\n",
        "    Invalid fn returns None or a nonfinite value. No extrapolation is performed.\n",
        "    Coverage is geometric/modelled response support, NOT statistical confidence.\n",
        "    \"\"\"\n",
        "    if not all(math.isfinite(v) for v in (centre, fwhm, step)) or fwhm <= 0 or step <= 0:\n",
        "        raise ValueError(\"Finite centre, positive FWHM and step required\")\n",
        "    sigma = fwhm/SQ\n",
        "    a,b = centre-4*sigma, centre+4*sigma\n",
        "    n = math.ceil((b-a)/step)\n",
        "    dx = (b-a)/n\n",
        "    total=known=num=0.\n",
        "    for i in range(n):\n",
        "        x=a+(i+.5)*dx\n",
        "        w=math.exp(-.5*((x-centre)/sigma)**2)*dx\n",
        "        y=fn(x)\n",
        "        total+=w\n",
        "        if y is not None and math.isfinite(y):\n",
        "            known+=w\n",
        "            num+=w*y\n",
        "    coverage=known/total\n",
        "    return {\"value\": num/known if coverage >= .999 else None, \"coverage\":coverage}\n",
        "\n",
        "def interpolate(xs, ys, x):\n",
        "    \"\"\"Piecewise linear reconstruction of source band AVERAGES, not point truth.\"\"\"\n",
        "    if not math.isfinite(x) or x<xs[0] or x>xs[-1]:\n",
        "        return None\n",
        "    i=0\n",
        "    while i<len(xs)-2 and xs[i+1]<x:\n",
        "        i+=1\n",
        "    if x==xs[i]: return ys[i]\n",
        "    if x==xs[i+1]: return ys[i+1]\n",
        "    if ys[i] is None or ys[i+1] is None: return None\n",
        "    t=(x-xs[i])/(xs[i+1]-xs[i])\n",
        "    return ys[i]*(1-t)+ys[i+1]*t\n",
        "\n",
        "def validate(state):\n",
        "    for k,lo,hi in [(\"sourceFwhm\",4,100),(\"targetFwhm\",8,120),(\"featureSigma\",3,40),(\"shift\",-20,20)]:\n",
        "        v=state[k]\n",
        "        if not isinstance(v,(int,float)) or not math.isfinite(v) or not lo<=v<=hi or v!=int(v):\n",
        "            raise ValueError(f\"Invalid {k}\")\n",
        "    if state['spacing'] not in [5,20,60] or state['coverage'] not in ['full','gap','short']:\n",
        "        raise ValueError('Invalid spacing or coverage')\n",
        "\n",
        "def run(**kwargs):\n",
        "    state={**DEFAULTS, **kwargs}; validate(state)\n",
        "    xs=list(range(1800,2501,state['spacing']))\n",
        "    ys=[]\n",
        "    for c in xs:\n",
        "        invalid=(state['coverage']=='gap' and 2150<=c<=2250) or (state['coverage']=='short' and c>2200)\n",
        "        ys.append(None if invalid else integrate(lambda x:truth(x,state['featureSigma']),c,state['sourceFwhm'])['value'])\n",
        "    recon=lambda x:interpolate(xs,ys,x)\n",
        "    rows=[]\n",
        "    for c in range(2000,2301,25):\n",
        "        ref=integrate(lambda x:truth(x,state['featureSigma']),c+state['shift'],state['targetFwhm'])\n",
        "        est=integrate(recon,c,state['targetFwhm'])\n",
        "        rows.append(dict(nominalNm=c,trueCentreNm=c+state['shift'],reference=ref['value'],idealCentre=truth(c,state['featureSigma']),centre=recon(c),resampled=est['value'],coverage=est['coverage']))\n",
        "    common=[r for r in rows if r['centre'] is not None and r['resampled'] is not None]\n",
        "    metrics={}\n",
        "    for key in ['centre','resampled','idealCentre']:\n",
        "        e=[abs(r[key]-r['reference']) for r in common]\n",
        "        metrics[key]=dict(n=len(e),mae=sum(e)/len(e) if e else None,max=max(e) if e else None)\n",
        "    return dict(state=state,xs=xs,ys=ys,rows=rows,metrics=metrics)\n",
        "\n",
        "def write_outputs(directory='outputs'):\n",
        "    p=Path(directory);p.mkdir(parents=True,exist_ok=True)\n",
        "    results={name:run(**settings) for name,settings in PRESETS.items()}\n",
        "    (p/'results.json').write_text(json.dumps(results,indent=2))\n",
        "    for name,result in results.items():\n",
        "        with (p/f'{name}.csv').open('w',newline='') as f:\n",
        "            writer=csv.DictWriter(f,fieldnames=list(result['rows'][0]));writer.writeheader();writer.writerows(result['rows'])\n",
        "    return results\n",
        "\n"
      ],
      "execution_count": null,
      "outputs": [],
      "id": "sensor-02"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 2. Check units, constants and support\n",
        "\n",
        "A constant spectrum must remain constant. Expressing the same wavelengths in micrometres must give the same integral after converting the function and integration step. Half-supported bands and gaps are rejected."
      ],
      "id": "sensor-03"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "assert abs(integrate(lambda x: 0.37, 2200, 80)[\"value\"] - 0.37) < 1e-12\n",
        "a = integrate(truth, 2200, 80)[\"value\"]\n",
        "b = integrate(lambda um: truth(1000*um), 2.2, 0.08, step=0.0005)[\"value\"]\n",
        "assert abs(a-b) < 1e-9\n",
        "assert integrate(lambda x: .4 if x < 2200 else None, 2200, 80)[\"value\"] is None\n",
        "assert interpolate([2100,2200,2300], [.4,None,.5], 2150) is None\n",
        "print(\"Units, constant preservation, missing-support checks passed\")"
      ],
      "execution_count": null,
      "outputs": [],
      "id": "sensor-04"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 3. Run all six controlled experiments\n",
        "\n",
        "MAE is measured in reflectance fraction against the known synthetic target. All metrics within a run use the same eligible target bands. Do not compare the two-band gap subset as if it were the full thirteen-band task."
      ],
      "id": "sensor-05"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "results = {name: run(**settings) for name,settings in PRESETS.items()}\n",
        "for name, result in results.items():\n",
        "    m=result[\"metrics\"]\n",
        "    print(name, \"retained:\",m[\"resampled\"][\"n\"], \"response MAE:\",m[\"resampled\"][\"mae\"], \"centre MAE:\",m[\"centre\"][\"mae\"])\n",
        "assert results[\"baseline\"][\"metrics\"][\"resampled\"][\"mae\"] < 0.0002\n",
        "assert results[\"broad\"][\"metrics\"][\"resampled\"][\"mae\"] > 0.02\n",
        "assert results[\"sparse\"][\"metrics\"][\"resampled\"][\"max\"] > 0.13\n",
        "assert results[\"gap\"][\"metrics\"][\"resampled\"][\"n\"] == 2"
      ],
      "execution_count": null,
      "outputs": [],
      "id": "sensor-06"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 4. Inspect the 2,200 nm counterexample\n",
        "\n",
        "The broad-source preset loses the narrow trough before target conversion. Integrating a reconstructed source signal is not inversion of the source sensor. The shift preset deliberately keeps nominal metadata while moving the true target centre."
      ],
      "id": "sensor-07"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "for name in [\"baseline\", \"broad\", \"sparse\", \"gap\", \"shift\"]:\n",
        "    row=next(r for r in results[name][\"rows\"] if r[\"nominalNm\"]==2200)\n",
        "    print(name, row)"
      ],
      "execution_count": null,
      "outputs": [],
      "id": "sensor-08"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 5. Sweep width and check numerical convergence\n",
        "\n",
        "This sweep changes target FWHM with everything else fixed. It is a sensitivity experiment, not a statistical interval. Refining quadrature checks numerical convergence under the assumed model; it does not validate sensor physics."
      ],
      "id": "sensor-09"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "for fwhm in [8,20,40,80,120]:\n",
        "    r=run(targetFwhm=fwhm)\n",
        "    print(\"Target FWHM:\",fwhm,\"MAE:\",r[\"metrics\"][\"resampled\"][\"mae\"])\n",
        "for width in [4,8,12,120]:\n",
        "    coarse=integrate(truth,2200,width)[\"value\"]\n",
        "    fine=integrate(truth,2200,width,step=0.1)[\"value\"]\n",
        "    assert abs(coarse-fine)<1e-6\n",
        "print(\"Quadrature convergence checks passed\")"
      ],
      "execution_count": null,
      "outputs": [],
      "id": "sensor-10"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 6. Save reproducible outputs\n",
        "\n",
        "Outputs include all settings, source measurements, target values and per-band modeled coverage. Rejected values are JSON null or empty CSV fields, never zero-filled. No measured sensor files or pretrained model weights are downloaded."
      ],
      "id": "sensor-11"
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "saved=write_outputs(\"outputs\")\n",
        "print(\"Saved\",len(saved),\"presets to outputs/results.json and six CSV files\")"
      ],
      "execution_count": null,
      "outputs": [],
      "id": "sensor-12"
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Where this ends\n",
        "\n",
        "This notebook demonstrates measurement operators only. It does not train an adapter or support a cross-sensor accuracy claim. For model adaptation: pin checkpoint/pretraining provenance, split independent scenes or sites before fitting, fit normalization and adapters only on training data, select on validation data and evaluate once on untouched test units. Frozen-backbone protocols must also specify batch-normalization buffer behavior.\n",
        "\n",
        "Primary resources:\n",
        "- [DOFA v3](https://arxiv.org/html/2403.15356v3)\n",
        "- [SpectralEarth v2](https://arxiv.org/html/2408.08447v2)\n",
        "- [EnMAP-Box custom sensor responses](https://enmap-box.readthedocs.io/en/latest/usr_section/usr_manual/processing_algorithms/spectral_resampling/spectral_resampling_to_custom_sensor.html)"
      ],
      "id": "sensor-13"
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "name": "python",
      "version": "3.10"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}