1 files changed, 91 insertions, 57 deletions

diff --git a/include/genetic.h b/include/genetic.h
index eadcab5..2ebf985 100755
--- a/include/genetic.h
+++ b/include/genetic.h
@@ -1,31 +1,35 @@
 // genetic.h - A simple genetic algorithm implementation
+//
+// Thread safety: Not thread-safe. External synchronization required.
+// Behavior: rate_selection must be in (0, 1). A value of 1.0 is rejected
+//   because it results in no evolution (all individuals survive).
+//   rate_mutation must be >= 0. Returns -1.0 on any parameter error.
 
 #ifndef GENETIC_H
 #define GENETIC_H
 
-#include <assert.h>
-#include <stdint.h>
+#include "type.h"
 
-struct ga_individual_t;
+// ga_individual_t must be typedef'd by the user before including this header.
 
 // ga_crossover_t
 //
 // Type definition for a crossover function
 //
-// uintptr_t * - The address for the child pointer
-// uintptr_t - Parent pointer
-// uintptr_t - Parent pointer
+// ga_individual_t * - The address for the child individual
+// ga_individual_t - Parent individual 0
+// ga_individual_t - Parent individual 1
 
 typedef void (*ga_crossover_t)(ga_individual_t *, ga_individual_t,
-                               ga_individual_t);
+                                ga_individual_t);
 
 // ga_mutation_t
 //
-// Type definition for a crossover function
+// Type definition for a mutation function
 //
-// uintptr_t - The individual being mutated
+// ga_individual_t * - The individual being mutated
 
-typedef void (*ga_mutation_t)(ga_individual_t);
+typedef void (*ga_mutation_t)(ga_individual_t *);
 
 // ga_fitness_t()
 //
@@ -34,12 +38,12 @@ typedef void (*ga_mutation_t)(ga_individual_t);
 // uintptr_t - The individual on which to run the fitness function
 // void const * - The user defined data pointer
 
-typedef double (*ga_fitness_t)(ga_individual_t, void const *);
+typedef f64 (*ga_fitness_t)(ga_individual_t, void const *);
 
-double ga_generation(ga_individual_t *, uint32_t, double, double, ga_fitness_t,
-                     void const *, ga_crossover_t, ga_mutation_t);
+f64 ga_generation(ga_individual_t *, u32, f64, f64, ga_fitness_t,
+                    void const *, ga_crossover_t, ga_mutation_t);
 
-#endif // ga_H
+#endif
 
 #ifdef GENETIC_IMPLEMENTATION // IMPLEMENTATION
 
@@ -48,7 +52,7 @@ double ga_generation(ga_individual_t *, uint32_t, double, double, ga_fitness_t,
 #include <string.h>
 #include <tgmath.h>
 
-// ga_finess_sort()
+// ga_fitness_sort()
 //
 // Run the fitness function on
 //
@@ -59,38 +63,61 @@ double ga_generation(ga_individual_t *, uint32_t, double, double, ga_fitness_t,
 //
 // return - The error for the best performing individual
 
-static double ga_fitness_sort(ga_individual_t *population,
-                              uint32_t population_n, ga_fitness_t fitness,
-                              void const *fitness_data) {
-  double *fitness_buffer = calloc(population_n, sizeof(double));
-  double fitness_current = 0;
-  ga_individual_t individual;
+static f64 ga_fitness_sort(ga_individual_t *population,
+                                 u32 population_n, ga_fitness_t fitness,
+                                 void const *fitness_data) {
+  if (!population || !fitness || population_n == 0)
+    return -1;
 
-  // Insertion sort
-  uint32_t i, j, k;
-  for (i = 0; i < population_n; i++) {
-    fitness_current = fabs(fitness(population[i], fitness_data));
-    individual = population[i];
-    // for (j = i; j > 0 && fitness_current > fitness_buffer[j - 1]; j--) {
-    for (j = i; j > 0 && fitness_current < fitness_buffer[j - 1]; j--) {
-    }
+  f64 *fitness_buffer = calloc((size_t)population_n, sizeof(f64));
+  if (!fitness_buffer)
+    return -1;
 
-    // Shift array back by one for sorting
-    for (k = i; k > j; k--) {
-      population[k] = population[k - 1];
-      fitness_buffer[k] = fitness_buffer[k - 1];
-    }
+  for (u32 i = 0; i < population_n; i++) {
+    fitness_buffer[i] = fabs(fitness(population[i], fitness_data));
+  }
 
-    population[j] = individual;
-    fitness_buffer[j] = fitness_current;
+  // Index array for indirect sort
+  u32 *indices = malloc((size_t)population_n * sizeof(u32));
+  if (!indices) {
+    free(fitness_buffer);
+    return -1;
   }
+  for (u32 i = 0; i < population_n; i++)
+    indices[i] = i;
 
-  // for (i = 0; i < population_n; i++) {
-  //   printf("%u - %lf\n", i, fitness_buffer[i]);
-  // }
+  // Shell sort on indices by fitness value (indirect sort, only rearrange indices)
+  {
+    u32 gap = population_n / 2;
+    while (gap > 0) {
+      for (u32 i = gap; i < population_n; i++) {
+        u32 temp = indices[i];
+        f64 temp_f = fitness_buffer[temp];
+        u32 j;
+        for (j = i; j >= gap && fitness_buffer[indices[j - gap]] > temp_f; j -= gap) {
+          indices[j] = indices[j - gap];
+        }
+        indices[j] = temp;
+      }
+      gap /= 2;
+    }
+  }
+
+  // Copy sorted individuals back
+  ga_individual_t *sorted = calloc((size_t)population_n, sizeof(ga_individual_t));
+  if (!sorted) {
+    free(fitness_buffer);
+    free(indices);
+    return -1;
+  }
+  for (u32 i = 0; i < population_n; i++)
+    sorted[i] = population[indices[i]];
+  memcpy(population, sorted, (size_t)population_n * sizeof(ga_individual_t));
 
-  double error = fitness_buffer[0];
+  f64 error = fitness_buffer[indices[0]];
+  free(sorted);
   free(fitness_buffer);
+  free(indices);
 
   return error;
 }
@@ -113,30 +140,37 @@ static double ga_fitness_sort(ga_individual_t *population,
 //
 // return - The error of the best performing individual
 
-double ga_generation(ga_individual_t *population, uint32_t population_n,
-                     double rate_selection, double rate_mutation,
-                     ga_fitness_t fitness, void const *fitness_data,
-                     ga_crossover_t crossover, ga_mutation_t mutation) {
-  //   assert(0 < rate_selection && rate_selection <= 1.0);
-  //   assert(0 < rate_mutation && rate_mutation <= 1.0);
-  assert(2 < population_n);
+f64 ga_generation(ga_individual_t *population, u32 population_n,
+                       f64 rate_selection, f64 rate_mutation,
+                       ga_fitness_t fitness, void const *fitness_data,
+                       ga_crossover_t crossover, ga_mutation_t mutation) {
+  if (!population || !fitness || !crossover || !mutation)
+    return -1;
+  if (rate_selection <= 0 || rate_selection >= 1.0)
+    return -1;
+  if (rate_mutation < 0)
+    return -1;
+  if (population_n < 3)
+    return -1;
 
-  uint32_t selection_n = rate_selection * population_n;
+  u32 selection_n = (u32)(rate_selection * (f64)population_n);
+  if (selection_n >= population_n)
+    return -1;
 
   // Crossover
-  for (uint32_t i = selection_n; i < population_n; i++) {
-    uint32_t parent_i = rand() % selection_n;
-    uint32_t parent_j = rand() % selection_n;
+  for (u32 i = selection_n; i < population_n; i++) {
+    u32 parent_i = (u32)((u64)rand() * (u64)selection_n / ((u64)RAND_MAX + 1));
+    u32 parent_j = (u32)((u64)rand() * (u64)selection_n / ((u64)RAND_MAX + 1));
 
-    // printf("%p - %p\n", &population[i], population[i]);
     crossover(&population[i], population[parent_i], population[parent_j]);
   }
 
-  // Mutation
-  uint32_t mutation_n = rate_mutation * (population_n - selection_n);
-  for (uint32_t i = 0; i < mutation_n; i++) {
-    uint32_t mutation_i = rand() % (population_n - selection_n) + selection_n;
-    mutation(population[mutation_i]);
+  // Mutation count
+  u32 mutation_n = (u32)(rate_mutation * (f64)(population_n - selection_n));
+  for (u32 i = 0; i < mutation_n; i++) {
+    u32 range = population_n - selection_n;
+    u32 mutation_i = selection_n + (u32)((u64)rand() * (u64)range / ((u64)RAND_MAX + 1));
+    mutation(&population[mutation_i]);
   }
 
   // Fitness