Trained software weights with ridge regression

1 files changed, 152 insertions, 15 deletions

diff --git a/src/main.c b/src/main.c
index 467b419..e8a6d92 100644
--- a/src/main.c
+++ b/src/main.c
@@ -12,8 +12,10 @@
 #define NUM_INPUTS 4
 #define NUM_OUTPUTS 2
 #define NUM_TRAINING_STEPS 500
+#define NUM_STEPS 2000
 
 #define SPIKE_THRESHOLD 0.5
+#define SPIKE_AMPLITUDE 0.1
 
 #define PI 3.14159265358979323846
 
@@ -23,6 +25,17 @@ static int plot_raster(
     double spike_threshold
 ); 
 
+
+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
+);
+
 int main(void) {
     //discrete LIF parameters for spires
     double lif_config[] = {
@@ -58,15 +71,24 @@ int main(void) {
         return -1;
     }
 
-    //create data set (sin wave time series prediction)
+    //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;
+        }
+    }
+
     Reservoir_State_Matrix state_matrix = {0};
     if (collect_reservoir_states(reservoir, training_inputs, NUM_TRAINING_STEPS,
                 &state_matrix) != 0) {
@@ -77,6 +99,37 @@ int main(void) {
     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;
+    }
+
+    print_software_outputs(
+        reservoir,
+        training_inputs,
+        NUM_TRAINING_STEPS,
+        NUM_INPUTS,
+        NUM_OUTPUTS,
+        10,
+        "Software outputs after training:"
+    );
+
     //generate raster plot for verification
     if (plot_raster(&state_matrix, NUM_NEURONS, SPIKE_THRESHOLD) != 0) {
         fprintf(stderr, "Failed to plot raster\n");
@@ -88,6 +141,8 @@ int main(void) {
         fprintf(stderr, "Failed to allocate memory for initial resistances");
         return -1;
     }
+
+    //resistances are inversely proportional to the software weigts
     for (size_t i = 0; i < NUM_NEURONS * NUM_OUTPUTS; i++) {
         initial_resistances[i] = 80000;
     }
@@ -96,12 +151,20 @@ int main(void) {
     double *spikes_voltages = malloc(state_matrix.num_features * 
             state_matrix.num_samples * sizeof(*spikes_voltages));
 
+    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;
+    }
+
     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; 
 
            spikes_voltages[index] = 
-               state_matrix.states [index] > SPIKE_THRESHOLD ? 0.1 : 0.0;
+               state_matrix.states [index] > SPIKE_THRESHOLD ? SPIKE_AMPLITUDE : 0.0;
         }
     }
 
@@ -126,6 +189,7 @@ int main(void) {
         spires_reservoir_destroy(reservoir);
         return -1;
     }
+    printf("Generated crossbar!!");
 
     //call ngspice for crossbar
     if (run_ngspice("crossbar.cir") < 0) {
@@ -135,7 +199,9 @@ int main(void) {
         spires_reservoir_destroy(reservoir);
         return -1;
     }
+    printf("ran ngspice!!");
 
+    //crossbar parameters needed for reading
     Crossbar_Output_Matrix crossbar_output = {
         .num_samples = NUM_TRAINING_STEPS, //this isnt right?
         .num_outputs = NUM_OUTPUTS,
@@ -151,20 +217,22 @@ int main(void) {
         free_crossbar_output_matrix(&crossbar_output);
         return -1;
     }
+    printf("read the crossbar outputs!!");
 
-    //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 * NUM_OUTPUTS + output]);
-        }
-        printf("\n");
-    }
+    //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");
+    // }
 
-    //clean up
     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);
@@ -223,8 +291,8 @@ static int plot_raster(
     }
 
     //output to png
-    plsdev("pngcairo");
-    plsfnam("reservoir_raster.png");
+    plsdev("svg");
+    plsfnam("reservoir_raster.svg");
     
     plsetopt("geometry", "1600x1200");
     plscolbg(255, 255, 255);
@@ -268,3 +336,72 @@ static int plot_raster(
     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;
+    }
+
+    double *output = malloc(num_outputs * sizeof(*output));
+    if (!output) {
+        return -1;
+    }
+
+    if (spires_reservoir_reset(reservoir) != SPIRES_OK) {
+        free(output);
+        return -1;
+    }
+
+    if (samples_to_print > series_length) {
+        samples_to_print = series_length;
+    }
+
+    printf("\n%s\n", label);
+
+    for (size_t timestep = 0;
+         timestep < series_length;
+         timestep++) {
+
+        const double *current_input =
+            &input_series[timestep * num_inputs];
+
+        if (spires_step(reservoir, current_input) != SPIRES_OK) {
+            free(output);
+            return -1;
+        }
+
+        if (spires_compute_output(reservoir, output) != SPIRES_OK) {
+            free(output);
+            return -1;
+        }
+
+        if (timestep < samples_to_print) {
+            printf("timestep %zu:", timestep);
+
+            for (size_t output_index = 0;
+                 output_index < num_outputs;
+                 output_index++) {
+
+                printf(
+                    " output[%zu]=%+.8e",
+                    output_index,
+                    output[output_index]
+                );
+            }
+
+            printf("\n");
+        }
+    }
+
+    free(output);
+    return 0;
+}