SPICE crossbar readout layer

1 files changed, 245 insertions, 344 deletions

diff --git a/src/main.c b/src/main.c
index e8a6d92..e238d96 100644
--- a/src/main.c
+++ b/src/main.c
@@ -1,407 +1,308 @@
-#include "spires_interface.h"
 #include "crossbar_generator.h"
 #include "read_crossbar.h"
+#include "spires_interface.h"
 
 #include <math.h>
+#include <plplot/plplot.h>
 #include <spires.h>
 #include <stdio.h>
 #include <stdlib.h>
-#include <plplot/plplot.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.5
+#define SPIKE_THRESHOLD 0.9
 #define SPIKE_AMPLITUDE 0.1
 
 #define PI 3.14159265358979323846
 
-static int plot_raster(
-    const Reservoir_State_Matrix *matrix,
-    size_t neurons_to_plot,
-    double spike_threshold
-); 
+static int plot_raster(const Reservoir_State_Matrix *matrix,
+		       size_t neurons_to_plot, double spike_threshold);
 
+int main(void)
+{
+	// discrete LIF parameters for spires
+	double lif_config[] = {
+	    0.0, // V_off
+	    1.0, // V_th
+	    0.2, // leak rate
+	    0.5, // bias
+	};
 
-static int print_software_outputs(
-    spires_reservoir *reservoir,
-    const double *input_series,
-    size_t series_length,
-    size_t num_inputs,
-    size_t num_outputs,
-    size_t samples_to_print,
-    const char *label
-);
+	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};
 
-int main(void) {
-    //discrete LIF parameters for spires
-    double lif_config[] = {
-        0.0, //V_off
-        1.0, //V_th
-        0.2, //leak rate
-        0.5, //bias
-    };
+	spires_reservoir *reservoir = NULL;
 
-    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_status status = spires_reservoir_create(&config, &reservoir);
 
-    spires_reservoir *reservoir = NULL;
-    
-    spires_status status = spires_reservoir_create(
-            &config,
-            &reservoir
-    );
+	if (status != SPIRES_OK) {
+		fprintf(stderr, "Failed to create reservoir");
+		return -1;
+	}
 
-    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++) {
+			training_inputs[timestep * NUM_INPUTS + input] =
+			    sin(2.0 * PI * (double)timestep / 50.0);
+		}
+	}
 
-    //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++) {
-            training_inputs[timestep * NUM_INPUTS + input] = sin(2.0 * PI * (double)timestep / 50.0);
-        }
-    }
-    
-    //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 = sin(2.0 * PI * (double)next_timestep / 50.0);
-        for (size_t output = 0; output < NUM_OUTPUTS; output++) {
-            target_outputs[timestep * NUM_OUTPUTS + output] = target;
-        }
-    }
+	// 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 = sin(2.0 * PI * (double)next_timestep / 50.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);
+	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);
 
-    print_software_outputs(
-        reservoir,
-        training_inputs,
-        NUM_TRAINING_STEPS,
-        NUM_INPUTS,
-        NUM_OUTPUTS,
-        10,
-        "Software outputs before training:"
-    );
+	// 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;
+	}
 
-    //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");
+	}
 
-    print_software_outputs(
-        reservoir,
-        training_inputs,
-        NUM_TRAINING_STEPS,
-        NUM_INPUTS,
-        NUM_OUTPUTS,
-        10,
-        "Software outputs after training:"
-    );
+	// copy readout weights and convert to conductances
+	double *initial_resistances = NULL;
+	conductance_mapping mapping;
 
-    //generate raster plot for verification
-    if (plot_raster(&state_matrix, NUM_NEURONS, SPIKE_THRESHOLD) != 0) {
-        fprintf(stderr, "Failed to plot raster\n");
-    }
+	if (convert_weights_to_resistances(
+		reservoir, NUM_NEURONS, NUM_OUTPUTS, 1000.0, 100000.0,
+		&initial_resistances, &mapping) != 0) {
+		return -1;
+	}
 
-    //fill out initial resistances
-    double *initial_resistances = malloc(NUM_NEURONS * NUM_OUTPUTS * sizeof(*initial_resistances));
-    if (!initial_resistances) {
-        fprintf(stderr, "Failed to allocate memory for initial resistances");
-        return -1;
-    }
+	double *row_voltages =
+	    malloc(state_matrix.num_features * state_matrix.num_samples *
+		   sizeof(double));
 
-    //resistances are inversely proportional to the software weigts
-    for (size_t i = 0; i < NUM_NEURONS * NUM_OUTPUTS; i++) {
-        initial_resistances[i] = 80000;
-    }
+	if (row_voltages == NULL) {
+		fprintf(stderr, "Failed to allocate spikes voltages");
+		free(initial_resistances);
+		free_reservoir_state_matrix(&state_matrix);
+		spires_reservoir_destroy(reservoir);
+		return -1;
+	}
 
-    //convert continous states to spikes
-    double *spikes_voltages = malloc(state_matrix.num_features * 
-            state_matrix.num_samples * sizeof(*spikes_voltages));
+	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;
 
-    if (spikes_voltages == NULL) {
-        fprintf(stderr, "Failed to allocate spikes voltages");
-        free(initial_resistances);
-        free_reservoir_state_matrix(&state_matrix);
-        spires_reservoir_destroy(reservoir);
-        return -1;
-    }
+			row_voltages[index] =
+			    SPIKE_AMPLITUDE * state_matrix.states[index];
+		}
+	}
 
-    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; 
+	const Crossbar_Config crossbar_config = {
+	    .rows = state_matrix.num_features,
+	    .columns = NUM_CROSSBAR_COLUMNS,
+	    .input_series = row_voltages,
+	    .num_samples = state_matrix.num_samples,
+	    .initial_resistance = initial_resistances,
+	    .model_path = "hp_memristor.cir",
+	    .subcircuit_name = "memristor",
+	    .load_resistance = 50.0,
+	    .time_step = 1e-6,
+	    .stop_time = state_matrix.num_samples * 1e-6,
+	    .print_state_nodes = 0};
 
-           spikes_voltages[index] = 
-               state_matrix.states [index] > SPIKE_THRESHOLD ? SPIKE_AMPLITUDE : 0.0;
-        }
-    }
+	if (generate_crossbar("crossbar.cir", &crossbar_config) < 0) {
+		fprintf(stderr, "failed to create crossbar config");
+		free(initial_resistances);
+		free_reservoir_state_matrix(&state_matrix);
+		spires_reservoir_destroy(reservoir);
+		return -1;
+	}
+	printf("Generated crossbar!!");
 
-    const Crossbar_Config crossbar_config = {
-        .rows = state_matrix.num_features,
-        .columns = NUM_OUTPUTS,
-        .input_series = spikes_voltages,
-        .num_samples = state_matrix.num_samples,
-        .initial_resistance = initial_resistances,
-        .model_path = "hp_memristor.cir",
-        .subcircuit_name = "memristor",
-        .load_resistance = 50.0,
-        .time_step = 1e-6,
-        .stop_time = state_matrix.num_samples * 1e-6,
-        .print_state_nodes = 0
-    };
+	// call ngspice for crossbar
+	if (run_ngspice("crossbar.cir") < 0) {
+		fprintf(stderr, "Failed to run_ngspice");
+		free(initial_resistances);
+		free_reservoir_state_matrix(&state_matrix);
+		spires_reservoir_destroy(reservoir);
+		return -1;
+	}
+	printf("ran ngspice!!");
 
-    if (generate_crossbar("crossbar.cir", &crossbar_config) < 0) {
-        fprintf(stderr, "failed to create crossbar config");
-        free(initial_resistances);
-        free_reservoir_state_matrix(&state_matrix);
-        spires_reservoir_destroy(reservoir);
-        return -1;
-    }
-    printf("Generated crossbar!!");
+	// crossbar parameters needed for reading
+	Crossbar_Output_Matrix crossbar_output = {
+	    .num_samples = NUM_TRAINING_STEPS,
+	    .num_outputs = NUM_CROSSBAR_COLUMNS,
+	    .time = NULL,
+	    .voltages = NULL};
 
-    //call ngspice for crossbar
-    if (run_ngspice("crossbar.cir") < 0) {
-        fprintf(stderr, "Failed to run_ngspice");
-        free(initial_resistances);
-        free_reservoir_state_matrix(&state_matrix);
-        spires_reservoir_destroy(reservoir);
-        return -1;
-    }
-    printf("ran ngspice!!");
+	if (read_crossbar("crossbar_output.dat", NUM_CROSSBAR_COLUMNS,
+			  &crossbar_output) < 0) {
+		fprintf(stderr, "Failed to read crossbar output file");
+		free(initial_resistances);
+		free_reservoir_state_matrix(&state_matrix);
+		spires_reservoir_destroy(reservoir);
+		free_crossbar_output_matrix(&crossbar_output);
+		return -1;
+	}
+	printf("read the crossbar outputs!!\n");
 
-    //crossbar parameters needed for reading
-    Crossbar_Output_Matrix crossbar_output = {
-        .num_samples = NUM_TRAINING_STEPS, //this isnt right?
-        .num_outputs = NUM_OUTPUTS,
-        .time = NULL,
-        .voltages = NULL
-    };
-    
-    if (read_crossbar("crossbar_output.dat", NUM_OUTPUTS, &crossbar_output) < 0) {
-        fprintf(stderr, "Failed to read crossbar output file");
-        free(initial_resistances);
-        free_reservoir_state_matrix(&state_matrix);
-        spires_reservoir_destroy(reservoir);
-        free_crossbar_output_matrix(&crossbar_output);
-        return -1;
-    }
-    printf("read the crossbar outputs!!");
+	double *decoded_outputs =
+	    malloc(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,
+		crossbar_output.voltages, initial_resistances,
+		crossbar_config.load_resistance, &mapping, row_voltages,
+		SPIKE_AMPLITUDE, decoded_outputs) < 0) {
+		fprintf(stderr,
+			"Failed to convert crossbar outputs back to software");
+		return -1;
+	}
 
-    //printing for testing purposes 
-    // printf("Read data:\n"); 
-    // for (size_t sample = 0; sample < crossbar_output.num_samples; sample++) {
-    //     printf("%f", crossbar_output.time[sample]);
-    //     for (size_t output = 0; output < crossbar_output.num_outputs; output++) {
-    //         printf(" %f", 
-    //                 crossbar_output.voltages[sample * crossbar_output.num_outputs + output]);
-    //     }
-    //     printf("\n");
-    // }
+	// comparing prediction
+	for (int i = 0; i < NUM_TRAINING_STEPS; i++) {
+		printf("real: %g , predicted: %g\n", target_outputs[i],
+		       decoded_outputs[i]);
+	}
 
-    printf("YAY IT WORKED!!! Cleaning up :)");
-    free(initial_resistances);
-    free(spikes_voltages);
-    free_reservoir_state_matrix(&state_matrix);
-    spires_reservoir_destroy(reservoir);
-    free_crossbar_output_matrix(&crossbar_output);
+	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;
+	return 0;
 }
 
-static 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;
-    }
-
-    if (neurons_to_plot > matrix->num_features) {
-        neurons_to_plot = matrix->num_features;
-    }
-
-    //count spikes
-    size_t spike_count = 0;
-    for (size_t t = 0; t < matrix->num_samples; t++) {
-        for (size_t n = 0; n < neurons_to_plot; n++) {
-            double value = matrix->states[t * matrix->num_features + n];
-            if (value > spike_threshold) {
-                spike_count++;
-            }
-        }
-    }
-
-    if (spike_count == 0) {
-        fprintf(stderr, "No spikes found above threshold %.3f\n", spike_threshold);
-        return -1;
-    }
-
-    PLFLT *x = malloc(spike_count * sizeof(*x));
-    PLFLT *y = malloc(spike_count * sizeof(*y));
-    if (!x || !y) {
-        free(x);
-        free(y);
-        return -1;
-    }
-
-    //fill spike coordinates
-    size_t k = 0;
-    for (size_t t = 0; t < matrix->num_samples; t++) {
-        for (size_t n = 0; n < neurons_to_plot; n++) {
-            double value = matrix->states[t * matrix->num_features + n];
-            if (value > spike_threshold) {
-                x[k] = (PLFLT)t;
-                y[k] = (PLFLT)n;
-                k++;
-            }
-        }
-    }
-
-    //output to png
-    plsdev("svg");
-    plsfnam("reservoir_raster.svg");
-    
-    plsetopt("geometry", "1600x1200");
-    plscolbg(255, 255, 255);
-
-    plinit();
-
-    plscol0(1, 40, 40, 40); //gray axis
-    plscol0(2, 0, 0, 0); //blue points
-
-    plcol0(1);
-    plwidth(1.0);
-
-    plenv(
-        0.0,
-        (PLFLT)(matrix->num_samples - 1),
-        0.0,
-        (PLFLT)(neurons_to_plot - 1),
-        0,
-        0
-    );
-
-    pllab(
-        "Timestep",
-        "Neuron index",
-        "SPIRES Reservoir Raster Plot"
-    );
+static 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;
+	}
 
-    plcol0(2);
-    plwidth(1.0);
+	if (neurons_to_plot > matrix->num_features) {
+		neurons_to_plot = matrix->num_features;
+	}
 
-    for (size_t i = 0; i < spike_count; i++) {
-        PLFLT xline[2] = {x[i], x[i]};
-        PLFLT yline[2] = {y[i] - 0.35, y[i] + 0.35};
+	// count spikes
+	size_t spike_count = 0;
+	for (size_t t = 0; t < matrix->num_samples; t++) {
+		for (size_t n = 0; n < neurons_to_plot; n++) {
+			double value =
+			    matrix->states[t * matrix->num_features + n];
+			if (value > spike_threshold) {
+				spike_count++;
+			}
+		}
+	}
 
-        plline(2, xline, yline);
-    }
-    
-    plend();
+	if (spike_count == 0) {
+		fprintf(stderr, "No spikes found above threshold %.3f\n",
+			spike_threshold);
+		return -1;
+	}
 
-    free(x);
-    free(y);
-    return 0;
-}
-
-static int print_software_outputs(
-    spires_reservoir *reservoir,
-    const double *input_series,
-    size_t series_length,
-    size_t num_inputs,
-    size_t num_outputs,
-    size_t samples_to_print,
-    const char *label
-)
-{
-    if (!reservoir || !input_series || !label) {
-        return -1;
-    }
+	PLFLT *x = malloc(spike_count * sizeof(*x));
+	PLFLT *y = malloc(spike_count * sizeof(*y));
+	if (!x || !y) {
+		free(x);
+		free(y);
+		return -1;
+	}
 
-    double *output = malloc(num_outputs * sizeof(*output));
-    if (!output) {
-        return -1;
-    }
+	// fill spike coordinates
+	size_t k = 0;
+	for (size_t t = 0; t < matrix->num_samples; t++) {
+		for (size_t n = 0; n < neurons_to_plot; n++) {
+			double value =
+			    matrix->states[t * matrix->num_features + n];
+			if (value > spike_threshold) {
+				x[k] = (PLFLT)t;
+				y[k] = (PLFLT)n;
+				k++;
+			}
+		}
+	}
 
-    if (spires_reservoir_reset(reservoir) != SPIRES_OK) {
-        free(output);
-        return -1;
-    }
+	// output to png
+	plsdev("svg");
+	plsfnam("reservoir_raster.svg");
 
-    if (samples_to_print > series_length) {
-        samples_to_print = series_length;
-    }
+	plsetopt("geometry", "1600x1200");
+	plscolbg(255, 255, 255);
 
-    printf("\n%s\n", label);
+	plinit();
 
-    for (size_t timestep = 0;
-         timestep < series_length;
-         timestep++) {
+	plscol0(1, 40, 40, 40); // gray axis
+	plscol0(2, 0, 0, 0);	// blue points
 
-        const double *current_input =
-            &input_series[timestep * num_inputs];
+	plcol0(1);
+	plwidth(1.0);
 
-        if (spires_step(reservoir, current_input) != SPIRES_OK) {
-            free(output);
-            return -1;
-        }
+	plenv(0.0, (PLFLT)(matrix->num_samples - 1), 0.0,
+	      (PLFLT)(neurons_to_plot - 1), 0, 0);
 
-        if (spires_compute_output(reservoir, output) != SPIRES_OK) {
-            free(output);
-            return -1;
-        }
+	pllab("Timestep", "Neuron index", "SPIRES Reservoir Raster Plot");
 
-        if (timestep < samples_to_print) {
-            printf("timestep %zu:", timestep);
+	plcol0(2);
+	plwidth(1.0);
 
-            for (size_t output_index = 0;
-                 output_index < num_outputs;
-                 output_index++) {
+	for (size_t i = 0; i < spike_count; i++) {
+		PLFLT xline[2] = {x[i], x[i]};
+		PLFLT yline[2] = {y[i] - 0.35, y[i] + 0.35};
 
-                printf(
-                    " output[%zu]=%+.8e",
-                    output_index,
-                    output[output_index]
-                );
-            }
+		plline(2, xline, yline);
+	}
 
-            printf("\n");
-        }
-    }
+	plend();
 
-    free(output);
-    return 0;
+	free(x);
+	free(y);
+	return 0;
 }