Demos

The examples elsewhere in this guide are short and meant to be read. The programs under demos/ are meant to be run. Each is a whole program you can build, execute and change, and most print something worth looking at.

Each entry below says what the demo is for, then shows the few lines where it actually calls clad. Those lines are pulled from the demo itself, so they are the real code and not a copy that can go out of date. The name links to the whole program in the repository, which is where to go once the excerpt has told you whether it is the one you want.

They come in two groups. The basics are short and each one shows a single thing clad does; read those first if you are new to it. The rest are programs with a purpose of their own that happen to need a derivative, and are worth reading once you know what you are looking at.

Building one

Nothing in the build system builds the demos. A demo is compiled the way any program using clad is, by loading the plugin:

clang++ -std=c++17 -I /path/to/clad/include \
        -fplugin=/path/to/clad.so demos/GradientDescent.cpp -o demo

Installation and usage explains the flags and the options reference lists what clad accepts. Demos needing more than that say so.

The basics

Each of these is one file that builds with clad and a compiler, and exists to show a single thing clad can do.

Taking a derivative

Gradient.cpp

The same trick in miniature and easier to read first: the direction a sphere faces at a point, from three derivatives of the sphere’s equation.

  auto sphere_implicit_func_dx = clad::differentiate(sphere_implicit_func, 0);
  auto sphere_implicit_func_dy = clad::differentiate(sphere_implicit_func, 1);
  auto sphere_implicit_func_dz = clad::differentiate(sphere_implicit_func, 2);

Using the second derivative

NewtonsMethod.cpp

Newton’s method walks to the bottom of the Rosenbrock function, a long curved valley whose floor is nearly flat. Following the slope alone crawls once you are in the valley; the second derivative says how the slope is itself changing, which is what lets a step cross the floor rather than inch along it. Clad writes both, so the method is a few lines of arithmetic.

double rosenbrock(double x, double y) {
  return (x - 1) * (x - 1) + 100 * (y - x * x) * (y - x * x);
}

  auto grad = clad::gradient(rosenbrock);
  auto hess = clad::hessian(rosenbrock);

Every derivative at once

CoordinateChange.cpp

A gradient is for a function with one output. This one has three, and the whole table of partial derivatives is the jacobian, which clad fills in a single pass. The demo takes its determinant and prints it beside r*r*sin(theta), the factor every integral in spherical coordinates carries, so you can see they agree.

void spherical_to_cartesian(double r, double theta, double phi, double p[]) {
  p[0] = r * std::sin(theta) * std::cos(phi);
  p[1] = r * std::sin(theta) * std::sin(phi);
  p[2] = r * std::cos(theta);
}

  auto jac = clad::jacobian(spherical_to_cartesian);

VectorForwardMode.cpp

Forward mode normally costs one pass per input you ask about. Vector mode does them together, which matters here because the inputs are two arrays whose length is only known while the program runs.

  auto weighted_sum_grad =
      clad::differentiate<clad::opts::vector_mode>(weighted_sum, "arr,weights");

Differentiating code that is not a formula

KeplerEquation.cpp

Where a body is on its orbit has no closed form: you iterate until the answer stops moving. The loop runs as many times as its arguments make it run and stops on a value computed inside it, so there is no formula to differentiate. Clad differentiates what the program does, and the demo prints its answer beside the one worked out by hand so you can see they agree.

double eccentric_anomaly(double M, double e) {
  double E = M;
  for (int i = 0; i < 40; ++i) {
    double step = (E - e * std::sin(E) - M) / (1 - e * std::cos(E));
    E -= step;
    if (std::fabs(step) < 1e-15)
      break;
  }
  return E;
}

  auto grad = clad::gradient(eccentric_anomaly);

Reaching past what clad can differentiate

CustomDerivative.cpp

Some code is not worth differentiating. This raises a number to a power by treating its bit pattern as a logarithm – fast, a few percent out, and with no derivative you would want. So you write the derivative down instead and clad uses it without reading the body, which is also what you do for a function from a library whose source you do not have. The demo prints the approximate value beside the exact derivatives to show which came from where.

namespace clad {
namespace custom_derivatives {
void fast_pow_pullback(float base, float exponent, float d_result,
                       float* d_base, float* d_exponent) {
  *d_base += d_result * exponent * ::std::pow(base, exponent - 1.f);
  *d_exponent += d_result * ::std::pow(base, exponent) * ::std::log(base);
}
} // namespace custom_derivatives
} // namespace clad

  auto grad = clad::gradient(model);

CustomTypeNumDiff.cpp

Some types clad cannot take apart, like a number stored as a scaled integer. It falls back to measuring the derivative instead of deriving it, by evaluating the function at nearby points.

  numerical_diff::central_difference(func, grad, /*printErrors=*/0, x, y);

Templates.cpp

A template whose long double version computes a different formula from the general one. Clad differentiates whichever version the compiler actually picked, not the one you wrote first.

  auto d_E_double = clad::differentiate(E_double, "i");
  auto d_E_long_double = clad::differentiate(E_long_double, "j");

Trusting the floating point

ErrorEstimation/FloatSum.cpp

Adds the same numbers twice, once plainly and once with a trick that compensates for rounding, and has clad estimate how much error each one built up. The gnuplot lines in the file plot the two against each other.

float vanillaSum(float x, unsigned int n) {
  float sum = 0.0;
  for (unsigned int i = 0; i < n; i++) {
    sum = sum + x;
  }
  return sum;
}

  auto df = clad::estimate_error(vanillaSum);

ErrorEstimation/CustomModel/

Clad’s built-in guess at the error is deliberately pessimistic: it reports the worst case. If you know more about your numbers you can say so by supplying your own model, which is what this does. It has a README.

namespace clad {
double getErrorVal(double dx, double x, const char* name) { return dx * x; }
} // namespace clad

ErrorEstimation/PrintModel/

A model that estimates nothing and simply reports every place clad would have accounted for error, which is the easiest way to see what the machinery is doing. It has a README.

namespace clad
{
    __attribute__((always_inline)) double getErrorVal(double dx, double x, const char *name)
    {
        double error = std::abs(dx * x * std::numeric_limits<float>::epsilon());
        std::cout << "Error in " << name << " : " << error << std::endl;
        return error;
    }
} // clad

Inside a real program

These are programs with a purpose of their own that happen to need a derivative. Some want a toolkit, and say so.

Learning from a derivative

GradientDescent.cpp

Fits a straight line to data. To know which way to nudge the line, you need the slope of the error with respect to each parameter, and clad supplies it. This is the shortest answer to why reverse mode exists: many parameters, one number to make smaller.

double cost(double theta_0, double theta_1, double x, double y) {
  double f_x = f(theta_0, theta_1, x);
  return (f_x - y) * (f_x - y);
}

  auto clad_grad = clad::gradient(cost);

HelixFit.cpp

A charged particle in a magnetic field moves along a helix, and a detector reports points it passed through; recovering the helix is how its momentum is measured. Each point gives two numbers that should be zero when the helix is right, and clad differentiates both. Levenberg-Marquardt turns those derivatives into a fit, which is how data is fitted in practice when the model is not a straight line. Reduced from a fitter contributed in #1202.

double radial_residual(Helix h, double x, double y, double z) {
  double dx = x - h.cx, dy = y - h.cy;
  return std::sqrt(dx * dx + dy * dy) - h.r;
}

double z_residual(Helix h, double x, double y, double z) {
  return z - (h.z0 + h.lam * std::atan2(y - h.cy, x - h.cx));
}

  auto d_radial = clad::gradient(radial_residual, "h");
  auto d_z = clad::gradient(z_residual, "h");

XorNetwork.cpp

A network with one hidden layer learning exclusive or, the standard first example of something no straight line can separate. Nine weights, four training cases, and one call to clad; everything else is arithmetic on the gradient it returns. Between this and the GPT-2 above, the only thing that changes is the size.

double loss(const double w[9]) {
  double total = 0;
  for (int i = 0; i < 4; ++i) {
    double miss = net(w, kInput[i][0], kInput[i][1]) - kWanted[i];
    total += miss * miss;
  }
  return total;
}

  auto grad = clad::gradient(loss);

cladtorch/

The same idea, grown up: a GPT-2 that trains and writes text. One call asks for the gradient of the loss, and the training loop does the rest. Needs libtorch; download_training_data.sh fetches the text it learns from.

static float gpt2_loss(const GPT2& model, const ITensor& input,
                       const ITensor& targets) {
  auto probs = model.forward(input);
  auto loss = cross_entropy_loss(probs, targets);
  return loss.scalar();
}

static float gpt2_loss(const GPT2& model, const ITensor& input,
                       const ITensor& targets) {
  auto probs = model.forward(input);
  auto loss = cross_entropy_loss(probs, targets);
  return loss.scalar();
}

Differentiating a whole program

ODESolverSensitivity.cpp

Asks how much the answer of a differential equation would move if you changed the numbers you started with. Clad differentiates the solver itself, loops and all, so changing the equation or swapping the integrator gives new derivatives without deriving anything by hand.

double solution(double a, double b, double c, double x) {
  return rungeKutta(0, 0, x, 0.001, a, b, c);
}

  auto h = clad::gradient(solution);

ComputerGraphics/smallpt/

A path tracer. Its shapes are described by a function that says how far away a surface is, and the direction a surface faces is the derivative of that function, so clad works out the shading rather than a formula written by hand for every shape.

double sphere_distance_func(const Vec& p, const Vec& p0, double r) {
  return sqrt((p.x - p0.x) * (p.x - p0.x) + (p.y - p0.y) * (p.y - p0.y) +
              (p.z - p0.z) * (p.z - p0.z)) -
         r;
}

    auto dist_grad = clad::gradient(sphere_distance_func, "p, p0");

Running on a GPU

CUDA/

Six programs: VectorAddition, ParticleSimulation, TensorContraction, LinearRegression, BoWLogisticRegression, and a port of NVIDIA’s BlackScholes that keeps the original CPU version and checks the gradients against it. Most differentiate host code that drives the GPU through Thrust; TensorContraction differentiates a function that launches a kernel, which is what the excerpt shows. Needs a CUDA toolkit.

  auto tensor_grad = clad::gradient(launchTensorContraction3D, "C, A, B");

Without installing anything

Clad is on Compiler Explorer, which is the quickest way to see what it does to a function: write one, ask for its derivative, and read the code clad generated beside it. The build there comes from clad’s own master, so what you try is what is current.

Jupyter/Intro.ipynb

Clad in a notebook, a cell at a time. Binder runs this one in a browser; locally it needs a C++ Jupyter kernel, and the notebook uses xeus-cpp.