Hot Swapping models

12 files changed, 247 insertions, 173 deletions

diff --git a/bin/benchmark b/bin/benchmark
index b4de420..6cbd4fb 100755
--- a/bin/benchmark
+++ b/bin/benchmark
diff --git a/build/application.o b/build/application.o
new file mode 100644
index 0000000..575b272
--- /dev/null
+++ b/build/application.o
diff --git a/build/benchmark.o b/build/benchmark.o
new file mode 100644
index 0000000..1dfd78b
--- /dev/null
+++ b/build/benchmark.o
diff --git a/build/main.o b/build/main.o
index e2959f0..a4df765 100644
--- a/build/main.o
+++ b/build/main.o
diff --git a/build/read_crossbar.o b/build/read_crossbar.o
index 8383e08..eb1bb7f 100644
--- a/build/read_crossbar.o
+++ b/build/read_crossbar.o
diff --git a/build/spires_interface.o b/build/spires_interface.o
index 758c963..5f10f0e 100644
--- a/build/spires_interface.o
+++ b/build/spires_interface.o
diff --git a/makefile b/makefile
index c62c247..21c8ca4 100644
--- a/makefile
+++ b/makefile
@@ -8,10 +8,11 @@ BIN_DIR = bin
 TARGET = $(BIN_DIR)/benchmark
 
 SOURCES = \
-	$(SRC_DIR)/main.c \
+	$(SRC_DIR)/application.c \
 	$(SRC_DIR)/spires_interface.c \
 	$(SRC_DIR)/crossbar_generator.c \
-	$(SRC_DIR)/read_crossbar.c
+	$(SRC_DIR)/read_crossbar.c \
+	$(SRC_DIR)/benchmark.c 
 
 OBJECTS = $(SOURCES:$(SRC_DIR)/%.c=$(BUILD_DIR)/%.o)
 
diff --git a/src/README.md b/src/README.md
index b26f45d..cff34af 100644
--- a/src/README.md
+++ b/src/README.md
@@ -1,29 +1,11 @@
 OH MY GOD SO MUCH HAS CHANGED NEED TO UPDATE
 
-# 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
-
+# memristor crossbar readout layer for reservoir computer
+##setup 
+### download 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
 
diff --git a/src/application.c b/src/application.c
new file mode 100644
index 0000000..bda6aa7
--- /dev/null
+++ b/src/application.c
@@ -0,0 +1,163 @@
+#include "benchmark.h"
+
+#include <dirent.h>
+#include <math.h>
+#include <spires.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+// spires reservoir parameters
+#define NUM_NEURONS 600
+#define NUM_INPUTS 1
+#define NUM_OUTPUTS 1
+#define SPECTRAL_RADIUS 0.95
+#define EI_RATIO 0.8
+#define INPUT_STRENGTH 0.1
+#define CONNECTIVITY 0.1
+#define DT 1.0
+
+// ridge regression training
+#define PI 3.14159265358979323846
+#define LAMBDA 1.0e-4
+
+#define NUM_TRAINING_STEPS 500
+#define NUM_CROSSBAR_COLUMNS (NUM_OUTPUTS * 2)
+// #define NUM_STEPS 2000
+
+int main(void)
+{
+	// get number of files in model directory
+	// This can be replaced later in for loop with number of models to test
+	// size_t model_count = 0;
+	// DIR *dirp;
+	// struct dirent *entry;
+	//
+	// dirp = opendir("models");
+	// while ((entry = readdir(dirp)) != NULL) {
+	// 	if (entry->d_type == DT_REG) {
+	// 		model_count++;
+	// 	}
+	// }
+	// closedir(dirp);
+
+	/* ---------- LIST ALL MODELS HERE ----------*/
+	MemModel models[] = {
+	    {
+		.model_path = "models/hp_memristor.cir",
+		.subcircuit_name = "memristor",
+	    },
+	    {
+		.model_path = "models/fixed_resistor.cir",
+		.subcircuit_name = "fixed_resistor",
+	    },
+	};
+
+	size_t model_count = sizeof(models) / sizeof(models[0]);
+
+	// MemModel hp_memristor = {.model_path = "models/hp_memristor.cir",
+	// 			 .subcircuit_name = "memristor"};
+	//
+	// MemModel normal_resistor = {.model_path =
+	// "models/fixed_resistor.cir", .subcircuit_name = "fixed_resistor"};
+	//
+	// MemModel *models = calloc(model_count, sizeof(MemModel));
+	// if (models == NULL) {
+	// 	fprintf(stderr, "Failed to allocate for model list");
+	// 	return -1;
+	// }
+
+	/* ---------- SPIRES SET UP ----------*/
+	// discrete LIF parameters for spires
+	double lif_config[] = {
+	    0.0, // V_off
+	    1.0, // V_th
+	    0.2, // leak rate
+	    0.5, // bias
+	};
+
+	const spires_reservoir_config config = {
+	    .num_neurons = NUM_NEURONS,
+	    .num_inputs = NUM_INPUTS,
+	    .num_outputs = NUM_OUTPUTS,
+	    .spectral_radius = SPECTRAL_RADIUS,
+	    .ei_ratio = EI_RATIO,
+	    .input_strength = INPUT_STRENGTH,
+	    .connectivity = CONNECTIVITY,
+	    .dt = DT,
+	    .connectivity_type = SPIRES_CONN_RANDOM,
+	    .neuron_type = SPIRES_NEURON_LIF_DISCRETE,
+	    .neuron_params = lif_config};
+
+	spires_reservoir *reservoir = NULL;
+	if (spires_reservoir_create(&config, &reservoir) != 0) {
+		fprintf(stderr, "Failed to create reservoir");
+		return -1;
+	}
+
+	/* ---------- Train spires readout layer ----------*/
+	/* ---------- Training inputs ----------*/
+	double training_inputs[NUM_TRAINING_STEPS * NUM_INPUTS];
+	for (size_t timestep = 0; timestep < NUM_TRAINING_STEPS; timestep++) {
+		for (size_t input = 0; input < NUM_INPUTS; input++) {
+			double signal =
+			    0.7 * sin(2.0 * PI * (double)timestep / 50.0) +
+			    0.3 * sin(2.0 * PI * (double)timestep / 17.0);
+			training_inputs[timestep * NUM_INPUTS + input] = signal;
+		}
+	}
+
+	/* ---------- Target outputs  ----------*/
+	double target_outputs[NUM_TRAINING_STEPS * NUM_OUTPUTS];
+	for (size_t timestep = 0; timestep < NUM_TRAINING_STEPS; timestep++) {
+		size_t next_timestep = (timestep + 1) % NUM_TRAINING_STEPS;
+
+		double target =
+		    0.7 * sin(2.0 * PI * (double)next_timestep / 50.0) +
+		    0.3 * sin(2.0 * PI * (double)next_timestep / 17.0);
+		for (size_t output = 0; output < NUM_OUTPUTS; output++) {
+			target_outputs[timestep * NUM_OUTPUTS + output] =
+			    target;
+		}
+	}
+
+	/* ---------- Ridge Regression ----------*/
+	if (train_reservoir(reservoir, training_inputs, target_outputs,
+			    NUM_TRAINING_STEPS, LAMBDA) < 0) {
+		fprintf(stderr, "Failed to train the spires reservoir");
+		spires_reservoir_destroy(reservoir);
+		return -1;
+	}
+
+	/* ---------- Collect reservoir states ----------*/
+	Reservoir_State_Matrix state_matrix = {0};
+	if (collect_reservoir_states(reservoir, training_inputs,
+				     NUM_TRAINING_STEPS, &state_matrix) != 0) {
+		fprintf(stderr, "Failed to collect reservoir states");
+		spires_reservoir_destroy(reservoir);
+		return -1;
+	}
+	printf("collected state matrix size: %zu x %zu\n",
+	       state_matrix.num_samples, state_matrix.num_features);
+
+	/* ---------- Generate raster plot ----------*/
+	if (plot_raster(&state_matrix, NUM_NEURONS, 0.5) != 0) {
+		fprintf(stderr, "Failed to plot raster");
+	}
+
+	/* ---------- Run Benchmark on each model ----------*/
+	for (size_t i = 0; i < model_count; i++) {
+		printf("Running benchmark on %s\n", models[i].model_path);
+		if (run_benchmark(&config, reservoir, &state_matrix,
+				  target_outputs, models[i].model_path,
+				  models[i].subcircuit_name) < 0) {
+			fprintf(stderr, "Failed to run benchmark");
+			spires_reservoir_destroy(reservoir);
+			return -1;
+		}
+	}
+
+	free_reservoir_state_matrix(&state_matrix);
+	spires_reservoir_destroy(reservoir);
+
+	return 0;
+}
diff --git a/src/main.c b/src/benchmark.c
index ab96823..4ced485 100644
--- a/src/main.c
+++ b/src/benchmark.c
@@ -8,160 +8,71 @@
 #include <stdio.h>
 #include <stdlib.h>
 
-#define NUM_NEURONS 400
-#define NUM_INPUTS 4
-#define NUM_OUTPUTS 2
-#define NUM_CROSSBAR_COLUMNS (NUM_OUTPUTS * 2)
 #define NUM_TRAINING_STEPS 500
-#define NUM_STEPS 2000
 
-#define SPIKE_THRESHOLD 0.9
-#define SPIKE_AMPLITUDE 0.1
+#define SPIKE_THRESHOLD 0.1
+#define SPIKE_AMPLITUDE 1
 
-#define PI 3.14159265358979323846
+int plot_reservoir_predictions(const double *expected, const double *predicted,
+			       size_t num_samples, size_t num_outputs,
+			       size_t output_to_plot);
 
-static int plot_raster(const Reservoir_State_Matrix *matrix,
-		       size_t neurons_to_plot, double spike_threshold);
-
-static int plot_reservoir_predictions(const double *expected,
-				      const double *predicted,
-				      size_t num_samples, size_t num_outputs,
-				      size_t output_to_plot);
-
-int main(void)
+int run_benchmark(const spires_reservoir_config *config,
+		  spires_reservoir *reservoir,
+		  Reservoir_State_Matrix *state_matrix,
+		  const double *target_outputs, const char *model_path,
+		  const char *subcircuit_name)
 {
-	// discrete LIF parameters for spires
-	double lif_config[] = {
-	    0.0, // V_off
-	    1.0, // V_th
-	    0.2, // leak rate
-	    0.5, // bias
-	};
-
-	const spires_reservoir_config config = {
-	    .num_neurons = NUM_NEURONS,
-	    .num_inputs = NUM_INPUTS,
-	    .num_outputs = NUM_OUTPUTS,
-	    .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_config};
-
-	spires_reservoir *reservoir = NULL;
-
-	spires_status status = spires_reservoir_create(&config, &reservoir);
-
-	if (status != SPIRES_OK) {
-		fprintf(stderr, "Failed to create reservoir");
-		return -1;
-	}
-
-	// create training inputs
-	double training_inputs[NUM_TRAINING_STEPS * NUM_INPUTS];
-	for (size_t timestep = 0; timestep < NUM_TRAINING_STEPS; timestep++) {
-		for (size_t input = 0; input < NUM_INPUTS; input++) {
-			double signal =
-			    0.7 * sin(2.0 * PI * (double)timestep / 50.0) +
-			    0.3 * sin(2.0 * PI * (double)timestep / 17.0);
-			training_inputs[timestep * NUM_INPUTS + input] = signal;
-		}
-	}
-
-	// create target outputs
-	double target_outputs[NUM_TRAINING_STEPS * NUM_OUTPUTS];
-	for (size_t timestep = 0; timestep < NUM_TRAINING_STEPS; timestep++) {
-		size_t next_timestep = (timestep + 1) % NUM_TRAINING_STEPS;
-
-		double target =
-		    0.7 * sin(2.0 * PI * (double)next_timestep / 50.0) +
-		    0.3 * sin(2.0 * PI * (double)next_timestep / 17.0);
-		for (size_t output = 0; output < NUM_OUTPUTS; output++) {
-			target_outputs[timestep * NUM_OUTPUTS + output] =
-			    target;
-		}
-	}
-
-	Reservoir_State_Matrix state_matrix = {0};
-	if (collect_reservoir_states(reservoir, training_inputs,
-				     NUM_TRAINING_STEPS, &state_matrix) != 0) {
-		fprintf(stderr, "Failed to collect reservoir states");
-		spires_reservoir_destroy(reservoir);
-		return -1;
-	}
-	printf("collected state matrix: %zu x %zu\n", state_matrix.num_samples,
-	       state_matrix.num_features);
-
-	// training the readout layer
-	const double lambda = 1.0e-4;
-	int training_status =
-	    train_reservoir(reservoir, training_inputs, target_outputs,
-			    NUM_TRAINING_STEPS, lambda);
-	if (training_status < 0) {
-		fprintf(stderr, "Failed to train the reservoir");
-		free_reservoir_state_matrix(&state_matrix);
-		spires_reservoir_destroy(reservoir);
-		return -1;
-	}
-
-	// generate raster plot for verification
-	if (plot_raster(&state_matrix, NUM_NEURONS, SPIKE_THRESHOLD) != 0) {
-		fprintf(stderr, "Failed to plot raster\n");
-	}
-
 	// copy readout weights and convert to conductances
 	double *initial_resistances = NULL;
 	conductance_mapping mapping;
 
 	if (convert_weights_to_resistances(
-		reservoir, NUM_NEURONS, NUM_OUTPUTS, 1000.0, 100000.0,
-		&initial_resistances, &mapping) != 0) {
+		reservoir, config->num_neurons, config->num_outputs, 1000.0,
+		100000.0, &initial_resistances, &mapping) != 0) {
 		return -1;
 	}
 
 	double *row_voltages =
-	    malloc(state_matrix.num_features * state_matrix.num_samples *
+	    malloc(state_matrix->num_features * state_matrix->num_samples *
 		   sizeof(double));
 
 	if (row_voltages == NULL) {
 		fprintf(stderr, "Failed to allocate spikes voltages");
 		free(initial_resistances);
-		free_reservoir_state_matrix(&state_matrix);
+		free_reservoir_state_matrix(state_matrix);
 		spires_reservoir_destroy(reservoir);
 		return -1;
 	}
 
-	for (size_t sample = 0; sample < state_matrix.num_samples; sample++) {
-		for (size_t neuron = 0; neuron < state_matrix.num_features;
+	for (size_t sample = 0; sample < state_matrix->num_samples; sample++) {
+		for (size_t neuron = 0; neuron < state_matrix->num_features;
 		     neuron++) {
 			size_t index =
-			    sample * state_matrix.num_features + neuron;
+			    sample * state_matrix->num_features + neuron;
 
 			row_voltages[index] =
-			    SPIKE_AMPLITUDE * state_matrix.states[index];
+			    SPIKE_AMPLITUDE * state_matrix->states[index];
 		}
 	}
 
 	const Crossbar_Config crossbar_config = {
-	    .rows = state_matrix.num_features,
-	    .columns = NUM_CROSSBAR_COLUMNS,
+	    .rows = state_matrix->num_features,
+	    .columns = config->num_outputs * 2,
 	    .input_series = row_voltages,
-	    .num_samples = state_matrix.num_samples,
+	    .num_samples = state_matrix->num_samples,
 	    .initial_resistance = initial_resistances,
-	    .model_path = "models/hp_memristor.cir",
-	    .subcircuit_name = "memristor",
 	    .load_resistance = 50.0,
+	    .model_path = model_path,
+	    .subcircuit_name = subcircuit_name,
 	    .time_step = 1e-6,
-	    .stop_time = state_matrix.num_samples * 1e-6,
-	    .print_state_nodes = 0};
+	    .stop_time = state_matrix->num_samples * 1e-6,
+	    .print_state_nodes = 0}; // state nodes is not acutally implemented
 
 	if (generate_crossbar("output/crossbar.cir", &crossbar_config) < 0) {
 		fprintf(stderr, "failed to create crossbar config");
 		free(initial_resistances);
-		free_reservoir_state_matrix(&state_matrix);
+		free_reservoir_state_matrix(state_matrix);
 		spires_reservoir_destroy(reservoir);
 		return -1;
 	}
@@ -171,7 +82,7 @@ int main(void)
 	if (run_ngspice("output/crossbar.cir") < 0) {
 		fprintf(stderr, "Failed to run_ngspice");
 		free(initial_resistances);
-		free_reservoir_state_matrix(&state_matrix);
+		free_reservoir_state_matrix(state_matrix);
 		spires_reservoir_destroy(reservoir);
 		return -1;
 	}
@@ -179,16 +90,16 @@ int main(void)
 
 	// crossbar parameters needed for reading
 	Crossbar_Output_Matrix crossbar_output = {
-	    .num_samples = NUM_TRAINING_STEPS,
-	    .num_outputs = NUM_CROSSBAR_COLUMNS,
+	    .num_samples = state_matrix->num_samples,
+	    .num_outputs = config->num_outputs * 2,
 	    .time = NULL,
 	    .voltages = NULL};
 
-	if (read_crossbar("output/crossbar_output.dat", NUM_CROSSBAR_COLUMNS,
+	if (read_crossbar("output/crossbar_output.dat", config->num_outputs * 2,
 			  &crossbar_output) < 0) {
 		fprintf(stderr, "Failed to read crossbar output file");
 		free(initial_resistances);
-		free_reservoir_state_matrix(&state_matrix);
+		free_reservoir_state_matrix(state_matrix);
 		spires_reservoir_destroy(reservoir);
 		free_crossbar_output_matrix(&crossbar_output);
 		return -1;
@@ -196,13 +107,13 @@ int main(void)
 	printf("read the crossbar outputs!!\n");
 
 	double *decoded_outputs =
-	    malloc(NUM_OUTPUTS * NUM_TRAINING_STEPS * sizeof(double));
+	    malloc(config->num_outputs * NUM_TRAINING_STEPS * sizeof(double));
 	if (decoded_outputs == NULL) {
 		fprintf(stderr,
 			"Failed to allocate memory for decoded outputs");
 	}
 	if (convert_output_to_software(
-		NUM_NEURONS, NUM_OUTPUTS, NUM_TRAINING_STEPS,
+		config->num_neurons, config->num_outputs, NUM_TRAINING_STEPS,
 		crossbar_output.voltages, initial_resistances,
 		crossbar_config.load_resistance, &mapping, row_voltages,
 		SPIKE_AMPLITUDE, decoded_outputs) < 0) {
@@ -219,22 +130,19 @@ int main(void)
 
 	// plotting expected vs. prediction
 	plot_reservoir_predictions(target_outputs, decoded_outputs,
-				   NUM_TRAINING_STEPS, NUM_OUTPUTS, 0);
+				   NUM_TRAINING_STEPS, config->num_outputs, 0);
 
 	printf("YAY IT WORKED!!! Cleaning up :)");
 	free(initial_resistances);
-	// free(spikes_voltages);
 	free(row_voltages);
 	free(decoded_outputs);
-	free_reservoir_state_matrix(&state_matrix);
-	spires_reservoir_destroy(reservoir);
 	free_crossbar_output_matrix(&crossbar_output);
 
 	return 0;
 }
 
-static int plot_raster(const Reservoir_State_Matrix *matrix,
-		       size_t neurons_to_plot, double spike_threshold)
+int plot_raster(const Reservoir_State_Matrix *matrix, size_t neurons_to_plot,
+		double spike_threshold)
 {
 	if (!matrix || !matrix->states || matrix->num_samples == 0) {
 		return -1;
@@ -321,10 +229,9 @@ static int plot_raster(const Reservoir_State_Matrix *matrix,
 	return 0;
 }
 
-static int plot_reservoir_predictions(const double *expected,
-				      const double *predicted,
-				      size_t num_samples, size_t num_outputs,
-				      size_t output_to_plot)
+int plot_reservoir_predictions(const double *expected, const double *predicted,
+			       size_t num_samples, size_t num_outputs,
+			       size_t output_to_plot)
 {
 	PLFLT *x;
 	PLFLT *y_expected;
diff --git a/src/benchmark.h b/src/benchmark.h
new file mode 100644
index 0000000..cd79d50
--- /dev/null
+++ b/src/benchmark.h
@@ -0,0 +1,27 @@
+#ifndef BENCHMARK_H
+#define BENCHMARK_H
+
+#include "crossbar_generator.h"
+#include "read_crossbar.h"
+#include "spires_interface.h"
+#include <stddef.h>
+
+typedef struct {
+	const char *model_path;
+	const char *subcircuit_name;
+} MemModel;
+
+int run_benchmark(const spires_reservoir_config *config,
+		  spires_reservoir *reservoir,
+		  Reservoir_State_Matrix *state_matrix,
+		  const double *target_outputs, const char *model_path,
+		  const char *subcircuit_name);
+
+int plot_raster(const Reservoir_State_Matrix *matrix, size_t neurons_to_plot,
+		double spike_threshold);
+
+int plot_reservoir_predictions(const double *expected, const double *predicted,
+			       size_t num_samples, size_t num_outputs,
+			       size_t output_to_plot);
+
+#endif
diff --git a/src/spires_interface.c b/src/spires_interface.c
index 9c37155..06e23bf 100644
--- a/src/spires_interface.c
+++ b/src/spires_interface.c
@@ -72,25 +72,17 @@ int convert_weights_to_resistances(const spires_reservoir *reservoir,
 				   double **resistances_out,
 				   conductance_mapping *mapping)
 {
-	double *readout;
-	double *resistances;
 	double max_abs_weight = 0.0;
-	size_t weight_count;
-	size_t weight_index;
-	size_t positive_index;
-	size_t negative_index;
-	size_t positive_column;
-	size_t negative_column;
-	double positive_conductance;
-	double negative_conductance;
 
 	*resistances_out = NULL;
-	weight_count = num_neurons * num_outputs;
+	size_t weight_count = num_neurons * num_outputs;
 	size_t num_physical_columns = num_outputs * 2;
 
-	readout = spires_copy_readout(reservoir);
+	double *readout = malloc(num_neurons * num_outputs * sizeof(double));
 
-	resistances =
+	spires_read_readout(reservoir, readout);
+
+	double *resistances =
 	    malloc(num_neurons * num_physical_columns * sizeof(double));
 	if (resistances == NULL) {
 		fprintf(stderr, "Failed to allocated crossbar resistances");
@@ -126,13 +118,15 @@ int convert_weights_to_resistances(const spires_reservoir *reservoir,
 	// differential pair mapping
 	for (size_t neuron = 0; neuron < num_neurons; neuron++) {
 		for (size_t output = 0; output < num_outputs; output++) {
-			weight_index = output * num_neurons + neuron;
-			positive_column = 2 * output;
-			negative_column = positive_column + 1;
+			size_t weight_index = output * num_neurons + neuron;
+			size_t positive_column = 2 * output;
+			size_t negative_column = positive_column + 1;
+			double positive_conductance;
+			double negative_conductance;
 
-			positive_index =
+			size_t positive_index =
 			    neuron * num_physical_columns + positive_column;
-			negative_index =
+			size_t negative_index =
 			    neuron * num_physical_columns + negative_column;
 
 			double weight = readout[weight_index];