Functions

Writing Good Functions #

Functions are wonderfully simple and we want to use them liberally. They’re excellent for extracting shared functionality and structuring the code into different layers of detail. When used well they’re highly effective at reducing the amount of code while increasing the readability of the remaining code.

There’s three aspects to writing good functions: they have few arguments, they express one layer of abstraction and they’re short.

Few Arguments #

Functions then to have about two or three arguments. Occasionally, more maybe four or five; but that’s likely already a sign that something else is missing.

If an algorithm requires 8 parameters, then those can be passed in as a single struct.

The reasons to keep the number of arguments down are:

  • Easier to test, since otherwise one often needs to cover too many combinations. For example if there’s 5 boolean flags that each toggle some feature. We must test 32 variations. Maybe not all of those are valid.
  • It’s harder to remember all arguments and the order they must be in. Hence, it’s easier to make a mistake when calling the function.
  • The callsite gets quite lengthy, which may makes it harder to understand what the function does. It also makes the function less attractive to use and there’s a higher chance that not all arguments are available in the requested form.

Single Layer of Abstraction #

This is about controlling the amount of things one needs to keep track of while reading a function. There’s the constant threat that some piece of code modifies a value one didn’t expect. There’s also the problem that separating two part of an operation by lots of unrelated code will require scrolling and makes it much harder to be confident. Finally, 15 lines of code implementing one cohesive idea are not a suitable or obvious description of that idea.

We’ll use an example to demonstrate the principle.

std::vector<double>
solve_pde(const std::vector<double>& x, double t_stop) {
  size_t n = x.size();
  std::vector<double> u0(n);
  std::vector<double> u1(n);

  // Initial conditions
  for(size_t i = 0; i < n; ++i) {
    u0[i] = std::sin(x);
  }

  // Advance PDE
  while(t < t_stop) {
    for(size_t iL = 0; iL < n-1; ++iL) {
      size_t iR = iL + 1;

      double f = 0.5*(flux(u0[iL]) + flux(u0[iR)) + c * (u0[iL] - u1[iR]);
      u1[iR] = u0[iR] + dt * f;
      u1[iL] -= dt * f;
    }

    // Boundary conditions
    u1[0] = u1[1];
    u1[n-1] = u1[n-2];

    // Prepare for next iteration.
    std::swap(u0, u1);
    t += dt;
  }

  return u0;
}

(Note that in a real example the details of picking c and dt would add further lines of complexity.)

The comment before the initial conditions has already been identified that block as a function. Inside the function solve_pde the details of how to pick the initial conditions don’t matter. Next the for-loop inside a while-loop feels like it’s crossing levels of detail. It’s exposing lots of details about how the numerical flux f is being computed. We don’t immediate need to know, and it’s unrelated to the details about how the initial conditions where computed. The boundary conditions are implemented as a block of code and again, computing the boundary conditions is independent of the initial conditions and the numerical flux.

Let see what happens if we split the function into multiple functions:

std::vector<double>
solve_pde(const std::vector<double>& x, double t_stop) {
  size_t n = x.size();

  auto u0 = initial_conditions(x);
  auto u1 = std::vector<double>(n);
  while(t < t_stop) {
    double dt = cfl_condition(u0);
    advance_solution(u1, u0, dt);

    // Prepare for next iteration
    std::swap(u0, u1);
    t += dt;
  }
}

It’s clear that to solve a PDE we must first evaluate the initial conditions. Then on each time-step we advance the solution from buffer u0 into buffer u1 and then swap the two buffers. Since there’s few distractions, we can easily keep track of which buffer is which. Additionally, we now have space to express how to pick dt and we can also see where we’d add the logic for writing snapshots.

While a good argument could be made for fusing cfl_condition into advance_solution for performance reasons, the important point is that clear code is easy to understand and high performance can only be achieved through profiling and a precise mental model of how the algorithm is implemented.

Let’s look at advance_solution. The advantage is that it does not need to keep track of how u1 and u0 are being toggled. Its only job is to advance from state u0 and save the result in u1. As a result it’s shorter and there’s less context that might affect it, but ultimately doesn’t:

void advance_solution(      std::vector<double>& u1,
                      const std::vector<double>& u0,
                            double               dt) {

  size_t n = u1.size();
  for(size_t iL = 0; iL < n-1; ++iL) {
    size_t iR = iL + 1;

    double f = 0.5*(flux(u0[iL]) + flux(u0[iR)) + c * (u0[iL] - u1[iR]);
    u1[iR] = u0[iR] + dt * f;
    u1[iL] -= dt * f;
  }

  apply_boundary_conditions(u1);
}

One might again take offence that advance_solution has mixed levels of abstraction: the details of computing the numerical flux are still independent of the boundary conditions. Same as how the precise formula used for f is independent of the looping logic.

Note, that there are performance related difficulties to rigorously subdividing and making things nice from a purely abstract point of view. For example we could have reorganized the algorithm into two phases. First compute the rate of change:

std::vector<double>
compute_rate_of_change(const std::vector<double>& u) {

  auto dudt = std::vector<double>(u.size());
  for(size_t iL = 0; iL < n-1; ++iL) {
    size_t iR = iL + 1;

    double f = 0.5*(flux(u0[iL]) + flux(u0[iR)) + c * u0[iL] - u1[iR];
    dudt[iR] = f;
    dudt[iL] -= f;
  }

  return dudt;
}

and extract the forward Euler integration into a separate function:

std::vector<double>
advance_ode(const std::vector<double>& u0,
                  double               dt,
            const std::vector<double>& dudt) {

  auto u1 = std::vector<double>(u.size());
  for(size_t i = 0; i < n; ++i) {
    u1[i] = u0[i] + dt*dudt[i];
  }

  reutrn u1;
}

Personally, I feel this further improves the clarity with which we express the different algorithms used to solve the PDE. However, at the same time we’ve split one long loop into two equally long loops. Note that if we choose the fuse the two loops again, as long advance_solution doesn’t change it’s behaviour, solve_pde can’t possibly be affected by how we chose to organize advance_solution internally.

Two more comments about the example. One, it uses offensively generic names for highly specific algorithms, e.g. initial_conditions for one precise choice of initial conditions; apply_boundary_conditions for outflow boundary conditions; and advance_ode for forward Euler. By organizing the code as shown we immediately recognize the need to be able to customize each of those components of solving a PDE. You’ll likely want more powerful tools than functions to inject these parts nicely. Two, the example is very liberal about allocating memory. This might not be acceptable from a performance point of view, but since programming techniques exist for eliminating the overhead, we consider this aspect a simplification for the sake of demonstrating a different point.

Few Lines of Code #

This is more of a helpful mindset than a hard rule. It’s also an effective way of spotting severe violations of the single layer of abstraction principle.

What’s short enough? Maybe 5-10 lines. Even 20 lines may not be a huge concern; but somewhere as of 30 lines and certainly when it stops fitting on one screen (in landscape mode), one would probably want to start looking for missing layers of abstraction.