A demonstration code for Hopfield Network to learn a sequence of bit patterns.

hopnet_1d.c

/* *****************************************************************************
A demonstration for Hopfield Network to learn a sequence of bit patterns.

Copyright (C) 2025 by Akira TAMAMORI

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
***************************************************************************** */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <unistd.h>

typedef struct {
    int num_neurons;          // Number of neurons
    int num_patterns;         // Number of patterns to be learned
    double **weights;         // Weights of network
    double neuron_threshold;  // Activation threshold for neurons
    double noise_strength;    // Strength of bit-flip noise to each neuron
    char update_rule[6];      // Update rule for state (sync/async)
    int max_iterations;       // Maximum number of iterations for recall
    double recall_threshold;  // Threshold for recall convergence
    int async_iterations;     // Number of iterations for asynchronous update
} HopfieldConfig;

void print_usage(const char *program_name) {
    /* Displays the usage information of the program.

    Prints the available command-line options and usage examples.
    This function is called when invalid arguments are provided
    or when the help option is requested.

    Args:
        program_name (const char*): The name of the executable program.
    */
    printf("Usage: %s [options]\n", program_name);
    printf("Options:\n");
    printf("  -n Number of neurons (default: 20)\n");
    printf("  -p Number of patterns (default: 3)\n");
    printf("  -t Neuron threshold (default: 0.0)\n");
    printf("  -s Strength of bit-flip noise (default: 0.2)\n");
    printf("  -u Update rule (sync/async) (default: sync)\n");
    printf("  -i Maximum iterations (default: 100)\n");
    printf("  -r Recall threshold (default: 1e-6)\n");
    printf("  -a Async iterations (default: 100)\n");
    printf("  -h Show this help message and exit\n");
    printf("Example:\n");
    printf("  %s -n 30 -t 0.1 -u async -p 5 -i 200 -r 1e-5 -s 0.3 -a 150\n",
           program_name);
}

void configure(int argc, char *argv[], HopfieldConfig *config) {
    /* Configures the Hopfield network based on command-line arguments.

    Initializes default configuration values and updates them based on the
    provided command-line options. If invalid options are detected, it
    displays an error message and exits the program.

    Args:
        argc (int): The number of command-line arguments.
        argv (char *): The array of command-line argument strings.
        config (HopfieldConfig *): Pointer to the HopfieldConfig structure to be
            initialized.
    */
    int opt;

    // Set up default configuration
    config->num_neurons = 20;
    config->num_patterns = 3;
    config->weights = NULL;
    config->neuron_threshold = 0.0;
    config->noise_strength = 0.2;
    strcpy(config->update_rule, "sync");
    config->max_iterations = 100;
    config->recall_threshold = 1e-6;
    config->async_iterations = 100;

    while ((opt = getopt(argc, argv, "n:p:t:s:u:i:r:a:h")) != -1) {
        switch (opt) {
            case 'n':  // Number of neurons
                config->num_neurons = atoi(optarg);
                break;
            case 'p':  // Number of patterns
                config->num_patterns = atoi(optarg);
                break;
            case 't':  // Neuron threshold
                config->neuron_threshold = atof(optarg);
                break;
            case 's':  // Noise strength
                config->noise_strength = atof(optarg);
                break;
            case 'u':  // Update rule (sync/async)
                strcpy(config->update_rule, optarg);
                break;
            case 'i':  // Maximum iterations
                config->max_iterations = atoi(optarg);
                break;
            case 'r':  // Recall threshold
                config->recall_threshold = atof(optarg);
                break;
            case 'a':  // Async iterations
                config->async_iterations = atoi(optarg);
                break;
            case 'h':
                print_usage(argv[0]);
                exit(EXIT_SUCCESS);
            default:
                fprintf(stderr, "Unknown option: %c\n", optopt);
                print_usage(argv[0]);
                exit(EXIT_FAILURE);
        }
    }
    if (optind < argc) {
        fprintf(stderr, "Unexpected arguments:\n");
        for (int i = optind; i < argc; i++) {
            fprintf(stderr, "%s ", argv[i]);
        }
        fprintf(stderr, "\n");
        print_usage(argv[0]);
        exit(EXIT_FAILURE);
    }
}

double **allocate_double_matrix(int rows, int cols) {
    /* Allocates memory for a 2D matrix of type double.

    Dynamically allocates memory for a matrix of the specified size
    and initializes all elements to 0.0. If memory allocation fails,
    an error message is displayed and the program exits.

    Args:
        rows (int): The number of rows in the matrix.
        cols (int): The number of columns in the matrix.

    Return:
        matrix (double **): A pointer to the allocated matrix.
    */
    double **matrix = (double **)malloc(rows * sizeof(double *));
    if (matrix == NULL) {
        perror("malloc failed");
        exit(EXIT_FAILURE);
    }
    for (int i = 0; i < rows; i++) {
        matrix[i] = (double *)malloc(cols * sizeof(double));
        if (matrix[i] == NULL) {
            perror("malloc failed");
            exit(EXIT_FAILURE);
        }
        for (int j = 0; j < cols; j++) {
            matrix[i][j] = 0.0;
        }
    }
    return matrix;
}

void free_double_matrix(double **matrix, int rows) {
    /* Frees the memory allocated for a 2D matrix of type double.

    Frees the memory allocated for each row and then frees the matrix itself.
    If a NULL pointer is passed, the function performs no operation.

    Args:
        matrix (double **): Pointer to the matrix to be freed.
        rows (int): The number of rows in the matrix.
    */
    if (matrix == NULL) {
        return;
    }
    for (int i = 0; i < rows; i++) {
        free(matrix[i]);
    }
    free(matrix);
}

int *allocate_int_vector(int size) {
    /* Allocates memory for an integer vector (1D array).

    Dynamically allocates memory for a vector of the specified size.
    If memory allocation fails, an error message is displayed and
    the program exits.

    Args:
        size (int): The number of elements in the vector.

    Return:
        vector (int *): A pointer to the allocated vector.
    */
    int *vector = (int *)malloc(size * sizeof(int));
    if (vector == NULL) {
        perror("malloc failed");
        exit(EXIT_FAILURE);
    }
    return vector;
}

int **allocate_int_matrix(int rows, int cols) {
    /* Allocates memory for a 2D matrix of type int.

    Dynamically allocates memory for a matrix of the specified size
    and initializes all elements to 0. If memory allocation fails,
    an error message is displayed and the program exits.

    Args:
        rows (int): The number of rows in the matrix.
        cols (int): The number of columns in the matrix.

    Return:
        matrix (int **): A pointer to the allocated matrix.
    */
    int **matrix = (int **)malloc(rows * sizeof(int *));
    if (matrix == NULL) {
        perror("malloc failed");
        exit(EXIT_FAILURE);
    }
    for (int i = 0; i < rows; i++) {
        matrix[i] = allocate_int_vector(cols);
        for (int j = 0; j < cols; j++) {
            matrix[i][j] = 0;
        }
    }
    return matrix;
}

void free_int_matrix(int **matrix, int rows) {
    /* Frees the memory allocated for a 2D matrix of type int.

    Frees the memory allocated for each row and then frees the matrix itself.
    If a NULL pointer is passed, the function performs no operation.

    Args:
        matrix (int **): Pointer to the matrix to be freed.
        rows (int): The number of rows in the matrix.
    */
    if (matrix == NULL) {
        return;
    }
    for (int i = 0; i < rows; i++) {
        free(matrix[i]);
    }
    free(matrix);
}

void learn(const HopfieldConfig *config, const int *const *patterns) {
    /* Learn the given patterns by Hebbian learning rule.

    Args:
        config (const HopfieldConfig *): The configuration of the network.
        pattern (int **): The bit patterns to be learned.
    */
    int num_neurons = config->num_neurons;
    int num_patterns = config->num_patterns;
    double **weights = config->weights;
    for (int p = 0; p < num_patterns; p++) {
        for (int i = 0; i < num_neurons; i++) {
            for (int j = 0; j < num_neurons; j++) {
                weights[i][j] += (double)patterns[p][i] * patterns[p][j];
            }
        }
    }
    for (int i = 0; i < num_neurons; i++) {
        for (int j = 0; j < num_neurons; j++) {
            weights[i][j] /= (double)num_patterns;
        }
    }
    for (int i = 0; i < num_neurons; i++) {
        weights[i][i] = 0.0;
    }
}

double energy(const HopfieldConfig *config, const int *state) {
    /* Calculate the energy of the network for a given state.

    Args:
        config (const HopfieldConfig *): The configuration of the network.
        state (const int *): The state of the network.

    Return:
        energy (double): The energy of the network.
    */
    int num_neurons = config->num_neurons;
    double **weights = config->weights;
    double neuron_threshold = config->neuron_threshold;
    double energy = 0.0;
    for (int i = 0; i < num_neurons; i++) {
        for (int j = 0; j < num_neurons; j++) {
            energy += -0.5 * state[i] * weights[i][j] * state[j];
        }
        energy += state[i] * neuron_threshold;
    }
    return energy;
}

void synchronous_update(const HopfieldConfig *config, int *recalled_state) {
    /* Perform synchronous update of the network state.

    Args:
        config (const HopfieldConfig *): The configuration of the network.
        recalled_state (int *): The current state to be updated.
    */
    int num_neurons = config->num_neurons;
    double **weights = config->weights;
    double neuron_threshold = config->neuron_threshold;
    int max_iterations = config->max_iterations;
    double recall_threshold = config->recall_threshold;
    double prev_energy = energy(config, recalled_state);
    int *new_state = allocate_int_vector(num_neurons);
    for (int iter = 0; iter < max_iterations; iter++) {
        for (int i = 0; i < num_neurons; i++) {
            double potential = 0.0;
            for (int j = 0; j < num_neurons; j++) {
                potential += weights[i][j] * recalled_state[j];
            }
            new_state[i] = (potential >= neuron_threshold) ? 1 : -1;
        }
        memcpy(recalled_state, new_state, num_neurons * sizeof(int));
        double current_energy = energy(config, recalled_state);
        if (prev_energy - current_energy < recall_threshold) {
            break;
        }
        prev_energy = current_energy;
    }
    free(new_state);
}

void asynchronous_update(const HopfieldConfig *config, int *recalled_state) {
    /* Perform asynchronous update of the network state.

    Args:
        config (const HopfieldConfig *): The configuration of the network.
        recalled_state (int *): The current state to be updated.
    */
    int num_neurons = config->num_neurons;
    double **weights = config->weights;
    double neuron_threshold = config->neuron_threshold;
    int async_iterations = config->async_iterations;
    for (int iter = 0; iter < async_iterations; iter++) {
        int i = rand() % num_neurons;
        double potential = 0.0;
        for (int j = 0; j < num_neurons; j++) {
            potential += weights[i][j] * recalled_state[j];
        }
        recalled_state[i] = (potential >= neuron_threshold) ? 1 : -1;
    }
}

void recall(const HopfieldConfig *config, const int *initial_state,
            int *recalled_state) {
    /* Recall a pattern from the network, starting from the initial state.

    Args:
        config (const HopfieldConfig *): The configuration of the network.
        initial_state (const int *): The initial state of the network.
        recalled_state (int *): The recalled state of the network.
    */
    int num_neurons = config->num_neurons;
    memcpy(recalled_state, initial_state, num_neurons * sizeof(int));
    if (strcmp(config->update_rule, "sync") == 0) {
        synchronous_update(config, recalled_state);
    } else {
        asynchronous_update(config, recalled_state);
    }
}

void define_patterns(const HopfieldConfig *config, int **patterns) {
    /* Define bit patterns to be learned.

    Args:
        config (const HopfieldConfig *): The configuration of the network.
        pattern (int **): A 2D array where each row represents a pattern.
    */
    int num_neurons = config->num_neurons;
    int num_patterns = config->num_patterns;
    for (int i = 0; i < num_patterns; i++) {
        for (int j = 0; j < num_neurons; j++) {
            patterns[i][j] = (rand() % 2 == 0) ? 1 : -1;
        }
    }
}

void noisy_initial_state(const HopfieldConfig *config, const int *pattern,
                         int *initial_state) {
    /* Create a noisy initial state for a given pattern.

    Args:
        config (const HopfieldConfig *): The configuration of the network.
        pattern (const int *): The original pattern.
        initial_state (int *): The initial state of the network.
    */
    int num_neurons = config->num_neurons;
    double noise_strength = config->noise_strength;
    for (int i = 0; i < num_neurons; i++) {
        double rand_val = (double)rand() / RAND_MAX;
        int mask = (rand_val < noise_strength) ? -1 : 1;
        initial_state[i] = mask * pattern[i];
    }
}

void print_pattern(const HopfieldConfig *config, const int *pattern) {
    /* Print a 1D pattern to the console.

    Args:
        config (const HopfieldConfig *): The configuration of the network.
        pattern (const int *): The bit pattern to print.
    */
    int num_neurons = config->num_neurons;
    for (int i = 0; i < num_neurons; i++) {
        if (pattern[i] == 1) {
            printf("%c", '+');
        } else {
            printf("%c", '-');
        }
    }
    printf("\n");
}

int hamming_distance(const HopfieldConfig *config, const int *pattern1,
                     const int *pattern2) {
    /* Calculate Hamming distance between two patterns.

    Args:
        config (const HopfieldConfig *): The configuration of the network.
        pattern1 (int *): The bit pattern.
        pattern2 (int *): The bit pattern.

    Return:
        distance (int): Hamming distance.
    */
    int num_neurons = config->num_neurons;
    int distance = 0;
    for (int i = 0; i < num_neurons; i++) {
        if (pattern1[i] != pattern2[i]) {
            distance++;
        }
    }
    return distance;
}

void search_nearest_pattern(const HopfieldConfig *config, int *min_distance,
                            int *min_index, const int *recalled_state,
                            const int *const *patterns) {
    /* Search the nearest pattern to the recalled state given a set of patterns.

    Args:
        config (const HopfieldConfig *): The configuration of the network.
        min_distance (int *): The minimum Hamming distance.
        min_index (int *): The index of the nearest pattern.
        recalled_state (const int *): The recalled state from the network.
        patterns (const int *const *): The patterns to test.
    */
    int num_neurons = config->num_neurons;
    int num_patterns = config->num_patterns;
    *min_distance = num_neurons + 1;
    *min_index = -1;
    for (int i = 0; i < num_patterns; i++) {
        int distance = hamming_distance(config, recalled_state, patterns[i]);
        if (distance < *min_distance) {
            *min_distance = distance;
            *min_index = i;
        }
    }
}

void recall_test(const HopfieldConfig *config, const int *const *patterns) {
    /* Run the recall test for the Hopfield network.

    Args:
        config (const HopfieldConfig *): The configuration of the network.
        patterns (const int *const *): The patterns to test.
    */
    int num_neurons = config->num_neurons;
    int num_patterns = config->num_patterns;
    int *initial_state = allocate_int_vector(num_neurons);
    int *recalled_state = allocate_int_vector(num_neurons);
    for (int i = 0; i < num_patterns; i++) {
        printf("Pattern %d: ", i);
        noisy_initial_state(config, patterns[i], initial_state);
        recall(config, initial_state, recalled_state);
        int min_distance;
        int min_index;
        search_nearest_pattern(config, &min_distance, &min_index,
                               recalled_state, patterns);
        if (min_distance == 0 && min_index == i) {
            printf("Recall successful!\n");
        } else {
            printf("Recall failed.\n");
        }
        printf("  Correct pattern: ");
        print_pattern(config, patterns[i]);
        printf("  Initial pattern: ");
        print_pattern(config, initial_state);
        printf("  Recalled state:  ");
        print_pattern(config, recalled_state);
        printf("  Closest pattern: %d\n", min_index);
        printf("  Minimum distance: %d\n", min_distance);
        printf("\n");
    }
    free(initial_state);
    free(recalled_state);
}

int main(int argc, char *argv[]) {
    HopfieldConfig config;
    srand((unsigned int)time(NULL));

    // 1. Configure the Hopfield network
    configure(argc, argv, &config);

    // 2. Define bit patterns to be learned
    int **patterns =
        allocate_int_matrix(config.num_patterns, config.num_neurons);
    define_patterns(&config, patterns);

    // 3. Learn patterns by Hebbian learning rule
    config.weights =
        allocate_double_matrix(config.num_neurons, config.num_neurons);
    learn(&config, (const int *const *)patterns);

    // 4. Run recall test (initial state is a noisy version of each pattern)
    recall_test(&config, (const int *const *)patterns);

    // 5. Cleaning up
    free_double_matrix(config.weights, config.num_neurons);
    free_int_matrix(patterns, config.num_patterns);

    return 0;
}