Tutorial

This tutorial builds a test in five levels, and each level adds one idea. A worked example at the end puts them together.

Level 0: one line

Give all_pairs a list of values for each argument of the function under test:

using UnitTestDesign

cases = all_pairs([1, 2, 3], ["a", "b"], [1.0, 2.0])
6 cases (minimal) · strength 2 · Auto: Construction() · 3 parameters · 12 combinations
    p1  p2   p3
 1  1   "a"  1.0
 2  2   "a"  2.0
 3  3   "a"  1.0
 4  1   "b"  2.0
 5  2   "b"  1.0
 6  3   "b"  2.0

The three lists have 3 × 2 × 2 = 12 combinations. The six cases hold every pair of values: each value of the first argument appears beside each value of the second, each value of the first beside each value of the third, and each value of the second beside each value of the third. The summary line says what was asked (strength 2, which means pairs), which engine built it (the default, Auto, and what it chose: here the catalog's array, Construction()), and how large the full product is. The columns are named p1, p2, p3 because the arguments have no names yet.

cases is a vector of tuples, so a test loops over it:

using Test

scaled(n, unit, x) = string(n * x, " ", unit)     # stands in for the code under test

@testset "scaled" begin
    for (n, unit, x) in cases
        number, suffix = split(scaled(n, unit, x))
        @test parse(Float64, number) == n * x
        @test suffix == unit
    end
end
Test Summary: | Pass  Total  Time
scaled        |   12     12  0.0s

Pairs matter because many bugs need two things at once: an if on one option inside a branch on another. The saving grows with the number of parameters. Four parameters of three values each have 81 combinations, and every pair of their values fits in 9 cases:

all_pairs([1, 2, 3], ["low", "mid", "high"], [1.0, 3.7, 4.9], [:greedy, :relax, :optim])
9 cases (minimal) · strength 2 · Auto: Construction() · 4 parameters · 81 combinations
    p1  p2      p3   p4
 1  1   "low"   1.0  :greedy
 2  2   "mid"   3.7  :greedy
 3  3   "high"  4.9  :greedy
 4  1   "mid"   4.9  :relax
 5  2   "high"  1.0  :relax
 6  3   "low"   3.7  :relax
 7  1   "high"  3.7  :optim
 8  2   "low"   4.9  :optim
 9  3   "mid"   1.0  :optim

Which values to list is the most important decision in this kind of testing, and no function makes it for you. Choosing values and oracles explains how to pick representatives, boundaries, and a check that knows the right answer.

Level 1: named parameters

Naming the parameters unlocks everything that follows: rules, partial must-include cases, groups of parameters that need a higher strength, and readable test names. Give name => values pairs, or a NamedTuple of lists:

all_pairs(:mode => [:fast, :exact], :solver => [:none, :lu, :qr], :tol => [1e-3, 1e-6])
6 cases (minimal) · strength 2 · Auto: Construction() · 3 parameters · 12 combinations
    mode    solver  tol
 1  :fast   :none   0.001
 2  :fast   :lu     1.0e-6
 3  :fast   :qr     0.001
 4  :exact  :none   1.0e-6
 5  :exact  :lu     0.001
 6  :exact  :qr     1.0e-6

The cases are now NamedTuples. The rest of this tutorial tests this function:

using LinearAlgebra

"""
Solve `A x = b`. In `:fast` mode, iterate until the residual is below `tol`,
ignoring `solver`. In `:exact` mode, meant for tight tolerances, factor `A`
with `solver`.
"""
function solve(A, b; mode = :exact, solver = :lu, tol = 1e-6)
    tol > 0 || throw(ArgumentError("tol must be positive, got $tol"))
    if mode == :fast                                  # Jacobi iteration
        x = zero(b)
        while norm(A * x - b) > tol * norm(b)
            x = x + (b - A * x) ./ diag(A)
        end
        return x
    end
    F = solver == :lu ? lu(A) : solver == :qr ? qr(A) : A
    return F \ b
end

When the same parameters feed several designs, or need rules, put them in a TestSpace:

space = TestSpace((mode = [:fast, :exact], solver = [:none, :lu, :qr], tol = [1e-3, 1e-6]))
TestSpace with 3 parameters, 12 combinations, 0 constraints
  mode:   :fast, :exact
  solver: :none, :lu, :qr
  tol:    0.001, 1.0e-6
cases = all_pairs(space)
6 cases (minimal) · strength 2 · Auto: Construction() · 3 parameters · 12 combinations
    mode    solver  tol
 1  :fast   :none   0.001
 2  :fast   :lu     1.0e-6
 3  :fast   :qr     0.001
 4  :exact  :none   1.0e-6
 5  :exact  :lu     0.001
 6  :exact  :qr     1.0e-6

A case splats into keyword arguments, and its fields read by name. The check here is an oracle that needs no reference answer: the residual must be below the tolerance.

A = [4.0 1.0; 1.0 3.0]
b = [1.0, 2.0]

@testset "solve" begin
    @testset "$case" for case in cases
        x = solve(A, b; case...)
        @test norm(A * x - b) <= case.tol * norm(b)
    end
end
Test Summary: | Pass  Total  Time
solve         |    6      6  0.1s

Level 2: constraints

Look at the cases above. solve ignores solver in fast mode, so the cases with mode = :fast and solver = :lu or :qr repeat what solver = :none already tests, while taking the place of combinations that matter. And exact mode is meant for tight tolerances, so mode = :exact with tol = 1e-3 is not a configuration anyone runs. Constraints say so, in three spellings. A pattern forbids one exact combination:

forbid((mode = :exact, tol = 1e-3); reason = "exact mode needs a tight tolerance")
exact mode needs a tight tolerance

An expression over bare parameter names either forbids the combinations that make it true (@forbid) or allows only those (@require):

@forbid mode == :fast && solver != :none
@forbid(mode == :fast && solver != :none)
@require mode == :exact || solver == :none
@require(mode == :exact || solver == :none)

@require e is the same rule as @forbid !(e); having both verbs keeps the intent visible. In a macro, a bare name is a parameter, and $x takes the value of a variable x from the surrounding code. When the logic is too big for one expression, list the names and write a function:

forbid(:mode, :solver) do mode, solver
    mode == :fast && solver != :none
end
forbid rule on (mode, solver)

A rule only ever sees complete values of the parameters it names, drawn from their lists. The solver rule is the common pattern for a dependent parameter: an option that applies only when another has some value gets a sentinel value (:none) and one @require.

Rules belong to the space:

space = TestSpace(
    (mode = [:fast, :exact], solver = [:none, :lu, :qr], tol = [1e-3, 1e-6]);
    constraints = [
        @require(mode == :exact || solver == :none),
        forbid((mode = :exact, tol = 1e-3); reason = "exact mode needs a tight tolerance"),
    ])

cases = all_pairs(space)
5 cases (lower bound 4) · strength 2 · Auto: IPOG() · 3 parameters · 12 combinations
excluded: 3 pairs forbidden, 2 impossible under the constraints; see report(cases)
    mode    solver  tol
 1  :fast   :none   0.001
 2  :exact  :none   1.0e-6
 3  :exact  :lu     1.0e-6
 4  :exact  :qr     1.0e-6
 5  :fast   :none   1.0e-6

Every case satisfies both rules, and every pair of values that some valid case can hold appears in a case. The second line says what was left out. Three pairs are forbidden by a rule that names them. Two more are impossible: no single rule covers both of their parameters, but together the rules leave no valid case that holds them. explain says why:

explain(space, (solver = :lu, tol = 1e-3))
infeasible: no valid case contains (solver = :lu, tol = 0.001); rules 1 and 2 together exclude it (rule 1: @require(mode == :exact || solver == :none); rule 2: exact mode needs a tight tolerance)

solver = :lu needs exact mode, and exact mode needs tol = 1e-6. You never write that third rule; the package finds it, and does not ask for a case that cannot exist. explain takes any partial or complete assignment:

explain(space, (mode = :fast, solver = :lu))
forbidden by rule 1 (@require(mode == :exact || solver == :none))
explain(space, (solver = :lu,))
completable, e.g. (mode = :exact, solver = :lu, tol = 1.0e-6)

A rule that names a parameter the space lacks is an error when the space is built, and the message names what went wrong:

try
    TestSpace((mode = [:fast, :exact], solver = [:none, :lu, :qr]);
              constraints = [@forbid(mode == :fast && solvr != :none)])
catch err
    showerror(stdout, err)
end
ArgumentError: rule 1 (@forbid(mode == :fast && solvr != :none)) names `solvr`, which is not a parameter of this space. The parameters are mode, solver. If `solvr` is a variable, write `$solvr` to use its value.

The test loop does not change, and it now runs only valid cases:

@testset "solve" begin
    @testset "$case" for case in cases
        x = solve(A, b; case...)
        @test norm(A * x - b) <= case.tol * norm(b)
    end
end
Test Summary: | Pass  Total  Time
solve         |    5      5  0.0s

Constraints describes exactly what a rule excludes and how the package decides that a combination is impossible.

Level 3: must-include cases and strength

Top up tests you already have

Suppose the suite already has two hand-written cases. coverage measures their interaction coverage, which pairs they hold (this is not line coverage):

handwritten = [(mode = :fast, solver = :none, tol = 1e-3),
               (mode = :exact, solver = :lu, tol = 1e-6)]
coverage(handwritten, space)
covers 6 of 11 feasible pairs, 5 missing: (mode = :exact, solver = :none), (mode = :exact, solver = :qr), (mode = :fast, tol = 1.0e-6), (solver = :none, tol = 1.0e-6), (solver = :qr, tol = 1.0e-6)
excluded: 3 pairs forbidden, 2 impossible under the constraints

Pass them as must_include. They come first, in the order given, and the design adds cases only for what they miss:

topped = all_pairs(space; must_include = handwritten)
5 cases (2 must-include, lower bound 4) · strength 2 · Auto: IPOG() · 3 parameters · 12 combinations
excluded: 3 pairs forbidden, 2 impossible under the constraints; see report(cases)
    mode    solver  tol
 1  :fast   :none   0.001
 2  :exact  :lu     1.0e-6
 3  :exact  :none   1.0e-6
 4  :exact  :qr     1.0e-6
 5  :fast   :none   1.0e-6

A must-include case may be partial; the package completes it with valid values:

all_pairs(space; must_include = [(mode = :fast,), (solver = :qr,)])
5 cases (2 must-include, lower bound 4) · strength 2 · Auto: IPOG() · 3 parameters · 12 combinations
excluded: 3 pairs forbidden, 2 impossible under the constraints; see report(cases)
    mode    solver  tol
 1  :fast   :none   0.001
 2  :exact  :qr     1.0e-6
 3  :exact  :none   1.0e-6
 4  :exact  :lu     1.0e-6
 5  :fast   :none   1.0e-6

A must-include case that breaks a rule is an error that names the case and the rule:

try
    all_pairs(space; must_include = [(mode = :fast, solver = :lu, tol = 1e-6)])
catch err
    showerror(stdout, err)
end
ArgumentError: must_include row 1, (mode = :fast, solver = :lu, tol = 1.0e-6), breaks rule 1 (@require(mode == :exact || solver == :none)) (contract §10.3)

Strength

all_pairs is strength 2. all_values is strength 1 (every value appears at least once), all_triples is strength 3, and covering takes any strength. With four parameters:

wide = (n = [1, 2, 3], level = ["low", "mid", "high"], tol = [1.0, 3.7, 4.9],
        kind = [:greedy, :relax, :optim])
all_triples(wide)
27 cases (minimal) · strength 3 · Auto: Construction() · 4 parameters · 81 combinations
     n  level   tol  kind
  1  1  "low"   1.0  :greedy
  2  2  "low"   1.0  :optim
  3  3  "low"   1.0  :relax
  4  1  "mid"   1.0  :optim
  5  2  "mid"   1.0  :relax
  6  3  "mid"   1.0  :greedy
  7  1  "high"  1.0  :relax
  8  2  "high"  1.0  :greedy
  9  3  "high"  1.0  :optim
 10  1  "low"   3.7  :optim
  ⋮  ⋮  ⋮       ⋮    ⋮
 23  2  "mid"   4.9  :optim
 24  3  "mid"   4.9  :relax
 25  1  "high"  4.9  :optim
 26  2  "high"  4.9  :relax
 27  3  "high"  4.9  :greedy

Often only a few parameters interact strongly. stronger raises the strength for a group of them and keeps pairs for the rest:

all_pairs(wide; stronger = [(:n, :tol, :kind) => 3])
27 cases (minimal) · strength 2, 3 within (n, tol, kind) · Auto: IPOG() · 4 parameters · 81 combinations
     n  level   tol  kind
  1  1  "low"   1.0  :greedy
  2  1  "mid"   3.7  :relax
  3  1  "high"  4.9  :optim
  4  2  "high"  1.0  :relax
  5  2  "low"   3.7  :greedy
  6  2  "mid"   4.9  :greedy
  7  3  "mid"   1.0  :optim
  8  3  "high"  3.7  :greedy
  9  3  "low"   4.9  :relax
 10  2  "low"   3.7  :optim
  ⋮  ⋮  ⋮       ⋮    ⋮
 23  3  "mid"   1.0  :relax
 24  3  "high"  3.7  :relax
 25  3  "low"   3.7  :optim
 26  3  "mid"   4.9  :greedy
 27  3  "high"  4.9  :optim

The group's three parameters have three values each, so its triples are all 27 of their combinations, and the pairs with level fit into those same cases. A positional call names the group by position, stronger = [(1, 3, 4) => 3], and several groups may be listed, such as stronger = [(3, 4, 5, 6) => 3, (25, 26, 27) => 3] for 30 parameters.

To go from pairs to triples without losing the cases you have, pass them as must_include:

triples = all_triples(space; must_include = cases)
coverage(triples)
covers 5 of 5 feasible triples
excluded: 7 triples forbidden

Excursions

When there is one known-good configuration, excursions varies it: every valid case within distance changed parameters of the base. This is not a covering design; it promises only that distance.

excursions(space; from = (mode = :exact, solver = :lu, tol = 1e-6), distance = 1)
3 cases · excursion, distance 1 from (mode = :exact, …) · 3 parameters · 12 combinations, 2 rows dropped
    mode    solver  tol
 1  :exact  :lu     1.0e-6
 2  :exact  :none   1.0e-6
 3  :exact  :qr     1.0e-6
Choosing an engine

Auto(), the default, builds cases with IPOG, one parameter at a time, and also builds the catalog's algebraic array where every parameter has the same number of values, keeping whichever design is smaller; neither uses randomness, and engine = IPOG() builds IPOG's design alone. Auto(goal = :compact) then removes rows with a reducer, for tests that are expensive to run. GND builds each case from random candidates drawn from a fixed seed, so it too gives the same cases on every run, and GND(seed = 7) gives a different design with the same guarantee. Every engine covers every feasible combination; none promises a minimum number of cases, but each result states a lower bound beside its count. See Engines.

all_pairs(space; engine = GND(seed = 7))
5 cases (lower bound 4) · strength 2 · GND seed 7 · 3 parameters · 12 combinations
excluded: 3 pairs forbidden, 2 impossible under the constraints; see report(cases)
    mode    solver  tol
 1  :exact  :none   1.0e-6
 2  :fast   :none   0.001
 3  :exact  :lu     1.0e-6
 4  :exact  :qr     1.0e-6
 5  :fast   :none   1.0e-6

Level 4: inspection

Before committing: design_sizes

design_sizes runs each strategy and says how many cases it gives and what it covers, so you can choose before writing the test:

design_sizes(wide)
strategy        cases   share  pairs  triples
full_factorial     81  100.0%  54/54  108/108  valid 81 of 81
covering(1)         3    3.7%  18/54   12/108
covering(2)         9   11.1%  54/54   36/108
covering(3)        27   33.3%  54/54  108/108
excursions(1)       9   11.1%  30/54   28/108
excursions(2)      33   40.7%  54/54   76/108
case counts are the rows each strategy produced with Auto, not lower bounds

The share column is each design's size as a share of the valid cases. The pairs and triples columns are what each design covers, so a strength-2 design also reports the triples it happens to hold.

Exporting cases

A result of named cases is a Tables.jl table:

using DataFrames
DataFrame(cases)
5×3 DataFrame
Rowmodesolvertol
SymbolSymbolFloat64
1fastnone0.001
2exactnone1.0e-6
3exactlu1.0e-6
4exactqr1.0e-6
5fastnone1.0e-6

github_matrix writes the cases as a GitHub Actions matrix, one job per case. See Plan a CI matrix.

github_matrix(cases)
{"include":[{"mode":"fast","solver":"none","tol":0.001},{"mode":"exact","solver":"none","tol":1.0e-6},{"mode":"exact","solver":"lu","tol":1.0e-6},{"mode":"exact","solver":"qr","tol":1.0e-6},{"mode":"fast","solver":"none","tol":1.0e-6}]}

Invalid inputs

Mark a value Invalid to test that the function rejects it. The design covers the ordinary values as before, then adds negative cases, each holding exactly one invalid value, marked !, so one error cannot hide another. A rule that reads tol does not apply to a case whose tol is invalid:

negative = TestSpace(
    (mode = [:fast, :exact], solver = [:none, :lu, :qr],
     tol = [1e-3, 1e-6, Invalid(0.0), Invalid(-1.0)]);
    constraints = [
        @require(mode == :exact || solver == :none),
        forbid((mode = :exact, tol = 1e-3); reason = "exact mode needs a tight tolerance"),
    ])
negcases = all_pairs(negative)
11 cases (lower bound 10) · strength 2 · Auto: IPOG() · 3 parameters · 24 combinations · 10 negative targets
excluded: 3 pairs forbidden, 2 impossible under the constraints; see report(cases)
      mode    solver  tol
  1   :fast   :none   0.001
  2   :exact  :none   1.0e-6
  3   :exact  :lu     1.0e-6
  4   :exact  :qr     1.0e-6
  5   :fast   :none   1.0e-6
  6!  :fast   :none   Invalid(0.0)
  7!  :exact  :lu     Invalid(0.0)
  8!  :exact  :qr     Invalid(0.0)
  9!  :fast   :none   Invalid(-1.0)
 10!  :exact  :lu     Invalid(-1.0)
 11!  :exact  :qr     Invalid(-1.0)

hasinvalid tells a test body which kind of case it has, and the wrapped value of an Invalid x is x.value:

unwrap(x) = x isa Invalid ? x.value : x

@testset "solve with bad tolerances" begin
    for case in negcases
        args = map(unwrap, case)
        if hasinvalid(case)
            @test_throws ArgumentError solve(A, b; args...)
        else
            x = solve(A, b; args...)
            @test norm(A * x - b) <= args.tol * norm(b)
        end
    end
end
Test Summary:             | Pass  Total  Time
solve with bad tolerances |   11     11  0.1s

See Test invalid inputs.

Classes of values

A Partition names a class of values and draws a concrete one when the test runs. The design, the rules and coverage all see the name; realize draws the values from a random number generator you pass:

using Random

sized = TestSpace((n = [Partition(:small, rng -> rand(rng, 2:5)),
                        Partition(:large, rng -> rand(rng, 200:300))],
                   solver = [:lu, :qr]))
labeled = all_pairs(sized)
4 cases (minimal) · strength 2 · Auto: Construction() · 2 parameters · 4 combinations, 4 valid
    n                  solver
 1  Partition(:small)  :lu
 2  Partition(:large)  :qr
 3  Partition(:large)  :lu
 4  Partition(:small)  :qr
realize(labeled; rng = Xoshiro(1))
4-element Vector{@NamedTuple{n::Int64, solver::Symbol}}:
 (n = 2, solver = :lu)
 (n = 235, solver = :qr)
 (n = 270, solver = :lu)
 (n = 4, solver = :qr)

See Combine with property-based testing.

After a failure

diagnose takes the cases and which of them passed, and ranks the combinations that appear only in failing cases. Here a bug fails every case with method = :newton on a sparse matrix:

suite = all_pairs((n = [10, 100, 1000, 10000], method = [:newton, :bicg, :gmres],
                   tol = [1e-3, 1e-6], sparse = [false, true]))
passed = [!(c.method == :newton && c.sparse) for c in suite]
diagnose(suite, passed)
3 failures of 12 cases; 6 suspects in 5 groups (hypotheses, not proof)
1. (method = :newton, sparse = true) — in 3 of 3 failures
2. (method = :newton, tol = 0.001) — in 2 of 3 failures
3. (n = 10, method = :newton) — in 1 of 3 failures
4. (n = 100, method = :newton) — in 1 of 3 failures
5. (n = 10000, method = :newton) — in 1 of 3 failures — same failures as (n = 10000, sparse = true)

The true cause ranks first. The ranking is a set of hypotheses, not a proof: the other suspects appeared only in failing cases, so nothing yet says whether they work. followups proposes a case to separate each one. See Diagnose a failure.

Check the claim: coverage and report

coverage of a result measures it again, from its rows alone, at the strength it was built for. iscomplete is the one-line check that every feasible combination is covered and none is unresolved:

coverage(cases)
covers 11 of 11 feasible pairs
excluded: 3 pairs forbidden, 2 impossible under the constraints
iscomplete(coverage(cases))
true

report is the full account: the guarantee, its coverage figures measured from the rows; each excluded combination with the rules that exclude it; the triples the pairwise design covers as a bonus; how much the first cases cover, for a suite that runs only some of them; and the seed.

report(cases)
5 cases cover all 11 feasible pairs of a 12-combination space (3 pairs forbidden, 2 impossible under the constraints)
excluded:
  (mode = :fast, solver = :lu): forbidden by rule 1 (@require(mode == :exact || solver == :none))
  (mode = :fast, solver = :qr): forbidden by rule 1 (@require(mode == :exact || solver == :none))
  (mode = :exact, tol = 0.001): forbidden by rule 2 (exact mode needs a tight tolerance)
  (solver = :lu, tol = 0.001): impossible because rules 1 and 2 combine (rule 1: @require(mode == :exact || solver == :none); rule 2: exact mode needs a tight tolerance)
  (solver = :qr, tol = 0.001): impossible because rules 1 and 2 combine (rule 1: @require(mode == :exact || solver == :none); rule 2: exact mode needs a tight tolerance)
size: 5 cases; lower bound 4: the 4 feasible combinations of mode and solver need a case each
bonus: 5 of 5 feasible triples covered
prefix curve:
  first 1 of 5 cover 27% (3 of 11)
  first 2 of 5 cover 54% (6 of 11)
  first 3 of 5 cover 72% (8 of 11)
  first 4 of 5 cover 90% (10 of 11)
  first 5 of 5 cover 100% (11 of 11)
seed: none (Auto uses no randomness)

A worked example: generic code across types

Julia code is often generic, so one function is really one function per combination of argument types, and the types are parameters to test. This example tests a Hilbert curve index, which maps a point (x, y) on a grid to its position z along a space-filling curve, and back. The code has many branches, and its author worried about signed and unsigned integers, integer sizes, and zero-based against one-based counting.

struct Simple2D{T} end

function encode_hilbert_zero(::Simple2D{T}, X::Vector{A})::T where {A, T}
    x = X[1]
    y = X[2]
    z = zero(T)
    if x == zero(A) && y == zero(A)
        return z
    end
    rmin = convert(Int, floor(log2(max(x, y))) + 1)
    w = one(A) << (rmin - 1)
    while rmin > 0
        z <<= 2
        if rmin & 1 == 1  # odd
            if x < w
                if y >= w
                    x, y = (y - w, x)
                    z += one(T)
                end
            else
                if y < w
                    x, y = ((w << 1) - x - one(w), w - y - one(w))
                    z += T(3)
                else
                    x, y = (y - w, x - w)
                    z += T(2)
                end
            end
        else  # even
            if x < w
                if y >= w
                    x, y = (w - x - one(w), (w << 1) - y - one(w))
                    z += T(3)
                end
            else
                if y < w
                    x, y = (y, x - w)
                    z += one(T)
                else
                    x, y = (y - w, x - w)
                    z += T(2)
                end
            end
        end
        rmin -= 1
        w >>= 1
    end
    z
end

function decode_hilbert_zero!(::Simple2D{T}, X::Vector{A}, z::T) where {A, T}
    r = z & T(3)
    x, y = [(zero(A), zero(A)), (zero(A), one(A)), (one(A), one(A)), (one(A), zero(A))][r + 1]
    z >>= 2
    rmin = 2
    w = one(A) << 1
    while z > zero(T)
        r = z & T(3)
        if rmin & 1 != 0
            if r == 1
                x, y = (y, x + w)
            elseif r == 2
                x, y = (y + w, x + w)
            elseif r == 3
                x, y = ((w << 1) - x - one(A), w - y - one(A))
            end
        else
            if r == 1
                x, y = (y + w, x)
            elseif r == 2
                x, y = (y + w, x + w)
            elseif r == 3
                x, y = (w - x - one(A), (w << 1) - y - one(A))
            end
        end
        z >>= 2
        rmin += 1
        w <<= 1
    end
    X[1] = x
    X[2] = y
end

# The same functions, counting from one.
encode_hilbert(gg::Simple2D{T}, X::Vector{A}) where {A, T} =
    encode_hilbert_zero(gg, X .- one(A)) + one(T)

function decode_hilbert!(gg::Simple2D{T}, X::Vector{A}, h::T) where {A, T}
    decode_hilbert_zero!(gg, X, h - one(T))
    X .+= one(A)
end

The parameters are the integer type of the coordinates, the integer type of the index, whether counting starts at zero or one, the length of the coordinate vector (the code accepts one longer than it needs), and the number of bits per coordinate. A domain can hold types. Two rules say that the index and the coordinates must have room for the bits the curve uses:

using UnitTestDesign, Test

ints = [Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Int128, UInt128]
space = TestSpace(
    (axis = ints, index = ints, offset = [0, 1], dims = [2, 3, 4], bits = [2, 3, 4, 5]);
    constraints = [
        @forbid(bits * dims > log2(typemax(index))),
        @forbid(bits * dims > log2(typemax(axis))),
    ])
cases = all_pairs(space)
100 cases (minimal) · strength 2 · Auto: IPOG() · 5 parameters · 2400 combinations
excluded: 12 pairs impossible under the constraints; see report(cases)
      axis     index    offset  dims  bits
   1  UInt128  UInt128  0       2     2
   2  Int128   UInt128  1       3     3
   3  UInt64   UInt128  0       4     4
   4  Int64    UInt128  1       2     5
   5  UInt32   UInt128  0       3     2
   6  Int32    UInt128  1       4     2
   7  UInt16   UInt128  0       2     2
   8  Int16    UInt128  0       2     2
   9  UInt8    UInt128  0       2     2
  10  Int8     UInt128  0       2     2
   ⋮  ⋮        ⋮        ⋮       ⋮     ⋮
  96  Int32    Int8     1       2     3
  97  UInt16   Int8     0       3     2
  98  Int16    Int8     1       2     3
  99  UInt8    Int8     0       3     2
 100  Int8     Int8     1       2     3

typemax and log2 are called, so the macro reads them as functions; index, axis, bits and dims are parameters. The test decodes the first and last few indices and encodes them again, which must give back the same index, of the same type:

@testset "Hilbert round trip" begin
    for (A, I, C, D, B) in cases
        gg = Simple2D{I}()
        last = (one(I) << (B * D)) - one(I) + I(C)
        mid = one(I) << (B * D - 1)
        X = zeros(A, D)
        for hl in vcat(C:min(mid, 5), max(mid + 1, last - 5):last)
            h = I(hl)
            if C == 0
                decode_hilbert_zero!(gg, X, h)
                @test encode_hilbert_zero(gg, X) === h
            else
                decode_hilbert!(gg, X, h)
                @test encode_hilbert(gg, X) === h
            end
        end
    end
end
Test Summary:      | Pass  Total  Time
Hilbert round trip | 1150   1150  5.1s

When this test was first written, there were no rules: it called all_pairs on the five lists and skipped a case with continue when the types were too small. coverage shows what that cost. Of the 100 cases, 64 ran, and they missed 62 pairs that a valid case could have held:

unconstrained = all_pairs(ints, ints, [0, 1], [2, 3, 4], [2, 3, 4, 5])
ran = [row for row in unconstrained if isallowed(space, row)]
coverage(ran, space)
covers 232 of 294 feasible pairs, 62 missing: (axis = UInt8, index = Int8), (axis = Int16, index = Int8), (axis = UInt16, index = Int8), (axis = Int8, index = UInt8), (axis = Int16, index = UInt8), (axis = UInt16, index = UInt8), (axis = Int32, index = UInt8), (axis = UInt32, index = UInt8), (axis = Int64, index = UInt8), (axis = UInt64, index = UInt8), and 52 more
excluded: 12 pairs impossible under the constraints

With the rules in the space, every one of those pairs is in a case that runs, and the report lists the pairs that no case can hold, each with its rule:

report(cases)
100 cases cover all 294 feasible pairs of a 2400-combination space (12 pairs impossible under the constraints)
excluded:
  (axis = Int8, dims = 4): impossible because of rule 2 (@forbid(bits * dims > log2(typemax(axis))))
  (axis = UInt8, dims = 4): impossible because of rule 2 (@forbid(bits * dims > log2(typemax(axis))))
  (axis = Int8, bits = 4): impossible because of rule 2 (@forbid(bits * dims > log2(typemax(axis))))
  (axis = UInt8, bits = 4): impossible because of rule 2 (@forbid(bits * dims > log2(typemax(axis))))
  (axis = Int8, bits = 5): impossible because of rule 2 (@forbid(bits * dims > log2(typemax(axis))))
  (axis = UInt8, bits = 5): impossible because of rule 2 (@forbid(bits * dims > log2(typemax(axis))))
  (index = Int8, dims = 4): impossible because of rule 1 (@forbid(bits * dims > log2(typemax(index))))
  (index = UInt8, dims = 4): impossible because of rule 1 (@forbid(bits * dims > log2(typemax(index))))
  (index = Int8, bits = 4): impossible because of rule 1 (@forbid(bits * dims > log2(typemax(index))))
  (index = UInt8, bits = 4): impossible because of rule 1 (@forbid(bits * dims > log2(typemax(index))))
  (index = Int8, bits = 5): impossible because of rule 1 (@forbid(bits * dims > log2(typemax(index))))
  (index = UInt8, bits = 5): impossible because of rule 1 (@forbid(bits * dims > log2(typemax(index))))
size: 100 cases, minimal: the 10 × 10 = 100 combinations of axis and index need a case each
bonus: 618 of 1266 feasible triples covered
prefix curve:
  first 20 of 100 cover 39% (117 of 294)
  first 40 of 100 cover 62% (184 of 294)
  first 60 of 100 cover 76% (224 of 294)
  first 80 of 100 cover 89% (262 of 294)
  first 100 of 100 cover 100% (294 of 294)
seed: none (Auto uses no randomness)