Trained software weights with ridge regression

3 files changed, 80 insertions, 134 deletions

diff --git a/references/kaevin_memristor.py b/references/kaevin_memristor.py
deleted file mode 100644
index 890eb9b..0000000
--- a/references/kaevin_memristor.py
+++ /dev/null
@@ -1,134 +0,0 @@
-import numpy as np
-from scipy.integrate import solve_ivp
-
-# TEAM (Threshold Adaptive Memristor) model as a reusable class
-class TEAMMemristor:
-
-    def __init__(
-        self,
-        
-        k_on=-1,          # switching rate for on-state
-        alpha_off=5,      # exponent controlling nonlinearity when switching off
-        alpha_on=5,       # exponent controlling nonlinearity when switching on
-        i_off=0.5e-3,     # threshold current to trigger off state switching
-        i_on=-0.5e-3,     # threshold current to trigger on state switching
-        g_on=1/1e3,       # maximum conductance (1 kohm = 1 ms)
-        g_off=1/10e3,     # minimum conductance (10 kohm = 0.1 ms)
-        w_init=0.5,       # initial state variable (0=off, 1=on)
-        p=2               # window function exponent
-    ):
-        self.k_off = k_off
-        self.k_on = k_on
-        self.alpha_off = alpha_off
-        self.alpha_on = alpha_on
-        self.i_off = i_off
-        self.i_on = i_on
-        self.g_on = g_on
-        self.g_off = g_off
-        self.w_init = w_init
-        self.p = p
-
-    def set_state(self, w):
-        # Update initial condition for next simulation
-        self.w_init = np.clip(w, 0, 1)
-
-    def window(self, w, i):
-        # Nonlinear window function: reduces switching rate near the boundaries
-        w = np.clip(w, 0.0, 1.0)
-        if i >= 0:
-            return 1 - w**(2*self.p)  # switching off: slower near w=1
-        else:
-            return 1 - (1-w)**(2*self.p)  # switching on: slower near w=0
-
-    def conductance(self, w):
-        # Linear interpolation between off and on conductance based on state w
-        w = np.clip(w, 0.0, 1.0)
-        return self.G_off + w*(self.G_on - self.G_off)
-
-    def dw_dt(self, w, i):
-        # TEAM state dynamics: dw/dt depends on current magnitude and direction
-        w = np.clip(w, 0.0, 1.0)
-        
-        if i >= self.i_off:  # positive current above threshold = switch off
-            dw = (
-                self.k_off
-                * ((i/self.i_off)-1)**self.alpha_off
-                * self.window(w, i)
-            )
-        elif i <= self.i_on:  # negative current below threshold = switch on
-            dw = (
-                self.k_on
-                * (((-i)/abs(self.i_on))-1)**self.alpha_on
-                * self.window(w, i)
-            )
-        else:  # between thresholds = no switching
-            dw = 0.0
-        
-        # Enforce physical bounds: prevent state from leaving [0,1]
-        if w <= 0 and dw < 0:
-            dw = 0
-        if w >= 1 and dw > 0:
-            dw = 0
-        
-        return dw
-
-    def simulate(self,
-                 freq=1,          # excitation frequency (Hz)
-                 V_amp=1.5,       # sinusoid amplitude (V)
-                 cycles=3):       # number of periods to simulate
-        # Solve the memristor ODE for given frequency and voltage amplitude
-        
-        def voltage(t):  # sinusoidal excitation signal
-            return V_amp*np.sin(2*np.pi*freq*t)
-        
-        def ode(t, y):  # dy/dt: current through memristor
-            w = y[0]
-            v = voltage(t)
-            G = self.conductance(w)
-            i = G*v  # Ohm's law: i = G*v
-            return [self.dw_dt(w, i)]
-        
-        T = 1/freq  # period
-        t_end = cycles*T
-        t_eval = np.linspace(0, t_end, 10000)  # dense time grid for smooth curves
-        
-        # Solve with RK45, small max step
-        sol = solve_ivp(
-            ode,
-            [0, t_end],
-            [self.w_init],
-            t_eval=t_eval,
-            method='RK45',
-            max_step=T/1000,  # max step keeps resolution within one period
-            rtol=1e-8,
-            atol=1e-10
-        )
-        
-        raw_w = sol.y[0]
-        
-        # Check if numerical solver violated physical bounds
-        eps = 1e-6
-        if np.any(raw_w < -eps) or np.any(raw_w > 1+eps):
-            print(
-                "WARNING: solver left bounds "
-                f"min={raw_w.min():.12f}, "
-                f"max={raw_w.max():.12f}"
-            )
-        
-        w = np.clip(raw_w, 0, 1)  # enforce bounds just in case
-        t = sol.t
-        v = voltage(t)
-        G = self.conductance(w)
-        i = G*v  # current throughout simulation
-        
-        return t, w, v, i
-
-    def resistance(self, w):
-        # Compute resistance as reciprocal of conductance
-        return 1/self.conductance(w)
-    
-    def reset(self):
-        # Reset state to default initial condition
-        self.w_init = 0.5
-
-
diff --git a/references/memristor_spires_proof.png b/references/memristor_spires_proof.png
new file mode 100644
index 0000000..ea294fc
--- /dev/null
+++ b/references/memristor_spires_proof.png
diff --git a/references/spires_tutorial.c b/references/spires_tutorial.c
new file mode 100644
index 0000000..f4eba10
--- /dev/null
+++ b/references/spires_tutorial.c
@@ -0,0 +1,80 @@
+#include <../src/neurons/lif_discrete.h>
+#include <math.h>
+#include <spires.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#define N_TRAIN 500
+#define N_TEST 100
+#define PI 3.14159265358979323846
+
+int main(void) {
+  /* 1. Configure the reservoir */
+  double lif_cfg[] = {0.0, 1.0, 0.2, 0.5};
+
+  spires_reservoir_config cfg = {
+      .num_neurons = 400,
+      .num_inputs = 1,
+      .num_outputs = 1,
+      .spectral_radius = 0.95,
+      .ei_ratio = 0.8,
+      .input_strength = 0.1,
+      .connectivity = 0.1,
+      .dt = 1.0,
+      .connectivity_type = SPIRES_CONN_RANDOM,
+      .neuron_type = SPIRES_NEURON_LIF_DISCRETE,
+      .neuron_params = lif_cfg,
+  };
+
+  /* 2. Create the reservoir */
+  spires_reservoir *r = NULL;
+  spires_status s = spires_reservoir_create(&cfg, &r);
+  if (s != SPIRES_OK) {
+    fprintf(stderr, "Failed to create reservoir: %d\n", s);
+    return 1;
+  }
+
+  /* 3. Generate training data — predict sin(t+1) from sin(t) */
+  double input_train[N_TRAIN];
+  double target_train[N_TRAIN];
+  for (int i = 0; i < N_TRAIN; i++) {
+    input_train[i] = sin(2.0 * PI * i / 50.0);
+    target_train[i] = sin(2.0 * PI * (i + 1) / 50.0);
+  }
+
+  /* 4. Train with ridge regression */
+  s = spires_train_ridge(r, input_train, target_train, N_TRAIN, 1e-6);
+
+  if (s != SPIRES_OK) {
+    fprintf(stderr, "Training failed: %d\n", s);
+    spires_reservoir_destroy(r);
+    return 1;
+  }
+
+  /* 5. Generate test input */
+  double input_test[N_TEST];
+  for (int i = 0; i < N_TEST; i++) {
+    input_test[i] = sin(2.0 * PI * (N_TRAIN + i) / 50.0);
+  }
+
+  /* 6. Run inference */
+  double *predictions = spires_run(r, input_test, N_TEST);
+  if (!predictions) {
+    fprintf(stderr, "Inference failed\n");
+    spires_reservoir_destroy(r);
+    return 1;
+  }
+
+  /* 7. Print a few predictions vs. expected values */
+  printf("Step | Predicted | Expected\n");
+  printf("-----+-----------+---------\n");
+  for (int i = 0; i < 10; i++) {
+    double expected = sin(2.0 * PI * (N_TRAIN + i + 1) / 50.0);
+    printf("%4d | %+.5f | %+.5f\n", i, predictions[i], expected);
+  }
+
+  /* 8. Clean up */
+  free(predictions);
+  spires_reservoir_destroy(r);
+  return 0;
+}