Evolutionary Solver

main.c

#include <stdlib.h>
#include <assert.h>
#include <stdio.h>
#include <stdint.h>
#include <stdbool.h>
#include <math.h>
#include <float.h>
#include <time.h>

typedef uint32_t U32;
typedef uint64_t U64;
typedef float F32;
typedef bool B32;

#define MAX_VALUES 128

typedef enum {
	GoalMinimize = 0,
	GoalMaximize = 1,
} Goal;

typedef struct {
	F32 values[MAX_VALUES];
	F32 fitness;
} Individual;

typedef struct {
	F32 min;
	F32 max;
} Slider;

typedef struct {
	U32 population_size;
	U32 generation_count;
	F32 crossover_rate;
	F32 mutation_rate;
	B32 record_history;
	B32 print_progress;
} Config;

F32 squared(F32 x) {
	return x * x;
}

// Returns the smaller of the two
F32 minf(F32 lhs, F32 rhs) {
	return lhs < rhs ? lhs : rhs;
}

// Returns the larger of the two
F32 maxf(F32 lhs, F32 rhs) {
	return lhs < rhs ? rhs : lhs;
}

// Ensures val is in range [lo, hi]
F32 clampf(F32 val, F32 lo, F32 hi) {
	return maxf(lo, minf(val, hi));
}

// Returns a random float in the range: [0.0, 1.0]
F32 rand_float(void) {
	return (F32)rand() / RAND_MAX;
}

// Returns a random index value in the range [0, length)
U64 rand_idx(U64 length) {
	return (U64)floorf(rand_float() * (length - 1));
}

// (a - x)^2 + b(y - x^2)^2
F32 rosenbrock(F32 *params) {
	F32 a = 1.0f;
	F32 b = 100.0f;
	F32 x = params[0];
	F32 y = params[1];
	return squared(a - x) + b * squared(y - squared(x));
}

// x^2 + y^2
//   Goal: Minimize
//   Expect: fitness=0, params={0,0}
F32 sphere(F32* params) {
	F32 x = params[0];
	F32 y = params[1];
	return squared(x) + squared(y);
}

// 20 + x^2 + y^2 - 10(cos(2pi x) + cos(2pi y))
//   Goal: Minimize
//   Expect: fitness=0, params={0,0}
F32 rastrigin(F32* params) {
	const F32 pi = 3.14159265359;
	F32 x = params[0];
	F32 y = params[1];
	return 20 + squared(x) + squared(y) - 10 * (cos(2 * pi * x) + cos(2 * pi * y));
}

Individual compete_two(Individual lhs, Individual rhs, Goal goal) {
	if (goal == GoalMaximize) {
		return lhs.fitness < rhs.fitness ? rhs : lhs;
	}
	assert(goal == GoalMinimize);
	return lhs.fitness < rhs.fitness ? lhs : rhs;
}

Individual evo_solve(F32 (f)(F32 *params), Slider* sliders, U32 num_sliders, Goal goal, const Config *config) {
	// Sanity check
	for (U64 i = 0; i < num_sliders; i += 1) {
		assert(sliders[i].min < sliders[i].max);
	}

	// Initialize the population
	Individual* pool = (Individual*)malloc(sizeof(Individual) * config->population_size);
	for (U64 i = 0; i < config->population_size; i += 1) {
		for (U64 j = 0; j < num_sliders; j += 1) {
			F32 val = rand_float() * (sliders[j].max - sliders[j].min) + sliders[j].min;
			assert(val >= sliders[j].min && val <= sliders[j].max);
			pool[i].values[j] = val;
		}
	}

	// Strongest individual across generations
	Individual global_winner = { 0 };
	global_winner.fitness = goal == GoalMaximize ? FLT_MIN : FLT_MAX;
	Individual* winners = (Individual*)malloc(sizeof(Individual) * config->population_size);
	for (U64 gen = 0; gen < config->generation_count; gen += 1) {
		// Strongest individual within its generation
		Individual local_winner = { 0 };
		local_winner.fitness = goal == GoalMaximize ? FLT_MIN : FLT_MAX;
		for (U64 i = 0; i < config->population_size; i += 1) {
			// Host a tournament
			//   - select two random contestants
			//   - winner = MIN/MAX for function
			Individual lhs = pool[rand_idx(config->population_size)];
			Individual rhs = pool[rand_idx(config->population_size)];

			lhs.fitness = f(lhs.values);
			rhs.fitness = f(rhs.values);

			Individual winner = compete_two(lhs, rhs, goal);
			memcpy(winners + i, &winner, sizeof(Individual));
			local_winner = compete_two(local_winner, winner, goal);
		}

		// Cross-over the winners
		//   - select midpoint
		//   - child gets mom's values left of mid point
		//   - child gets dad's values right of mid point
		for (U64 i = 0; i < config->population_size; i += 2) {
			if (rand_float() < config->crossover_rate) {
				Individual mom = winners[i + 0];
				Individual dad = winners[i + 1];
				U64 pivot = rand_idx(num_sliders);
				for (U64 j = 0; j < num_sliders; j += 1) {
					winners[i + 0].values[j] = j < pivot ? mom.values[j] : dad.values[j];
					winners[i + 1].values[j] = j < pivot ? dad.values[j] : mom.values[j];
				}
			}
			
		}

		// Mutate the winners
		for (U64 i = 0; i < config->population_size; i += 1) {
			for (U64 j = 0; j < num_sliders; j += 1) {
				if (rand_float() < config->mutation_rate) {
					// nudge +/- 10%
					winners[i].values[j] *= rand_float() < 0.5 ? 0.95 : 1.05;
					winners[i].values[j] = clampf(winners[i].values[j], sliders[j].min, sliders[j].max);
				}
			}
		}
		
		// Repeat with (candidates = children)
		pool = winners;

		// NOTE: Remember to store the most fit individual we've ever seen!
		global_winner = compete_two(global_winner, local_winner, goal);

		if (config->print_progress) {
			printf("Best Candidate: fitness=%0.2f\n", global_winner.fitness);
		}
	}

	return global_winner;
}

// Convenience function
void slider_init(Slider* slider, U64 i, F32 min, F32 max) {
	slider[i].min = min;
	slider[i].max = max;
}

int main(void) {
	srand(time(NULL)); // ensure different random values each run

	U64 num_sliders = 2;
	Slider* sliders = (Slider *)malloc(sizeof(Slider) * num_sliders);
	slider_init(sliders, 0, -5.12f, +5.12f);
	slider_init(sliders, 1, -5.12f, +5.12f);

	Config config = {
		.population_size = 2000,
		.generation_count = 5000,
		.crossover_rate = 0.70f,
		.mutation_rate = 0.10f,
		.record_history = false,
		.print_progress = true,
	};
	Individual solution = evo_solve(sphere, sliders, num_sliders, GoalMinimize, &config);
	printf("Found solution (pop=%d, gen=%d):\n", config.population_size, config.generation_count);
	for (U64 i = 0; i < num_sliders; i += 1) {
		printf("  param[%d] = %0.2f [%0.2f,%0.2f]\n", i, solution.values[i], sliders[i].min, sliders[i].max);
	}
	printf("  fitness = %0.2f\n", solution.fitness);

	return 0;
}