Cogs and Levers A blog full of technical stuff

Backpropagation from Scratch

Introduction

One of the most powerful ideas behind deep learning is backpropagation—the algorithm that lets a neural network learn from its mistakes. But while modern tools like PyTorch and TensorFlow make it easy to use backprop, they also hide the magic.

In this post, we’ll strip things down to the fundamentals and implement a neural network from scratch in NumPy to solve the XOR problem.

Along the way, we’ll dig into what backprop really is, how it works, and why it matters.

What Is Backpropagation?

Backpropagation is a method for computing how to adjust the weights—the tunable parameters of a neural network—so that it improves its predictions. It does this by minimizing a loss function, which measures how far off the network’s outputs are from the correct answers. To do that, it calculates gradients, which tell us how much each weight contributes to the overall error and how to adjust it to reduce that error.

Think of it like this:

  • In calculus, we use derivatives to understand how one variable changes with respect to another.
  • In neural networks, we want to know: How much does this weight affect the final error?
  • Enter the chain rule—a calculus technique that lets us break down complex derivatives into manageable parts.

The Chain Rule

Mathematically, if

\[z = f(g(x))\]

then:

\[\frac{dz}{dx} = \frac{df}{dg} \cdot \frac{dg}{dx}\]

Backpropagation applies the chain rule across all the layers in a network, allowing us to efficiently compute the gradient of the loss function for every weight.

Neural Network Flow

graph TD A[Input Layer] --> B[Hidden Layer] B --> C[Output Layer] C --> D[Loss Function] D -->|Backpropagate| C C -->|Backpropagate| B B -->|Backpropagate| A

We push inputs forward through the network to get predictions (forward pass), then pull error gradients backward to adjust the weights (backward pass).

Solving XOR with a Neural Network

The XOR problem is a classic test for neural networks. It looks like this:

Input Output
[0, 0] 0
[0, 1] 1
[1, 0] 1
[1, 1] 0

A simple linear model can’t solve XOR because it’s not linearly separable. But with a small neural network—just one hidden layer—we can crack it.

We’ll walk through our implementation step by step.

Activation Functions

import numpy as np

def sigmoid(x):
    return 1 / (1 + np.exp(-x))

def sigmoid_derivative(x):
    return x * (1 - x)

def mse_loss(y_true, y_pred):
    return np.mean((y_true - y_pred) ** 2)

We’re using the sigmoid function for both hidden and output layers.

The sigmoid activation function is defined as:

\[\sigma(x) = \frac{1}{1 + e^{-x}}\]

Its smooth curve is perfect for computing gradients.

Its derivative, used during backpropagation, is:

\[\sigma'(x) = \sigma(x) \cdot (1 - \sigma(x))\]

The mse_loss function computes the mean squared error between the network’s predictions and the known correct values (y).

Mathematically, the mean squared error is given by:

\[\text{MSE}(y, \hat{y}) = \frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2\]

Where:

  • \(y_i\) is the actual target value (y_true),
  • \(\hat{y}_i\) is the network’s predicted output (y_pred),
  • \(n\) is the number of training samples.

Data and Network Setup

X = np.array([
    [0, 0],
    [0, 1],
    [1, 0],
    [1, 1]
])

y = np.array([
    [0],
    [1],
    [1],
    [0]
])

The x matrix defines all of our inputs. You can see these as the bit pairs that you’d normally pass through an xor operation. The y matrix then defines the “well known” outputs.

np.random.seed(42)
input_size = 2
hidden_size = 2
output_size = 1
learning_rate = 0.1

The input_size is the number of input features. We have two values going in as an input here.

The hidden_size is the number of “neurons” in the hidden layer. Hidden layers are where the network transforms input into internal features. XOR requires non-linear transformation, so at least one hidden layer is essential. Setting this to 2 keeps the network small, but expressive enough to learn XOR.

output_size is the number of output neurons. XOR is a binary classification problem so we only need a single output.

Finally, learning_rate controls how fast the network learns. This value scales the size of the weight updates during training. By increasing this value, we get the network to learn faster but we risk overshooting optimal values. Lower values are safer, but slower.

W1 = np.random.randn(input_size, hidden_size)
b1 = np.zeros((1, hidden_size))
W2 = np.random.randn(hidden_size, output_size)
b2 = np.zeros((1, output_size))

We initialize weights randomly and biases to zero. The small network has two hidden units.

Training Loop

We run a “forward pass” and a “backward pass” many times (we refer to these as epochs).

Forward pass

The forward pass takes the input X, feeds it through the network layer by layer, and computes the output a2. Then it calculates how far off the prediction is using a loss function.

# Forward pass
z1 = np.dot(X, W1) + b1
a1 = sigmoid(z1)

z2 = np.dot(a1, W2) + b2
a2 = sigmoid(z2)

loss = mse_loss(y, a2)

In this step, we are calculating the loss for the current set of weights.

This loss is a measure of how “wrong” the network is, and it’s what drives the learning process in the backward pass.

Backward pass

The backward pass is how the network learns—by adjusting the weights based on how much they contributed to the final error. This is done by applying the chain rule in reverse across the network.

# Step 1: Derivative of loss with respect to output (a2)
d_loss_a2 = 2 * (a2 - y) / y.size

This computes the gradient of the mean squared error loss with respect to the output. It answers: How much does a small change in the output affect the loss?

\[\frac{\partial \text{Loss}}{\partial \hat{y}} = \frac{2}{n} (\hat{y} - y)\]
# Step 2: Derivative of sigmoid at output layer
d_a2_z2 = sigmoid_derivative(a2)
d_z2 = d_loss_a2 * d_a2_z2

Now we apply the chain rule. Since the output passed through a sigmoid function, we compute the derivative of the sigmoid to see how a change in the pre-activation \(z_2\) affects the output.

# Step 3: Gradients for W2 and b2
d_W2 = np.dot(a1.T, d_z2)
d_b2 = np.sum(d_z2, axis=0, keepdims=True)
  • a1.T is the transposed output from the hidden layer.
  • d_z2 is the error signal coming back from the output.
  • The dot product calculates how much each weight in W2 contributed to the error.
  • The bias gradient is simply the sum across all samples.
# Step 4: Propagate error back to hidden layer
d_a1 = np.dot(d_z2, W2.T)
d_z1 = d_a1 * sigmoid_derivative(a1)

Now we move the error back to the hidden layer:

  • d_a1 is the effect of the output error on the hidden layer output.
  • We multiply by the derivative of the hidden layer activation to get the true gradient of the hidden pre-activations.
# Step 5: Gradients for W1 and b1
d_W1 = np.dot(X.T, d_z1)
d_b1 = np.sum(d_z1, axis=0, keepdims=True)
  • X.T is the input data, transposed.
  • We compute how each input feature contributed to the hidden layer error.

This entire sequence completes one application of backpropagation—moving from output to hidden to input layer, using the chain rule and computing gradients at each step.

The final gradients (d_W1, d_W2, d_b1, d_b2) are then used in the weight update step:

# Apply the gradients to update the weights
W2 -= learning_rate * d_W2
b2 -= learning_rate * d_b2
W1 -= learning_rate * d_W1
b1 -= learning_rate * d_b1

This updates the model just a little bit—nudging the weights toward values that reduce the overall loss.

Final Predictions

print("\nFinal predictions:")
print(a2)

When we ran this code, we saw:

Epoch 0, Loss: 0.2558
...
Epoch 9000, Loss: 0.1438

Final predictions:
[[0.1241]
[0.4808]
[0.8914]
[0.5080]]

Interpreting the Results

The network is getting better, but not perfect. Let’s look at what these predictions mean:

Input Expected Predicted Interpreted
[0, 0] 0 0.1241 0
[0, 1] 1 0.4808 ~0.5
[1, 0] 1 0.8914 1
[1, 1] 0 0.5080 ~0.5

It’s nailed [1, 0] and is close on [0, 0], but it’s uncertain about [0, 1] and [1, 1]. That’s okay—XOR is a tough problem when learning from scratch with minimal capacity.

This ambiguity is actually a great teaching point: neural networks don’t just “flip a switch” to get things right. They learn gradually, and sometimes unevenly, especially when training conditions (like architecture or learning rate) are modest.

You can tweak the hidden layer size, activation functions, or even the optimizer to get better results—but the core algorithm stays the same: forward pass, loss computation, backpropagation, weight update.

Conclusion

As it stands, this tiny XOR network is a full demonstration of what makes neural networks learn.

You’ve now seen backpropagation from the inside.

A full version of this program can be found as a gist.

Building Lazy Composition in Rust

Introduction

Rust programmers encounter combinators all the time: map(), and_then(), filter(). They’re everywhere in Option, Result, Iterator, and of course, Future. But if you’re coming from a functional programming background — or just curious how these things work — you might ask:

What actually is a combinator?

Let’s strip it down to the bare minimum: a value, a function, and some deferred execution.

A Lazy Computation

We’ll start with a structure called Thunk. It wraps a closure that does some work, and it defers that work until we explicitly ask for it via .run().

pub struct Thunk<F> {
    f: Option<F>,
}

impl<F, R> Thunk<F>
where
    F: FnOnce() -> R,
{
    pub fn new(f: F) -> Self {
        Self { f: Some(f) }
    }

    pub fn run(mut self) -> R {
        let f = self.f.take().expect("already run");
        f()
    }
}

It’s essentially a one-shot deferred computation. We stash a closure inside, and we invoke it only when we’re ready.

Here, F is the type of the closure (the function) we’re wrapping, and R is the result it will produce once called. This lets Thunk be generic over any one-shot computation.

The work here is really wrapped up by self.f.take() which will force the value.

Simple.

Example

Here’s what this looks like in practice:

fn main() {
    let add_one = Thunk::new(|| 3 + 1);
    let result = add_one.run();
    println!("Result: {}", result); // prints 4
}

No magic. No threading. No async. Just a delayed function call.

Composing Thunks

The real value in combinators is that they compose. We can make more complex computations out of simpler ones — without immediately evaluating them.

Here’s how we can build on top of multiple Thunks:

fn main() {
    let m = Thunk::new(|| 3 + 1); // 4
    let n = Thunk::new(|| 9 + 1); // 10

    let o = Thunk::new(|| m.run() + n.run()); // 14
    let result = o.run();

    println!("Result: {}", result);
}

We’ve built a new computation (o) that depends on two others (m and n). They won’t run until o.run() is called — and then they run in the correct order, and just once.

Look Familiar?

If you’ve spent time in Haskell, this structure might look suspiciously familiar:

fmap :: Functor f => (a -> b) -> f a -> f b

This is a form of fmap. We’re not building a full trait implementation here, but the shape is the same. We can even imagine extending Thunk with a map() method:

impl<F, R> Thunk<F>
where
    F: FnOnce() -> R,
{
    pub fn map<G, S>(self, g: G) -> Thunk<impl FnOnce() -> S>
    where
        G: FnOnce(R) -> S,
    {
        Thunk::new(|| g(self.run()))
    }
}

And now:

let t = Thunk::new(|| 42);
let u = t.map(|x| x * 2);
assert_eq!(u.run(), 84);

No typeclasses, no lifetimes — just combinator building blocks.

From Lazy to Async

Now here’s the twist. What if our .run() method couldn’t give us a value right away? What if it needed to register a waker, yield, and be polled later?

That’s exactly what happens in Rust’s async system. The structure is the same — a value and a function bundled together — but the execution context changes. Instead of calling .run(), we implement Future and respond to .poll().

Here’s what that looks like for a simple async Map combinator:

use std::future::Future;
use std::pin::Pin;
use std::task::{Context, Poll};
use pin_project::pin_project;

// Our Map combinator
#[pin_project]
pub struct Map<Fut, F> {
    #[pin]
    future: Fut,

    f: Option<F>, // Option to allow taking ownership in poll
}

impl<Fut, F> Map<Fut, F> {
    pub fn new(future: Fut, f: F) -> Self {
        Self { future, f: Some(f) }
    }
}

impl<Fut, F, T, U> Future for Map<Fut, F>
where
    Fut: Future<Output = T>,
    F: FnOnce(T) -> U,
{
    type Output = U;

    fn poll(self: Pin<&mut Self>, cx: &mut Context<'_>) -> Poll<Self::Output> {
        let mut this = self.project();

        match this.future.poll(cx) {
            Poll::Pending => Poll::Pending,
            Poll::Ready(val) => {
                let f = this.f.take().expect("polled Map after completion");
                Poll::Ready(f(val))
            }
        }
    }
}

// Helper function to use it ergonomically
pub fn map<Fut, F, T, U>(future: Fut, f: F) -> Map<Fut, F>
where
    Fut: Future<Output = T>,
    F: FnOnce(T) -> U,
{
    Map::new(future, f)
}

Let’s take a step back and notice something: this structure is almost identical to Thunk. We’re still storing a value (future) and a function (f), and the combinator (Map) still controls when that function is applied. The difference is that we now interact with the asynchronous task system via poll(), instead of calling .run() ourselves.

This is how Future combinators in futures and tokio work under the hood — by carefully pinning, polling, and composing smaller futures into larger ones.

This is essentially a hand-rolled version of what futures::FutureExt::map() gives you for free.

As a simple example, we can use this as follows:

#[tokio::main]
async fn main() {
    let fut = async { 21 };
    let mapped = map(fut, |x| x * 2);
    let result = mapped.await;
    println!("Result: {}", result); // Should print 42
}

Conclusion

We often think of combinators as “just utility functions.” But they’re really more than that: they’re a way of thinking. Package a value and a transformation together. Delay the work. Compose more when you’re ready.

So the next time you write .map(), remember — it’s just a Thunk waiting to happen.

Create your own Filesystem with FUSE

Introduction

FUSE is a powerful Linux kernel module that lets you implement your own filesystems entirely in user space. No kernel hacking required. With it, building your own virtual filesystem becomes surprisingly achievable and even… fun.

In today’s article, we’ll build a filesystem that’s powered entirely by HTTP. Every file operation — reading a file, listing a directory, even getting file metadata — will be handled by a REST API. On the client side, we’ll use libcurl to perform HTTP calls from C, and on the server side, a simple Python Flask app will serve as our in-memory file store.

Along the way, you’ll learn how to:

  • Use FUSE to handle filesystem operations in user space
  • Make REST calls from C using libcurl
  • Create a minimal RESTful backend for serving file content
  • Mount and interact with your filesystem like any other directory

Up in my github repository I have added this project if you’d like to pull it down and try it. It’s called restfs.

Let’s get into it.

Defining a FUSE Filesystem

Every FUSE-based filesystem starts with a fuse_operations struct. This is essentially a table of function pointers — you provide implementations for the operations you want your filesystem to support.

Here’s the one used in restfs:

static struct fuse_operations restfs_ops = {
    .getattr = restfs_getattr,
    .readdir = restfs_readdir,
    .open    = restfs_open,
    .read    = restfs_read
};

This tells FUSE: “When someone calls stat() on a file, use restfs_getattr. When they list a directory, use restfs_readdir, and so on.”

Let’s break these down:

  • getattr: Fills in a struct stat with metadata about a file or directory — size, mode, timestamps, etc. It’s the equivalent of stat(2).
  • readdir: Lists the contents of a directory. It’s how ls knows what to show.
  • open: Verifies that a file can be opened. You don’t need to return a file descriptor — just confirm the file exists and is readable.
  • read: Reads data from a file into a buffer. This is where the real I/O happens.

Each function corresponds to a familiar POSIX operation. For this demo, we’re implementing just the basics — enough to mount the FS, ls it, and cat a file.

If you leave an operation out, FUSE assumes it’s unsupported — for example, we haven’t implemented write, mkdir, or unlink, so the filesystem will be effectively read-only.

Making REST Calls from C with libcurl

To interact with our HTTP-based server, we use libcurl, a powerful and flexible HTTP client library for C. In restfs, we wrap libcurl in a helper function called http_io() that performs an HTTP request and returns a parsed response object.

Here’s the core of the function:

struct _rest_response* http_io(const char *url, const char *body, const char *type) {
   CURL *curl = NULL;
   CURLcode res;
   long status = 0L;

   struct _http_write_buffer buf;
   buf.data = malloc(1);
   buf.size = 0;

   curl = curl_easy_init();

   if (curl) {
      curl_easy_setopt(curl, CURLOPT_URL, url);
      curl_easy_setopt(curl, CURLOPT_CUSTOMREQUEST, type);

      if (body) {
         curl_easy_setopt(curl, CURLOPT_POSTFIELDS, body);
         curl_easy_setopt(curl, CURLOPT_POSTFIELDSIZE, strlen(body));
      }

      curl_easy_setopt(curl, CURLOPT_WRITEFUNCTION, http_write_callback);
      curl_easy_setopt(curl, CURLOPT_WRITEDATA, (void *)&buf);

      curl_easy_setopt(curl, CURLOPT_USERAGENT, _http_user_agent);

      struct curl_slist *headers = NULL;
      headers = curl_slist_append(headers, "Content-Type: application/json");
      curl_easy_setopt(curl, CURLOPT_HTTPHEADER, headers);

      res = curl_easy_perform(curl);
      curl_easy_getinfo(curl, CURLINFO_RESPONSE_CODE, &status);
      curl_easy_cleanup(curl);
      curl_slist_free_all(headers);

      if (res != CURLE_OK) {
         fprintf(stderr, "error: %s\n", curl_easy_strerror(res));
         if (buf.data) free(buf.data);
         return NULL;
      }
   }

   return rest_make_response(buf.data, buf.size, status);
}

Let’s break it down:

  • curl_easy_init() creates a new easy handle.
  • CURLOPT_URL sets the request URL.
  • CURLOPT_CUSTOMREQUEST lets us specify GET, POST, PUT, DELETE, etc.
  • If a body is provided (e.g. for POST/PUT), we pass it in using CURLOPT_POSTFIELDS.
  • CURLOPT_WRITEFUNCTION and CURLOPT_WRITEDATA capture the server’s response into a buffer.
  • Headers are added manually to indicate we’re sending/expecting JSON.
  • After the call, we extract the HTTP status code and clean up.

The result is returned as a _rest_response struct:

struct _rest_response {
   int status;
   json_object *json;
   char *data;     // raw response body
   size_t length;  // response size in bytes
};

This makes it easy to access either the full raw data or a parsed JSON object depending on the use case.

To parse the JSON responses from the server, we use the json-c library — a lightweight and widely used C library for working with JSON data. This allows us to easily extract fields like st_mode, st_size, or timestamps directly from the server’s responses.

To simplify calling common HTTP methods, we define a few handy macros:

#define rest_get(uri)         http_io(uri, NULL, "GET")
#define rest_delete(uri)      http_io(uri, NULL, "DELETE")
#define rest_post(uri, body)  http_io(uri, body, "POST")
#define rest_put(uri, body)   http_io(uri, body, "PUT")

With these in place, calling a REST endpoint is as simple as:

struct _rest_response *res = rest_get("/getattr?path=/hello.txt");

This layer abstracts away the curl boilerplate so each FUSE handler can focus on interpreting the result.

The Backend

So far we’ve focused on the FUSE client — how file operations are translated into HTTP requests. But for the system to work, we need something on the other side of the wire to respond.

Enter: a minimal Python server built with Flask.

This server acts as a fake in-memory filesystem. It knows nothing about actual disk files — it just stores a few predefined paths and returns metadata and file contents in response to requests.

Let’s look at the key parts:

  • A Python dictionary (fs) holds a small set of files and their byte contents.
  • The /getattr endpoint returns a JSON version of struct stat for a given file path.
  • The /readdir endpoint lists all available files (we only support the root directory).
  • The /read endpoint returns a slice of the file contents, based on offset and size.

Here’s a simplified version of the server:

from flask import Flask, request, jsonify
from urllib.parse import unquote
import os, stat, time

app = Flask(__name__)
fs = { '/hello.txt': b"Hello, RESTFS!\\n" }

def now(): return { "tv_sec": int(time.time()), "tv_nsec": 0 }

@app.route('/getattr')
def getattr():
    path = unquote(request.args.get('path', ''))
    if path == "/":
        return jsonify({ "st_mode": stat.S_IFDIR | 0o755, ... })
    if path in fs:
        return jsonify({ "st_mode": stat.S_IFREG | 0o644, "st_size": len(fs[path]), ... })
    return ('Not Found', 404)

@app.route('/readdir')
def readdir():
    return jsonify([name[1:] for name in fs.keys()])  # ['hello.txt']

@app.route('/read')
def read():
    path = request.args.get('path')
    offset = int(request.args.get('offset', 0))
    size = int(request.args.get('size', 4096))
    return fs[path][offset:offset+size]

This is enough to make ls and cat work on the mounted filesystem. The client calls getattr and readdir to explore the directory, and uses read to pull down bytes from the file.

End to End

With the server running and the client compiled, we can now bring it all together.

Start the Flask server in one terminal:

python server.py

Then, in another terminal, create a mountpoint and run the restfs client:

mkdir /tmp/restmnt
./restfs --base http://localhost:5000/ /tmp/restmnt -f

Now try interacting with your mounted filesystem just like any other directory:

➜  restmnt ls -l
total 1
-rw-r--r-- 1 michael michael  6 Jan  1  1970 data.bin
-rw-r--r-- 1 michael michael 15 Jan  1  1970 hello.txt

➜  restmnt cat hello.txt
Hello, RESTFS!

You should see logs from the server indicating incoming requests:

[GETATTR] path=/
127.0.0.1 - - [18/Aug/2025 21:29:46] "GET /getattr?path=/ HTTP/1.1" 200 -
[READDIR] path=/
127.0.0.1 - - [18/Aug/2025 21:29:46] "GET /readdir?path=/ HTTP/1.1" 200 -
[GETATTR] path=/hello.txt
127.0.0.1 - - [18/Aug/2025 21:29:46] "GET /getattr?path=/hello.txt HTTP/1.1" 200 -
127.0.0.1 - - [18/Aug/2025 21:29:47] "GET /open?path=/hello.txt HTTP/1.1" 200 -
127.0.0.1 - - [18/Aug/2025 21:29:47] "GET /read?path=/hello.txt&offset=0&size=4096 HTTP/1.1" 200 -
[GETATTR] path=/
127.0.0.1 - - [18/Aug/2025 21:29:47] "GET /getattr?path=/ HTTP/1.1" 200 -

Under the hood, every file operation is being translated into a REST call, logged by the Flask server, and fulfilled by your in-memory dictionary.

This is where the whole thing becomes delightfully real — you’ve mounted an HTTP API as if it were a native part of your filesystem.

Conclusion

restfs is a fun and minimal example of what FUSE can unlock — filesystems that aren’t really filesystems at all. Instead of reading from disk, we’re routing every file operation over HTTP, backed by a tiny REST server.

While this project is intentionally lightweight and a bit absurd, the underlying ideas are surprisingly practical. FUSE is widely used for things like encrypted filesystems, network mounts, and user-space views over application state. And libcurl remains a workhorse for robust HTTP communication in C programs.

What you’ve seen here is just the start. You could extend restfs to support writing files, persisting data to disk, mounting a remote object store, or even representing entirely virtual data (like logs, metrics, or debug views).

Sometimes the best way to understand a system is to reinvent it — badly, on purpose.

Time Integration in Physics Simulations

Introduction

When simulating physical systems—whether it’s a bouncing ball, orbiting planets, or particles under gravity—accurately updating positions and velocities over time is crucial. This process is known as time integration, and it’s the backbone of most game physics and real-time simulations.

In this post, we’ll explore two fundamental methods for time integration: Euler’s method and Runge-Kutta 4 (RK4).

We’ll go through how each of these methods is represented mathemtically, and then we’ll translate that into code. We’ll build a small visual simulation in Python using pygame to see how the two methods behave differently when applied to the same system.

The Simulation

Our simulation consists of a central massive object (a “sun”) and several orbiting bodies, similar to a simplified solar system. Each body is influenced by the gravitational pull of the others, and we update their positions and velocities in each frame of the simulation loop.

At the heart of this simulation lies a decision: how should we advance these objects forward in time? This is where the integration method comes in.

Euler’s Method

Euler’s method is the simplest way to update motion over time. It uses the current velocity to update position, and the current acceleration to update velocity:

\[\begin{aligned} \vec{x}_{t+\Delta t} &= \vec{x}_t + \vec{v}_t \cdot \Delta t \\\\ \vec{v}_{t+\Delta t} &= \vec{v}_t + \vec{a}_t \cdot \Delta t \end{aligned}\]

This translates down into the following python code:

def step_euler(bodies):
    accs = [compute_acc(bodies, i) for i in range(len(bodies))]
    for b, a in zip(bodies, accs):
        b.pos += b.vel * DT
        b.vel += a * DT

This is easy to implement, but has a major downside: error accumulates quickly, especially in systems with strong forces or rapidly changing directions.

Here’s an example of it running:

RK4

Runge-Kutta 4 (RK4) improves on Euler by sampling the system at multiple points within a single timestep. It estimates what will happen halfway through the step, not just at the beginning. This gives a much better approximation of curved motion and reduces numerical instability.

Runge-Kutta 4 samples the derivative at four points:

\[\begin{aligned} \vec{k}_1 &= f(t, \vec{y}) \\\\ \vec{k}_2 &= f\left(t + \frac{\Delta t}{2}, \vec{y} + \frac{\vec{k}_1 \cdot \Delta t}{2}\right) \\\\ \vec{k}_3 &= f\left(t + \frac{\Delta t}{2}, \vec{y} + \frac{\vec{k}_2 \cdot \Delta t}{2}\right) \\\\ \vec{k}_4 &= f\left(t + \Delta t, \vec{y} + \vec{k}_3 \cdot \Delta t\right) \\\\ \vec{y}_{t+\Delta t} &= \vec{y}_t + \frac{\Delta t}{6}(\vec{k}_1 + 2\vec{k}_2 + 2\vec{k}_3 + \vec{k}_4) \end{aligned}\]
Combined State! The y vector here represents both position and velocity as a combined state vector.

This translates down into the following python code:

def step_rk4(bodies):
    n = len(bodies)
    pos0 = [b.pos.copy() for b in bodies]
    vel0 = [b.vel.copy() for b in bodies]

    a1 = [compute_acc(bodies, i) for i in range(n)]

    for i, b in enumerate(bodies):
        b.pos = pos0[i] + vel0[i] * (DT / 2)
        b.vel = vel0[i] + a1[i] * (DT / 2)
    a2 = [compute_acc(bodies, i) for i in range(n)]

    for i, b in enumerate(bodies):
        b.pos = pos0[i] + b.vel * (DT / 2)
        b.vel = vel0[i] + a2[i] * (DT / 2)
    a3 = [compute_acc(bodies, i) for i in range(n)]

    for i, b in enumerate(bodies):
        b.pos = pos0[i] + b.vel * DT
        b.vel = vel0[i] + a3[i] * DT
    a4 = [compute_acc(bodies, i) for i in range(n)]

    for i, b in enumerate(bodies):
        b.pos = pos0[i] + vel0[i] * DT + (DT**2 / 6) * (a1[i] + 2*a2[i] + 2*a3[i] + a4[i])
        b.vel = vel0[i] + (DT / 6) * (a1[i] + 2*a2[i] + 2*a3[i] + a4[i])

RK4 requires more code and computation, but the visual payoff is immediately clear: smoother orbits, fewer explosions, and longer-lasting simulations.

Here’s an example of it running:

Trade-offs

Euler is fast and simple. It’s great for prototyping, simple games, or systems where precision isn’t critical.

RK4 is more accurate and stable, especially in chaotic or sensitive systems—but it’s computationally more expensive. In real-time applications (like games), you’ll need to weigh performance vs. quality.

Also, both methods depend heavily on the size of the timestep. Larger steps amplify error; smaller ones improve accuracy at the cost of performance.

Conclusion

Switching from Euler to RK4 doesn’t just mean writing more code—it fundamentally changes how your simulation evolves over time. If you’re seeing odd behaviors like spiraling orbits, exploding systems, or jittery motion, trying a higher-order integrator like RK4 might fix it.

Or, it might inspire a deeper dive into the world of numerical simulation—welcome to the rabbit hole!

You can find the full code listing here as a gist, so you can tweak and run it for yourself.

Getting Started with ClojureScript

Introduction

I recently decided to dip my toes into ClojureScript. As someone who enjoys exploring different language ecosystems, I figured getting a basic “Hello, World!” running in the browser would be a fun starting point. It turns out that even this small journey taught me quite a bit about how ClojureScript projects are wired together.

This post captures my first successful setup: a minimal ClojureScript app compiled with lein-cljsbuild, rendering output in the browser console.

A Rough Start

I began with the following command to create a new, blank project:

lein new cljtest

First job from here is to organise dependencies, and configure the build system for the project.

project.clj

There’s a few things to understand in the configuration of the project:

  • We add org.clojure/clojurescript "1.11.132" as a dependency
  • To assist with our builds, we add the plugin lein-cljsbuild "1.1.8"
  • The source path is normally src, but we change this for ClojureScript to src-cljs
  • The output will be javascript output for a website, and all of our web assets go into resources/public
(defproject cljtest "0.1.0-SNAPSHOT"
  :min-lein-version "2.9.1"
  :description "Minimal ClojureScript Hello World"
  :dependencies [[org.clojure/clojure "1.11.1"]
                 [org.clojure/clojurescript "1.11.132"]]
  :plugins [[lein-cljsbuild "1.1.8"]]
  :source-paths ["src-cljs"]
  :clean-targets ^{:protect false} ["resources/public/js" "target"]

  :cljsbuild
  {:builds
   {:dev
    {:source-paths ["src-cljs"]
     :compiler {:main cljtest.core
                :output-to "resources/public/js/main.js"
                :output-dir "resources/public/js/out"
                :asset-path "js/out"
                :optimizations :none
                :source-map true
                :pretty-print true}}

    :prod
    {:source-paths ["src-cljs"]
     :compiler {:main cljtest.core
                :output-to "resources/public/js/main.js"
                :optimizations :advanced
                :pretty-print false}}}})

We have two different build configurations here: dev and prod.

The dev configuration focuses on being much quicker to build so that the change / update cycle during development is quicker. Source maps, pretty printing, and no optimisations provide the verbose output appropriate for debugging.

The prod configuration applies all the optimisations. This build is slower, but produces one single output file: main.js. This is the configuration that you use to “ship” your application.

Your First ClojureScript File

Place this in src-cljs/cljtest/core.cljs:

(ns cljtest.core)

(enable-console-print!)
(println "Hello from ClojureScript!")

HTML Page to Load It

Create a file at resources/public/index.html:

<!doctype html>
<html>
  <head><meta charset="utf-8"><title>cljtest</title></head>
  <body>
    <h1>cljtest</h1>
    <script src="js/out/goog/base.js"></script>
    <script src="js/main.js"></script>
    <script>goog.require('cljtest.core');</script>
  </body>
</html>

Build & Run

Compile your dev build:

lein clean
lein cljsbuild once dev

Then open resources/public/index.html in your browser, and check the developer console — you should see your message.

If you want to iterate while coding:

lein cljsbuild auto dev

When you’re ready to build a production bundle:

lein cljsbuild once prod

Then you can simplify the HTML:

<script src="js/main.js"></script>

No goog.require needed — it all gets bundled.

Step it up

Next, we’ll step up to something a little more useful. We’ll put together a table of names that we can add, edit, delete, etc. Just a really simple CRUD style application.

In order to do this, we’re going to rely on a pretty cool library called reagent.

We add the following dependency to project.clj:

[reagent "1.0.0"]

State

Our little application requires some state:

(defonce names (r/atom [{:id 1 :name "Alice"}
                        {:id 2 :name "Bob"}]))

(defonce next-id (r/atom 3))
(defonce editing-id (r/atom nil))
(defonce edit-text (r/atom ""))

names is the currentl list of names. next-id gives us the next value that we’ll use an ID when adding a new record. editing-id and edit-text manage the state for updates.

Table

We can now render our table using a simple function:

(defn name-table []
  [:div
   [:h2 "Name Table"]
   [:table
    [:thead
     [:tr [:th "Name"] [:th "Edit"] [:th "Delete"]]]
    [:tbody
     (for [n @names]
       ^{:key (:id n)} [name-row n])]]
   [:div
    [:input {:placeholder "New name"
             :value @edit-text
             :on-change #(reset! edit-text (.. % -target -value))}]
    [:button {:on-click
              #(when-not (clojure.string/blank? @edit-text)
                 (swap! names conj {:id @next-id :name @edit-text})
                 (swap! next-id inc)
                 (reset! edit-text ""))}
     "Add"]]])

The table renders all of the names, as well and handles the create case. The edit case is a little more complex and requires a function of its own. The name-row function manages this complexity for us.

(defn name-row [{:keys [id name]}]
  [:tr
   [:td name]
   [:td
    (if (= id @editing-id)
      [:<>
       [:input {:value @edit-text
                :on-change #(reset! edit-text (.. % -target -value))}]
       [:button {:on-click
                 (fn []
                   (swap! names (fn [ns]
                                  (mapv (fn [n]
                                          (if (= (:id n) id)
                                            (assoc n :name @edit-text)
                                            n))
                                        ns)))
                   (reset! editing-id nil))}
        "Save"]]
      [:<>
       [:button {:on-click #(do (reset! editing-id id)
                                (reset! edit-text name))}
        "Edit"]])]
   [:td
    [:button {:on-click
              (fn []
                (swap! names (fn [ns]
                               (vec (remove (fn [n] (= (:id n) id)) ns)))))} ;; FIX
     "Delete"]]])

Mounting!

Now we’re going to make sure that these functions end up on our web page.

(defn mount-root []
  (dom/render [name-table] (.getElementById js/document "app")))

(defn init []
  (enable-console-print!)
  (mount-root))

We need an app element in our HTML page.

<!doctype html>
<html>
  <head>
    <meta charset="utf-8">
    <title>cljtest</title>
  </head>
  <body>
    <h1>cljtest</h1>

    <!-- This is our new element! -->
    <div id="app"></div>

    <script src="js/out/goog/base.js"></script>
    <script src="js/main.js"></script>
    <script>goog.require('cljtest.core'); cljtest.core.init();</script>

  </body>
</html>

Conclusion

This journey started with a humble goal: get a simple ClojureScript app running in the browser. Along the way, I tripped over version mismatches, namespace assumptions, and nested anonymous functions — but I also discovered the elegance of Reagent and the power of functional UIs in ClojureScript.

While the setup using lein-cljsbuild and Reagent 1.0.0 may feel a bit dated, it’s still a solid way to learn the fundamentals. From here, I’m looking forward to exploring more advanced tooling like Shadow CLJS, integrating external JavaScript libraries, and building more interactive UIs.

This was my first real toe-dip into ClojureScript, and already I’m hooked. Stay tuned — there’s more to come.