1 files changed, 286 insertions, 318 deletions

diff --git a/include/ann.h b/include/ann.h
index 087daa5..e34176a 100755
--- a/include/ann.h
+++ b/include/ann.h
@@ -1,10 +1,14 @@
-// ann.h - Artificial Neural Nework Libray
+// ann.h - Artificial Neural Network Library
+//
+// Thread safety: Not thread-safe. External synchronization required.
+// Memory: Single-allocation design — all neurons, weights, and deltas
+//   stored in one contiguous buffer for easy persistence and cloning.
+//   Use ann_copy for deep clone; ann_free to release.
 
 #ifndef ANN_H
 #define ANN_H
 
-#include <stdint.h>
-#include "./type.h"
+#include "type.h"
 
 #ifdef __cplusplus
 extern "C" {
@@ -21,64 +25,66 @@ typedef enum {
 } ann_activation_t;
 
 typedef struct {
+  void *buffer;
+
   // The full size of the allocated structure
-  uint_t n;
+  u64 n;
 
   // The number of layers in the neural network
-  uint_t layer_n;
+  u64 layer_n;
 
   // The neuron count ( hidden )
-  uint_t neuron_n;
+  u64 neuron_n;
 
   // The total number of weights and biases
-  uint_t weight_n;
+  u64 weight_n;
 
   // The number of neurons in each layer
-  uint_t *layer_neuron_n;
+  u64 *layer_neuron_n;
 
   // The pointer to the neurons
-  fp_t *neuron;
+  f64 *neuron;
 
   // The weights and biases for each neuron
   //   - w_lji denotes the connection between the jth neuron in layer l, and
   //     the ith neuron in the previous layer
   //   - b_lj denotes the bias for the jth neuron in layer l
   // { w_000, w_001, w_002, .... , b_00, ...., w_ijk, b_ij }`
-  fp_t *weight;
+  f64 *weight;
 
   // The delta between between the actual and the cost function
-  fp_t *delta;
+  f64 *delta;
 
   // The activation function used in the hidden layer neurons
-  fp_t (*activation_hidden)(fp_t);
+  f64 (*activation_hidden)(f64);
 
   // The partial derivative of the activation function used in the hidden
   // layer neurons
-  fp_t (*activation_hidden_partial)(fp_t);
+  f64 (*activation_hidden_partial)(f64);
 
   // The activation function used in the ouput neurons
-  fp_t (*activation_output)(fp_t);
+  f64 (*activation_output)(f64);
 
   // The partial derivative of the activation function used in the output
   // layer neurons
-  fp_t (*activation_output_partial)(fp_t);
+  f64 (*activation_output_partial)(f64);
 } ann_t;
 
-ann_t *ann_init(uint_t, uint_t *);
-ann_t *ann_copy(ann_t const *);
+int ann_init(ann_t *, u64, u64 *);
+void ann_copy(ann_t *, ann_t const *);
 void ann_free(ann_t *);
 void ann_random(ann_t *);
 
-void ann_propagation_forward(ann_t *, fp_t const *const, fp_t *);
-void ann_propagation_backward(ann_t *, fp_t const *, fp_t *, fp_t const *,
-                              fp_t);
-void ann_train_numeric(ann_t *, fp_t const *, fp_t const *, fp_t);
+void ann_propagation_forward(ann_t *, f64 const *const, f64 *);
+void ann_propagation_backward(ann_t *, f64 const *, f64 *, f64 const *,
+                              f64);
+void ann_train_numeric(ann_t *, f64 const *, f64 const *, u64, f64);
 
-fp_t ann_error_total(fp_t const *, fp_t const *, uint_t);
+f64 ann_error_total(f64 const *, f64 const *, u64);
 void ann_set_activation(ann_t *, ann_activation_t, ann_activation_t);
 
 void ann_print_weight(ann_t *);
-void ann_print_neuron(ann_t *, fp_t const *const, fp_t const *const);
+void ann_print_neuron(ann_t *, f64 const *const, f64 const *const);
 
 #ifdef __cplusplus
 }
@@ -87,7 +93,7 @@ void ann_print_neuron(ann_t *, fp_t const *const, fp_t const *const);
 #endif // ANN_H
 
 ////////////////////////////////////////////////////////////////////////////////
-// IMPELEMENTATION SECTION
+// IMPLEMENTATION SECTION
 ////////////////////////////////////////////////////////////////////////////////
 
 #ifdef ANN_IMPLEMENTATION
@@ -108,27 +114,27 @@ void ann_print_neuron(ann_t *, fp_t const *const, fp_t const *const);
 #define ELU_ALPHA 0.2
 #define LRELU_ALPHA 0.2
 
-static fp_t ann_random_range(fp_t, fp_t);
-
-static fp_t ann_error(fp_t, fp_t);
-static fp_t ann_error_partial(fp_t, fp_t);
+static void ann_zero_weights(ann_t *);
+static f64 ann_random_range(f64, f64);
+static f64 ann_error(f64, f64);
+static f64 ann_error_partial(f64, f64);
 
-static fp_t ann_activation_identity(fp_t);
-static fp_t ann_activation_identity_partial(fp_t);
-static fp_t ann_activation_binary(fp_t);
-static fp_t ann_activation_binary_partial(fp_t);
-static fp_t ann_activation_sigmoid(fp_t);
-static fp_t ann_activation_sigmoid_partial(fp_t);
-static fp_t ann_activation_relu(fp_t);
-static fp_t ann_activation_relu_partial(fp_t);
-static fp_t ann_activation_elu(fp_t);
-static fp_t ann_activation_elu_partial(fp_t);
-static fp_t ann_activation_lrelu(fp_t);
-static fp_t ann_activation_lrelu_partial(fp_t);
-static fp_t ann_activation_tanh(fp_t);
-static fp_t ann_activation_tanh_partial(fp_t);
+static f64 ann_activation_identity(f64);
+static f64 ann_activation_identity_partial(f64);
+static f64 ann_activation_binary(f64);
+static f64 ann_activation_binary_partial(f64);
+static f64 ann_activation_sigmoid(f64);
+static f64 ann_activation_sigmoid_partial(f64);
+static f64 ann_activation_relu(f64);
+static f64 ann_activation_relu_partial(f64);
+static f64 ann_activation_elu(f64);
+static f64 ann_activation_elu_partial(f64);
+static f64 ann_activation_lrelu(f64);
+static f64 ann_activation_lrelu_partial(f64);
+static f64 ann_activation_tanh(f64);
+static f64 ann_activation_tanh_partial(f64);
 
-static fp_t (*ACTIVATION[][2])(fp_t) = {
+static f64 (*ACTIVATION[][2])(f64) = {
     {ann_activation_identity, ann_activation_identity_partial}, // IDENTITY
     {ann_activation_binary, ann_activation_binary_partial},     // BINARY
     {ann_activation_sigmoid, ann_activation_sigmoid_partial},   // SIGMOID
@@ -138,103 +144,135 @@ static fp_t (*ACTIVATION[][2])(fp_t) = {
     {ann_activation_tanh, ann_activation_tanh_partial},         // TANH
 };
 
-// ann_init()
-//
-// Initialize the neural network
-//
-// layer_n - The layer count for the layer count
-// layer_neuron_n - An array of neuron counts for each layer
-//
-// return - The created ann_t instance
-
-ann_t *ann_init(uint_t layer_n, uint_t *layer_neuron_n) {
-  assert(layer_n >= 2);
+int ann_init(ann_t *ann, u64 layer_n, u64 *layer_neuron_n) {
+  if (!ann || !layer_neuron_n || layer_n < 2)
+    return 0;
 
-  uint_t neuron_n = 0;
-  uint_t weight_n = 0;
+  u64 neuron_n = 0;
+  u64 weight_n = 0;
 
-  // Calculate the number of neurons and weights / biases
-  uint_t l = 1;
+  u64 l = 1;
   for (; l < layer_n - 1; l++) {
+    if (neuron_n > UINT64_MAX - layer_neuron_n[l])
+      return 0;
     neuron_n += layer_neuron_n[l];
-    weight_n += layer_neuron_n[l] * (layer_neuron_n[l - 1] + 1);
+    if (layer_neuron_n[l - 1] == UINT64_MAX)
+      return 0;
+    {
+      u64 rhs = layer_neuron_n[l - 1] + 1;
+      if (layer_neuron_n[l] > UINT64_MAX / rhs)
+        return 0;
+      u64 prod = layer_neuron_n[l] * rhs;
+      if (weight_n > UINT64_MAX - prod)
+        return 0;
+      weight_n += prod;
+    }
+  }
+
+  if (layer_neuron_n[l - 1] == UINT64_MAX)
+    return 0;
+  {
+    u64 rhs = layer_neuron_n[l - 1] + 1;
+    if (layer_neuron_n[l] > UINT64_MAX / rhs)
+      return 0;
+    u64 prod = layer_neuron_n[l] * rhs;
+    if (weight_n > UINT64_MAX - prod)
+      return 0;
+    weight_n += prod;
   }
 
-  weight_n += layer_neuron_n[l] * (layer_neuron_n[l - 1] + 1);
+  if (layer_n > SIZE_MAX / sizeof(u64))
+    return 0;
+  size_t layer_neuron_sz = sizeof(u64) * layer_n;
+  size_t padding = (sizeof(f64) - (layer_neuron_sz % sizeof(f64))) % sizeof(f64);
 
-  uint_t n =
-      sizeof(ann_t) +              // ANN
-      (sizeof(uint_t) * layer_n) + // layer_neuron_n[]
-      (sizeof(fp_t) * (neuron_n +  // neuron[]
-                       weight_n +  // weight[]
-                       neuron_n + layer_neuron_n[layer_n - 1])); // delta[]
+  u64 total_fp = neuron_n + weight_n;
+  if (total_fp < neuron_n || total_fp < weight_n)
+    return 0;
+  if (total_fp > UINT64_MAX - neuron_n)
+    return 0;
+  total_fp += neuron_n;
+  if (total_fp > UINT64_MAX - layer_neuron_n[layer_n - 1])
+    return 0;
+  total_fp += layer_neuron_n[layer_n - 1];
+  if (total_fp > SIZE_MAX / sizeof(f64))
+    return 0;
+  u64 n = (u64)(layer_neuron_sz + padding + sizeof(f64) * total_fp);
 
-  // Allocate everything as one structure
-  ann_t *ann = (ann_t *)malloc(n);
+  ann->buffer = malloc(n);
+  if (!ann->buffer)
+    return 0;
 
-  // ann_t | layer_neuron_n[] | neuron[] | weight[] | delta[]
   ann->n = n;
   ann->layer_n = layer_n;
   ann->weight_n = weight_n;
   ann->neuron_n = neuron_n;
-  ann->layer_neuron_n = (uint_t *)((uint8_t *)ann + sizeof(ann_t));
-  memcpy(ann->layer_neuron_n, layer_neuron_n, sizeof(uint_t) * layer_n);
-  ann->neuron = (fp_t *)(ann->layer_neuron_n + ann->layer_n);
+  ann->layer_neuron_n = (u64 *)ann->buffer;
+  memcpy(ann->layer_neuron_n, layer_neuron_n, sizeof(u64) * layer_n);
+  ann->neuron = (f64 *)((char *)ann->buffer + layer_neuron_sz + padding);
   ann->weight = ann->neuron + ann->neuron_n;
   ann->delta = ann->weight + ann->weight_n;
 
   ann_set_activation(ann, SIGMOID, SIGMOID);
+  ann_zero_weights(ann);
+  return 1;
+}
 
-  return ann;
+void ann_free(ann_t *ann) {
+  if (ann) {
+    free(ann->buffer);
+    ann->buffer = NULL;
+    ann->layer_neuron_n = NULL;
+    ann->neuron = NULL;
+    ann->weight = NULL;
+    ann->delta = NULL;
+  }
 }
 
-// ann_free()
-//
-// Free the neural network's memory
-//
-// ann - The instance of ann_t to free
+void ann_copy(ann_t *dest, ann_t const *source) {
+  if (!dest || !source || !source->buffer)
+    return;
+  void *new_buffer = malloc(source->n);
+  if (!new_buffer)
+    return;
 
-void ann_free(ann_t *ann) { free(ann); }
+  memcpy(new_buffer, source->buffer, source->n);
 
-// ann_copy()
-//
-// Make a copy of the given ann_t instance
-//
-// ann - The instance of ann_t to copy
-//
-// return - The copied ann_t instance
+  if (dest->buffer)
+    free(dest->buffer);
 
-ann_t *ann_copy(ann_t const *ann) {
-  ann_t *copy = (ann_t *)malloc(ann->n);
-  memcpy(copy, ann, ann->n);
+  size_t layer_neuron_sz = sizeof(u64) * source->layer_n;
+  size_t padding = (sizeof(f64) - (layer_neuron_sz % sizeof(f64))) % sizeof(f64);
 
-  return copy;
+  dest->buffer = new_buffer;
+  dest->n = source->n;
+  dest->layer_n = source->layer_n;
+  dest->weight_n = source->weight_n;
+  dest->neuron_n = source->neuron_n;
+  dest->layer_neuron_n = (u64 *)dest->buffer;
+  dest->neuron = (f64 *)((char *)dest->buffer + layer_neuron_sz + padding);
+  dest->weight = dest->neuron + source->neuron_n;
+  dest->delta = dest->weight + source->weight_n;
+  dest->activation_hidden = source->activation_hidden;
+  dest->activation_hidden_partial = source->activation_hidden_partial;
+  dest->activation_output = source->activation_output;
+  dest->activation_output_partial = source->activation_output_partial;
 }
 
-// ann_propagation_forward()
-//
-// Perform forward propagation on the ann_t instance
-//
-// ann - The ann_t instance to perform the propagation on
-// input - An array containing the input vector
-// output - The destination array for the output vector
-
-void ann_propagation_forward(ann_t *ann, fp_t const *const input,
-                             fp_t *output) {
-  // printf("ann:\t\t%p\n", ann);
-  // printf("output:\t\t%p\n", output);
-  fp_t *w_ij = ann->weight;
-  // printf("ann->weight:\t%p\n", w_ij);
-  fp_t const *x = input; // Input neuron into y
-  fp_t *y = ann->neuron; // The current neuron being calculated
+void ann_propagation_forward(ann_t *ann, f64 const *const input,
+                              f64 *output) {
+  if (!ann || !input || !output)
+    return;
+  f64 *w_ij = ann->weight;
+  f64 const *x = input;
+  f64 *y = ann->neuron;
 
-  uint_t l = 1;
+  u64 l = 1;
 
   for (; l < ann->layer_n - 1; l++) {
-    for (uint_t j = 0; j < ann->layer_neuron_n[l]; j++) {
-      // printf("neuron[%u]:\t%p\n", j, &(y[j]));
+    for (u64 j = 0; j < ann->layer_neuron_n[l]; j++) {
       y[j] = 0;
-      for (uint_t i = 0; i < ann->layer_neuron_n[l - 1]; i++) {
+      for (u64 i = 0; i < ann->layer_neuron_n[l - 1]; i++) {
         y[j] += x[i] * *w_ij++;
       }
 
@@ -246,12 +284,10 @@ void ann_propagation_forward(ann_t *ann, fp_t const *const input,
     y += ann->layer_neuron_n[l];
   }
 
-  // Last layer
-  for (uint_t j = 0; j < ann->layer_neuron_n[l]; j++) {
-    // printf("output[%u]:\t%p\n", j, &(output[j]));
+  for (u64 j = 0; j < ann->layer_neuron_n[l]; j++) {
     output[j] = 0;
 
-    for (uint_t i = 0; i < ann->layer_neuron_n[l - 1]; i++) {
+    for (u64 i = 0; i < ann->layer_neuron_n[l - 1]; i++) {
       output[j] += x[i] * *w_ij++;
     }
 
@@ -260,321 +296,253 @@ void ann_propagation_forward(ann_t *ann, fp_t const *const input,
   }
 }
 
-// ann_propagation_backward()
-//
-// Perform backpropagation on the ann_t instance
-//
-// ann - The ann_t instance to perform backpropagation upon
-// input - Input vector array
-// output - Output vector array
-// target - Target output vector array
-// rate - Learning rate
-
-void ann_propagation_backward(ann_t *ann, fp_t const *input, fp_t *output,
-                              fp_t const *target, fp_t rate) {
-  int_t l = ann->layer_n - 1;
-  uint_t i, j, q;
+void ann_propagation_backward(ann_t *ann, f64 const *input, f64 *output,
+                               f64 const *target, f64 rate) {
+  if (!ann || !input || !output || !target)
+    return;
+  i64 l = (i64)ann->layer_n - 1;
+  u64 i = 0, j = 0, q = 0;
 
-  // First output layer delta
-  fp_t *d_j = ann->delta + ann->neuron_n;
+  f64 *d_j = ann->delta + ann->neuron_n;
 
-  // Output Deltas
   for (j = 0; j < ann->layer_neuron_n[l]; j++) {
-    d_j[j] = ann->activation_output_partial(output[j]) *
-             ann_error_partial(output[j], target[j]);
+    d_j[j] = ann_error_partial(output[j], target[j]) * ann->activation_output_partial(output[j]);
   }
 
-  // First weight in the set between the last layer and the current
-  fp_t *w_jq = ann->weight + ann->weight_n -
-               ann->layer_neuron_n[l] * ann->layer_neuron_n[l - 1] -
-               ann->layer_neuron_n[l];
+ {
+    u64 last_layer_n = ann->layer_neuron_n[l];
+    u64 prev_layer_plus_bias = ann->layer_neuron_n[l - 1] + 1;
+    if (last_layer_n > UINT64_MAX / prev_layer_plus_bias)
+      return;
+    f64 *w_jq = ann->weight + ann->weight_n - last_layer_n * prev_layer_plus_bias;
 
-  fp_t *o_j = ann->neuron + ann->neuron_n;
-  fp_t *d_q;
+    f64 *o_j = ann->neuron + ann->neuron_n;
+    f64 *d_q;
 
-  l--;
+    l--;
 
-  // Hidden Deltas
-  for (; l > 0; --l) {
-    d_q = d_j;
-    d_j -= ann->layer_neuron_n[l];
-    o_j -= ann->layer_neuron_n[l];
+    for (; l > 0; --l) {
+      d_q = d_j;
+      d_j -= ann->layer_neuron_n[l];
+      o_j -= ann->layer_neuron_n[l];
 
-    for (j = 0; j < ann->layer_neuron_n[l]; j++) {
-      d_j[j] = 0;
-
-      for (q = 0; q < ann->layer_neuron_n[l + 1]; q++) {
-        d_j[j] += w_jq[q * (ann->layer_neuron_n[l] + 1) + j] * d_q[q];
+      for (j = 0; j < ann->layer_neuron_n[l]; j++) {
+        d_j[j] = 0;
+        for (q = 0; q < ann->layer_neuron_n[l + 1]; q++) {
+          d_j[j] += w_jq[q * (ann->layer_neuron_n[l] + 1) + j] * d_q[q];
+        }
+        d_j[j] *= ann->activation_hidden_partial(o_j[j]);
       }
 
-      d_j[j] *= ann->activation_hidden_partial(o_j[j]);
+      if (l > 1) {
+        u64 shift = ann->layer_neuron_n[l];
+        u64 prev_plus_bias = ann->layer_neuron_n[l - 1] + 1;
+        if (shift <= UINT64_MAX / prev_plus_bias) {
+          w_jq -= shift * prev_plus_bias;
+        } else {
+          return;
+        }
+      }
     }
-
-    w_jq -= (ann->layer_neuron_n[l - 1] * ann->layer_neuron_n[l] +
-             ann->layer_neuron_n[l - 1]);
   }
 
-  fp_t *w_ij = ann->weight;
+  f64 *w_ij = ann->weight;
   d_j = ann->delta;
 
   l = 1;
 
-  // Input training
   for (j = 0; j < ann->layer_neuron_n[l]; j++) {
     for (i = 0; i < ann->layer_neuron_n[l - 1]; i++) {
       *w_ij -= rate * input[i] * d_j[j];
       w_ij++;
     }
-
     *w_ij -= rate * d_j[j];
     w_ij++;
   }
 
   l++;
-  fp_t *i_i = ann->neuron;
+  f64 *i_i = ann->neuron;
 
-  // Hidden training
-  for (; l < (int_t)ann->layer_n; l++) {
+  for (; l < (i64)ann->layer_n; l++) {
     d_j += ann->layer_neuron_n[l - 1];
-
     for (j = 0; j < ann->layer_neuron_n[l]; j++) {
       for (i = 0; i < ann->layer_neuron_n[l - 1]; i++) {
         *w_ij -= rate * i_i[i] * d_j[j];
         w_ij++;
       }
-
       *w_ij -= rate * d_j[j];
       w_ij++;
     }
-
     i_i += ann->layer_neuron_n[l - 1];
   }
 }
 
-// ann_activation_identity()
-//
-// y = x
-// dy/dx = 1
-
-static fp_t ann_activation_identity(fp_t x) { return x; }
-static fp_t ann_activation_identity_partial(fp_t x) { return 1; }
-
-// ann_activation_binary()
-//
-// y = { x < 0 -> 0, x > 0 -> 1 }
-// dy/dx = 0
+static f64 ann_activation_identity(f64 x) { return x; }
+static f64 ann_activation_identity_partial(f64 x) { (void)x; return 1.0; }
 
-static fp_t ann_activation_binary(fp_t x) {
-  return (x > 0.0) ? 1.0 : (x < 0.0) ? -1.0 : 0.0;
+static f64 ann_activation_binary(f64 x) {
+  return (x >= 0.0) ? 1.0 : 0.0;
 }
-static fp_t ann_activation_binary_partial(fp_t x) { return 0.0; }
-
-// ann_activation_sign()
-//
-// y = { x < 0 -> -1, x = 0 -> 0, x > 0 -> 1 }
-// dy/dx = 0
-
-static fp_t ann_activation_sign(fp_t x) { return (x > 0.0) ? 1.0 : 0.0; }
-static fp_t ann_activation_sign_partial(fp_t x) { return 0.0; }
-
-// ann_activation_sigmoid()
-//
-// y = 1.0 / ( 1.0 + e^-x )
-// dy/dx = f(x) * ( 1.0 - f(x) )
+static f64 ann_activation_binary_partial(f64 x) { (void)x; return 0.0; }
 
-static fp_t ann_activation_sigmoid(fp_t x) { return 1.0 / (1.0 + exp(-x)); }
-static fp_t ann_activation_sigmoid_partial(fp_t x) { return x * (1.0 - x); }
-
-// ann_activation_relu()
-//
-// y = { x < 0 -> 0, x >= 0 -> x }
-// dy/dx = { x < 0 -> 0, x >= 0 -> 1 }
+static f64 ann_activation_sigmoid(f64 x) {
+  return (x >= 0) ? (1.0 / (1.0 + exp(-x))) : (exp(x) / (1.0 + exp(x)));
+}
+static f64 ann_activation_sigmoid_partial(f64 x) {
+  // x is the already-activated sigmoid output σ(z), so σ'(z) = σ(z)(1-σ(z)) = x(1-x)
+  return x * (1.0 - x);
+}
 
-static fp_t ann_activation_relu(fp_t x) { return (x > 0.0) ? x : 0.0; }
-static fp_t ann_activation_relu_partial(fp_t x) {
+static f64 ann_activation_relu(f64 x) { return (x > 0.0) ? x : 0.0; }
+static f64 ann_activation_relu_partial(f64 x) {
   return (x > 0.0) ? 1.0 : 0.0;
 }
 
-// ann_activation_elu()
-//
-// y = { x < 0 -> a * ( e^x - 1 ), x >= 0 -> x }
-// dy/dx = { x < 0 -> x + a, x >= 0 -> 1 }
-
-static fp_t ann_activation_elu(fp_t x) {
+static f64 ann_activation_elu(f64 x) {
   return (x > 0.0) ? x : ELU_ALPHA * (expm1(x));
 }
-static fp_t ann_activation_elu_partial(fp_t x) {
+static f64 ann_activation_elu_partial(f64 x) {
   return (x > 0.0) ? 1.0 : x + ELU_ALPHA;
 }
 
-// ann_activation_lelu()
-//
-// y = { x < 0 -> a * x, x >= 0 -> x }
-// dy/dx = { x < 0 -> a, x >= 0 -> 1 }
-
-static fp_t ann_activation_lrelu(fp_t x) {
+static f64 ann_activation_lrelu(f64 x) {
   return (x > 0.0) ? x : LRELU_ALPHA * x;
 }
-static fp_t ann_activation_lrelu_partial(fp_t x) {
+static f64 ann_activation_lrelu_partial(f64 x) {
   return (x > 0.0) ? 1.0 : LRELU_ALPHA;
 }
 
-// ann_activation_tanh()
-//
-// y = sinh( x ) / cosh( x ) = ( e^2x - 1 ) / ( e^2x + 1 )
-// dy/dx = 1.0 - f(x)^2
-
-static fp_t ann_activation_tanh(fp_t x) { return tanh(x); }
-static fp_t ann_activation_tanh_partial(fp_t x) { return 1.0 - (x * x); }
+static f64 ann_activation_tanh(f64 x) { return tanh(x); }
+static f64 ann_activation_tanh_partial(f64 x) { return 1.0 - (x * x); }
 
-// ann_error()
-//
-// Implementation of mean squared error for calculating the difference between
-// the current and target outputs for the ann_t instance
-//
-// output - The current ann_t output
-// target - The target ann_t output
+static f64 ann_error(f64 output, f64 target) {
+  f64 error = output - target;
+  return 0.5 * error * error;
+}
 
-fp_t ann_error(fp_t output, fp_t target) {
-  fp_t error = output - target;
-  error = 0.5 * error * error;
+static f64 ann_error_partial(f64 output, f64 target) { return (output - target); }
 
+f64 ann_error_total(f64 const *output, f64 const *target, u64 n) {
+  if (!output || !target)
+    return 0;
+  f64 error = 0;
+  for (u64 i = 0; i < n; i++) {
+    error += ann_error(output[i], target[i]);
+  }
   return error;
 }
 
-// ann_error_partial()
-//
-// The partial derivative of mean squared error where n = 2
-//
-// output - The current ann_t output
-// target - The target ann_t output
-
-fp_t ann_error_partial(fp_t output, fp_t target) { return (output - target); }
-
-// ann_error_total()
-//
-//
-//
-// output - The array of current outputs
-// target - The array of target outputs
-// n - The size of the output and target arrays
+void ann_train_numeric(ann_t *ann, f64 const *input, f64 const *target,
+                            u64 n, f64 rate) {
+  if (!ann || !input || !target || rate == 0)
+    return;
+  if (n > SIZE_MAX / sizeof(f64))
+    return;
+  f64 *output = malloc(n * sizeof(f64));
+  if (!output)
+    return;
 
-fp_t ann_error_total(fp_t const *output, fp_t const *target, uint_t n) {
-  fp_t error = 0;
-  for (uint_t i = 0; i < n; i++) {
-    error += ann_error(output[i], target[i]);
+  for (u64 w = 0; w < ann->weight_n; w++) {
+    f64 original = ann->weight[w];
+    ann->weight[w] = original + rate;
+    ann_propagation_forward(ann, input, output);
+    f64 error_p = ann_error_total(output, target, n);
+    ann->weight[w] = original - rate;
+    ann_propagation_forward(ann, input, output);
+    f64 error_m = ann_error_total(output, target, n);
+    f64 gradient = (error_p - error_m) / (2.0 * rate);
+    ann->weight[w] = original - rate * gradient;
   }
 
-  return error;
+  free(output);
 }
 
-// ann_set_activation()
-//
-// Set the activation functions for the hidden and output layers
-//
-// ann - The current ann_t instance
-// activation_hidden - The hidden layer activation function
-// activation_output - The output layer activation function
-
 void ann_set_activation(ann_t *ann, ann_activation_t activation_hidden,
-                        ann_activation_t activation_output) {
+                          ann_activation_t activation_output) {
+  if (!ann)
+    return;
+  if ((int)activation_hidden < 0 || (int)activation_output < 0)
+    return;
+  if (activation_hidden > TANH || activation_output > TANH)
+    return;
   ann->activation_hidden = ACTIVATION[activation_hidden][0];
   ann->activation_hidden_partial = ACTIVATION[activation_hidden][1];
   ann->activation_output = ACTIVATION[activation_output][0];
   ann->activation_output_partial = ACTIVATION[activation_output][1];
 }
 
-// ann_random()
-//
-// Set the weights of the ann_t instance to random numbers in the [-1, 1] range
-//
-// ann - The current ann_t instance
+static void ann_zero_weights(ann_t *ann) {
+  if (!ann)
+    return;
+  memset(ann->weight, 0, ann->weight_n * sizeof(f64));
+}
 
 void ann_random(ann_t *ann) {
-  uint_t l, j, i;
-  fp_t *wb = ann->weight;
+  if (!ann)
+    return;
+  u64 l, j, i;
+  f64 *wb = ann->weight;
   for (l = 1; l < ann->layer_n; l++) {
+    if (ann->layer_neuron_n[l - 1] == 0)
+      continue;
+    if (ann->layer_neuron_n[l] == 0)
+      continue;
+    f64 fan_in = (f64)ann->layer_neuron_n[l - 1];
+    f64 fan_out = (f64)ann->layer_neuron_n[l];
+    f64 range = sqrt(6.0 / (fan_in + fan_out));
     for (j = 0; j < ann->layer_neuron_n[l]; j++) {
-      i = 0;
-      for (; i < ann->layer_neuron_n[l - 1]; i++) {
-        *wb++ = ann_random_range(-1.0, 1.0);
+      for (i = 0; i < ann->layer_neuron_n[l - 1]; i++) {
+        *wb++ = ann_random_range(-range, range);
       }
-
       *wb++ = 0;
     }
   }
 }
 
-// ann_random_range()
-//
-// Remap a random value to a given range
-//
-// low - The lower bounds of the range
-// high - The upper bounds of the range
-
-static fp_t ann_random_range(fp_t low, fp_t high) {
-  return (low + ((fp_t)rand()) * (high - low) / RAND_MAX);
+static f64 ann_random_range(f64 low, f64 high) {
+  return (low + ((f64)rand()) * (high - low) / (f64)RAND_MAX);
 }
 
-// ann_print_neuron()
-//
-// Print a text representation of the neurons in a given ann_t instance
-//
-// ann - The current ann_t instance
-// input - The input used for the ann_t propagation being printed
-// output - The output used for the ann_t propagation being printed
-
-void ann_print_neuron(ann_t *ann, fp_t const *const input,
-                      fp_t const *const output) {
+void ann_print_neuron(ann_t *ann, f64 const *const input,
+                       f64 const *const output) {
+  if (!ann || !input || !output)
+    return;
   fputs("-> ", stderr);
-  for (uint_t i = 0; i < ann->layer_neuron_n[0]; i++) {
+  for (u64 i = 0; i < ann->layer_neuron_n[0]; i++) {
     fprintf(stderr, "  %+.*f", PRINT_PRECISION, input[i]);
   }
-
   fputs("\n", stderr);
 
-  fp_t *n = ann->neuron;
-
-  for (uint_t l = 1; l < ann->layer_n - 1; l++) {
-    fprintf(stderr, "   ", l);
-    for (uint_t i = 0; i < ann->layer_neuron_n[l]; i++) {
+  f64 *n = ann->neuron;
+  for (u64 l = 1; l < ann->layer_n - 1; l++) {
+    fprintf(stderr, "   ");
+    for (u64 i = 0; i < ann->layer_neuron_n[l]; i++) {
       fprintf(stderr, "  %+.*f", PRINT_PRECISION, *n);
       n++;
     }
-
     fputs("\n", stderr);
   }
 
   fputs("=  ", stderr);
-  for (uint_t i = 0; i < ann->layer_neuron_n[ann->layer_n - 1]; i++) {
+  for (u64 i = 0; i < ann->layer_neuron_n[ann->layer_n - 1]; i++) {
     fprintf(stderr, "  %+.*f", PRINT_PRECISION, output[i]);
   }
-
   fputs("\n\n", stderr);
 }
 
-// ann_print_weight()
-//
-// Print a text representation of weights in the given ann_t instance
-//
-// ann - The current ann_t instance
-
 void ann_print_weight(ann_t *ann) {
-  fp_t *weight = ann->weight;
-
-  for (uint_t i = 1; i < ann->layer_n; i++) {
-    for (uint_t j = 0; j < ann->layer_neuron_n[i - 1] + 1; j++) {
-      for (uint_t k = 0; k < ann->layer_neuron_n[i]; k++) {
+  if (!ann)
+    return;
+  f64 *weight = ann->weight;
+  for (u64 i = 1; i < ann->layer_n; i++) {
+    for (u64 j = 0; j < ann->layer_neuron_n[i - 1] + 1; j++) {
+      for (u64 k = 0; k < ann->layer_neuron_n[i]; k++) {
         fprintf(stderr, "  %+.*f", PRINT_PRECISION,
-                *(weight + j + k * (ann->layer_neuron_n[i - 1] + 1)));
+                *(weight + k * (ann->layer_neuron_n[i - 1] + 1) + j));
       }
-
       fputs("\n", stderr);
     }
-
-    weight += ann->layer_neuron_n[i - 1] * ann->layer_neuron_n[i] +
-              ann->layer_neuron_n[i];
-
+    weight += ann->layer_neuron_n[i] * (ann->layer_neuron_n[i - 1] + 1);
     fputs("\n", stderr);
   }
 }