{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# HSI dimensionality reduction: runnable CPU experiment\n",
        "\n",
        "Original synthetic spectra, not measured HSI. Python 3.12 was used for verification. Install the pinned requirements from the accompanying source bundle in your own environment. No GPU, API key or external dataset is needed.\n",
        "\n",
        "This notebook contains the generator and full NumPy autoencoder implementation, so it does not depend on a local Python helper file. PCA and regularised Fisher LDA use training rows only. t-SNE and UMAP below fit all 120 rows and are explicitly exploratory. Class labels annotate unsupervised plots; only LDA receives training labels. Repeatedly examining held-out MSE turns that set into development data."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "execution_count": null,
      "outputs": [],
      "source": [
        "# Install in a fresh environment before running; uncomment only if needed.\n",
        "# %pip install numpy==2.3.5 scipy==1.17.0 scikit-learn==1.8.0 umap-learn==0.5.9.post2 numba==0.63.1 llvmlite==0.46.0 pynndescent==0.5.13 tqdm==4.67.1 matplotlib==3.10.8\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "execution_count": null,
      "outputs": [],
      "source": [
        "#!/usr/bin/env python3\n",
        "\"\"\"Original synthetic spectral example. CPU only. Regenerates real embeddings/fixtures.\n",
        "python python/reproduce.py --output learning --presets\n",
        "\"\"\"\n",
        "import argparse, json, time, platform\n",
        "from pathlib import Path\n",
        "import numpy as np\n",
        "from scipy.linalg import eigh\n",
        "from sklearn.decomposition import PCA\n",
        "from sklearn.manifold import TSNE\n",
        "import sklearn\n",
        "\n",
        "def dataset(seed=314159):\n",
        "    rng=np.random.default_rng(seed); wave=np.linspace(450,1000,48); X=[]; y=[]; train=[]; test=[]\n",
        "    for c in range(4):\n",
        "        for k in range(30):\n",
        "            t=rng.uniform(-1,1); bright=rng.normal(0,.075); slope=rng.normal(0,.025)\n",
        "            base=.42+.07*np.sin((wave-450)/550*np.pi)+bright+slope*(wave-725)/275\n",
        "            trough=-.085*np.exp(-.5*((wave-[550,670,810,930][c])/28)**2)\n",
        "            latent=.035*t*np.cos(wave/91)+.023*(t*t-.33)*np.sin(wave/43)\n",
        "            X.append(np.clip(base+trough+latent+rng.normal(0,.006,48),.02,.98));y.append(c)\n",
        "            (train if k<20 else test).append(len(y)-1)\n",
        "    return {'X':np.asarray(X),'labels':np.asarray(y),'train':np.asarray(train),'test':np.asarray(test),'wavelengths':wave,'classes':['A \u00b7 550 nm dip','B \u00b7 670 nm dip','C \u00b7 810 nm dip','D \u00b7 930 nm dip'],'seed':seed}\n",
        "\n",
        "def prepare(X,train,scaling):\n",
        "    mean=X[train].mean(axis=0); scale=X[train].std(axis=0) if scaling=='standard' else np.ones(X.shape[1]);scale[scale==0]=1\n",
        "    return (X-mean)/scale,mean,scale\n",
        "\n",
        "def fisher(X,labels,train,dims=2,ridge=.1):\n",
        "    \"\"\"Regularised Fisher LDA: Sb v = lambda (Sw + ridge*tr(Sw)/p I) v.\n",
        "    Biased covariance (divide by training N); labels outside train are never read.\n",
        "    \"\"\"\n",
        "    Xtr=X[train]; ytr=labels[train]; mean=Xtr.mean(0); p=X.shape[1]; sw=np.zeros((p,p));sb=np.zeros((p,p))\n",
        "    classes=np.unique(ytr)\n",
        "    if not 1<=dims<=min(p,len(classes)-1):raise ValueError('LDA dimension exceeds C-1')\n",
        "    for c in classes:\n",
        "        group=Xtr[ytr==c]; mu=group.mean(0); centered=group-mu\n",
        "        sw+=centered.T@centered/len(Xtr);d=mu-mean;sb+=len(group)/len(Xtr)*np.outer(d,d)\n",
        "    sw+=ridge*max(np.trace(sw)/p,1e-12)*np.eye(p)\n",
        "    vals,vecs=eigh(sb,sw); order=np.argsort(vals)[::-1];basis=vecs[:,order[:dims]]\n",
        "    return (X-mean)@basis, vals[order[:len(classes)-1]]\n",
        "\n",
        "def mulberry(seed):\n",
        "    a=seed\n",
        "    def rand():\n",
        "        nonlocal a\n",
        "        a=(a+0x6D2B79F5)&0xffffffff;t=a\n",
        "        t=((t^(t>>15))*(t|1))&0xffffffff\n",
        "        t^=(t+(((t^(t>>7))*(t|61))&0xffffffff))&0xffffffff\n",
        "        return ((t^(t>>14))&0xffffffff)/4294967296\n",
        "    return rand\n",
        "\n",
        "def autoencoder(X,train,dims=2,seed=7,epochs=250,rate=.01):\n",
        "    \"\"\"Independent NumPy implementation of the browser's 48-12-d-12-48 tanh AE.\n",
        "    Adam, full batch, 2/(N*P) MSE derivative, linear output, train only.\n",
        "    \"\"\"\n",
        "    rand=mulberry(seed); widths=[X.shape[1],12,dims,12,X.shape[1]];net=[]\n",
        "    for i,o in zip(widths[:-1],widths[1:]):\n",
        "        limit=np.sqrt(6/(i+o));net.append([np.array([(rand()*2-1)*limit for _ in range(i*o)]).reshape(i,o),np.zeros(o)])\n",
        "    m=[[np.zeros_like(p) for p in layer] for layer in net];v=[[np.zeros_like(p) for p in layer] for layer in net]\n",
        "    def forward(x):\n",
        "        acts=[x]\n",
        "        for l,(w,b) in enumerate(net):\n",
        "            z=acts[-1]@w+b;acts.append(np.tanh(z) if l<3 else z)\n",
        "        return acts\n",
        "    xt=X[train];hist=[{'epoch':0,'loss':float(np.mean((forward(xt)[-1]-xt)**2))}]\n",
        "    for epoch in range(1,epochs+1):\n",
        "        acts=forward(xt); delta=2*(acts[-1]-xt)/xt.size; grads=[None]*4\n",
        "        for l in reversed(range(4)):\n",
        "            grads[l]=[acts[l].T@delta,delta.sum(0)]\n",
        "            if l:delta=(delta@net[l][0].T)*(1-acts[l]**2)\n",
        "        for l in range(4):\n",
        "            for k in range(2):\n",
        "                m[l][k]=.9*m[l][k]+.1*grads[l][k];v[l][k]=.999*v[l][k]+.001*grads[l][k]**2\n",
        "                net[l][k]-=rate*(m[l][k]/(1-.9**epoch))/(np.sqrt(v[l][k]/(1-.999**epoch))+1e-8)\n",
        "        if epoch%10==0 or epoch==epochs:hist.append({'epoch':epoch,'loss':float(np.mean((forward(xt)[-1]-xt)**2))})\n",
        "    a=forward(X);return {'z':a[2], 'reconstruction':a[-1], 'history':hist}\n",
        "\n",
        "def native(x):\n",
        "    if isinstance(x,np.ndarray):return x.tolist()\n",
        "    if isinstance(x,np.generic):return x.item()\n",
        "    raise TypeError(type(x))\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 1. Fix the experimental split before fitting preprocessing"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "execution_count": null,
      "outputs": [],
      "source": [
        "d = dataset()\n",
        "X, mean, scale = prepare(d['X'], d['train'], 'standard')\n",
        "print('Original spectra:', d['X'].shape, 'train:', len(d['train']), 'held out:', len(d['test']))\n",
        "assert set(d['train']).isdisjoint(d['test'])\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 2. Fit the three inductive methods\n",
        "\n",
        "PCA is checked using sklearn; regularised Fisher LDA uses SciPy\u2019s generalised eigensolver. The autoencoder is trained from scratch by the complete NumPy code above."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "execution_count": null,
      "outputs": [],
      "source": [
        "pca_model = PCA(n_components=2, svd_solver='full').fit(X[d['train']])\n",
        "z_pca = pca_model.transform(X)\n",
        "z_lda, lda_eigenvalues = fisher(X, d['labels'], d['train'], dims=2, ridge=.1)\n",
        "ae = autoencoder(X, d['train'], dims=2, seed=7, epochs=250, rate=.01)\n",
        "print('PCA retained training variance:', pca_model.explained_variance_ratio_.sum())\n",
        "print('AE training scaled MSE, initial/final:', ae['history'][0]['loss'], ae['history'][-1]['loss'])\n",
        "for name, reconstructed in [('PCA', pca_model.inverse_transform(z_pca)), ('AE', ae['reconstruction'])]:\n",
        "    raw = reconstructed * scale + mean\n",
        "    print(name, 'train raw MSE', np.mean((d['X'][d['train']] - raw[d['train']])**2),\n",
        "          'held-out raw MSE', np.mean((d['X'][d['test']] - raw[d['test']])**2))\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 3. Fit exploratory full-data maps\n",
        "\n",
        "Every point participates in these fits. They are not out-of-sample projections. UMAP may take longer on its first call due to CPU compilation. Change the seed and parameters before comparing local neighbourhood stability; do not read island gaps as physical distances."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "execution_count": null,
      "outputs": [],
      "source": [
        "import umap\n",
        "z_tsne = TSNE(n_components=2, perplexity=20, random_state=7, init='random',\n",
        "              learning_rate='auto', max_iter=750, method='exact').fit_transform(X)\n",
        "z_umap = umap.UMAP(n_neighbors=15, min_dist=.1, random_state=7, n_jobs=1,\n",
        "                   n_components=2, n_epochs=300, metric='euclidean').fit_transform(X)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 4. Plot actual coordinates and reconstruction\n",
        "\n",
        "Each scatter independently scales its axes. Colours annotate toy classes; the methods were not all trained on the same information, so this is not a predictive ranking."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "execution_count": null,
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\n",
        "colours = ['#146a62', '#805a00', '#5847ad', '#a23956']\n",
        "fig, axes = plt.subplots(1, 5, figsize=(18, 3.5))\n",
        "for ax, name, z in zip(axes, ['PCA (train fit)', 'LDA (train labels)', 't-SNE (all rows)', 'UMAP (all rows)', 'AE (train fit)'],\n",
        "                     [z_pca, z_lda, z_tsne, z_umap, ae['z']]):\n",
        "    for c, colour in enumerate(colours):\n",
        "        ids = np.flatnonzero(d['labels'] == c)\n",
        "        ax.scatter(z[ids,0], z[ids,1], color=colour, s=16, label='ABCD'[c])\n",
        "    ax.set_title(name); ax.set_xlabel('Coordinate 1'); ax.set_ylabel('Coordinate 2')\n",
        "fig.tight_layout(); plt.show()\n",
        "selected = int(d['test'][0]); reconstructed = ae['reconstruction'][selected]*scale+mean\n",
        "plt.plot(d['wavelengths'], d['X'][selected], label='Original')\n",
        "plt.plot(d['wavelengths'], reconstructed, '--', label='Autoencoder reconstruction')\n",
        "plt.xlabel('Wavelength (nm)'); plt.ylabel('Synthetic unitless value'); plt.legend(); plt.show()\n",
        "plt.plot([h['epoch'] for h in ae['history']], [h['loss'] for h in ae['history']])\n",
        "plt.xlabel('Epoch'); plt.ylabel('Scaled training MSE'); plt.show()\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 5. Check leakage resistance and parameter constraints"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "execution_count": null,
      "outputs": [],
      "source": [
        "altered = d['labels'].copy(); altered[d['test']] = (altered[d['test']] + 1) % 4\n",
        "z_altered, _ = fisher(X, altered, d['train'], dims=2, ridge=.1)\n",
        "assert np.allclose(z_lda, z_altered)\n",
        "try:\n",
        "    fisher(X, d['labels'], d['train'], dims=4)\n",
        "except ValueError:\n",
        "    print('Correctly rejected LDA dimension > C\u22121')\n",
        "else:\n",
        "    raise AssertionError('Invalid dimension was accepted')\n",
        "assert ae['history'][-1]['loss'] < ae['history'][0]['loss']\n",
        "print('Notebook checks passed')\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Reproduce the browser presets and run the numerical tests\n",
        "\n",
        "From the extracted source directory:\n",
        "\n",
        "    python python/reproduce.py --presets\n",
        "    node --test tests/core.test.mjs\n",
        "    python -m http.server 8000\n",
        "\n",
        "Then open http://localhost:8000 in your browser. A web server is needed for module workers and JSON loading; double-clicking the HTML file is insufficient. See README.md for exact scope, references and verification limits."
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "name": "python",
      "version": "3.12"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}