First example of Zig code

[ndrean]

I wanted to test Zig for simple computational algorithms after going through the "ziglings" exercises.

I looked into genetic algorithms for the simplicity of these kind of algorithms. It is an optimisation algorithm: instead of a "standard" say gradient descent, you perform a sorting, an elitism and mutation of the population.

What is this challenge?

You are given a phrase (a list of characters) and the program must find it. 🤔

Why?

You can use genetic algorithms for this kind of task, and it is conceptually simple to program and is performant because this problem is 1) simple (it is a kind of sorting problem) 2) has a "global minimum" (only one possible solution).

Genetic algorithms are not silver bullet. Read the "limitations" paragraph on wikipedia.

How?

You start with a list (or "population") of random phrases (or "chromosomes"). Each phrase is made of letters (or "genes"). The length of a phrase is of course the same length as your target phrase.

You will then produce a new generation from the old one. Concretely, you modify your list of genes at each step. This is a convergente process.

Convergence measure

We need of course a metric. The measure of convergence - aka as "fitness" - is simply the number of letters ("genes") in your guess (phrases) that are equal and at the right place (in your target): if "I love you" is your target, you produce a chromosome of random genes of length 10.

Note that the implementation of a chromosome depends upon your language. In Node.js, I used a string, a concatenation of characters. With Zig, an ArrayList, a dynamic array with memory allocation.

The algorithm

• Your first generation is an initial random "population" of guesses. This is done by taking a random letter (or "gene") from a given pool of letters (all a..zA..Z0..9 plus special characters), and building a "chromosome" (a string in the Node.js implementation below) of the same size as the length of your target string

• You sort your population of "chromosomes" (strings in fact here) with your previously defined "fitness" measure.

• You run your genetic operator. You then produce a genetic admixture from this population for your next generation:

• you keep the e% best candidates, the "elitism" parameter,

• for the rest of your population, you take 45% of one candidate (or "parent") and 45% of another one and 10% random and rebuild a new candidate from it.

Node.js reference code

Node.js: Find my phrase

// Number of individuals in each generation
const POPULATION_SIZE = 100;

// Valid Genes
const GENES =
  "abcdefghijklmnopqrstuvwxyz" + "ABCDEFGHIJKLMNOPQRSTUVWXYZ " + "1234567890" + ", .-;:_!#%&/()=?@${[]}";

const GLEN = GENES.length;

// Target string to be generated
const TARGET = "I love you baby!";
const TLen = TARGET.length;

const ELISTIM = 10;

// Function to generate random numbers in given range
function RandomNum(start, end) {
  return Math.floor(Math.random() * (end - start + 1)) + start;
}

// Create random genes for mutation
function TakeRandomGene() {
  let r = RandomNum(0, GLEN - 1);
  return GENES.charAt(r);
}

// Create chromosome or string of genes
function CreateGnome() {
  let gnome = "";
  for (let i = 0; i < TLen; i++) {
    gnome += TakeRandomGene();
  }
  return gnome;
}

// Class representing individual in population
class Individual {
  constructor(chromosome) {
    this.Chromosome = chromosome;
    this.Fitness = this.CalculateFitness();
  }

  // Calculate fitness score, it is the number of
  // characters in string which differ from target string.
  CalculateFitness() {
    let fitness = 0;
    for (let i = 0; i < this.Chromosome.length; i++) {
      if (this.Chromosome[i] !== TARGET[i]) {
        fitness++;
      }
    }
    return fitness;
  }

  // Perform mating and produce new offspring:
  // take 45% of the genes from the parent 1 and
  // 45% of the genes from the parent 2, and the
  // rest from a random gene
  Mate(parent2) {
    let childChromosome = "";
    for (let i = 0; i < this.Chromosome.length; i++) {
      // get a random number between 0 and 1
      let p = Math.random();
      if (p < 0.45) childChromosome += this.Chromosome[i];
      else if (p < 0.9) childChromosome += parent2.Chromosome[i];
      else childChromosome += TakeRandomGene();
    }
    return new Individual(childChromosome);
  }
}

// Overloading < operator
class FitnessComparer {
  static Compare(ind1, ind2) {
    return ind1.Fitness - ind2.Fitness;
  }
}

// Driver code
function Main() {
  // current generation
  let generation = 0;

  let population = [];
  let found = false;

  // create initial population
  for (let i = 0; i < POPULATION_SIZE; i++) {
    let gnome = CreateGnome();
    population.push(new Individual(gnome));
  }

  while (!found) {
    // mutate and sort the population in increasing order of fitness score
    population.sort((a, b) => FitnessComparer.Compare(a, b));

    // if the individual having lowest fitness score ie.
    // 0 then we know that we have reached the target
    // and break the loop
    if (population[0].Fitness === 0) {
      found = true;
      break;
    }

    // Otherwise generate new offsprings for new generation
    let newGeneration = [];

    // Perform Elitism, that means 10% of fittest population
    // goes to the next generation
    let s = Math.floor((ELISTIM * POPULATION_SIZE) / 100);
    for (let i = 0; i < s; i++) newGeneration.push(population[i]);

    // From 50% of fittest population, Individuals
    // will mate to produce offspring
    // s = Math.floor((90 * POPULATION_SIZE) / 100);
    s = POPULATION_SIZE - s;
    for (let i = 0; i < s; i++) {
      let r = RandomNum(0, 50);
      let parent1 = population[r];
      r = RandomNum(0, 50);
      let parent2 = population[r];
      let offspring = parent1.Mate(parent2);
      newGeneration.push(offspring);
    }
    population = newGeneration;
    console.log(
      "Generation: " +
        generation +
        "\t" +
        "String: " +
        population[0].Chromosome +
        "\t" +
        "Fitness: " +
        population[0].Fitness
    );

    generation++;
  }

  console.log(
    "Generation: " +
      generation +
      "\t" +
      "String: " +
      population[0].Chromosome +
      "\t" +
      "Fitness: " +
      population[0].Fitness
  );
}

// Execute the main function
Main();

Results: plotting the influence of "elitism" on the number of runs to find the solution

For each elitism parameter, I ran 1000 times the algorithm, so for each elitism parameter, I have 1000 number of generations needed to find the guess. From this, I can extract a mean and a standard deviation and the min and the Max.

The axis is the "elitism" factor: how much of the best population do you keep at each generation.

The upper points are the mean. In average, the minimum runs is 80 at approximatively 20% of elitism. Then it is surprising to see that beyond 25% of elitism, it goes quickly up.
The middle points are the standard deviation. Above 25% of elitism, it is really increasing. You see that the variance is 50%, (stddev/avg) which is high.
The lower curve is the minimum generations needed to find the target phrase. It is rather flat between 0 and 40%.

An Elixir Livebook to explore a CSV file produced by the Zig code

This graph has been produced by importing the CSV file into a Livebook, and using Explorer to convert the CSV into a dataframe, and using Vegalite to plot the data.

The Livebook code

The Notebook code

# ExPlore your Genetic Algorithm

```elixir
Mix.install([
  {:kino_explorer, "~> 0.1.20"},
  {:explorer, "~> 0.9.2"},
  {:kino_vega_lite, "~> 0.1.11"}
])

Section

require Explorer.DataFrame, as: DF

{:ok, data} = data = DF.from_csv("/Users/nevendrean/code/zig/exercises/exercises/genetic_algo/genetic_algorithm_results.csv")

Elitism,Avg,StdDev,Min,Max
2,101.03,51.36,25,400
3,96.89,52.05,22,466
4,98.15,54.23,23,499
5,93.62,49.65,23,441
...

data

VegaLite.new(title: "Influence of Elitism", width: "500", height: "500")
|> VegaLite.data_from_values(data, only: ["Elitism", "Avg", "StdDev", "Min"])
|> VegaLite.layers([
  VegaLite.new()
  |> VegaLite.mark(:point)
  |> VegaLite.encode_field(:x, "Elitism", type: :quantitative)
  |> VegaLite.encode_field(:y, "Avg", type: :quantitative)
  |> VegaLite.encode_field(:color, "Avg", type: :quantitative),
  VegaLite.new()
  |> VegaLite.mark(:point)
  |> VegaLite.encode_field(:x, "Elitism", type: :quantitative)
  |> VegaLite.encode_field(:y, "StdDev", type: :quantitative)
  |> VegaLite.encode_field(:color, "Min", type: :quantitative),
  VegaLite.new()
  |> VegaLite.mark(:line)
  |> VegaLite.encode_field(:x, "Elitism", type: :quantitative)
  |> VegaLite.encode_field(:y, "Min", type: :quantitative)
])

The Zig code

I reproduced this code in Zig but sophisticated it a bit to get the data aboveo:

• run the algorithm say X times and compute the mean, standard deviation, min and Max.
• and then I increment the "elitism" parameter,
• and finally save this into a CSV file
• so that I can use Explorer in a Livebook and chart the data.
  The Zig code. It takes quite a long to run (with the parameters below) and consumes at peak 2MB of memory (really low). This is because I write to disk and output to stdout.

it uses a number of hardcoded parameters. In particular, there is no check on them; for example, I used GENE_MIX=0.45 as the amount of genes a chromosome gives up to the next generation, but anything above 0.50 will make the program crash.

const std = @import("std");
const math = @import("std").math;
const print = std.debug.print;
const ArrayList = std.ArrayList;
const Random = std.Random;
const ArenaAllocator = std.heap.ArenaAllocator;

const POPULATION_SIZE: comptime_int = 100;
const GENES = "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ 1234567890, .-;:_!\"#%&/()=?@${[]}";
const MAX_TARGET_LENGTH: comptime_int = 100;
const MAX_GENERATIONS: comptime_int = 10_000; // Prevent infinite loops
const MAX_TRIALS: comptime_int = 2_000;
const MAX_ELITISM: comptime_int = 75;
const ELITISM_STEP: comptime_int = 1;
const ELITISM_START: comptime_int = 2;
const GENE_MIX: comptime_float = 0.45;

const TARGET = "I love you baby!";

const Individual = struct {
    string: []u8,
    fitness: usize,

    fn init(allocator: std.mem.Allocator, string: []const u8) !Individual {
        var self = Individual{
            .string = try allocator.dupe(u8, string),
            .fitness = 0,
        };
        self.calculateFitness();
        return self;
    }

    // mutate the Individual struct and compute the fitness:
    // it is the number of differences between the genes and the target
    fn calculateFitness(self: *Individual) void {
        self.fitness = 0;
        for (self.string, 0..) |char, i| {
            if (char != TARGET[i]) {
                self.fitness += 1;
            }
        }
    }

    fn mate(self: *const Individual, allocator: std.mem.Allocator, parent2: *const Individual, rng: Random) !Individual {
        const child_string = try allocator.alloc(u8, TARGET.len);
        for (child_string, 0..) |*gene, i| {
            const p = rng.float(f32);
            if (p < GENE_MIX) {
                gene.* = self.string[i];
            } else if (p < GENE_MIX*2) {
                gene.* = parent2.string[i];
            } else {
                gene.* = GENES[rng.intRangeAtMost(usize, 0, GENES.len - 1)];
            }
        }
        return try Individual.init(allocator, child_string);
    }
};

fn createGnome(allocator: std.mem.Allocator, ran: Random) ![]u8 {
    const gnome = try allocator.alloc(u8, TARGET.len);
    for (gnome) |*gene| {
        gene.* = GENES[ran.intRangeAtMost(usize, 0, GENES.len - 1)];
    }
    return gnome;
}

fn individualLessThan(context: void, a: Individual, b: Individual) bool {
    _ = context;
    return a.fitness < b.fitness;
}

fn runGeneticAlgorithm(allocator: std.mem.Allocator, elitism: usize, ran: Random) !usize {
    var population = ArrayList(Individual).init(allocator);
    defer population.deinit();

    // Create initial population
    for (0..POPULATION_SIZE) |_| {
        const gnome = try createGnome(allocator, ran);
        const individual = try Individual.init(allocator, gnome);
        try population.append(individual);
    }

    var generation: usize = 0;
    while (generation < MAX_GENERATIONS) : (generation += 1) {
        // Sort population by fitness
        std.mem.sort(Individual, population.items, {}, individualLessThan);

        // Check if we've found the target
        if (population.items[0].fitness == 0) {
            return generation;
        }

        // Generate new population
        var new_generation = ArrayList(Individual).init(allocator);

        // Elitism
        const elitism_count = @as(usize, (elitism * POPULATION_SIZE) / 100);
        for (population.items[0..elitism_count]) |individual| {
            try new_generation.append(try Individual.init(allocator, individual.string));
        }

        // Mating
        const mating_count = POPULATION_SIZE - elitism_count;
        for (0..mating_count) |_| {
            const parent1 = &population.items[ran.intRangeAtMost(usize, 0, 50)];
            const parent2 = &population.items[ran.intRangeAtMost(usize, 0, 50)];
            const offspring = try parent1.mate(allocator, parent2, ran);
            try new_generation.append(offspring);
        }

        // Replace old population
        population.deinit();
        population = new_generation;
    }

    return MAX_GENERATIONS; // If solution not found
}

const RunningStats = struct {
    count: usize,
    mean: f64,
    m2: f64,
    min: f64,
    max: f64,

    fn init() RunningStats {
        return RunningStats{
            .count = 0,
            .mean = 0,
            .m2 = 0,
            .min = 0,
            .max = 0,
        };
    }

    // mutate the Individual struct
    fn update(self: *RunningStats, value: usize) void {
        self.count += 1;
        const delta = @as(f64, @floatFromInt(value)) - self.mean;
        self.mean += delta / @as(f64, @floatFromInt(self.count));
        const delta2 = @as(f64, @floatFromInt(value)) - self.mean;
        self.m2 += delta * delta2;
        self.mixMax(value);
    }

    fn getStandardDeviation(self: *const RunningStats) f64 {
        if (self.count < 2) return 0;
        return std.math.sqrt(self.m2 / @as(f64, @floatFromInt(self.count - 1)));
    }

    fn mixMax(self: *RunningStats, value: usize) void {
        const v = @as(f64, @floatFromInt(value));
        if (self.min == 0) self.min = v;
        if (v > self.max) {
            self.max = v;
        }
        if (v < self.min) {
            self.min = v;
        }
    }
};

pub fn main() !void {
    var arena = ArenaAllocator.init(std.heap.page_allocator);
    defer arena.deinit();
    const allocator = arena.allocator();

    var rand = Random.DefaultPrng.init(@intCast(std.time.timestamp()));
    const ran = rand.random();

    const file = try std.fs.cwd().createFile("genetic_algorithm_results2.csv", .{});
    defer file.close();
    const writer = file.writer();

    // Write CSV header
    try writer.writeAll("Elitism,Avg,StdDev,Min,Max\n");

    var elitism: usize = ELITISM_START;

    while (elitism <= MAX_ELITISM) : (elitism += ELITISM_STEP) {
        var stats = RunningStats.init();

        for (0..MAX_TRIALS) |_| {
            const generations = try runGeneticAlgorithm(allocator, elitism, ran);
            stats.update(generations);

            // Free memory after each trial
            arena.deinit();
            arena = ArenaAllocator.init(std.heap.page_allocator);
        }

        // Write summary statistics to CSV
        try writer.print("{d},{d:.2},{d:.2},{d},{d}\n", .{ elitism, stats.mean, stats.getStandardDeviation(), stats.min, stats.max });

        // Print to console for immediate feedback
        std.debug.print("Elitism: {d}, Avg: {d:.2}, StdDev: {d:.2}, min: {d}, max: {d}\n", .{ elitism, stats.mean, stats.getStandardDeviation(), stats.min, stats.max });
    }

    // Ensure all data is written to the file
    try file.sync();

    std.debug.print("The results have been saved to 'genetic_algorithm_results.csv'\n", .{});
}

First conclusion:
I used a lot the ziglings exercises, this helped a lot, but it was quite challenging (for me). I had plenty of errors. In particular, you have to be careful with the types (ints vs floats and cast quite a few times as functions take ints or floats), you must be careful between runtime code and compiled time code, manage the memory.
I also cheated and I used Claude.ai: how to export a CSV file, how to generate random numbers in Zig, how to sort a list with the sorting function.

Next step: use Elixir and its easy parallelism to check if I can speed up the computations.

[ndrean]

I just discovered GA in Elixir by S Moriarity. Of course, parallelism is a good choice. Good news.

Although the Zig code above explores more deeply how the parameters influences the results, some shorter Elixir code:

• a gist from Sean Moriarity: https://gist.github.com/seanmor5/feff05877f3c5cb04b539bf24dda0c96
  and the companion blog he made

• another application: https://github.com/b0a04gl/np-hard-solver-elixir/tree/master