A demonstration of Hopfield Network to evaluate dynamics of recall process.

hopnet_dynamics.c

/* *****************************************************************************
A demonstration of Hopfield Network to evaluate dynamics of recall process.

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 <stdbool.h>
#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 threshold;             // Activation threshold for neurons
    int max_iterations;           // Maximum number of iterations for recall
    double similarity;            // Similarity of initial pattern
    bool self_connection;         // Flag to allow self-connection in weights
    char filename[FILENAME_MAX];  // Filename to save direction cosines
} HopfieldConfig;

// Function prototype: argument parser
void print_usage(const char *program_name);
void configure(int argc, char *argv[], HopfieldConfig *config);

// Function prototype: memory allocator
double **allocate_double_matrix(int rows, int cols);
double *allocate_double_vector(int size);
int *allocate_int_vector(int size);
int **allocate_int_matrix(int rows, int cols);
void free_double_matrix(double **matrix, int rows);
void free_int_matrix(int **matrix, int rows);

// Fuction prototype: utility
FILE *open_file(const char *filename, const char *mode);
void shuffle(int *array, size_t n);
double direction_cosine(const HopfieldConfig *config, const int *pattern1,
                        const int *pattern2);
void define_patterns(const HopfieldConfig *config, int **patterns);
void init_state(const HopfieldConfig *config, const int *pattern,
                int *initial_state);

// Function prototype: Hopfield network
void learn(const HopfieldConfig *config, const int *const *patterns);
void recall(const HopfieldConfig *config, const int *initial_state,
            const int *reference_state, double *dircos_array);

// Function prototype: evaluation
void eval_dynamics(const HopfieldConfig *config, const int *const *patterns);

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("  -i Maximum iterations (default: 30)\n");
    printf("  -s Similarity of initial pattern (default: 0.5)\n");
    printf("  -c Flag to allow self-connection in weights (default: false)\n");
    printf("  -f Filename to save direction cosines (default: log.txt)\n");
    printf("  -h Show this help message and exit\n");
    printf("Example:\n");
    printf("  %s -n 30 -p 5 -t 0.1 -i 50 -s 0.3 -c true -f mylog.txt\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->threshold = 0.0;
    config->max_iterations = 30;
    config->similarity = 0.5;
    config->self_connection = false;
    strcpy(config->filename, "log.txt");

    while ((opt = getopt(argc, argv, "n:p:t:i:s:c:f: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->threshold = atof(optarg);
                break;
            case 'i':  // Maximum iterations
                config->max_iterations = atoi(optarg);
                break;
            case 's':  // Similarity
                config->similarity = atof(optarg);
                break;
            case 'c':  // Flag to allow self-connection in weights
                if (strcmp(optarg, "true") == 0) {
                    config->self_connection = true;
                } else if (strcmp(optarg, "false") == 0) {
                    config->self_connection = false;
                } else {
                    fprintf(stderr, "Invalid argument: %s\n", optarg);
                    print_usage(argv[0]);
                    exit(EXIT_FAILURE);
                }
                break;
            case 'f':  // Filename to save direction cosines
                strcpy(config->filename, 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_vector(int size) {
    /* Allocates memory for an double 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 (double *): A pointer to the allocated vector.
    */
    double *vector = (double *)malloc(size * sizeof(double));
    if (vector == NULL) {
        perror("malloc failed");
        exit(EXIT_FAILURE);
    }
    for (int i = 0; i < size; i++) {
        vector[i] = 0.0;
    }
    return vector;
}

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] = allocate_double_vector(cols);
    }
    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);
    }
    for (int i = 0; i < size; i++) {
        vector[i] = 0;
    }
    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);
    }
    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);
}

FILE *open_file(const char *filename, const char *mode) {
    /* A wrapper function for fopen that provides basic error handling.

    This function attempts to open the specified file with the given mode.
    If fopen fails, it prints an error message to stderr using perror and exits
    the program.

    Args:
        filename (const char*): The name of the file to open.
        mode (const char*): The mode in which to open the file.

    Return:
        fp (FILE *): A file pointer on success. If the open fails, the function
            prints an error message and terminates the program, so it will never
            return NULL.
    */
    FILE *fp = fopen(filename, mode);
    if (fp == NULL) {
        perror("fopen failed");  // Print error message to stderr
        exit(EXIT_FAILURE);      // Exit the program with a failure status
    }
    return fp;
}

void shuffle(int *array, size_t n) {
    /* Shuffles the elements of an array using the Fisher-Yates shuffle
    algorithm.

    This function randomly reorders the elements of the given array.
    It utilizes the Fisher-Yates shuffle algorithm for efficient shuffling.

    Args:
        array (int *): A pointer to the integer array to be shuffled.
        n (int): The number of elements in the array.
    */
    if (n > 1) {
        size_t i;
        for (i = 0; i < n - 1; i++) {
            size_t j = i + rand() / (RAND_MAX / (n - i) + 1);
            int t = array[j];
            array[j] = array[i];
            array[i] = t;
        }
    }
}

double direction_cosine(const HopfieldConfig *config, const int *pattern1,
                        const int *pattern2) {
    /* Calculate a direction cosine (similarity) between two patterns.

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

    Return:
        dircos (double): Direction cosine (similarity).
    */
    int num_neurons = config->num_neurons;
    int distance = 0;
    for (int i = 0; i < num_neurons; i++) {
        distance += pattern1[i] * pattern2[i];
    }
    double dircos = (double)distance / num_neurons;
    return dircos;
}

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 init_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 hamming_dist =
        (1 - config->similarity) / 2;  // a normalized Hamming distance
    int num_to_extract = (int)(num_neurons * hamming_dist);
    int numbers[num_neurons];
    for (int i = 0; i < num_neurons; i++) {
        numbers[i] = i;
    }
    shuffle(numbers, num_neurons);
    memcpy(initial_state, pattern, num_neurons * sizeof(int));
    for (int i = 0; i < num_to_extract; i++) {
        initial_state[numbers[i]] *= -1;
    }
}

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_neurons;
        }
    }

    if (config->self_connection == false) {
        for (int i = 0; i < num_neurons; i++) {
            weights[i][i] = 0.0;
        }
    }
}

void recall(const HopfieldConfig *config, const int *initial_state,
            const int *reference_state, double *dircos_array) {
    /* Perform synchronous update to recall the state.

    This function also calculates the direcition cosine (similarity) between the
    recalled state and the reference state during the recall process.

    Args:
        config (const HopfieldConfig *): The configuration of the network.
        initial_state (const int *): The initial state.
        reference_state (const int*): The reference state.
        dircos_array (double *): The array of direction cosines.
    */
    int num_neurons = config->num_neurons;
    double **weights = config->weights;
    double threshold = config->threshold;
    int max_iterations = config->max_iterations;
    int recalled_state[num_neurons];
    int new_state[num_neurons];
    memcpy(recalled_state, initial_state, num_neurons * sizeof(int));
    memset(new_state, 0, num_neurons * sizeof(int));
    dircos_array[0] = direction_cosine(config, recalled_state, reference_state);
    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 >= threshold) ? 1 : -1;
        }
        memcpy(recalled_state, new_state, num_neurons * sizeof(int));
        dircos_array[iter + 1] =
            direction_cosine(config, recalled_state, reference_state);
    }
}

void eval_dynamics(const HopfieldConfig *config, const int *const *patterns) {
    /* Evaluate dynamics of recall process in the 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 max_iterations = config->max_iterations;
    FILE *fp = open_file(config->filename, "w");
    int initial_state[num_neurons];
    double dircos_array[max_iterations + 1];
    int index_eval = rand() % num_patterns;
    init_state(config, patterns[index_eval], initial_state);
    recall(config, initial_state, patterns[index_eval], dircos_array);
    for (int i = 0; i < max_iterations + 1; i++) {
        fprintf(fp, "%f\n", dircos_array[i]);
    }
    fclose(fp);
}

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. Evaluate dynamics of recall process
    eval_dynamics(&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;
}