readout layer trained

verified with time-series forecasting task and plotted

6 files changed, 334 insertions, 0 deletions

diff --git a/src/README.md b/src/README.md
new file mode 100644
index 0000000..0d58f19
--- /dev/null
+++ b/src/README.md
@@ -0,0 +1,36 @@
+# Spires reservoir simulated on memristor crossbar using memTorch
+## setup
+### general
+python3 -m .venv venv
+
+pip install -r requirements.txt
+
+### memTorch
+
+git clone --recursive https://github.com/coreylammie/MemTorch
+
+cd memTorch
+
+python setup.py install
+
+### spires
+
+
+## overview
+### simulation
+Here we are using both memTorch and the Spiresrc libraries working together to simulate
+a spires reservoir on a memristor crossbar. memTorch is acting as the weights, simulating
+at a physics level including device to device differences, conductance drag, etc....
+Spires is acting as the actual reservoir neurons and determining the spikes.
+
+### benchmark
+
+## Goal
+The goal is to able to simulate a spires reservoir on a memristor crossbar using
+a custom memristor spice model(HIL), thus you can run the same simulation with 
+multiple models to compare how each performs using benchmark results.
+
+## task
+
+## metrics
+
diff --git a/src/custom_memristor.sub b/src/custom_memristor.sub
new file mode 100644
index 0000000..bd909c8
--- /dev/null
+++ b/src/custom_memristor.sub
@@ -0,0 +1,14 @@
+* This is completely ai generated to be honest
+* Custom Behavioral Memristor Model
+.subckt custom_memristor plus minus state_node
+* Internal state capacitor to store value 'w'
+Cstate state_node 0 1IC=0.1
+
+* Example behavioral state derivative derivative equation (dw/dt)
+* Gbar represents the rate of change based on voltage across plus and minus
+Gstate 0 state_node value={V(plus,minus) * 1e3}
+
+* Behavioral current equation (I = V / R(w))
+* Resistance varies inversely with the voltage at 'state_node'
+Gmem plus minus value={V(plus,minus) / (100 + (10000 * (1 - V(state_node))))}
+.ends custom_memristor
diff --git a/src/requirements.txt b/src/requirements.txt
new file mode 100644
index 0000000..cf8aac2
--- /dev/null
+++ b/src/requirements.txt
@@ -0,0 +1,10 @@
+numpy
+torch
+pandas
+torchvision
+ipython
+lmfit
+scipy
+seaborn
+scikit-learn
+pyspice
diff --git a/src/spice_memristor.py b/src/spice_memristor.py
new file mode 100644
index 0000000..87dcdc3
--- /dev/null
+++ b/src/spice_memristor.py
@@ -0,0 +1,7 @@
+class SpiceMemristor:
+    def __init__(self):
+        self.build_circuit()
+        
+        self.build_simulator()
+
+        self.time = 0
diff --git a/src/spires_interface.py b/src/spires_interface.py
new file mode 100644
index 0000000..3d307ad
--- /dev/null
+++ b/src/spires_interface.py
@@ -0,0 +1,116 @@
+import ctypes
+import numpy
+import torch
+import random
+
+LIF_DISCRETE = 0
+
+#load the spires library
+spires_lib = ctypes.CDLL("../spires/lib/libspires.so")
+
+# C signautes for creating reservoir
+spires_lib.create_reservoir.argtypes = [
+        ctypes.c_size_t,
+        ctypes.c_size_t,
+        ctypes.c_size_t,
+        ctypes.c_double,
+        ctypes.c_double,
+        ctypes.c_double,
+        ctypes.c_double,
+        ctypes.c_double,
+        ctypes.c_int,
+        ctypes.c_int,
+        ctypes.POINTER(ctypes.c_double)
+]
+spires_lib.create_reservoir.restype = ctypes.c_void_p
+
+# C signatures for reservoir destruction
+spires_lib.free_reservoir.argtypes = [ctypes.c_void_p]
+spires_lib.free_reservoir.restype = None
+
+# C signatures for step reservoir function
+spires_lib.update_neuron.argtypes = [
+        ctypes.c_void_p,                #Pointer to spires reservoir struct
+        ctypes.c_int,
+        ctypes.c_double,
+        ctypes.c_double
+]
+spires_lib.update_neuron.restype = None
+
+# C signatures for get neuron spike function
+spires_lib.read_reservoir_spikes.argtypes = [
+        ctypes.c_void_p,
+        ctypes.POINTER(ctypes.c_float)
+]
+spires_lib.read_reservoir_spikes.restype = None
+
+def init_spires_reservoir(reservoir_size):
+    print("Creating spires reservoir")
+
+    neuron_parameters = (ctypes.c_double * 4)(0.0, 1.0, 0.2, 0.5)
+    c_neuron_parameters = ctypes.cast(neuron_parameters, ctypes.POINTER(ctypes.c_double))
+
+    reservoir_ptr = spires_lib.create_reservoir(
+            ctypes.c_size_t(int(reservoir_size)), # num neurons
+            ctypes.c_size_t(int(reservoir_size)), # num neurons
+            ctypes.c_size_t(2),                  # num_outputs
+            ctypes.c_double(0),                # spectral radius
+            ctypes.c_double(0.8),                # ei_ratio
+            ctypes.c_double(1.0),                # input_strength
+            ctypes.c_double(0.0),                # connectivity
+            ctypes.c_double(1.0),                # dt
+            ctypes.c_int(1),                  # connectivity type ( 1 = sparse)
+            ctypes.c_int(LIF_DISCRETE),       # neuron type
+            c_neuron_parameters # neuron params
+    )
+    
+    #allocate empty void pointer, (reservoir will go here)
+    if not reservoir_ptr: 
+        raise RuntimeError(f"Spires faile to initialize the reservoir: {status_code}")
+    else:
+        print("Spires reservoir initialized")
+
+    # reservoir_ptr._keep_alive = neuron_parameters
+
+    return reservoir_ptr
+
+def free_spires_reservoir(reservoir_ptr):
+    print("Freeing the spires reservoir")
+    spires_lib.free_reservoir(reservoir_ptr)
+    return 0
+    
+
+#change currents from tensor to a flat C pointer array for spires to read
+def send_currents_to_spires(reservoir_ptr, currents_tensor):
+    # print("Sending currents to spires")
+    # ----- extract from pytorhc graph -----
+    # .detach() removes it from auto gradient tracking
+    # .cpu() make sure data is in RAM, not VRAM
+    # .numpy() maps it to numpy array
+    # ----- makes suren layout matches 64 bit double C-array -----
+    # .astype(npfloat64) forces standard CC float precistion
+    # .flatten() makes sure the memory is a 1D block
+    numpy_array = currents_tensor.detach().cpu().numpy().astype(numpy.float64).flatten()
+
+    some_ptr = ctypes.cast(reservoir_ptr, ctypes.POINTER(ctypes.POINTER(ctypes.c_void_p)))
+
+    neurons_array = some_ptr[0]
+
+    # print("updating the neurons")
+    for i in range(len(numpy_array)):
+        neuron_ptr = neurons_array[i]
+        input_current = numpy_array[i]
+
+        spires_lib.update_neuron(neuron_ptr, LIF_DISCRETE, input_current, 1.0)
+
+    return numpy_array
+
+def read_spikes_from_spires(reservoir_ptr, size=0):
+    # print("Recieved spikes from spires")
+    #didnt do any safety checking womp womp
+    returned_spikes = numpy.zeros(size, dtype=numpy.float64)
+    c_spike_ptr = returned_spikes.ctypes.data_as(ctypes.POINTER(ctypes.c_float))
+
+    spires_lib.read_reservoir_spikes(reservoir_ptr, c_spike_ptr)
+
+    return returned_spikes
diff --git a/src/torch_reservoir.py b/src/torch_reservoir.py
new file mode 100644
index 0000000..6e7ceb7
--- /dev/null
+++ b/src/torch_reservoir.py
@@ -0,0 +1,151 @@
+import numpy as np
+import torch 
+import memtorch
+from spires_interface import (
+        free_spires_reservoir,
+        init_spires_reservoir,
+        send_currents_to_spires,
+        read_spikes_from_spires,
+)
+import matplotlib.pyplot as plt
+from sklearn.linear_model import Ridge
+import time
+
+# ----- PARAMETERS -----
+NUM_INPUTS = 1
+NUM_OUTPUTS = 1
+NUM_NEURONS = 2000
+time_steps = 2000
+
+input_layer = torch.nn.Linear(NUM_INPUTS, NUM_NEURONS, bias=False)
+#this is a recurrent layer sort of???
+reservoir_layer = torch.nn.Linear(NUM_NEURONS, NUM_NEURONS, bias=False)
+readout_layer = torch.nn.Linear(NUM_NEURONS, NUM_OUTPUTS, bias=False)
+
+spires_reservoir = init_spires_reservoir(NUM_NEURONS)
+
+#make the reservoir sparse and random
+with torch.no_grad():
+    reservoir_layer.weight.data.normal_(0.0, 0.5) #random weight
+
+    #mask so 10% of connections exist
+    mask = (torch.rand(NUM_NEURONS, NUM_NEURONS) < 0.10).float()
+    reservoir_layer.weight.data *= mask
+
+# patch layers into memristor crossbars with memTorch
+memristor_model = memtorch.bh.memristor.VTEAM
+memristor_model_params = {
+        'time_series_resolution': 1e-10,
+        'r_on': 50,
+        'r_off': 1000,
+}
+
+mem_input_layer = memtorch.mn.Module.patch_model(
+        input_layer, 
+        memristor_model,
+        memristor_model_params,
+)
+
+mem_reservoir_layer = memtorch.mn.Module.patch_model(
+        reservoir_layer,
+        memristor_model,
+        memristor_model_params,
+)
+
+start_time = time.perf_counter()
+
+# ----- run spires & collect spikes-----
+previous_spikes = torch.zeros(1, NUM_NEURONS)
+spike_history = []
+
+input_signal = []
+t = np.linspace(0, 20, time_steps)
+
+input_signal = np.sin(-2.5 * t) + np.sin(5 * t)
+target_signal = np.roll(input_signal, -5) #predicting 5 time steps in the future
+
+for step in range(time_steps):
+
+    input_tensor = torch.tensor([input_signal[step]]).float().unsqueeze(0)
+
+    currents_in = mem_input_layer(input_tensor)
+    currents_recv = mem_reservoir_layer(previous_spikes)
+
+    noise_multiplier = 0.08
+    current_scaler = 0.05
+
+    total_currents = (currents_in + currents_recv)
+    total_currents = total_currents - total_currents.mean() # remove DC bias
+    total_currents = total_currents + (torch.randn_like(total_currents) * noise_multiplier )
+    total_currents = total_currents * current_scaler #scale to fit threshhold
+
+    send_currents_to_spires(spires_reservoir, total_currents)
+    current_spikes = read_spikes_from_spires(spires_reservoir, NUM_NEURONS)
+
+    spike_history.append(current_spikes.copy())
+
+    # convert spikes back to tensor
+    previous_spikes = torch.from_numpy(current_spikes).float().unsqueeze(0)
+    
+#low pass filter
+spike_matrix = np.array(spike_history)
+
+decay_rate = 0.85
+filtered_spikes = np.zeros_like(spike_matrix, dtype=float)
+
+current_trace = np.zeros(NUM_NEURONS)
+for i in range(len(spike_matrix)):
+    current_trace = current_trace * decay_rate + spike_matrix[i]
+    filtered_spikes[i] = current_trace
+
+#training readout layer
+print("training the readout layer")
+x_train = filtered_spikes[100:800]
+y_train = target_signal[100:800]
+
+x_test = filtered_spikes[800:]
+y_test = target_signal[800:]
+
+ridge = Ridge(alpha=5)
+ridge.fit(x_train, y_train)
+
+#predict and plot results
+print("Generating predictions...")
+predictions = ridge.predict(x_test)
+
+end_time = time.perf_counter()
+execution_time = end_time - start_time
+print(f"execution time: {execution_time}")
+
+test_steps = np.arange(800, time_steps)
+
+plt.figure(figsize=(10, 5))
+plt.plot(test_steps, y_test, label="True Future Wave", color="black", linestyle="dashed")
+plt.plot(test_steps, predictions, label="Reservoir Prediction", color="blue", alpha=0.8)
+plt.title("spires memristor Time-Series Forecasting")
+plt.xlabel("time steps")
+plt.ylabel("Amplitude")
+plt.legend()
+plt.tight_layout()
+plt.show()
+
+free_spires_reservoir(spires_reservoir)
+
+# # ---------- Plotting ----------
+# #This block is all AI generated to be transparent
+# #plotting to verify neurons are firing randomly
+#
+# print("Simulation complete, plotting results") #why? cuz its fun and almost 5
+# spike_matrix = np.array(spike_history)
+# # Plotting
+# plt.figure(figsize=(12, 6))
+# # Transpose (.T) so Time is the X-axis and Neurons are the Y-axis
+# plt.imshow(spike_matrix.T, aspect='auto', cmap='binary', interpolation='nearest')
+#
+# plt.title("Reservoir Spiking Activity (Raster Plot)")
+# plt.xlabel("Time Step")
+# plt.ylabel("Neuron ID (0 to 799)")
+# plt.colorbar(label="Spike (0 or 1)")
+# plt.tight_layout()
+# plt.show()
+#