Optimizing Pointless things

So, we’ve been playing rimworld a lot recently. One of the mods we use is called Rim Factory, which is quite nice overall for building automated bases(yes this breaks the game’s economy, what of it?). One of the mechanics it adds is a buildable object called a sprinkler, which, instead of watering your plants(this isn’t necessary in vanilla rimworld), it causes planets within a certain radius(the second tier is within twelve tiles), to grow by two hours. This is, completely unbalanced and, when you get down to it, just another item in the pursuit of excess, but that’s how we like to play the game.

Anyway, we were wondering about figuring out an optimal pattern for the usage, and since the more rigorous tools we have to measure such things are largely forgotten, we built a little simulation script in python to try and get a feel for various optimization strategies, and the most fun part, of course, is how we went about building this kind of tool.

This tool is remarkably sketchy, it does not attempt to be rigorous, or even follow the watering pattern other than cycling through all of the tiles in range.

Since we’ve been building stuff with python lately, this is in python, which, as always, never ceases to amaze how easy the ecosystem makes it to do this kind of task.

Anyway, to get an idea of where this is going, let’s start with the import block at the top:

from typing import *
import numpy as np
import itertools
import functools
import math
import statistics
import matplotlib
import matplotlib.pyplot as plt
import random
from enum import IntEnum, auto

The building blocks

The core of the script is the iterators that provide the information about where the focus is, and since this is a script without much intention of going particularly far, we can use globals and constants, so let’s start with those!

# The size of the map being kept track of.
#
# In this case, this means that the map consists of the area from (-24,-24) to (24,24)
MAP_SIZE = 24
# The radius of a sprinkler mkii is 12 blocks
RADIUS = 12

Iterators

The heart and soul of this program is in its iterators, so let’s start by defining one that will make the rest of our work easier.

def cartesian(start: int, end: int) -> Iterator[Tuple[int,int]]:
    """
    Turn a single range iterator as in range(start,end) into a cartesian product
    """
    return itertools.product(range(start,end),repeat=2)

That’ll be useful for iterating over the map. The next are a bit less sensible, but, they work, and this isn’t meant to be perfect.

def map_tiles():
    r = RADIUS+1
    z = itertools.product(range(-MAP_SIZE,MAP_SIZE), repeat=2)
    z = list(z)
    random.shuffle(z)
    for x in z:
        yield x

def sprinkler(position:Tuple[int,int]) -> Iterator[Tuple[int,int]]:
    for i, j in cartesian(0, RADIUS+1):
        if i==0 and j == 0:
            continue
        s = i**2+j**2
        if math.sqrt(s) <= RADIUS:
            yield (i+position[0], j+position[1])
            yield (position[0]+i,position[1]-j)
            yield (position[0]-i,j+position[1])
            yield (position[0]-i,position[1]-j)

Fitness

Next comes the fitness function, which, to start off with, will just be a simple mean of how many times each tile has been hit by a sprinkler.

def evaluate_fitness(field: Dict[Tuple[int,int], int]) -> float:
    return sum(field.values()) / (2*MAP_SIZE)**2

Each round in the simulation consists of iterating across a given number of sprinklers n times, and is called for each candidate to figure out the best one.

Putting that together

def simulate(positions: Iterable[Tuple[int,int]], rounds: int)\
             -> Dict[Tuple[int,int],int]:
    field = {(x,y):0 for x,y in cartesian(-MAP_SIZE, MAP_SIZE)
             if (x,y) not in positions}
    for i in positions:
        if i in field:
            del field[i]
    sprinklers = list(map(lambda x: itertools.cycle(sprinkler(x)), positions))
    for _ in range(rounds):
        for s in sprinklers:
            pos = next(s)
            if pos in field:
                field[pos] += 1
    return field

So now that we can evaluate each choice we might make, we have to make it possible to put these together for any number of sprinklers, so we write another function called evaluate_alternatives. Here we decided to limit the breadth of each search step to a given number, here it’s called to_try.

def evaluate_alternatives(positions: Iterable[Tuple[int,int]],
                          to_try: int,
                          rounds: int) -> List[Tuple[int,int]]:
    # Eliminate already extant positions
    options = set(map_tiles()) - set(positions)
    options = list(options)
    random.shuffle(options)
    to_try = min(to_try, len(options)-1)
    options = options[0:to_try]
    best_fitness = -2**48
    best_found = None
    for i in options:
        f = simulate(positions + [i], rounds)
        fitness = evaluate_fitness(f)
        if fitness > best_fitness:
            best_found = i
            best_fitness = fitness
    print(f"best fitness found this round {best_fitness}")
    return positions + [best_found]

Next we put together the final pieces in a function we called optimize.

def optimize(n: int, candidates:int, rounds: int):
    positions = []
    while len(positions) < n:
        positions = evaluate_alternatives(positions, candidates, rounds)
    fields = simulate(positions, rounds)
    heatmap = np.array([[(fields[x,y] if (x,y) in fields else -5)
                             for x in range(-MAP_SIZE, MAP_SIZE)]
                        for y in range(-MAP_SIZE,MAP_SIZE)])
    fig, ax = plt.subplots()
    im = ax.imshow(heatmap)
    # You can uncomment these if you want to have a label on each cell in the heatmap,
    # for 48x48, it is rather overwhelming
##    for (x,y), v in fields.items():
##        ax.text(x+MAP_SIZE,y+MAP_SIZE, math.floor(fields[x,y]), ha='center', va='center', color='w')
    fig.tight_layout()
    plt.show()

and to top it off, we invoke it at the end of the file:

optimize(10, 80, 1500)

The output should then look something like this

A Heatmap

Okay, that’s fine, but it’s not necessarily what you might want from such a tool. What if you’d rather have as much coverage as possible with a given number of sprinklers?

Well, we can do that.

Alternative fitness metrics

Okay, so what we need to optimize for to get the best coverage possible is to simply count the fields which are touched. In all likelihood, we might end up wanting further metrics to measure by, so we might as well make it a bit… well… Configurable, if you’re going to use a charitable term. So let’s introduce an Enum right before the fitness function.

class FitnessMetric(IntEnum):
    MEAN = auto()
    COVERAGE = auto()

fitness_metric = FitnessMetric.COVERAGE

Then we need to modify the evaluate_fitness function.

def evaluate_fitness(field: Dict[Tuple[int,int],int]) -> float:
    if fitness_metric == FitnessMetric.COVERAGE:
        r = len(list(filter(lambda x:x>0, field.values())))
    elif fitness_metric == FitnessMetric.MEAN:
        r = sum(field.values())/(2*MAP_SIZE)**2
    return r

Let’s run the script again and see what falls out.

Coverage That certainly seems to have the right effect. The majority of the map is covered by the sprinklers.

Complications

So, if you’re familiar with this mod, you might realize there are some things this isn’t even trying to account for, namely any infrastructure. Well, that should be easy to add. 🤞

Let’s add a set of positions that we don’t want filled. For now, let’s leave it to the center(we’re assuming that maybe you have a harvester station set up there), so let’s add another constant in front of evaluate_alternatives (yeah, this is getting pretty messy).

excluded_positions = frozenset([(0,0)])

Then we just have to modify evaluate_alternatives to avoid these positions.

def evaluate_alternatives(positions: Iterable[Tuple[int,int]],
                          to_try: int,
                          rounds: int) -> List[Tuple[int,int]]:
    # Eliminate already extant positions
    options = set(map_tiles()) - set(positions) - excluded_positions
#    --- SNIP ---
    heatmap = np.array([[(fields[x,y] if (x,y) in fields and (x,y) not in excluded_positions else -5)
                         for x in range(-MAP_SIZE, MAP_SIZE)]
                    for y in range(-MAP_SIZE,MAP_SIZE)])
#    --- SNIP ---

Alright, what does that make it do?

Alright, that’s a further improvement, we suppose, if nothing else it permits more flexibility in terms of describing a given scenario you want to optimize for.

Okay, so let’s say that you have a specific area that you want to optimize coverage for. How would we go about adding that?

Well, at this point it’s pretty clear that all this nonsense with enums is going to be an increasingly weak abstraction for dealing with additional fitness functions, but let’s keep it for now, if only because we’re not sure how we want to go about fixing this ugly feeling design as of writing this very paragraph.

Okay, so let’s assume that you want to optimize for an automatic harvester with an area between (-5,-5) and (5,5).

special_region = set(cartesian(-5,5))

And we can alter evaluate_fitness to this:

def evaluate_fitness(field: Dict[Tuple[int,int], int]) -> float:
    r=0
    if fitness_metric == Fitness.MEAN:
        r = sum(field.values())/(MAP_SIZE*2)**2
    elif fitness_metric == Fitness.COVERAGE:
        r = len(list(filter(lambda x:x>0, field.values())))
    elif fitness_metric == Fitness.SPECIAL_REGION:
        r = sum(map(lambda x:x[1], filter(lambda x: x[0] in special_region,field.items())))
    return r

We end up with output like this:

You might not like this output because, among other things, it doesn’t penalize the sprinklers for showing up inside of our growing zone. So, we need to add another line to our much suffering evaluate_fitness function.

def evaluate_fitness(field: Dict[Tuple[int,int], int]) -> float:
    r=0
    if fitness_metric == Fitness.MEAN:
        r = sum(field.values())/(MAP_SIZE*2)**2
    elif fitness_metric == Fitness.COVERAGE:
        r = len(list(filter(lambda x:x>0, field.values())))
    elif fitness_metric == Fitness.SPECIAL_REGION:
        r = sum(map(lambda x:x[1], filter(lambda x: x[0] in special_region,field.items())))
        r /= max(1, len([(x,y) for x,y in special_region
                         if (x,y) not in field.keys()
                         if (x,y) not in excluded_positions]))
    return r

And this helps a bit: Penalizing the region Let’s label a similar run to see if it has excluded the region altogether, if a square ends up labelled ‘Se’ it indicates that it’s a sprinkler in an excluded area.

Well, damn. Our code has a mistake in it that permitted the single offending sprinkler. We can fix this by changing the function, to this:

def evaluate_fitness(field: Dict[Tuple[int,int], int]) -> float:
    r=0
    if fitness_metric == Fitness.MEAN:
        r = sum(field.values())/(MAP_SIZE*2)**2
    elif fitness_metric == Fitness.COVERAGE:
        r = len(list(filter(lambda x:x>0, field.values())))
    elif fitness_metric == Fitness.SPECIAL_REGION:
        r = sum(map(lambda x:x[1], filter(lambda x: x[0] in special_region,field.items())))
        penalty_squares = len([(x,y) for x,y in special_region
                              if (x,y) not in field.keys()
                              if (x,y) not in excluded_positions])
        
        r/= penalty_squares+1
    return r

Alright, and this fixes that issue(and with another fix applied it labels the spots with sprinklers correctly as ‘S’).

Summary

Okay, so what we want you to be able to take away from this post:

  1. When you’re processing data for something pointless, don’t bother being clean if doing it in a dirty way lets you figure out what matters faster
  2. Separate it into concerns that model, however roughly, the components that matter
  3. Build a metric for evaluating how good a result is
  4. Descend along randomly chosen positions based on how optimal they are

Feel free to ask questions or comment. The full, considerably messier code can be found here.

Trying Out Zig

So, we have heard good things about Zig. These boil down to the following things:

  • Good speed
  • Fast compilation
  • Decent type system
  • Simple Syntax

So far, a lot of these things seem to be born out by the experience we've had, though, we have some criticism that is probably more along the lines of our taste in programming, rather than issues we expect to be universal. For fairness sake, let's start with things we like, once again expressing our preferences rather than universals.

  • Low Level
    Low level languages feel liberating to us, because we get to make the decisions, rather than the compiler or virtual machine.
  • Trusts you, mostly.
    This is mostly in comparison to Rust, whose borrow checker we have wrestled with a great deal and whose Turing complete type system has, in the past, left us mystified when we have run into errors.
    And if you really, really want to throw caution to the wind, you can mostly just tell zig that you don't give a damn about whether or not the pointer is to a single element or an array.
  • Compile Time computation, Generics
    This is how generics are implemented, and for the most part, it provides a good experience, especially since it compares exceedingly favorably with pre-c++17 C++.
  • Not afraid of Bare Pointers
    Bare pointers aren't something to worry about too much, especially if they're typed. Bare pointers are quite a bit faster and lighter than other (tagged) references.

And honestly, we're impressed at how tight the language feels. There are some facilities we're not sure about, and some that we wish it had, but overall, this isn't a bad language. We like it more than Rust(which, admittedly, really isn't saying much for us), though, from familiarity alone we're likely to revert to using C or C++, but we're not nearly as skeptical as we once were about this little language.

The syntax is familiar enough, pretty similar to C likes, though, with a few elements that we can't really attribute to any language we know(which doesn't mean very much). So, let's see what something simple looks like.

const std = @import("std");
pub fn main() !void {
    const stdout = std.io.getStdOut().writer();
    var i: usize = 1;
    while (i < 10) {
        try stdout.print("{}\n", .{i});
        i += 1;
    }
}

You can see that Zig is careful about errors, requiring that you either mark that main can return an error (designated with !void), or handle the error somehow. Anything that can fail in Zig has a type denoted by !<type>, which requires you to think about what errors can happen, or so communities as Go users and Rust users insist.

Let's see something that differs even more significantly than C.

const std = @import("std");
fn arraylist_user(alloc: *std.mem.Allocator) !std.ArrayList(i32) {
    var al = std.ArrayList(i32).init(alloc);
    var i: i32 = 0;
    while (i < 100) {
        var item = try al.addOne();
        item.* = i;
        i += 1;
    }
    return al;
}
pub fn main() !void {
    var gpalloc = std.heap.GeneralPurposeAllocator(.{}){};
    defer std.debug.assert(!gpalloc.deinit());
    const gpa = &gpalloc.allocator;
    const stdout = std.io.getStdOut().writer();
    var array_list = try arraylist_user(gpa);
    defer array_list.deinit();
    for (array_list.items) |*item| {
        try stdout.print("{}\n", .{item.*});
    }
}

So, in this little snippet we have used several new things, Generic types, Allocator, and the defer keyword(which go users will immediately recognize). Zig does not have a default allocator, and it wants you to pass the allocator you want to use to the function that allocates memory. This is rather different than C, where you would probably just use malloc/realloc/calloc/free to manage memory(Or define a macro that evaluates to that by default if you're building that kind of library). The reason that the documentation has an assert in gpalloc.deinit() is because this particular allocator can track memory leaks, so this causes it to report such occurrences. It also shows one of the syntactic divergences from C, for(array_list.items) |*item|. Pipes don't get much use in C other than as logical operators, but here it tells the loop to unwrap the value(a pointer to a position in the list), and calls it item. In Zig, you dereference this pointer with item.*. Another point of comparison to C++, the addOne method doesn't add the item itself, but returns a pointer to it, so that you can assign it afterwards.

Interestingly, for such a low level language, array_list.items isn't a pointer, it's a slice, which is a pointer+length, so in C++, we would say that it has most in common with a string_view.

Okay, so what if we want to go further? What if we want to do a reprisal of our C N-body simulator written with Raylib? In fact, that's not too bad. It's almost even easy. In fact, we can even use raylib without too much trouble.

To start off, we need to import raylib. Luckily, as it turns out, Zig's compiler can handle C interop quite nicely, it even has builtins for it. (Though, if you do this then you need to call the compiler with -lc -lraylib to make sure that it has all the c libraries linked into it).

const std = @import("std");
const ray = @cImport({
    @cInclude("raylib.h");
});

Okay, so the first thing we should do here is define a vector type, not the resizable array type of vector, but the directional vector. Since we prize reuse here, we're going to make it generic even though we probably don't really need to.

pub fn vector2(comptime value: type) type {
    return struct {
        const This = @This();
        x: value,
        y: value,
        //Write methods here!
    };
}

Okay, so that's an easy enough definition. But let's add some vector math oriented methods to make using it easier and more pleasant.

    return struct {
        const This = @This();
        x: value,
        y: value,
        pub fn add(this: This, other: This) This {
            var t: This = v2{ .x = 0, .y = 0 };
            t = .{ .x = this.x + other.x, .y = this.y + other.y };
            return t;
        }
        pub fn sub(this: This, other: This) This {
            return this.add(.{ .x = -other.x, .y = -other.y });
        }
        pub fn scale(this: This, v: value) This {
            return .{ .x = this.x * v, .y = this.y * v };
        }
        pub fn distance(this: This, other: This) value {
            const sqrt = std.math.sqrt;
            const pow = std.math.pow;
            return sqrt(pow(value, this.x - other.y, 2) + pow(value, this.y - other.y, 2));
        }
        //Wrap the value around, like when pacman reaches the edge of the map
        pub fn wrap(c: This, a: This, b: This) This {
            var r: This = .{ .x = c.x, .y = c.y };
            if (r.x < a.x) {
                r.x = b.x;
            } else if (r.x > b.x) {
                r.x = a.x;
            }
            if (r.y < a.y) {
                r.y = b.y;
            } else if (r.y > b.y) {
                r.y = a.y;
            }
            return r;
        }
        pub fn magnitude(a: This) value {
            const sqrt = std.math.sqrt;
            return sqrt(a.x * a.x + a.y * a.y);
        }
        pub fn normalize(a: This) This {
            return a.scale(1.0 / a.magnitude());
        }
    };
}

Alright, the main things to note here are how we need to call std.math.pow, namely that we need to call it with a compile time type, in this case value. Later on we'll see it called with f32.

Now we need to define the type we use for particles, and while we're at it, we're going to make a shortcut to the kind of vector we're using here.

const v2 = vector2(f32);
const particle = struct {
    const This = @This();
    position: vector2(f32) = v2{ .x = 0, .y = 0 },
    velocity: vector2(f32) = v2{ .x = 0.0, .y = 0.0 },
    acceleration: vector2(f32) = v2{ .x = 0.0, .y = 0.0 },
    mass: f32,
    //Methods here!
};

We also need a radius property, but since it's derived from the mass, and later on that can change as the bodies absorb each other, so it needs to be a method. We should also write methods to determine if the particles overlap and just moving the position along by the velocity, as well as calculating the attraction between two particles.

const particle = struct {
    const This = @This();
    position: vector2(f32) = v2{ .x = 0, .y = 0 },
    velocity: vector2(f32) = v2{ .x = 0.0, .y = 0.0 },
    acceleration: vector2(f32) = v2{ .x = 0.0, .y = 0.0 },
    mass: f32,
    pub fn radius(this: *const This) f32 {
        var result: f32 = std.math.ln(this.mass) / std.math.ln(3.0);
        if (result > 40) {
            return 40 - std.math.ln(this.mass) / std.math.ln(5.0);
        }
        return result;
    }
    //Returns true if the two particles overlap
    pub fn overlaps(this: *const This, other: *const this) bool {
        var r1 = this.radius();
        var r2 = other.radius();
        var dist = this.position.distance(other.position);
        return (r1 + r2) > dist;
    }
    //Handles the base movement
    pub fn motion_step(this: *This, timestep: f32) void {
        this.position = this.position.add(this.velocity.scale(timestep));
        this.velocity = this.velocity.add(this.acceleration.scale(timestep));
    }
    pub fn attraction(this: *const This, other: *const This, g: f32) vector2(f32) {
        var dist = this.position.distance(other.position);
        var vector_to = other.position.sub(this.position).normalize();
        return vector_to.scale(g * (this.mass * other.mass) / std.math.pow(f32, dist, 2));
    }
};

Now we should build a container for many particles, and the properties necessary to simulate them that aren't appropriate to put in the particle structures themselves. We need a method to initialize the particles to random positions, a method to handle the gravity simulation and such, and a method to actually draw the particles.

const ParticleCollection = struct {
    const This = @This();
    particles: [100]particle,
    window_start: vector2(f32) = v2{ .x = 0.0, .y = 0.0 },
    window_end: vector2(f32) = v2{ .x = 100, .y = 100 },
    timestep: f32 = 0.01,
    gravitational_constant: f32 = 1e-3,
    pub fn init_particles(this: *This, rand: *std.rand.Random) void {
        for (this.particles) |*p| {
            p.mass = @intToFloat(f32, rand.intRangeLessThan(i32, 1, 100));
            p.position = .{
                .x = @intToFloat(f32, rand.intRangeLessThan(i32, @floatToInt(i32, this.window_start.x), @floatToInt(i32, this.window_end.x))),
                .y = @intToFloat(f32, rand.intRangeLessThan(i32, @floatToInt(i32, this.window_start.y), @floatToInt(i32, this.window_end.y))),
            };
            p.acceleration = .{ .x = 0.0, .y = 0.0 };
            p.velocity = .{ .x = 0.0, .y = 0.0 };
        }
    }
    pub fn step_world(this: *This) void {
        for (this.particles) |*p| {
            p.motion_step(this.timestep);
            p.position = p.position.wrap(this.window_start, this.window_end);
        }
        for (this.particles) |*a| {
            a.acceleration = .{ .x = 0.0, .y = 0.0 };
            for (this.particles) |*b| {
                //No self attraction please, allowing that would result in division by zero
                if (a == b)
                    continue;
                a.acceleration = a.acceleration.add(a.attraction(b, this.gravitational_constant));
            }
        }
    }
    pub fn drawSystem(this: *const This) void {
        for (this.particles) |p| {
            ray.DrawCircle(@floatToInt(c_int, p.position.x), @floatToInt(c_int, p.position.y), p.radius(), ray.BLACK);
        }
    }
};

Note how we have to pass a random number generator into init_particles, this is inline with how Zig also requires that you pass the allocators into functions that require memory allocations to be done. You also see some of the somewhat jagged interaction between Zig and C, namely that Zig doesn't specify that its i32 type is equivalent to C's int type(which on many architectures it might not be), it also requires explicit conversions between floating point numbers and integers.

The main function here is the simplest part yet.

pub fn main() !void {
    const width = 800;
    const height = 450;
    ray.InitWindow(width, height, "Nbody");
    ray.SetTargetFPS(60);
    //This is very much *not* good practice, but it's the easiest way to start this
    var rand = std.rand.Xoroshiro128.init(0);
    //Don't initialize the particles yet.
    var p: ParticleCollection = .{ .particles = undefined };
    p.window_end = .{ .x = width, .y = height };
    p.init_particles(&rand.random);
    while (!ray.WindowShouldClose()) {
        p.step_world();
        ray.BeginDrawing();
        defer ray.EndDrawing();
        ray.ClearBackground(ray.RAYWHITE);
        p.drawSystem();
    }
}

And so we have a working prototype for an nbody simulator, considerably shorter than the C version of the same program.

Interestingly, it appears to be smaller compiled than the original program in C, with just zig build-exe nbody.zig -lc -lraylib, we get an executable of around 784Kb. With zig build-exe nbody.zig -lc -lraylib -OReleaseFast, we can get it down to 92Kb, and with the -OReleaseSmall option, we can get down to 84Kb.

All in all, we'd definitely watch Zig carefully, it's very well thought out, and if they get build their package manager right, then their ecosystem might become competitive with Rust's quite quickly. The langauge already quite nice to use, and it might not a bad choice for your next project you might consider doing in C or Rust if you're looking for a new language to mess around in.