From 411f66a2540fa17c736116d865e0ceb0cfe5623b Mon Sep 17 00:00:00 2001 From: jeanne Date: Wed, 11 May 2022 09:54:38 -0700 Subject: Initial commit. --- src/bin/CMakeLists.txt | 3 + src/bin/mnist/CMakeLists.txt | 11 + src/bin/mnist/src/main.c | 473 +++++++++++++++++++++ src/lib/CMakeLists.txt | 37 ++ src/lib/include/neuralnet/matrix.h | 111 +++++ src/lib/include/neuralnet/neuralnet.h | 64 +++ src/lib/include/neuralnet/train.h | 42 ++ src/lib/include/neuralnet/types.h | 3 + src/lib/src/activation.h | 21 + src/lib/src/matrix.c | 298 +++++++++++++ src/lib/src/neuralnet.c | 228 ++++++++++ src/lib/src/neuralnet_impl.h | 36 ++ src/lib/src/train.c | 346 +++++++++++++++ src/lib/test/matrix_test.c | 350 +++++++++++++++ src/lib/test/neuralnet_test.c | 92 ++++ src/lib/test/test.h | 185 ++++++++ src/lib/test/test_main.c | 3 + src/lib/test/test_util.h | 22 + .../test/train_linear_perceptron_non_origin_test.c | 67 +++ src/lib/test/train_linear_perceptron_test.c | 62 +++ src/lib/test/train_sigmoid_test.c | 66 +++ src/lib/test/train_xor_test.c | 66 +++ 22 files changed, 2586 insertions(+) create mode 100644 src/bin/CMakeLists.txt create mode 100644 src/bin/mnist/CMakeLists.txt create mode 100644 src/bin/mnist/src/main.c create mode 100644 src/lib/CMakeLists.txt create mode 100644 src/lib/include/neuralnet/matrix.h create mode 100644 src/lib/include/neuralnet/neuralnet.h create mode 100644 src/lib/include/neuralnet/train.h create mode 100644 src/lib/include/neuralnet/types.h create mode 100644 src/lib/src/activation.h create mode 100644 src/lib/src/matrix.c create mode 100644 src/lib/src/neuralnet.c create mode 100644 src/lib/src/neuralnet_impl.h create mode 100644 src/lib/src/train.c create mode 100644 src/lib/test/matrix_test.c create mode 100644 src/lib/test/neuralnet_test.c create mode 100644 src/lib/test/test.h create mode 100644 src/lib/test/test_main.c create mode 100644 src/lib/test/test_util.h create mode 100644 src/lib/test/train_linear_perceptron_non_origin_test.c create mode 100644 src/lib/test/train_linear_perceptron_test.c create mode 100644 src/lib/test/train_sigmoid_test.c create mode 100644 src/lib/test/train_xor_test.c (limited to 'src') diff --git a/src/bin/CMakeLists.txt b/src/bin/CMakeLists.txt new file mode 100644 index 0000000..051a56f --- /dev/null +++ b/src/bin/CMakeLists.txt @@ -0,0 +1,3 @@ +cmake_minimum_required(VERSION 3.0) + +add_subdirectory(mnist) diff --git a/src/bin/mnist/CMakeLists.txt b/src/bin/mnist/CMakeLists.txt new file mode 100644 index 0000000..a6c54f2 --- /dev/null +++ b/src/bin/mnist/CMakeLists.txt @@ -0,0 +1,11 @@ +cmake_minimum_required(VERSION 3.0) + +add_executable(mnist + src/main.c) + +target_link_libraries(mnist PRIVATE + neuralnet + bsd + z) + +target_compile_options(mnist PRIVATE -Wall -Wextra) diff --git a/src/bin/mnist/src/main.c b/src/bin/mnist/src/main.c new file mode 100644 index 0000000..4d268ac --- /dev/null +++ b/src/bin/mnist/src/main.c @@ -0,0 +1,473 @@ +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +static const int TRAIN_ITERATIONS = 100; + +static const int32_t IMAGE_FILE_MAGIC = 0x00000803; +static const int32_t LABEL_FILE_MAGIC = 0x00000801; + +// Inputs of 0 cancel weights during training. This value is used to rescale the +// input pixels from [0,255] to [PIXEL_LOWER_BOUND, 1.0]. +static const double PIXEL_LOWER_BOUND = 0.01; + +// Scale the outputs to (0,1) since the sigmoid cannot produce 0 or 1. +static const double LABEL_LOWER_BOUND = 0.01; +static const double LABEL_UPPER_BOUND = 0.99; + +// Epsilon used to compare R values. +static const double EPS = 1e-10; + +#define min(a,b) ((a) < (b) ? (a) : (b)) + +typedef struct ImageSet { + nnMatrix images; // Images flattened into row vectors of the matrix. + nnMatrix labels; // One-hot-encoded labels. + int count; // Number of images and labels. + int rows; // Rows in an image. + int cols; // Columns in an image. +} ImageSet; + +static void usage(const char* argv0) { + fprintf(stderr, "Usage: %s [num images]\n", argv0); + fprintf(stderr, "\n"); + fprintf(stderr, " Use -1 for [num images] to use all the images in the data set\n"); +} + +static bool R_eq(R a, R b) { + return fabs(a-b) <= EPS; +} + +static void PrintImage(const nnMatrix* images, int rows, int cols, int image_index) { + assert(images); + assert((0 <= image_index) && (image_index < images->rows)); + + // Top line. + for (int j = 0; j < cols/2; ++j) { + printf(" -"); + } + printf("\n"); + + // Image. + const R* value = nnMatrixRow(images, image_index); + for (int i = 0; i < rows; ++i) { + printf("|"); + for (int j = 0; j < cols; ++j) { + if (*value > 0.8) { + printf("#"); + } else if (*value > 0.5) { + printf("*"); + } + else if (*value > PIXEL_LOWER_BOUND) { + printf(":"); + } else if (*value == 0.0) { + // Values should not be exactly 0, otherwise they cancel out weights + // during training. + printf("X"); + } else { + printf(" "); + } + value++; + } + printf("|\n"); + } + + // Bottom line. + for (int j = 0; j < cols/2; ++j) { + printf(" -"); + } + printf("\n"); +} + +static void PrintLabel(const nnMatrix* labels, int label_index) { + assert(labels); + assert((0 <= label_index) && (label_index < labels->rows)); + + // Compute the label from the one-hot encoding. + const R* value = nnMatrixRow(labels, label_index); + int label = -1; + for (int i = 0; i < 10; ++i) { + if (R_eq(*value++, LABEL_UPPER_BOUND)) { + label = i; + break; + } + } + assert((0 <= label) && (label <= 9)); + + printf("Label: %d ( ", label); + value = nnMatrixRow(labels, label_index); + for (int i = 0; i < 10; ++i) { + printf("%.3f ", *value++); + } + printf(")\n"); +} + +static R lerp(R a, R b, R t) { + return a + t*(b-a); +} + +/// Rescales a pixel from [0,255] to [PIXEL_LOWER_BOUND, 1.0]. +static R FormatPixel(uint8_t pixel) { + const R value = (R)(pixel) / 255.0 * (1.0 - PIXEL_LOWER_BOUND) + PIXEL_LOWER_BOUND; + assert(value >= PIXEL_LOWER_BOUND); + assert(value <= 1.0); + return value; +} + +/// Rescales a one-hot-encoded label value to (0,1). +static R FormatLabel(R label) { + const R value = lerp(LABEL_LOWER_BOUND, LABEL_UPPER_BOUND, label); + assert(value > 0.0); + assert(value < 1.0); + return value; +} + +static int32_t ReverseEndian32(int32_t x) { + const int32_t x0 = x & 0xff; + const int32_t x1 = (x >> 8) & 0xff; + const int32_t x2 = (x >> 16) & 0xff; + const int32_t x3 = (x >> 24) & 0xff; + return (x0 << 24) | (x1 << 16) | (x2 << 8) | x3; +} + +static void ImageToMatrix( + const uint8_t* pixels, int num_pixels, int row, nnMatrix* images) { + assert(pixels); + assert(images); + + for (int i = 0; i < num_pixels; ++i) { + const R pixel = FormatPixel(pixels[i]); + nnMatrixSet(images, row, i, pixel); + } +} + +static bool ReadImages(gzFile images_file, int max_num_images, ImageSet* image_set) { + assert(images_file != Z_NULL); + assert(image_set); + + bool success = false; + + uint8_t* pixels = 0; + + int32_t magic, total_images, rows, cols; + if ( (gzread(images_file, (char*)&magic, sizeof(int32_t)) != sizeof(int32_t)) || + (gzread(images_file, (char*)&total_images, sizeof(int32_t)) != sizeof(int32_t)) || + (gzread(images_file, (char*)&rows, sizeof(int32_t)) != sizeof(int32_t)) || + (gzread(images_file, (char*)&cols, sizeof(int32_t)) != sizeof(int32_t)) ) { + fprintf(stderr, "Failed to read header\n"); + goto cleanup; + } + + magic = ReverseEndian32(magic); + total_images = ReverseEndian32(total_images); + rows = ReverseEndian32(rows); + cols = ReverseEndian32(cols); + + if (magic != IMAGE_FILE_MAGIC) { + fprintf(stderr, "Magic number mismatch. Got %x, expected: %x\n", + magic, IMAGE_FILE_MAGIC); + goto cleanup; + } + + printf("Magic: %.8x\nTotal images: %d\nRows: %d\nCols: %d\n", + magic, total_images, rows, cols); + + total_images = max_num_images >= 0 ? min(total_images, max_num_images) : total_images; + + // Images are flattened into single row vectors. + const int num_pixels = rows * cols; + image_set->images = nnMatrixMake(total_images, num_pixels); + image_set->count = total_images; + image_set->rows = rows; + image_set->cols = cols; + + pixels = calloc(1, num_pixels); + if (!pixels) { + fprintf(stderr, "Failed to allocate image buffer\n"); + goto cleanup; + } + + for (int i = 0; i < total_images; ++i) { + const int bytes_read = gzread(images_file, pixels, num_pixels); + if (bytes_read < num_pixels) { + fprintf(stderr, "Failed to read image %d\n", i); + goto cleanup; + } + ImageToMatrix(pixels, num_pixels, i, &image_set->images); + } + + success = true; + +cleanup: + if (pixels) { + free(pixels); + } + if (!success) { + nnMatrixDel(&image_set->images); + } + return success; +} + +static void OneHotEncode(const uint8_t* labels_bytes, int num_labels, nnMatrix* labels) { + assert(labels_bytes); + assert(labels); + assert(labels->rows == num_labels); + assert(labels->cols == 10); + + static const R one_hot[10][10] = { + { 1, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 0, 1, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 0, 0, 1, 0, 0, 0, 0, 0, 0, 0 }, + { 0, 0, 0, 1, 0, 0, 0, 0, 0, 0 }, + { 0, 0, 0, 0, 1, 0, 0, 0, 0, 0 }, + { 0, 0, 0, 0, 0, 1, 0, 0, 0, 0 }, + { 0, 0, 0, 0, 0, 0, 1, 0, 0, 0 }, + { 0, 0, 0, 0, 0, 0, 0, 1, 0, 0 }, + { 0, 0, 0, 0, 0, 0, 0, 0, 1, 0 }, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 }, + }; + + R* value = labels->values; + + for (int i = 0; i < num_labels; ++i) { + const uint8_t label = labels_bytes[i]; + const R* one_hot_value = one_hot[label]; + + for (int j = 0; j < 10; ++j) { + *value++ = FormatLabel(*one_hot_value++); + } + } +} + +static int OneHotDecode(const nnMatrix* label_matrix) { + assert(label_matrix); + assert(label_matrix->cols == 1); + assert(label_matrix->rows == 10); + + R max_value = 0; + int pos_max = 0; + for (int i = 0; i < 10; ++i) { + const R value = nnMatrixAt(label_matrix, 0, i); + if (value > max_value) { + max_value = value; + pos_max = i; + } + } + assert(pos_max >= 0); + assert(pos_max <= 10); + return pos_max; +} + +static bool ReadLabels(gzFile labels_file, int max_num_labels, ImageSet* image_set) { + assert(labels_file != Z_NULL); + assert(image_set != 0); + + bool success = false; + + uint8_t* labels = 0; + + int32_t magic, total_labels; + if ( (gzread(labels_file, (char*)&magic, sizeof(int32_t)) != sizeof(int32_t)) || + (gzread(labels_file, (char*)&total_labels, sizeof(int32_t)) != sizeof(int32_t)) ) { + fprintf(stderr, "Failed to read header\n"); + goto cleanup; + } + + magic = ReverseEndian32(magic); + total_labels = ReverseEndian32(total_labels); + + if (magic != LABEL_FILE_MAGIC) { + fprintf(stderr, "Magic number mismatch. Got %x, expected: %x\n", + magic, LABEL_FILE_MAGIC); + goto cleanup; + } + + printf("Magic: %.8x\nTotal labels: %d\n", magic, total_labels); + + total_labels = max_num_labels >= 0 ? min(total_labels, max_num_labels) : total_labels; + + assert(image_set->count == total_labels); + + // One-hot encoding of labels, 10 values (digits) per label. + image_set->labels = nnMatrixMake(total_labels, 10); + + labels = calloc(total_labels, sizeof(uint8_t)); + if (!labels) { + fprintf(stderr, "Failed to allocate labels buffer\n"); + goto cleanup; + } + + if (gzread(labels_file, labels, total_labels * sizeof(uint8_t)) != total_labels) { + fprintf(stderr, "Failed to read labels\n"); + goto cleanup; + } + + OneHotEncode(labels, total_labels, &image_set->labels); + + success = true; + +cleanup: + if (labels) { + free(labels); + } + if (!success) { + nnMatrixDel(&image_set->labels); + } + return success; +} + +int main(int argc, const char** argv) { + if (argc < 2) { + usage(argv[0]); + return 1; + } + + bool success = false; + + gzFile train_images_file = Z_NULL; + gzFile train_labels_file = Z_NULL; + gzFile test_images_file = Z_NULL; + gzFile test_labels_file = Z_NULL; + ImageSet train_set = { 0 }; + ImageSet test_set = { 0 }; + nnNeuralNetwork* net = 0; + nnQueryObject* query = 0; + + const char* mnist_files_dir = argv[1]; + const int max_num_images = argc > 2 ? atoi(argv[2]) : -1; + + char train_labels_path[PATH_MAX]; + char train_images_path[PATH_MAX]; + char test_labels_path[PATH_MAX]; + char test_images_path[PATH_MAX]; + strlcpy(train_labels_path, mnist_files_dir, PATH_MAX); + strlcpy(train_images_path, mnist_files_dir, PATH_MAX); + strlcpy(test_labels_path, mnist_files_dir, PATH_MAX); + strlcpy(test_images_path, mnist_files_dir, PATH_MAX); + strlcat(train_labels_path, "/train-labels-idx1-ubyte.gz", PATH_MAX); + strlcat(train_images_path, "/train-images-idx3-ubyte.gz", PATH_MAX); + strlcat(test_labels_path, "/t10k-labels-idx1-ubyte.gz", PATH_MAX); + strlcat(test_images_path, "/t10k-images-idx3-ubyte.gz", PATH_MAX); + + train_images_file = gzopen(train_images_path, "r"); + if (train_images_file == Z_NULL) { + fprintf(stderr, "Failed to open file: %s\n", train_images_path); + goto cleanup; + } + + train_labels_file = gzopen(train_labels_path, "r"); + if (train_labels_file == Z_NULL) { + fprintf(stderr, "Failed to open file: %s\n", train_labels_path); + goto cleanup; + } + + test_images_file = gzopen(test_images_path, "r"); + if (test_images_file == Z_NULL) { + fprintf(stderr, "Failed to open file: %s\n", test_images_path); + goto cleanup; + } + + test_labels_file = gzopen(test_labels_path, "r"); + if (test_labels_file == Z_NULL) { + fprintf(stderr, "Failed to open file: %s\n", test_labels_path); + goto cleanup; + } + + if (!ReadImages(train_images_file, max_num_images, &train_set)) { + goto cleanup; + } + if (!ReadLabels(train_labels_file, max_num_images, &train_set)) { + goto cleanup; + } + + if (!ReadImages(test_images_file, max_num_images, &test_set)) { + goto cleanup; + } + if (!ReadLabels(test_labels_file, max_num_images, &test_set)) { + goto cleanup; + } + + printf("\nTraining image/label pair examples:\n"); + for (int i = 0; i < min(3, train_set.images.rows); ++i) { + PrintImage(&train_set.images, train_set.rows, train_set.cols, i); + PrintLabel(&train_set.labels, i); + printf("\n"); + } + + // Network definition. + const int image_size_pixels = train_set.rows * train_set.cols; + const int num_layers = 2; + const int layer_sizes[3] = { image_size_pixels, 100, 10 }; + const nnActivation layer_activations[2] = { nnSigmoid, nnSigmoid }; + if (!(net = nnMakeNet(num_layers, layer_sizes, layer_activations))) { + fprintf(stderr, "Failed to create neural network\n"); + goto cleanup; + } + + // Train. + printf("Training with up to %d images from the data set\n\n", max_num_images); + const nnTrainingParams training_params = { + .learning_rate = 0.1, + .max_iterations = TRAIN_ITERATIONS, + .seed = 0, + .weight_init = nnWeightInitNormal, + .debug = true, + }; + nnTrain(net, &train_set.images, &train_set.labels, &training_params); + + // Test. + int hits = 0; + query = nnMakeQueryObject(net, /*num_inputs=*/1); + for (int i = 0; i < test_set.count; ++i) { + const nnMatrix test_image = nnMatrixBorrowRows(&test_set.images, i, 1); + const nnMatrix test_label = nnMatrixBorrowRows(&test_set.labels, i, 1); + + nnQuery(net, query, &test_image); + + const int test_label_expected = OneHotDecode(&test_label); + const int test_label_actual = OneHotDecode(nnNetOutputs(query)); + + if (test_label_actual == test_label_expected) { + ++hits; + } + } + const R hit_ratio = (R)hits / (R)test_set.count; + printf("Test images: %d\n", test_set.count); + printf("Hits: %d/%d (%.3f%%)\n", hits, test_set.count, hit_ratio*100); + + success = true; + +cleanup: + if (query) { + nnDeleteQueryObject(&query); + } + if (net) { + nnDeleteNet(&net); + } + nnMatrixDel(&train_set.images); + nnMatrixDel(&test_set.images); + if (train_images_file != Z_NULL) { + gzclose(train_images_file); + } + if (train_labels_file != Z_NULL) { + gzclose(train_labels_file); + } + if (test_images_file != Z_NULL) { + gzclose(test_images_file); + } + if (test_labels_file != Z_NULL) { + gzclose(test_labels_file); + } + return success ? 0 : 1; +} diff --git a/src/lib/CMakeLists.txt b/src/lib/CMakeLists.txt new file mode 100644 index 0000000..9e0e924 --- /dev/null +++ b/src/lib/CMakeLists.txt @@ -0,0 +1,37 @@ +cmake_minimum_required(VERSION 3.0) + +# Library + +add_library(neuralnet + src/matrix.c + src/neuralnet.c + src/train.c) + +target_include_directories(neuralnet PUBLIC + include) + +target_link_libraries(neuralnet PRIVATE + math # System math library. + random) + +target_compile_options(neuralnet PRIVATE -Wall -Wextra) + +# Test + +add_executable(neuralnet-test + test/matrix_test.c + test/neuralnet_test.c + test/test_main.c + test/train_linear_perceptron_test.c + test/train_linear_perceptron_non_origin_test.c + test/train_sigmoid_test.c + test/train_xor_test.c) + +target_link_libraries(neuralnet-test PRIVATE + neuralnet) + +# So that we can include header files from the private implementation. +target_include_directories(neuralnet-test PRIVATE + src) + +target_compile_options(neuralnet-test PRIVATE -DUNIT_TEST -Wall -Wextra) diff --git a/src/lib/include/neuralnet/matrix.h b/src/lib/include/neuralnet/matrix.h new file mode 100644 index 0000000..9816b81 --- /dev/null +++ b/src/lib/include/neuralnet/matrix.h @@ -0,0 +1,111 @@ +#pragma once + +#include + +#include + +/// NxM matrix. +typedef struct nnMatrix { + int rows; + int cols; + R* values; +} nnMatrix; + +/// Construct a matrix. +nnMatrix nnMatrixMake(int rows, int cols); + +/// Delete a matrix and free its internal memory. +void nnMatrixDel(nnMatrix*); + +/// Move a matrix. +/// +/// |in| is an empty matrix after the move. +/// |out| is a matrix like |in| before the move. +void nnMatrixMove(nnMatrix* in, nnMatrix* out); + +/// Deep-copy a matrix. +void nnMatrixCopy(const nnMatrix* in, nnMatrix* out); + +/// Write the matrix values into an array in a row-major fashion. +void nnMatrixToArray(const nnMatrix* in, R* out); + +/// Write the given row of a matrix into an array. +void nnMatrixRowToArray(const nnMatrix* in, int row, R* out); + +/// Copy a column from a source to a target matrix. +void nnMatrixCopyCol(const nnMatrix* in, nnMatrix* out, int col_in, int col_out); + +/// Mutable borrow of a matrix. +nnMatrix nnMatrixBorrow(nnMatrix* in); + +/// Mutable borrow of a subrange of rows of a matrix. +nnMatrix nnMatrixBorrowRows(nnMatrix* in, int row_start, int num_rows); + +/// Initialize the matrix from an array of values. +/// +/// The array must hold values in a row-major fashion. +void nnMatrixInit(nnMatrix*, const R* values); + +/// Initialize all matrix values to a given constant. +void nnMatrixInitConstant(nnMatrix*, R value); + +/// Multiply two matrices. +void nnMatrixMul(const nnMatrix* left, const nnMatrix* right, nnMatrix* out); + +/// Matrix multiply-add. +/// +/// out = left + (right * scale) +void nnMatrixMulAdd(const nnMatrix* left, const nnMatrix* right, R scale, nnMatrix* out); + +/// Matrix multiply-subtract. +/// +/// out = left - (right * scale) +void nnMatrixMulSub(const nnMatrix* left, const nnMatrix* right, R scale, nnMatrix* out); + +/// Hadamard product of two matrices. +void nnMatrixMulPairs(const nnMatrix* left, const nnMatrix* right, nnMatrix* out); + +/// Add two matrices. +void nnMatrixAdd(const nnMatrix* left, const nnMatrix* right, nnMatrix* out); + +/// Subtract two matrices. +void nnMatrixSub(const nnMatrix* left, const nnMatrix* right, nnMatrix* out); + +/// Adds a row vector to all rows of the matrix. +void nnMatrixAddRow(const nnMatrix* matrix, const nnMatrix* row, nnMatrix* out); + +/// Scale a matrix. +void nnMatrixScale(nnMatrix*, R scale); + +/// Transpose a matrix. +/// |in| must be different than |out|. +void nnMatrixTranspose(const nnMatrix* in, nnMatrix* out); + +/// Threshold the values of a matrix using a greater-than operator. +/// +/// out[x,y] = 1 if in[x,y] > threshold else 0 +void nnMatrixGt(const nnMatrix* in, R threshold, nnMatrix* out); + +/// Return the matrix value at the given row and column. +static inline R nnMatrixAt(const nnMatrix* matrix, int row, int col) { + assert(matrix); + return matrix->values[row * matrix->cols + col]; +} + +/// Set the matrix value at the given row and column. +static inline void nnMatrixSet(nnMatrix* matrix, int row, int col, R value) { + assert(matrix); + matrix->values[row * matrix->cols + col] = value; +} + +/// Return a pointer to the given row in the matrix. +static inline const R* nnMatrixRow(const nnMatrix* matrix, int row) { + assert(matrix); + return &matrix->values[row * matrix->cols]; +} + +/// Return a mutable pointer to the given row in the matrix. +static inline R* nnMatrixRow_mut(nnMatrix* matrix, int row) { + assert(matrix); + return &matrix->values[row * matrix->cols]; +} diff --git a/src/lib/include/neuralnet/neuralnet.h b/src/lib/include/neuralnet/neuralnet.h new file mode 100644 index 0000000..1cf1c53 --- /dev/null +++ b/src/lib/include/neuralnet/neuralnet.h @@ -0,0 +1,64 @@ +#pragma once + +#include + +typedef struct nnMatrix nnMatrix; + +typedef struct nnNeuralNetwork nnNeuralNetwork; +typedef struct nnQueryObject nnQueryObject; + +/// Neuron activation. +typedef enum nnActivation { + nnIdentity, + nnSigmoid, + nnRelu, +} nnActivation; + +/// Create a network. +nnNeuralNetwork* nnMakeNet(int num_layers, const int* layer_sizes, const nnActivation* activations); + +/// Delete the network and free its internal memory. +void nnDeleteNet(nnNeuralNetwork**); + +/// Set the network's weights. +void nnSetWeights(nnNeuralNetwork*, const R* weights); + +/// Set the network's biases. +void nnSetBiases(nnNeuralNetwork*, const R* biases); + +/// Query the network. +/// +/// |input| is a matrix of inputs, one row per input and as many columns as the +/// input's dimension. +/// +/// The query object's output matrix (see nnQueryOutputs()) is a matrix of +/// outputs, one row per output and as many columns as the output's dimension. +void nnQuery(const nnNeuralNetwork*, nnQueryObject*, const nnMatrix* input); + +/// Query the network, array version. +void nnQueryArray(const nnNeuralNetwork*, nnQueryObject*, const R* input, R* output); + +/// Create a query object. +/// +/// The query object holds all the internal memory required to query a network. +/// Query objects allocate all memory up front so that network queries can run +/// without additional memory allocation. +nnQueryObject* nnMakeQueryObject(const nnNeuralNetwork*, int num_inputs); + +/// Delete the query object and free its internal memory. +void nnDeleteQueryObject(nnQueryObject**); + +/// Return the outputs of the query. +const nnMatrix* nnNetOutputs(const nnQueryObject*); + +/// Return the network's input size. +int nnNetInputSize(const nnNeuralNetwork*); + +/// Return the network's output size. +int nnNetOutputSize(const nnNeuralNetwork*); + +/// Return the layer's input size. +int nnLayerInputSize(const nnMatrix* weights); + +/// Return the layer's output size. +int nnLayerOutputSize(const nnMatrix* weights); diff --git a/src/lib/include/neuralnet/train.h b/src/lib/include/neuralnet/train.h new file mode 100644 index 0000000..79f8e7b --- /dev/null +++ b/src/lib/include/neuralnet/train.h @@ -0,0 +1,42 @@ +#pragma once + +#include + +#include +#include + +typedef struct nnMatrix nnMatrix; + +/// Weight initialization strategy. +/// +/// Note that regardless of strategy, a layer's weights are scaled by the +/// layer's size. This is to avoid saturation when, e.g., using a sigmoid +/// activation with many inputs. Thus, a (0,1) initialization is really +/// (0,scale), for example. +typedef enum nnWeightInitStrategy { + nnWeightInit01, // (0,1) range. + nnWeightInit11, // (-1,+1) range. + nnWeightInitNormal, // Normal distribution. +} nnWeightInitStrategy; + +/// Network training parameters. +typedef struct nnTrainingParams { + R learning_rate; + int max_iterations; + uint64_t seed; + nnWeightInitStrategy weight_init; + bool debug; +} nnTrainingParams; + +/// Train the network. +/// +/// |inputs| is a matrix of inputs, one row per input and as many columns as +/// the input's dimension. +/// +/// |targets| is a matrix of targets, one row per target and as many columns as +/// the target's dimension. +void nnTrain( + nnNeuralNetwork*, + const nnMatrix* inputs, + const nnMatrix* targets, + const nnTrainingParams*); diff --git a/src/lib/include/neuralnet/types.h b/src/lib/include/neuralnet/types.h new file mode 100644 index 0000000..e8d3942 --- /dev/null +++ b/src/lib/include/neuralnet/types.h @@ -0,0 +1,3 @@ +#pragma once + +typedef double R; diff --git a/src/lib/src/activation.h b/src/lib/src/activation.h new file mode 100644 index 0000000..42ab73f --- /dev/null +++ b/src/lib/src/activation.h @@ -0,0 +1,21 @@ +#pragma once + +#include + +#include + +static inline R sigmoid(R x) { + return 1. / (1. + exp(-x)); +} + +static inline R relu(R x) { + return fmax(0, x); +} + +#define NN_MAP_ARRAY(f, in, out, size) \ + for (int i = 0; i < size; ++i) { \ + out[i] = f(in[i]); \ + } + +#define sigmoid_array(in, out, size) NN_MAP_ARRAY(sigmoid, in, out, size) +#define relu_array(in, out, size) NN_MAP_ARRAY(relu, in, out, size) diff --git a/src/lib/src/matrix.c b/src/lib/src/matrix.c new file mode 100644 index 0000000..a7a4ce6 --- /dev/null +++ b/src/lib/src/matrix.c @@ -0,0 +1,298 @@ +#include + +#include +#include +#include + +nnMatrix nnMatrixMake(int rows, int cols) { + R* values = calloc(rows * cols, sizeof(R)); + assert(values != 0); + + return (nnMatrix) { + .rows = rows, + .cols = cols, + .values = values, + }; +} + +void nnMatrixDel(nnMatrix* matrix) { + assert(matrix != 0); + + if (matrix->values != 0) { + free(matrix->values); + matrix->values = 0; + matrix->rows = 0; + matrix->cols = 0; + } +} + +void nnMatrixMove(nnMatrix* in, nnMatrix* out) { + assert(in); + assert(out); + + out->rows = in->rows; + out->cols = in->cols; + out->values = in->values; + + in->rows = 0; + in->cols = 0; + in->values = 0; +} + +void nnMatrixCopy(const nnMatrix* in, nnMatrix* out) { + assert(in); + assert(out); + assert(in->rows == out->rows); + assert(in->cols == out->cols); + + const R* in_value = in->values; + R* out_value = out->values; + + for (int i = 0; i < in->rows * in->cols; ++i) { + *out_value++ = *in_value++; + } +} + +void nnMatrixToArray(const nnMatrix* in, R* out) { + assert(in); + assert(out); + + const R* values = in->values; + for (int i = 0; i < in->rows * in->cols; ++i) { + *out++ = *values++; + } +} + +void nnMatrixRowToArray(const nnMatrix* in, int row, R* out) { + assert(in); + assert(out); + + const R* values = in->values + row * in->cols; + for (int i = 0; i < in->cols; ++i) { + *out++ = *values++; + } +} + +void nnMatrixCopyCol(const nnMatrix* in, nnMatrix* out, int col_in, int col_out) { + assert(in); + assert(out); + assert(in->rows == out->rows); + assert(col_in < in->cols); + assert(col_out < out->cols); + + for (int row = 0; row < in->rows; ++row) { + nnMatrixSet(out, row, col_out, nnMatrixAt(in, row, col_in)); + } +} + +nnMatrix nnMatrixBorrow(nnMatrix* in) { + assert(in); + + nnMatrix out; + out.rows = in->rows; + out.cols = in->cols; + out.values = in->values; + return out; +} + +nnMatrix nnMatrixBorrowRows(nnMatrix* in, int row_start, int num_rows) { + assert(in); + assert(row_start < in->rows); + assert(row_start + num_rows <= in->rows); + + nnMatrix out; + out.rows = num_rows; + out.cols = in->cols; + out.values = nnMatrixRow_mut(in, row_start); + return out; +} + +void nnMatrixInit(nnMatrix* matrix, const R* values) { + assert(matrix); + assert(values); + memcpy(matrix->values, values, matrix->rows * matrix->cols * sizeof(R)); +} + +void nnMatrixInitConstant(nnMatrix* matrix, R value) { + assert(matrix); + for (int i = 0; i < matrix->rows * matrix->cols; ++i) { + matrix->values[i] = value; + } +} + +void nnMatrixMul(const nnMatrix* left, const nnMatrix* right, nnMatrix* out) { + assert(left != 0); + assert(right != 0); + assert(out != 0); + assert(out != left); + assert(out != right); + assert(left->cols == right->rows); + assert(out->rows == left->rows); + assert(out->cols == right->cols); + + R* out_value = out->values; + + for (int i = 0; i < left->rows; ++i) { + const R* left_row = &left->values[i * left->cols]; + + for (int j = 0; j < right->cols; ++j) { + const R* right_col = &right->values[j]; + *out_value = 0; + + // Vector dot product. + for (int k = 0; k < left->cols; ++k) { + *out_value += left_row[k] * right_col[0]; + right_col += right->cols; // Next row in the column. + } + + out_value++; + } + } +} + +void nnMatrixMulAdd(const nnMatrix* left, const nnMatrix* right, R scale, nnMatrix* out) { + assert(left); + assert(right); + assert(out); + assert(left->rows == right->rows); + assert(left->cols == right->cols); + assert(left->rows == out->rows); + assert(left->cols == out->cols); + + const R* left_value = left->values; + const R* right_value = right->values; + R* out_value = out->values; + + for (int i = 0; i < left->rows * left->cols; ++i) { + *out_value++ = *left_value++ + *right_value++ * scale; + } +} + +void nnMatrixMulSub(const nnMatrix* left, const nnMatrix* right, R scale, nnMatrix* out) { + assert(left); + assert(right); + assert(out); + assert(left->rows == right->rows); + assert(left->cols == right->cols); + assert(left->rows == out->rows); + assert(left->cols == out->cols); + + const R* left_value = left->values; + const R* right_value = right->values; + R* out_value = out->values; + + for (int i = 0; i < left->rows * left->cols; ++i) { + *out_value++ = *left_value++ - *right_value++ * scale; + } +} + +void nnMatrixMulPairs(const nnMatrix* left, const nnMatrix* right, nnMatrix* out) { + assert(left != 0); + assert(right != 0); + assert(out != 0); + assert(left->rows == right->rows); + assert(left->cols == right->cols); + assert(left->rows == out->rows); + assert(left->cols == out->cols); + + R* left_value = left->values; + R* right_value = right->values; + R* out_value = out->values; + + for (int i = 0; i < left->rows * left->cols; ++i) { + *out_value++ = *left_value++ * *right_value++; + } +} + +void nnMatrixAdd(const nnMatrix* left, const nnMatrix* right, nnMatrix* out) { + assert(left); + assert(right); + assert(out); + assert(left->rows == right->rows); + assert(left->cols == right->cols); + assert(left->rows == out->rows); + assert(left->cols == out->cols); + + const R* left_value = left->values; + const R* right_value = right->values; + R* out_value = out->values; + + for (int i = 0; i < left->rows * left->cols; ++i) { + *out_value++ = *left_value++ + *right_value++; + } +} + +void nnMatrixSub(const nnMatrix* left, const nnMatrix* right, nnMatrix* out) { + assert(left); + assert(right); + assert(out); + assert(left->rows == right->rows); + assert(left->cols == right->cols); + assert(left->rows == out->rows); + assert(left->cols == out->cols); + + const R* left_value = left->values; + const R* right_value = right->values; + R* out_value = out->values; + + for (int i = 0; i < left->rows * left->cols; ++i) { + *out_value++ = *left_value++ - *right_value++; + } +} + +void nnMatrixAddRow(const nnMatrix* matrix, const nnMatrix* row, nnMatrix* out) { + assert(matrix); + assert(row); + assert(out); + assert(row->rows == 1); + assert(matrix->cols == row->cols); + assert(matrix->rows == out->rows); + assert(matrix->cols == out->cols); + + const R* matrix_value = matrix->values; + R* out_value = out->values; + + for (int i = 0; i < matrix->rows; ++i) { + const R* row_value = row->values; + for (int j = 0; j < row->cols; ++j) { + *out_value++ = *matrix_value++ + *row_value++; + } + } +} + +void nnMatrixScale(nnMatrix* matrix, R scale) { + assert(matrix); + + R* value = matrix->values; + for (int i = 0; i < matrix->rows * matrix->cols; ++i) { + *value++ *= scale; + } +} + +void nnMatrixTranspose(const nnMatrix* in, nnMatrix* out) { + assert(in); + assert(out); + assert(in != out); + assert(in->rows == out->cols); + assert(in->cols == out->rows); + + for (int i = 0; i < in->rows; ++i) { + for (int j = 0; j < in->cols; ++j) { + nnMatrixSet(out, j, i, nnMatrixAt(in, i, j)); + } + } +} + +void nnMatrixGt(const nnMatrix* in, R threshold, nnMatrix* out) { + assert(in); + assert(out); + assert(in->rows == out->rows); + assert(in->cols == out->cols); + + const R* in_value = in->values; + R* out_value = out->values; + + for (int i = 0; i < in->rows * in->cols; ++i) { + *out_value++ = (*in_value++) > threshold ? 1 : 0; + } +} diff --git a/src/lib/src/neuralnet.c b/src/lib/src/neuralnet.c new file mode 100644 index 0000000..cac611a --- /dev/null +++ b/src/lib/src/neuralnet.c @@ -0,0 +1,228 @@ +#include + +#include +#include "activation.h" +#include "neuralnet_impl.h" + +#include +#include + +nnNeuralNetwork* nnMakeNet(int num_layers, const int* layer_sizes, const nnActivation* activations) { + assert(num_layers > 0); + assert(layer_sizes); + assert(activations); + + nnNeuralNetwork* net = calloc(1, sizeof(nnNeuralNetwork)); + if (net == 0) { + return 0; + } + + net->num_layers = num_layers; + + net->weights = calloc(num_layers, sizeof(nnMatrix)); + net->biases = calloc(num_layers, sizeof(nnMatrix)); + net->activations = calloc(num_layers, sizeof(nnActivation)); + if ( (net->weights == 0) || (net->biases == 0) || (net->activations == 0) ) { + nnDeleteNet(&net); + return 0; + } + + for (int l = 0; l < num_layers; ++l) { + // layer_sizes = { input layer size, first hidden layer size, ...} + const int layer_input_size = layer_sizes[l]; + const int layer_output_size = layer_sizes[l+1]; + + // We store the transpose of the weight matrix as written in textbooks. + // Our vectors are row vectors and the matrices row-major. + const int rows = layer_input_size; + const int cols = layer_output_size; + + net->weights[l] = nnMatrixMake(rows, cols); + net->biases[l] = nnMatrixMake(1, cols); + net->activations[l] = activations[l]; + } + + return net; +} + +void nnDeleteNet(nnNeuralNetwork** net) { + if ( (!net) || (!(*net)) ) { + return; + } + if ((*net)->weights != 0) { + for (int l = 0; l < (*net)->num_layers; ++l) { + nnMatrixDel(&(*net)->weights[l]); + } + free((*net)->weights); + (*net)->weights = 0; + } + if ((*net)->biases != 0) { + for (int l = 0; l < (*net)->num_layers; ++l) { + nnMatrixDel(&(*net)->biases[l]); + } + free((*net)->biases); + (*net)->biases = 0; + } + if ((*net)->activations) { + free((*net)->activations); + (*net)->activations = 0; + } + free(*net); + *net = 0; +} + +void nnSetWeights(nnNeuralNetwork* net, const R* weights) { + assert(net); + assert(weights); + + for (int l = 0; l < net->num_layers; ++l) { + nnMatrix* layer_weights = &net->weights[l]; + R* layer_values = layer_weights->values; + + for (int j = 0; j < layer_weights->rows * layer_weights->cols; ++j) { + *layer_values++ = *weights++; + } + } +} + +void nnSetBiases(nnNeuralNetwork* net, const R* biases) { + assert(net); + assert(biases); + + for (int l = 0; l < net->num_layers; ++l) { + nnMatrix* layer_biases = &net->biases[l]; + R* layer_values = layer_biases->values; + + for (int j = 0; j < layer_biases->rows * layer_biases->cols; ++j) { + *layer_values++ = *biases++; + } + } +} + +void nnQuery(const nnNeuralNetwork* net, nnQueryObject* query, const nnMatrix* input) { + assert(net); + assert(query); + assert(input); + assert(net->num_layers == query->num_layers); + assert(input->rows <= query->network_outputs->rows); + assert(input->cols == nnNetInputSize(net)); + + for (int i = 0; i < input->rows; ++i) { + // Not mutating the input, but we need the cast to borrow. + nnMatrix input_vector = nnMatrixBorrowRows((nnMatrix*)input, i, 1); + + for (int l = 0; l < net->num_layers; ++l) { + const nnMatrix* layer_weights = &net->weights[l]; + const nnMatrix* layer_biases = &net->biases[l]; + // Y^T = (W*X)^T = X^T*W^T + // + // TODO: If we had a row-row matrix multiplication, we could compute: + // Y^T = W ** X^T + // The row-row multiplication could be more cache-friendly. We just need + // to store W as is, without transposing. + // We could also rewrite the original Mul function to go row x row, + // decomposing the multiplication. Preserving the original meaning of Mul + // makes everything clearer. + nnMatrix output_vector = nnMatrixBorrowRows(&query->layer_outputs[l], i, 1); + nnMatrixMul(&input_vector, layer_weights, &output_vector); + nnMatrixAddRow(&output_vector, layer_biases, &output_vector); + + switch (net->activations[l]) { + case nnIdentity: + break; // Nothing to do for the identity function. + case nnSigmoid: + sigmoid_array(output_vector.values, output_vector.values, output_vector.cols); + break; + case nnRelu: + relu_array(output_vector.values, output_vector.values, output_vector.cols); + break; + default: + assert(0); + } + + input_vector = output_vector; // Borrow. + } + } +} + +void nnQueryArray(const nnNeuralNetwork* net, nnQueryObject* query, const R* input, R* output) { + assert(net); + assert(query); + assert(input); + assert(output); + assert(net->num_layers > 0); + + nnMatrix input_vector = nnMatrixMake(net->weights[0].cols, 1); + nnMatrixInit(&input_vector, input); + nnQuery(net, query, &input_vector); + nnMatrixRowToArray(query->network_outputs, 0, output); +} + +nnQueryObject* nnMakeQueryObject(const nnNeuralNetwork* net, int num_inputs) { + assert(net); + assert(num_inputs > 0); + assert(net->num_layers > 0); + + nnQueryObject* query = calloc(1, sizeof(nnQueryObject)); + if (!query) { + return 0; + } + + query->num_layers = net->num_layers; + + // Allocate the intermediate layer output matrices. + query->layer_outputs = calloc(net->num_layers, sizeof(nnMatrix)); + if (!query->layer_outputs) { + free(query); + return 0; + } + for (int l = 0; l < net->num_layers; ++l) { + const nnMatrix* layer_weights = &net->weights[l]; + const int layer_output_size = nnLayerOutputSize(layer_weights); + query->layer_outputs[l] = nnMatrixMake(num_inputs, layer_output_size); + } + query->network_outputs = &query->layer_outputs[net->num_layers - 1]; + + return query; +} + +void nnDeleteQueryObject(nnQueryObject** query) { + if ( (!query) || (!(*query)) ) { + return; + } + if ((*query)->layer_outputs != 0) { + for (int l = 0; l < (*query)->num_layers; ++l) { + nnMatrixDel(&(*query)->layer_outputs[l]); + } + } + free((*query)->layer_outputs); + free(*query); + *query = 0; +} + +const nnMatrix* nnNetOutputs(const nnQueryObject* query) { + assert(query); + return query->network_outputs; +} + +int nnNetInputSize(const nnNeuralNetwork* net) { + assert(net); + assert(net->num_layers > 0); + return net->weights[0].rows; +} + +int nnNetOutputSize(const nnNeuralNetwork* net) { + assert(net); + assert(net->num_layers > 0); + return net->weights[net->num_layers - 1].cols; +} + +int nnLayerInputSize(const nnMatrix* weights) { + assert(weights); + return weights->rows; +} + +int nnLayerOutputSize(const nnMatrix* weights) { + assert(weights); + return weights->cols; +} diff --git a/src/lib/src/neuralnet_impl.h b/src/lib/src/neuralnet_impl.h new file mode 100644 index 0000000..26107b5 --- /dev/null +++ b/src/lib/src/neuralnet_impl.h @@ -0,0 +1,36 @@ +#pragma once + +#include + +/// Neural network object. +/// +/// We store the transposes of the weight matrices so that we can do forward +/// passes with a minimal amount of work. That is, if in paper we write: +/// +/// [w11 w21] +/// [w12 w22] +/// +/// then the weight matrix in memory is stored as the following array: +/// +/// w11 w12 w21 w22 +typedef struct nnNeuralNetwork { + int num_layers; // Number of non-input layers (hidden + output). + nnMatrix* weights; // One matrix per non-input layer. + nnMatrix* biases; // One vector per non-input layer. + nnActivation* activations; // One per non-input layer. +} nnNeuralNetwork; + +/// A query object that holds all the memory necessary to query a network. +/// +/// |layer_outputs| is an array of matrices of intermediate layer outputs. There +/// is one matrix per intermediate layer. Each matrix holds the layer's output, +/// with one row per input, and as many columns as the layer's output size (the +/// output vector is transposed.) +/// +/// |network_outputs| points to the last output matrix in |layer_outputs| for +/// convenience. +typedef struct nnQueryObject { + int num_layers; + nnMatrix* layer_outputs; // Output matrices, one output per layer. + nnMatrix* network_outputs; // Points to the last output matrix. +} nnTrainingQueryObject; diff --git a/src/lib/src/train.c b/src/lib/src/train.c new file mode 100644 index 0000000..027de66 --- /dev/null +++ b/src/lib/src/train.c @@ -0,0 +1,346 @@ +#include + +#include +#include "neuralnet_impl.h" + +#include +#include + +#include +#include +#include + +#include +#define LOGD printf + +// If debug mode is requested, we will show progress every this many iterations. +static const int PROGRESS_THRESHOLD = 5; // % + +/// Computes the total MSE from the output error matrix. +R ComputeMSE(const nnMatrix* errors) { + R sum_sq = 0; + const int N = errors->rows * errors->cols; + const R* value = errors->values; + for (int i = 0; i < N; ++i) { + sum_sq += *value * *value; + value++; + } + return sum_sq / (R)N; +} + +/// Holds the bits required to compute a sigmoid gradient. +typedef struct nnSigmoidGradientElements { + nnMatrix ones; // A vector of just ones, same size as the layer. +} nnSigmoidGradientElements; + +/// Holds the various elements required to compute gradients. These depend on +/// what activation function are used, so they'll potentially be different for +/// each layer. A data type is defined for these because we allocate all the +/// required memory up front before entering the training loop. +typedef struct nnGradientElements { + nnActivation type; + // Gradient vector, same size as the layer. + // This will contain the gradient expression except for the output value of + // the previous layer. + nnMatrix gradient; + union { + nnSigmoidGradientElements sigmoid; + }; +} nnGradientElements; + +// Initialize the network's weights randomly and set their biases to 0. +void nnInitNet(nnNeuralNetwork* net, uint64_t seed, const nnWeightInitStrategy strategy) { + assert(net); + + mt19937_64 rng = mt19937_64_make(); + mt19937_64_init(&rng, seed); + + for (int l = 0; l < net->num_layers; ++l) { + nnMatrix* weights = &net->weights[l]; + nnMatrix* biases = &net->biases[l]; + + const R layer_size = (R)nnLayerInputSize(weights); + const R scale = 1. / layer_size; + const R stdev = 1. / sqrt((R)layer_size); + const R sigma = stdev * stdev; + + R* value = weights->values; + for (int k = 0; k < weights->rows * weights->cols; ++k) { + switch (strategy) { + case nnWeightInit01: { + const R x01 = mt19937_64_gen_real3(&rng); // (0, +1) interval. + *value++ = scale * x01; + break; + } + case nnWeightInit11: { + const R x11 = mt19937_64_gen_real4(&rng); // (-1, +1) interval. + *value++ = scale * x11; + break; + } + case nnWeightInitNormal: + // Using initialization with a normal distribution of standard + // deviation 1 / sqrt(num_layer_weights) to prevent saturation when + // multiplying inputs. + const R u01 = mt19937_64_gen_real3(&rng); // (0, +1) interval. + const R v01 = mt19937_64_gen_real3(&rng); // (0, +1) interval. + R z0, z1; + normal2(u01, v01, &z0, &z1); + z0 = normal_transform(z0, /*mu=*/0, sigma); + z1 = normal_transform(z1, /*mu=*/0, sigma); + *value++ = z0; + if (k < weights->rows * weights->cols - 1) { + *value++ = z1; + ++k; + } + break; + default: + assert(false); + } + } + + // Initialize biases. + // 0 is used so that functions originally go through the origin. + value = biases->values; + for (int k = 0; k < biases->rows * biases->cols; ++k, ++value) { + *value = 0; + } + } +} + +// |inputs| has one row vector per sample. +// |targets| has one row vector per sample. +// +// For now, each iteration trains with one sample (row) at a time. +void nnTrain( + nnNeuralNetwork* net, + const nnMatrix* inputs, + const nnMatrix* targets, + const nnTrainingParams* params) { + assert(net); + assert(inputs); + assert(targets); + assert(params); + assert(nnNetOutputSize(net) == targets->cols); + assert(net->num_layers > 0); + + // Allocate error vectors to hold the backpropagated error values. + // For now, these are one row vector per layer, meaning that we will train + // with one sample at a time. + nnMatrix* errors = calloc(net->num_layers, sizeof(nnMatrix)); + + // Allocate the weight transpose matrices up front for backpropagation. + nnMatrix* weights_T = calloc(net->num_layers, sizeof(nnMatrix)); + + // Allocate the weight delta matrices. + nnMatrix* weight_deltas = calloc(net->num_layers, sizeof(nnMatrix)); + + // Allocate the data structures required to compute gradients. + // This depends on each layer's activation type. + nnGradientElements* gradient_elems = calloc(net->num_layers, sizeof(nnGradientElements)); + + // Allocate the output transpose vectors for weight delta calculation. + // This is one column vector per layer. + nnMatrix* outputs_T = calloc(net->num_layers, sizeof(nnMatrix)); + + assert(errors != 0); + assert(weights_T != 0); + assert(weight_deltas != 0); + assert(gradient_elems); + assert(outputs_T); + + for (int l = 0; l < net->num_layers; ++l) { + const nnMatrix* layer_weights = &net->weights[l]; + const int layer_output_size = net->weights[l].cols; + const nnActivation activation = net->activations[l]; + + errors[l] = nnMatrixMake(1, layer_weights->cols); + + weights_T[l] = nnMatrixMake(layer_weights->cols, layer_weights->rows); + nnMatrixTranspose(layer_weights, &weights_T[l]); + + weight_deltas[l] = nnMatrixMake(layer_weights->rows, layer_weights->cols); + + outputs_T[l] = nnMatrixMake(layer_output_size, 1); + + // Allocate the gradient elements and vectors for weight delta calculation. + nnGradientElements* elems = &gradient_elems[l]; + elems->type = activation; + switch (activation) { + case nnIdentity: + break; // Gradient vector will be borrowed, no need to allocate. + + case nnSigmoid: + elems->gradient = nnMatrixMake(1, layer_output_size); + // Allocate the 1s vectors. + elems->sigmoid.ones = nnMatrixMake(1, layer_output_size); + nnMatrixInitConstant(&elems->sigmoid.ones, 1); + break; + + case nnRelu: + elems->gradient = nnMatrixMake(1, layer_output_size); + break; + } + } + + // Construct the query object with a size of 1 since we are training with one + // sample at a time. + nnQueryObject* query = nnMakeQueryObject(net, 1); + + // Network outputs are given by the query object. Every network query updates + // the outputs. + const nnMatrix* const training_outputs = query->network_outputs; + + // A vector to store the training input transposed. + nnMatrix training_inputs_T = nnMatrixMake(inputs->cols, 1); + + // If debug mode is requested, we will show progress every Nth iteration. + const int progress_frame = + (params->max_iterations < PROGRESS_THRESHOLD) + ? 1 + : (params->max_iterations * PROGRESS_THRESHOLD / 100); + + // --- TRAIN + + nnInitNet(net, params->seed, params->weight_init); + + for (int iter = 0; iter < params->max_iterations; ++iter) { + + // For now, we train with one sample at a time. + for (int sample = 0; sample < inputs->rows; ++sample) { + // Slice the input and target matrices with the batch size. + // We are not mutating the inputs, but we need the cast to borrow. + nnMatrix training_inputs = nnMatrixBorrowRows((nnMatrix*)inputs, sample, 1); + nnMatrix training_targets = nnMatrixBorrowRows((nnMatrix*)targets, sample, 1); + + // Will need the input transposed for backpropagation. + // Assuming one training input per iteration for now. + nnMatrixTranspose(&training_inputs, &training_inputs_T); + + // Run a forward pass and compute the output layer error. + // We don't square the error here; instead, we just compute t-o, which is + // part of the derivative, -2(t-o). Also, we compute o-t instead to + // remove that outer negative sign. + nnQuery(net, query, &training_inputs); + //nnMatrixSub(&training_targets, training_outputs, &errors[net->num_layers - 1]); + nnMatrixSub(training_outputs, &training_targets, &errors[net->num_layers - 1]); + + // Update outputs_T, which we need during weight updates. + for (int l = 0; l < net->num_layers; ++l) { + nnMatrixTranspose(&query->layer_outputs[l], &outputs_T[l]); + } + + // Update weights and biases for each internal layer, backpropagating + // errors along the way. + for (int l = net->num_layers - 1; l >= 0; --l) { + const nnMatrix* layer_output = &query->layer_outputs[l]; + nnMatrix* layer_weights = &net->weights[l]; + nnMatrix* layer_biases = &net->biases[l]; + nnGradientElements* elems = &gradient_elems[l]; + nnMatrix* gradient = &elems->gradient; + const nnActivation activation = net->activations[l]; + + // Compute the gradient (the part of the expression that does not + // contain the output of the previous layer). + // + // Identity: G = error_k + // Sigmoid: G = error_k * output_k * (1 - output_k). + // Relu: G = error_k * (output_k > 0 ? 1 : 0) + switch (activation) { + case nnIdentity: + // TODO: Just copy the pointer? + *gradient = nnMatrixBorrow(&errors[l]); + break; + case nnSigmoid: + nnMatrixSub(&elems->sigmoid.ones, layer_output, gradient); + nnMatrixMulPairs(layer_output, gradient, gradient); + nnMatrixMulPairs(&errors[l], gradient, gradient); + break; + case nnRelu: + nnMatrixGt(layer_output, 0, gradient); + nnMatrixMulPairs(&errors[l], gradient, gradient); + break; + } + + // Outer product to compute the weight deltas. + const nnMatrix* output_T = (l == 0) ? &training_inputs_T : &outputs_T[l-1]; + nnMatrixMul(output_T, gradient, &weight_deltas[l]); + + // Backpropagate the error before updating weights. + if (l > 0) { + nnMatrixMul(gradient, &weights_T[l], &errors[l-1]); + } + + // Update weights. + nnMatrixScale(&weight_deltas[l], params->learning_rate); + // The gradient has a negative sign from -(t - o), but we have computed + // e = o - t instead, so we can subtract directly. + //nnMatrixAdd(layer_weights, &weight_deltas[l], layer_weights); + nnMatrixSub(layer_weights, &weight_deltas[l], layer_weights); + + // Update weight transpose matrix for the next training iteration. + nnMatrixTranspose(layer_weights, &weights_T[l]); + + // Update biases. + // This is the same formula as for weights, except that the o_j term is + // just 1. We can simply re-use the gradient that we have already + // computed for the weight update. + //nnMatrixMulAdd(layer_biases, gradient, params->learning_rate, layer_biases); + nnMatrixMulSub(layer_biases, gradient, params->learning_rate, layer_biases); + } + + // TODO: Add this under a verbose debugging mode. + // if (params->debug) { + // LOGD("Iter: %d, Sample: %d, Error: %f\n", iter, sample, ComputeMSE(&errors[net->num_layers - 1])); + // LOGD("TGT: "); + // for (int i = 0; i < training_targets.cols; ++i) { + // printf("%.3f ", training_targets.values[i]); + // } + // printf("\n"); + // LOGD("OUT: "); + // for (int i = 0; i < training_outputs->cols; ++i) { + // printf("%.3f ", training_outputs->values[i]); + // } + // printf("\n"); + // } + } + + if (params->debug && ((iter % progress_frame) == 0)) { + LOGD("Iter: %d/%d, Error: %f\n", + iter, params->max_iterations, ComputeMSE(&errors[net->num_layers - 1])); + } + } + + // Print the final error. + if (params->debug) { + LOGD("Iter: %d/%d, Error: %f\n", + params->max_iterations, params->max_iterations, ComputeMSE(&errors[net->num_layers - 1])); + } + + for (int l = 0; l < net->num_layers; ++l) { + nnMatrixDel(&errors[l]); + nnMatrixDel(&outputs_T[l]); + nnMatrixDel(&weights_T[l]); + nnMatrixDel(&weight_deltas[l]); + + nnGradientElements* elems = &gradient_elems[l]; + switch (elems->type) { + case nnIdentity: + break; // Gradient vector is borrowed, no need to deallocate. + + case nnSigmoid: + nnMatrixDel(&elems->gradient); + nnMatrixDel(&elems->sigmoid.ones); + break; + + case nnRelu: + nnMatrixDel(&elems->gradient); + break; + } + } + nnMatrixDel(&training_inputs_T); + free(errors); + free(outputs_T); + free(weights_T); + free(weight_deltas); + free(gradient_elems); +} diff --git a/src/lib/test/matrix_test.c b/src/lib/test/matrix_test.c new file mode 100644 index 0000000..8191c97 --- /dev/null +++ b/src/lib/test/matrix_test.c @@ -0,0 +1,350 @@ +#include + +#include "test.h" +#include "test_util.h" + +#include +#include + +// static void PrintMatrix(const nnMatrix* matrix) { +// assert(matrix); + +// for (int i = 0; i < matrix->rows; ++i) { +// for (int j = 0; j < matrix->cols; ++j) { +// printf("%f ", nnMatrixAt(matrix, i, j)); +// } +// printf("\n"); +// } +// } + +TEST_CASE(nnMatrixMake_1x1) { + nnMatrix A = nnMatrixMake(1, 1); + TEST_EQUAL(A.rows, 1); + TEST_EQUAL(A.cols, 1); +} + +TEST_CASE(nnMatrixMake_3x1) { + nnMatrix A = nnMatrixMake(3, 1); + TEST_EQUAL(A.rows, 3); + TEST_EQUAL(A.cols, 1); +} + +TEST_CASE(nnMatrixInit_3x1) { + nnMatrix A = nnMatrixMake(3, 1); + nnMatrixInit(&A, (R[]) { 1, 2, 3 }); + TEST_EQUAL(A.values[0], 1); + TEST_EQUAL(A.values[1], 2); + TEST_EQUAL(A.values[2], 3); +} + +TEST_CASE(nnMatrixCopyCol_test) { + nnMatrix A = nnMatrixMake(3, 2); + nnMatrix B = nnMatrixMake(3, 1); + + nnMatrixInit(&A, (R[]) { + 1, 2, + 3, 4, + 5, 6, + }); + + nnMatrixCopyCol(&A, &B, 1, 0); + + TEST_EQUAL(nnMatrixAt(&B, 0, 0), 2); + TEST_EQUAL(nnMatrixAt(&B, 1, 0), 4); + TEST_EQUAL(nnMatrixAt(&B, 2, 0), 6); + + nnMatrixDel(&A); + nnMatrixDel(&B); +} + +TEST_CASE(nnMatrixMul_square_3x3) { + nnMatrix A = nnMatrixMake(3, 3); + nnMatrix B = nnMatrixMake(3, 3); + nnMatrix O = nnMatrixMake(3, 3); + + nnMatrixInit(&A, (const R[]){ + 1, 2, 3, + 4, 5, 6, + 7, 8, 9, + }); + nnMatrixInit(&B, (const R[]){ + 2, 4, 3, + 6, 8, 5, + 1, 7, 9, + }); + nnMatrixMul(&A, &B, &O); + + const R expected[3][3] = { + { 17, 41, 40 }, + { 44, 98, 91 }, + { 71, 155, 142 }, + }; + for (int i = 0; i < O.rows; ++i) { + for (int j = 0; j < O.cols; ++j) { + TEST_TRUE(double_eq(nnMatrixAt(&O, i, j), expected[i][j], EPS)); + } + } + + nnMatrixDel(&A); + nnMatrixDel(&B); + nnMatrixDel(&O); +} + +TEST_CASE(nnMatrixMul_non_square_2x3_3x1) { + nnMatrix A = nnMatrixMake(2, 3); + nnMatrix B = nnMatrixMake(3, 1); + nnMatrix O = nnMatrixMake(2, 1); + + nnMatrixInit(&A, (const R[]){ + 1, 2, 3, + 4, 5, 6, + }); + nnMatrixInit(&B, (const R[]){ + 2, + 6, + 1, + }); + nnMatrixMul(&A, &B, &O); + + const R expected[2][1] = { + { 17 }, + { 44 }, + }; + for (int i = 0; i < O.rows; ++i) { + for (int j = 0; j < O.cols; ++j) { + TEST_TRUE(double_eq(nnMatrixAt(&O, i, j), expected[i][j], EPS)); + } + } + + nnMatrixDel(&A); + nnMatrixDel(&B); + nnMatrixDel(&O); +} + +TEST_CASE(nnMatrixMulAdd_test) { + nnMatrix A = nnMatrixMake(2, 3); + nnMatrix B = nnMatrixMake(2, 3); + nnMatrix O = nnMatrixMake(2, 3); + const R scale = 2; + + nnMatrixInit(&A, (const R[]){ + 1, 2, 3, + 4, 5, 6, + }); + nnMatrixInit(&B, (const R[]){ + 2, 3, 1, + 7, 4, 3 + }); + nnMatrixMulAdd(&A, &B, scale, &O); // O = A + B * scale + + const R expected[2][3] = { + { 5, 8, 5 }, + { 18, 13, 12 }, + }; + for (int i = 0; i < O.rows; ++i) { + for (int j = 0; j < O.cols; ++j) { + TEST_TRUE(double_eq(nnMatrixAt(&O, i, j), expected[i][j], EPS)); + } + } + + nnMatrixDel(&A); + nnMatrixDel(&B); + nnMatrixDel(&O); +} + +TEST_CASE(nnMatrixMulSub_test) { + nnMatrix A = nnMatrixMake(2, 3); + nnMatrix B = nnMatrixMake(2, 3); + nnMatrix O = nnMatrixMake(2, 3); + const R scale = 2; + + nnMatrixInit(&A, (const R[]){ + 1, 2, 3, + 4, 5, 6, + }); + nnMatrixInit(&B, (const R[]){ + 2, 3, 1, + 7, 4, 3 + }); + nnMatrixMulSub(&A, &B, scale, &O); // O = A - B * scale + + const R expected[2][3] = { + { -3, -4, 1 }, + { -10, -3, 0 }, + }; + for (int i = 0; i < O.rows; ++i) { + for (int j = 0; j < O.cols; ++j) { + TEST_TRUE(double_eq(nnMatrixAt(&O, i, j), expected[i][j], EPS)); + } + } + + nnMatrixDel(&A); + nnMatrixDel(&B); + nnMatrixDel(&O); +} + +TEST_CASE(nnMatrixMulPairs_2x3) { + nnMatrix A = nnMatrixMake(2, 3); + nnMatrix B = nnMatrixMake(2, 3); + nnMatrix O = nnMatrixMake(2, 3); + + nnMatrixInit(&A, (const R[]){ + 1, 2, 3, + 4, 5, 6, + }); + nnMatrixInit(&B, (const R[]){ + 2, 3, 1, + 7, 4, 3 + }); + nnMatrixMulPairs(&A, &B, &O); + + const R expected[2][3] = { + { 2, 6, 3 }, + { 28, 20, 18 }, + }; + for (int i = 0; i < O.rows; ++i) { + for (int j = 0; j < O.cols; ++j) { + TEST_TRUE(double_eq(nnMatrixAt(&O, i, j), expected[i][j], EPS)); + } + } + + nnMatrixDel(&A); + nnMatrixDel(&B); + nnMatrixDel(&O); +} + +TEST_CASE(nnMatrixAdd_square_2x2) { + nnMatrix A = nnMatrixMake(2, 2); + nnMatrix B = nnMatrixMake(2, 2); + nnMatrix C = nnMatrixMake(2, 2); + + nnMatrixInit(&A, (R[]) { + 1, 2, + 3, 4, + }); + nnMatrixInit(&B, (R[]) { + 2, 1, + 5, 3, + }); + + nnMatrixAdd(&A, &B, &C); + + TEST_TRUE(double_eq(nnMatrixAt(&C, 0, 0), 3, EPS)); + TEST_TRUE(double_eq(nnMatrixAt(&C, 0, 1), 3, EPS)); + TEST_TRUE(double_eq(nnMatrixAt(&C, 1, 0), 8, EPS)); + TEST_TRUE(double_eq(nnMatrixAt(&C, 1, 1), 7, EPS)); + + nnMatrixDel(&A); + nnMatrixDel(&B); + nnMatrixDel(&C); +} + +TEST_CASE(nnMatrixSub_square_2x2) { + nnMatrix A = nnMatrixMake(2, 2); + nnMatrix B = nnMatrixMake(2, 2); + nnMatrix C = nnMatrixMake(2, 2); + + nnMatrixInit(&A, (R[]) { + 1, 2, + 3, 4, + }); + nnMatrixInit(&B, (R[]) { + 2, 1, + 5, 3, + }); + + nnMatrixSub(&A, &B, &C); + + TEST_TRUE(double_eq(nnMatrixAt(&C, 0, 0), -1, EPS)); + TEST_TRUE(double_eq(nnMatrixAt(&C, 0, 1), +1, EPS)); + TEST_TRUE(double_eq(nnMatrixAt(&C, 1, 0), -2, EPS)); + TEST_TRUE(double_eq(nnMatrixAt(&C, 1, 1), +1, EPS)); + + nnMatrixDel(&A); + nnMatrixDel(&B); + nnMatrixDel(&C); +} + +TEST_CASE(nnMatrixAddRow_test) { + nnMatrix A = nnMatrixMake(2, 3); + nnMatrix B = nnMatrixMake(1, 3); + nnMatrix C = nnMatrixMake(2, 3); + + nnMatrixInit(&A, (R[]) { + 1, 2, 3, + 4, 5, 6, + }); + nnMatrixInit(&B, (R[]) { + 2, 1, 3, + }); + + nnMatrixAddRow(&A, &B, &C); + + TEST_TRUE(double_eq(nnMatrixAt(&C, 0, 0), 3, EPS)); + TEST_TRUE(double_eq(nnMatrixAt(&C, 0, 1), 3, EPS)); + TEST_TRUE(double_eq(nnMatrixAt(&C, 0, 2), 6, EPS)); + TEST_TRUE(double_eq(nnMatrixAt(&C, 1, 0), 6, EPS)); + TEST_TRUE(double_eq(nnMatrixAt(&C, 1, 1), 6, EPS)); + TEST_TRUE(double_eq(nnMatrixAt(&C, 1, 2), 9, EPS)); + + nnMatrixDel(&A); + nnMatrixDel(&B); + nnMatrixDel(&C); +} + +TEST_CASE(nnMatrixTranspose_square_2x2) { + nnMatrix A = nnMatrixMake(2, 2); + nnMatrix B = nnMatrixMake(2, 2); + + nnMatrixInit(&A, (R[]) { + 1, 2, + 3, 4 + }); + + nnMatrixTranspose(&A, &B); + TEST_TRUE(double_eq(nnMatrixAt(&B, 0, 0), 1, EPS)); + TEST_TRUE(double_eq(nnMatrixAt(&B, 0, 1), 3, EPS)); + TEST_TRUE(double_eq(nnMatrixAt(&B, 1, 0), 2, EPS)); + TEST_TRUE(double_eq(nnMatrixAt(&B, 1, 1), 4, EPS)); + + nnMatrixDel(&A); + nnMatrixDel(&B); +} + +TEST_CASE(nnMatrixTranspose_non_square_2x1) { + nnMatrix A = nnMatrixMake(2, 1); + nnMatrix B = nnMatrixMake(1, 2); + + nnMatrixInit(&A, (R[]) { + 1, + 3, + }); + + nnMatrixTranspose(&A, &B); + TEST_TRUE(double_eq(nnMatrixAt(&B, 0, 0), 1, EPS)); + TEST_TRUE(double_eq(nnMatrixAt(&B, 0, 1), 3, EPS)); + + nnMatrixDel(&A); + nnMatrixDel(&B); +} + +TEST_CASE(nnMatrixGt_test) { + nnMatrix A = nnMatrixMake(2, 3); + nnMatrix B = nnMatrixMake(2, 3); + + nnMatrixInit(&A, (R[]) { + -3, 2, 0, + 4, -1, 5 + }); + + nnMatrixGt(&A, 0, &B); + TEST_TRUE(double_eq(nnMatrixAt(&B, 0, 0), 0, EPS)); + TEST_TRUE(double_eq(nnMatrixAt(&B, 0, 1), 1, EPS)); + TEST_TRUE(double_eq(nnMatrixAt(&B, 0, 2), 0, EPS)); + TEST_TRUE(double_eq(nnMatrixAt(&B, 1, 0), 1, EPS)); + TEST_TRUE(double_eq(nnMatrixAt(&B, 1, 1), 0, EPS)); + TEST_TRUE(double_eq(nnMatrixAt(&B, 1, 2), 1, EPS)); + + nnMatrixDel(&A); + nnMatrixDel(&B); +} diff --git a/src/lib/test/neuralnet_test.c b/src/lib/test/neuralnet_test.c new file mode 100644 index 0000000..14d9438 --- /dev/null +++ b/src/lib/test/neuralnet_test.c @@ -0,0 +1,92 @@ +#include + +#include +#include "activation.h" +#include "neuralnet_impl.h" + +#include "test.h" +#include "test_util.h" + +#include + +TEST_CASE(neuralnet_perceptron_test) { + const int num_layers = 1; + const int layer_sizes[] = { 1, 1 }; + const nnActivation layer_activations[] = { nnSigmoid }; + const R weights[] = { 0.3 }; + + nnNeuralNetwork* net = nnMakeNet(num_layers, layer_sizes, layer_activations); + assert(net); + nnSetWeights(net, weights); + + nnQueryObject* query = nnMakeQueryObject(net, /*num_inputs=*/1); + + const R input[] = { 0.9 }; + R output[1]; + nnQueryArray(net, query, input, output); + + const R expected_output = sigmoid(input[0] * weights[0]); + printf("\nOutput: %f, Expected: %f\n", output[0], expected_output); + TEST_TRUE(double_eq(output[0], expected_output, EPS)); + + nnDeleteQueryObject(&query); + nnDeleteNet(&net); +} + +TEST_CASE(neuralnet_xor_test) { + const int num_layers = 2; + const int layer_sizes[] = { 2, 2, 1 }; + const nnActivation layer_activations[] = { nnRelu, nnIdentity }; + const R weights[] = { + 1, 1, 1, 1, // First (hidden) layer. + 1, -2 // Second (output) layer. + }; + const R biases[] = { + 0, -1, // First (hidden) layer. + 0 // Second (output) layer. + }; + + nnNeuralNetwork* net = nnMakeNet(num_layers, layer_sizes, layer_activations); + assert(net); + nnSetWeights(net, weights); + nnSetBiases(net, biases); + + // First layer weights. + TEST_EQUAL(nnMatrixAt(&net->weights[0], 0, 0), 1); + TEST_EQUAL(nnMatrixAt(&net->weights[0], 0, 1), 1); + TEST_EQUAL(nnMatrixAt(&net->weights[0], 0, 2), 1); + TEST_EQUAL(nnMatrixAt(&net->weights[0], 0, 3), 1); + // Second layer weights. + TEST_EQUAL(nnMatrixAt(&net->weights[1], 0, 0), 1); + TEST_EQUAL(nnMatrixAt(&net->weights[1], 0, 1), -2); + // First layer biases. + TEST_EQUAL(nnMatrixAt(&net->biases[0], 0, 0), 0); + TEST_EQUAL(nnMatrixAt(&net->biases[0], 0, 1), -1); + // Second layer biases. + TEST_EQUAL(nnMatrixAt(&net->biases[1], 0, 0), 0); + + // Test. + + #define M 4 + + nnQueryObject* query = nnMakeQueryObject(net, /*num_inputs=*/M); + + const R test_inputs[M][2] = { { 0., 0. }, { 1., 0. }, { 0., 1. }, { 1., 1. } }; + nnMatrix test_inputs_matrix = nnMatrixMake(M, 2); + nnMatrixInit(&test_inputs_matrix, (const R*)test_inputs); + nnQuery(net, query, &test_inputs_matrix); + + const R expected_outputs[M] = { 0., 1., 1., 0. }; + for (int i = 0; i < M; ++i) { + const R test_output = nnMatrixAt(nnNetOutputs(query), i, 0); + printf("\nInput: (%f, %f), Output: %f, Expected: %f\n", + test_inputs[i][0], test_inputs[i][1], test_output, expected_outputs[i]); + } + for (int i = 0; i < M; ++i) { + const R test_output = nnMatrixAt(nnNetOutputs(query), i, 0); + TEST_TRUE(double_eq(test_output, expected_outputs[i], OUTPUT_EPS)); + } + + nnDeleteQueryObject(&query); + nnDeleteNet(&net); +} diff --git a/src/lib/test/test.h b/src/lib/test/test.h new file mode 100644 index 0000000..fd8dc22 --- /dev/null +++ b/src/lib/test/test.h @@ -0,0 +1,185 @@ +// SPDX-License-Identifier: MIT +#pragma once + +#ifdef UNIT_TEST + +#include +#include +#include +#include + +#if defined(__DragonFly__) || defined(__FreeBSD__) || defined(__FreeBSD_kernel__) || \ + defined(__NetBSD__) || defined(__OpenBSD__) +#define USE_SYSCTL_FOR_ARGS 1 +// clang-format off +#include +#include +// clang-format on +#include // getpid +#endif + +struct test_file_metadata; + +struct test_failure { + bool present; + const char *message; + const char *file; + int line; +}; + +struct test_case_metadata { + void (*fn)(struct test_case_metadata *, struct test_file_metadata *); + struct test_failure failure; + const char *name; + struct test_case_metadata *next; +}; + +struct test_file_metadata { + bool registered; + const char *name; + struct test_file_metadata *next; + struct test_case_metadata *tests; +}; + +struct test_file_metadata __attribute__((weak)) * test_file_head; + +#define SET_FAILURE(_message) \ + metadata->failure = (struct test_failure) { \ + .message = _message, .file = __FILE__, .line = __LINE__, .present = true, \ + } + +#define TEST_EQUAL(a, b) \ + do { \ + if ((a) != (b)) { \ + SET_FAILURE(#a " != " #b); \ + return; \ + } \ + } while (0) + +#define TEST_TRUE(a) \ + do { \ + if (!(a)) { \ + SET_FAILURE(#a " is not true"); \ + return; \ + } \ + } while (0) + +#define TEST_STREQUAL(a, b) \ + do { \ + if (strcmp(a, b) != 0) { \ + SET_FAILURE(#a " != " #b); \ + return; \ + } \ + } while (0) + +#define TEST_CASE(_name) \ + static void __test_h_##_name(struct test_case_metadata *, \ + struct test_file_metadata *); \ + static struct test_file_metadata __test_h_file; \ + static struct test_case_metadata __test_h_meta_##_name = { \ + .name = #_name, \ + .fn = __test_h_##_name, \ + }; \ + static void __attribute__((constructor(101))) __test_h_##_name##_register(void) { \ + __test_h_meta_##_name.next = __test_h_file.tests; \ + __test_h_file.tests = &__test_h_meta_##_name; \ + if (!__test_h_file.registered) { \ + __test_h_file.name = __FILE__; \ + __test_h_file.next = test_file_head; \ + test_file_head = &__test_h_file; \ + __test_h_file.registered = true; \ + } \ + } \ + static void __test_h_##_name( \ + struct test_case_metadata *metadata __attribute__((unused)), \ + struct test_file_metadata *file_metadata __attribute__((unused))) + +extern void __attribute__((weak)) (*test_h_unittest_setup)(void); +/// Run defined tests, return true if all tests succeeds +/// @param[out] tests_run if not NULL, set to whether tests were run +static inline void __attribute__((constructor(102))) run_tests(void) { + bool should_run = false; +#ifdef USE_SYSCTL_FOR_ARGS + int mib[] = { + CTL_KERN, +#if defined(__NetBSD__) || defined(__OpenBSD__) + KERN_PROC_ARGS, + getpid(), + KERN_PROC_ARGV, +#else + KERN_PROC, + KERN_PROC_ARGS, + getpid(), +#endif + }; + char *arg = NULL; + size_t arglen; + sysctl(mib, sizeof(mib) / sizeof(mib[0]), NULL, &arglen, NULL, 0); + arg = malloc(arglen); + sysctl(mib, sizeof(mib) / sizeof(mib[0]), arg, &arglen, NULL, 0); +#else + FILE *cmdlinef = fopen("/proc/self/cmdline", "r"); + char *arg = NULL; + int arglen; + fscanf(cmdlinef, "%ms%n", &arg, &arglen); + fclose(cmdlinef); +#endif + for (char *pos = arg; pos < arg + arglen; pos += strlen(pos) + 1) { + if (strcmp(pos, "--unittest") == 0) { + should_run = true; + break; + } + } + free(arg); + + if (!should_run) { + return; + } + + if (&test_h_unittest_setup) { + test_h_unittest_setup(); + } + + struct test_file_metadata *i = test_file_head; + int failed = 0, success = 0; + while (i) { + fprintf(stderr, "Running tests from %s:\n", i->name); + struct test_case_metadata *j = i->tests; + while (j) { + fprintf(stderr, "\t%s ... ", j->name); + j->failure.present = false; + j->fn(j, i); + if (j->failure.present) { + fprintf(stderr, "failed (%s at %s:%d)\n", j->failure.message, + j->failure.file, j->failure.line); + failed++; + } else { + fprintf(stderr, "passed\n"); + success++; + } + j = j->next; + } + fprintf(stderr, "\n"); + i = i->next; + } + int total = failed + success; + fprintf(stderr, "Test results: passed %d/%d, failed %d/%d\n", success, total, + failed, total); + exit(failed == 0 ? EXIT_SUCCESS : EXIT_FAILURE); +} + +#else + +#include + +#define TEST_CASE(name) static void __attribute__((unused)) __test_h_##name(void) + +#define TEST_EQUAL(a, b) \ + (void)(a); \ + (void)(b) +#define TEST_TRUE(a) (void)(a) +#define TEST_STREQUAL(a, b) \ + (void)(a); \ + (void)(b) + +#endif diff --git a/src/lib/test/test_main.c b/src/lib/test/test_main.c new file mode 100644 index 0000000..4cce7f6 --- /dev/null +++ b/src/lib/test/test_main.c @@ -0,0 +1,3 @@ +int main() { + return 0; +} diff --git a/src/lib/test/test_util.h b/src/lib/test/test_util.h new file mode 100644 index 0000000..8abb99a --- /dev/null +++ b/src/lib/test/test_util.h @@ -0,0 +1,22 @@ +#pragma once + +#include + +#include + +// General epsilon for comparing values. +static const R EPS = 1e-10; + +// Epsilon for comparing network weights after training. +static const R WEIGHT_EPS = 0.01; + +// Epsilon for comparing network outputs after training. +static const R OUTPUT_EPS = 0.01; + +static inline bool double_eq(double a, double b, double eps) { + return fabs(a - b) <= eps; +} + +static inline R lerp(R a, R b, R t) { + return a + t*(b-a); +} diff --git a/src/lib/test/train_linear_perceptron_non_origin_test.c b/src/lib/test/train_linear_perceptron_non_origin_test.c new file mode 100644 index 0000000..5a320ac --- /dev/null +++ b/src/lib/test/train_linear_perceptron_non_origin_test.c @@ -0,0 +1,67 @@ +#include + +#include +#include +#include "activation.h" +#include "neuralnet_impl.h" + +#include "test.h" +#include "test_util.h" + +#include + +TEST_CASE(neuralnet_train_linear_perceptron_non_origin_test) { + const int num_layers = 1; + const int layer_sizes[] = { 1, 1 }; + const nnActivation layer_activations[] = { nnIdentity }; + + nnNeuralNetwork* net = nnMakeNet(num_layers, layer_sizes, layer_activations); + assert(net); + + // Train. + + // Try to learn the Y = 2X + 1 line. + #define N 2 + const R inputs[N] = { 0., 1. }; + const R targets[N] = { 1., 3. }; + + nnMatrix inputs_matrix = nnMatrixMake(N, 1); + nnMatrix targets_matrix = nnMatrixMake(N, 1); + nnMatrixInit(&inputs_matrix, inputs); + nnMatrixInit(&targets_matrix, targets); + + nnTrainingParams params = { + .learning_rate = 0.7, + .max_iterations = 20, + .seed = 0, + .weight_init = nnWeightInit01, + .debug = false, + }; + + nnTrain(net, &inputs_matrix, &targets_matrix, ¶ms); + + const R weight = nnMatrixAt(&net->weights[0], 0, 0); + const R expected_weight = 2.0; + printf("\nTrained network weight: %f, Expected: %f\n", weight, expected_weight); + TEST_TRUE(double_eq(weight, expected_weight, WEIGHT_EPS)); + + const R bias = nnMatrixAt(&net->biases[0], 0, 0); + const R expected_bias = 1.0; + printf("Trained network bias: %f, Expected: %f\n", bias, expected_bias); + TEST_TRUE(double_eq(bias, expected_bias, WEIGHT_EPS)); + + // Test. + + nnQueryObject* query = nnMakeQueryObject(net, /*num_inputs=*/1); + + const R test_input[] = { 2.3 }; + R test_output[1]; + nnQueryArray(net, query, test_input, test_output); + + const R expected_output = test_input[0] * expected_weight + expected_bias; + printf("Output: %f, Expected: %f\n", test_output[0], expected_output); + TEST_TRUE(double_eq(test_output[0], expected_output, OUTPUT_EPS)); + + nnDeleteQueryObject(&query); + nnDeleteNet(&net); +} diff --git a/src/lib/test/train_linear_perceptron_test.c b/src/lib/test/train_linear_perceptron_test.c new file mode 100644 index 0000000..2b1336d --- /dev/null +++ b/src/lib/test/train_linear_perceptron_test.c @@ -0,0 +1,62 @@ +#include + +#include +#include +#include "activation.h" +#include "neuralnet_impl.h" + +#include "test.h" +#include "test_util.h" + +#include + +TEST_CASE(neuralnet_train_linear_perceptron_test) { + const int num_layers = 1; + const int layer_sizes[] = { 1, 1 }; + const nnActivation layer_activations[] = { nnIdentity }; + + nnNeuralNetwork* net = nnMakeNet(num_layers, layer_sizes, layer_activations); + assert(net); + + // Train. + + // Try to learn the Y=X line. + #define N 2 + const R inputs[N] = { 0., 1. }; + const R targets[N] = { 0., 1. }; + + nnMatrix inputs_matrix = nnMatrixMake(N, 1); + nnMatrix targets_matrix = nnMatrixMake(N, 1); + nnMatrixInit(&inputs_matrix, inputs); + nnMatrixInit(&targets_matrix, targets); + + nnTrainingParams params = { + .learning_rate = 0.7, + .max_iterations = 10, + .seed = 0, + .weight_init = nnWeightInit01, + .debug = false, + }; + + nnTrain(net, &inputs_matrix, &targets_matrix, ¶ms); + + const R weight = nnMatrixAt(&net->weights[0], 0, 0); + const R expected_weight = 1.0; + printf("\nTrained network weight: %f, Expected: %f\n", weight, expected_weight); + TEST_TRUE(double_eq(weight, expected_weight, WEIGHT_EPS)); + + // Test. + + nnQueryObject* query = nnMakeQueryObject(net, /*num_inputs=*/1); + + const R test_input[] = { 2.3 }; + R test_output[1]; + nnQueryArray(net, query, test_input, test_output); + + const R expected_output = test_input[0]; + printf("Output: %f, Expected: %f\n", test_output[0], expected_output); + TEST_TRUE(double_eq(test_output[0], expected_output, OUTPUT_EPS)); + + nnDeleteQueryObject(&query); + nnDeleteNet(&net); +} diff --git a/src/lib/test/train_sigmoid_test.c b/src/lib/test/train_sigmoid_test.c new file mode 100644 index 0000000..588e7ca --- /dev/null +++ b/src/lib/test/train_sigmoid_test.c @@ -0,0 +1,66 @@ +#include + +#include +#include +#include "activation.h" +#include "neuralnet_impl.h" + +#include "test.h" +#include "test_util.h" + +#include + +TEST_CASE(neuralnet_train_sigmoid_test) { + const int num_layers = 1; + const int layer_sizes[] = { 1, 1 }; + const nnActivation layer_activations[] = { nnSigmoid }; + + nnNeuralNetwork* net = nnMakeNet(num_layers, layer_sizes, layer_activations); + assert(net); + + // Train. + + // Try to learn the sigmoid function. + #define N 3 + R inputs[N]; + R targets[N]; + for (int i = 0; i < N; ++i) { + inputs[i] = lerp(-1, +1, (R)i / (R)(N-1)); + targets[i] = sigmoid(inputs[i]); + } + + nnMatrix inputs_matrix = nnMatrixMake(N, 1); + nnMatrix targets_matrix = nnMatrixMake(N, 1); + nnMatrixInit(&inputs_matrix, inputs); + nnMatrixInit(&targets_matrix, targets); + + nnTrainingParams params = { + .learning_rate = 0.9, + .max_iterations = 100, + .seed = 0, + .weight_init = nnWeightInit01, + .debug = false, + }; + + nnTrain(net, &inputs_matrix, &targets_matrix, ¶ms); + + const R weight = nnMatrixAt(&net->weights[0], 0, 0); + const R expected_weight = 1.0; + printf("\nTrained network weight: %f, Expected: %f\n", weight, expected_weight); + TEST_TRUE(double_eq(weight, expected_weight, WEIGHT_EPS)); + + // Test. + + nnQueryObject* query = nnMakeQueryObject(net, /*num_inputs=*/1); + + const R test_input[] = { 0.3 }; + R test_output[1]; + nnQueryArray(net, query, test_input, test_output); + + const R expected_output = 0.574442516811659; // sigmoid(0.3) + printf("Output: %f, Expected: %f\n", test_output[0], expected_output); + TEST_TRUE(double_eq(test_output[0], expected_output, OUTPUT_EPS)); + + nnDeleteQueryObject(&query); + nnDeleteNet(&net); +} diff --git a/src/lib/test/train_xor_test.c b/src/lib/test/train_xor_test.c new file mode 100644 index 0000000..6ddc6e0 --- /dev/null +++ b/src/lib/test/train_xor_test.c @@ -0,0 +1,66 @@ +#include + +#include +#include +#include "activation.h" +#include "neuralnet_impl.h" + +#include "test.h" +#include "test_util.h" + +#include + +TEST_CASE(neuralnet_train_xor_test) { + const int num_layers = 2; + const int layer_sizes[] = { 2, 2, 1 }; + const nnActivation layer_activations[] = { nnRelu, nnIdentity }; + + nnNeuralNetwork* net = nnMakeNet(num_layers, layer_sizes, layer_activations); + assert(net); + + // Train. + + #define N 4 + const R inputs[N][2] = { { 0., 0. }, { 0., 1. }, { 1., 0. }, { 1., 1. } }; + const R targets[N] = { 0., 1., 1., 0. }; + + nnMatrix inputs_matrix = nnMatrixMake(N, 2); + nnMatrix targets_matrix = nnMatrixMake(N, 1); + nnMatrixInit(&inputs_matrix, (const R*)inputs); + nnMatrixInit(&targets_matrix, targets); + + nnTrainingParams params = { + .learning_rate = 0.1, + .max_iterations = 500, + .seed = 0, + .weight_init = nnWeightInit01, + .debug = false, + }; + + nnTrain(net, &inputs_matrix, &targets_matrix, ¶ms); + + // Test. + + #define M 4 + + nnQueryObject* query = nnMakeQueryObject(net, /*num_inputs=*/M); + + const R test_inputs[M][2] = { { 0., 0. }, { 1., 0. }, { 0., 1. }, { 1., 1. } }; + nnMatrix test_inputs_matrix = nnMatrixMake(M, 2); + nnMatrixInit(&test_inputs_matrix, (const R*)test_inputs); + nnQuery(net, query, &test_inputs_matrix); + + const R expected_outputs[M] = { 0., 1., 1., 0. }; + for (int i = 0; i < M; ++i) { + const R test_output = nnMatrixAt(nnNetOutputs(query), i, 0); + printf("\nInput: (%f, %f), Output: %f, Expected: %f\n", + test_inputs[i][0], test_inputs[i][1], test_output, expected_outputs[i]); + } + for (int i = 0; i < M; ++i) { + const R test_output = nnMatrixAt(nnNetOutputs(query), i, 0); + TEST_TRUE(double_eq(test_output, expected_outputs[i], OUTPUT_EPS)); + } + + nnDeleteQueryObject(&query); + nnDeleteNet(&net); +} -- cgit v1.2.3