loopty

loopy, with types: the algorithm's theorem, the kernel and its transforms in one Python file, and one ledger that says how each of their claims is known.

If you write loopy kernels, this is what changes. A kernel is a typed Python function, which runs on numpy as it stands. A schedule is a chain of the transforms you already use, split, tag, prioritize, fuse, with contract to keep what passes between fused loops in a window. The theorem the algorithm rests on is a lanky claim in the same file, written as a Python annotation. Then one command prints every claim the file makes and what settled it:

• isl decides that every access is in bounds, that no two instances write one cell, and that every transform keeps every dependence, and refuses an illegal transform with the two instances it would reorder;
• Lean checks: each fact isl decided is stated again as a Lean theorem, from loopty's own description of the kernel and not from what isl computed, and proved there, with isl's hints checked rather than trusted; with Mathlib, Lean also proves the algorithm's theorems, and checks each where a kernel cites it, against what the kernel stores;
• the runs are compared: the compiled run with the native one, and the native one with the upstream code the kernel was written after.

You write no Lean, no isl and no loopy kernel. --show-lean prints the Lean, and loopty run --emit-code the generated code, for whoever wants to read them.

One file: sumpy's L2P

examples/l2p.py is sumpy's L2P for the 2-D Laplace kernel (E2PFromSingleBox), abridged here to its three layers: the theorem sumpy's compressed Taylor expansion rests on (lanky#86), the kernel that rebuilds each box's triangle of derivatives and cites it, and the schedules.

@theorem
def compressed_taylor(
    D: Fn[Nat, Nat, Real],
    R: Fn[Nat, Nat, Real],
    pde: all(D[a + 2, b] + D[a, b + 2] == 0 for a in Nat for b in Nat),
    stored: all(R[a, b] == D[a, b] for a in Nat for b in Nat if a < 2),
    recurrence: all(R[a + 2, b] == -R[a, b + 2] for a in Nat for b in Nat),
) -> all(R[a, b] == D[a, b] for a in Nat for b in Nat):
    """What sumpy stores of a harmonic table, reconstructed, is the table."""


@kernel(uses=[on_every_triangle])  # compressed_taylor, box by box on a + b <= P
def reconstruct(
    C: Arr[Where[box: Fin[nb], a: Fin[2], b: Fin[Q], a + b <= P], Real],
    R: Arr[Where[box: Fin[nb], a: Fin[Q], b: Fin[Q], a + b <= P], Real],
) -> ...:  # the copies and the recurrence: two of the theorem's hypotheses
    for box, a, b in R.dom:
        if a < 2:
            R[box, a, b] = C[box, a, b]
        else:
            R[box, a, b] = -R[box, a - 2, b + 2]


# ... unpack, which reads sumpy's coefficients through its storage index,
# evaluate, the Taylor sum at each target, and the program l2p that runs the three

fused = Schedule(l2p, target="c").fuse(  # the three calls' loops over boxes
    l2p.unpack.box, l2p.reconstruct.box, l2p.evaluate.box
)
contracted = fused.contract(l2p.R).contract(l2p.C)  # one box's triangle at a time
on_a_gpu = fused.retarget("opencl").on_groups(l2p.unpack.box)  # sumpy's own

A schedule names what it moves by where the file defines it: l2p.R is the table the program l2p makes, l2p.unpack.box the loop over boxes of the kernel unpack where l2p calls it, and reconstruct.b the loop over b in reconstruct itself. A name is checked when it is written, against the traced kernel, and every message, schedule and fact prints the same paths.

python examples/l2p.py runs it on numpy, refuses the interchange of the triangle's loops with the pair of instances it would reorder, refuses a window of one row with the read it would break, finds the window that keeps one box's triangle, and runs the schedule compiled:

$ python examples/l2p.py
...
refused: interchange(reconstruct.b, reconstruct.a) illegal: R[box, a, b] = C[box, a, b] (l2p.py:271) at box=0, a=1, b=9 writes R[0, 1, 9], read by R[box, a, b] = -R[box, a - 2, b + 2] (l2p.py:273) at box=0, a=3, b=7, scheduled earlier (at nb=1, isl's choice)

refused: contract(l2p.R, lambda box, a, b: b) illegal: R[box, a, b] = -R[box, a - 2, b + 2] (l2p.py:273) at box=0, a=9, b=0 reads R[0, 7, 2], which the same statement at box=0, a=7, b=2 stored in slot 2, after the same statement at box=0, a=8, b=2 stored R[0, 8, 2] there (at nb=1, isl's choice); contract(l2p.R) finds the smallest window that keeps every read, lambda box, a, b: (a, b)

schedule: Schedule(l2p, target='c').fuse(l2p.unpack.box, l2p.reconstruct.box, l2p.evaluate.box).contract(l2p.R, lambda box, a, b: (a, b)).contract(l2p.C, lambda box, a, b: (a, b))
...
  decided  isl  no read of l2p.R sees its slot overwritten after contract(l2p.R, lambda box, a, b: (a, b)): l2p.R kept in 121 cells, R[box, a, b] in slot (a, b)
...
  pot: difference 0 within 1.1e-06 (approx) -> tested
...

lanky check examples/l2p.py, on a machine with Lean and Mathlib, says how each claim is known and who settled it. Six of its 56 rows:

$ lanky check examples/l2p.py
STATUS   EFFECTIVE  BY           WHERE       OWNER              STATEMENT
-------  ---------  -----------  ----------  -----------------  ----------------------------------------------------------------------------------------------------
proved   tested     lean         l2p.py:129  compressed_taylor  D : Fn[Nat, Nat, Real], R : Fn[Nat, Nat, Real]
                                                                | pde: forall a in Nat, b in Nat. D[a + 2, b] + D[a, b + 2] == 0
                                                                | stored: forall a in Nat, b in Nat where a < 2. R[a, b] == D[a, b]
                                                                | recurrence: forall a in Nat, b in Nat. R[a + 2, b] == -R[a, b + 2]
                                                                |- forall a in Nat, b in Nat. R[a, b] == D[a, b]
...
proved   proved     lean         l2p.py:198  storage_index      storage_index stores a cell at every address of Fin(21)
...
proved   tested     lean         l2p.py:246  reconstruct        by on_every_triangle, for every D with pde and stored:
                                                                  forall i in Fin(nb), a in Nat, b in Nat where a + b <= 10. R[i, a, b] == D[i, a, b]
...
proved   proved     lean         l2p.py:334  l2p                no read of l2p.R sees its slot overwritten after contract(l2p.R, lambda box, a, b: (a, b)):
                                                                  l2p.R kept in 121 cells, R[box, a, b] in slot (a, b)
...
proved   proved     lean         l2p.py:342  l2p                the order after on_groups(l2p.unpack.box) runs every dependence of l2p forward, within one work item
tested   tested     upstream     l2p.py:464  l2p                the native run of l2p agrees with sumpy_l2p in pot within 1e-14 of the largest value,
                                                                  on its example inputs

56 facts: 3 decided, 46 proved, 7 tested

LEAN proved all 44 facts isl decided. Lean's statements are loopty's, from its own description, not isl's results, and isl's hints were checked, not trusted. --show-lean prints the Lean.

STATUS is how a claim was settled, and EFFECTIVE what it is worth given what it rests on. The theorem is proved for every order, by an induction whose steps sympy found and Lean checked, and is worth tested because lanky reads its annotations by running them. The application rests on the theorem and on reconstruct's postcondition, which isl decided from what the kernel stores and Lean proved again. The window's fact and the GPU schedule's cast rest on nothing weaker, and are worth proved. The last row is the agreement with sumpy's own kernel. The three rows left decided are the class of the Taylor sum, by type, and two postconditions restated after their calls, by the call.

Without Lean the same file reads decided isl where this reads proved lean, and the theorems read assumed, with a line under the table for each piece that is missing and the command that installs it. That ledger, which CI checks, is in examples/README.md.

The demos

Five files, each written after a kernel of sumpy, pytential, grudge or firedrake, codes whose kernels loopy generates. Each runs three ways, python FILE, lanky check FILE and loopty run FILE, and examples/README.md has what each prints. "Proved" below is with Lean and Mathlib.

The timings are benchmarks/'s drivers, run on a GPU and on CPUs beside the upstream code; their numbers describe the machines, so they are not kept here (benchmarks/README.md).

Run it

The demos run on the 0.4.0 releases of loopty and lanky on PyPI. The examples are in the repository and not in the package, so clone the release's tag for them. On Python 3.12 or 3.13 with a C compiler:

python3.13 -m venv .venv && . .venv/bin/activate
pip install loopty "lanky[lean,cas]"
git clone --branch v0.4.0 https://github.com/xywei/loopty.git && cd loopty
python examples/l2p.py
lanky check examples/l2p.py
loopty run examples/l2p.py

That much runs without Lean, and the ledger says what is missing, with the command for each piece. The other demos run the same three ways, and examples/README.md has what each prints. Its commands start with uv run, which is for a checkout with lanky cloned beside it (Install); in this clone, run them without it. Lean comes from elan, at the version the Lean REPL supports, and Mathlib, for the theorems about real numbers, is one command that fetches its compiled files:

elan default leanprover/lean4:v4.29.1
python -m lanky.mathlib ~/mathlib    # about 7 GB
export LANKY_LEAN_MATHLIB=~/mathlib

The first check with Lean builds Lean's REPL, which takes about a minute and the network, and the first check of a file asks Lean about everything it states. Checking the file again reads Lean's answers from a cache and starts no Lean, also for a question the first check ran out of time on. Install has the rest.

Where the limits are

• loopty checks the kernels it traces. A loopy kernel built elsewhere, sumpy's or TSFC's, is not checked, nor a transform applied to one. What ties a demo to the upstream code is a tested agreement row, not a proof.
• Some steps stay trusted: the printer from loopty's description to Lean, the Lean kernel and omega, and the step from the checked model to the code loopy and the compiler generate, which the comparison of the compiled run with the native one covers by sampling (docs/lean-checks.md). lanky checks its reading of a theorem's annotations by running them at drawn points, so a proved theorem is worth tested in the EFFECTIVE column.
• A proof about reals is about the reals. A theorem, or a postcondition Lean proves with Mathlib, is about the kernel's term read over the real or complex numbers, not about the floats it computes, and the row says so (#282). Without Mathlib, about 7 GB, those rows read assumed or tested.
• On a device, no local memory or barriers yet (#256), so a kernel that stages passes in a work group's local memory, as the fastest sum-factorized and libParanumal kernels do, cannot be written. A window each work item keeps is decided by isl and not yet proved by Lean (#295). CI runs no device.
• loopy 2025.2. loopty is held to the loopy release it is tested with. Its port to loopy main, on namedisl, is a draft (#180) until loopy releases it.
• Known limits lists the open issues you are most likely to meet, each with what to do about it.

A first kernel, row by row

The ledgers from here on are what CI prints, without Lean. Where Lean is installed, the rows below that read decided isl read proved lean, but for the few facts that stay isl's alone (Known limits).

A loopty kernel is a Python function. Run it under plain python and it computes on numpy. Trace it and you get a typed term whose obligations go into a ledger. Transform it and every step is checked before it is applied.

@kernel
def spmv(
    cnt: Arr[Fin[n], Nat],
    col: Arr[Fin[n], Fin[cnt], Fin[m]],
    val: Arr[Fin[n], Fin[cnt], Real],
    x: Arr[Fin[m], Real],
    y: Arr[Fin[n], Real],
):
    for r in y.dom:
        y[r] = reduce_sum(val[r, j] * x[col[r, j]] for j in val.dom[r])

reduce_sum is Loopty's reduction frontend: it binds j over the symbolic fiber and traces the operation as an ISL-backed loopty.term.Reduction.

val: Arr[Fin[n], Fin[cnt], Real] is a ragged shape: for each of the n rows, cnt[r] entries. The second axis names the counts array, which is what makes it a dependent sum rather than a rectangle, and val.dom[r] is the fiber over that row. Nothing in the kernel mentions offsets or flat storage; that is layout, and layout is the compiler's business.

Its three commands:

$ python examples/spmv.py            # runs natively on numpy; the theorem is tested
$ lanky check examples/spmv.py       # the ledger: every obligation and who decided it
$ loopty run examples/spmv.py        # lowers through loopy, compiles, runs, compares

and the ledger lanky check prints, abridged to the rows discussed here (each row verbatim, and checked by scripts/refresh_example_outputs.py):

$ lanky check examples/spmv.py
STATUS   BY             WHERE        OWNER          STATEMENT
-------  -------------  -----------  -------------  ---------------------------------------------------------------------------------------------------
decided  isl            spmv.py:81   scan           off[0] is in bounds for every instance of off[0] = 0 (spmv.py:81)
decided  isl            spmv.py:83   scan           off[r + 1] is in bounds for every instance of off[r + 1] = off[r] + cnt[r] (spmv.py:83)
decided  isl            spmv.py:81   scan           distinct instances of off[0] = 0 (spmv.py:81) write distinct cells of off
decided  isl            spmv.py:81   scan           no cell of off is stored twice, or read before its store
decided  isl            spmv.py:71   scan           off[0] == 0 and (forall r in Fin(n). off[r + 1] == off[r] + cnt[r])
tested   property-test  spmv.py:86   scan_monotone  n : Nat, cnt : Fn[Fin(n), Nat], off : Fn[Fin(n + 1), Nat]
                                                    | h0: off(0) == 0
                                                    | hs: forall r in Fin(n). off(r + 1) == off(r) + cnt(r)
                                                    |- forall a in Fin(n + 1), b in Fin(n + 1) where a <= b. off(a) <= off(b)
decided  isl            spmv.py:114  spmv           y[r] is in bounds for every instance of y[r] = reduce_sum(...) (spmv.py:114)
decided  type           spmv.py:114  spmv           x[col[r, j]] is in bounds by type (col[r, j] : Fin(m))
decided  isl            spmv.py:114  spmv           distinct instances of y[r] = reduce_sum(...) (spmv.py:114) write distinct cells of y
decided  type           spmv.py:114  spmv           the accumulation into y[r] over j is approx
tested   interpreter    spmv.py:104  spmv           the traced term computes what the body computes
decided  isl            spmv.py:195  spmv           the order after split(spmv.j, 2) runs every dependence of spmv forward
decided  call           spmv.py:117  solve          after scan(...) in solve: off[0] == 0 and (forall r in Fin(n). off[r + 1] == off[r] + cnt[r])
...
27 facts: 22 decided, 5 tested

WITHOUT LEAN: scan_monotone and 18 more rows are not proved. To have Lean prove what it can: pip install "lanky[lean]", and Lean v4.29.1 from elan (elan default leanprover/lean4:v4.29.1).

Look at the x[col[r, j]] row, and at what decided it. That indirection is the one obligation in a sparse product a polyhedral checker cannot settle on its own, and here it was settled by the type: the entries of col are points of Fin[m] and x has m cells, so the shape of the data discharges it and isl is never called.

And look at the interpreter row, the one fact that is about the trace rather than about the term. Every other row is a claim about what tracing recorded; this one checks that the record is the body. The term is run by an interpreter of its own, with numpy semantics, and compared with the body run natively, on the file's example inputs and on inputs drawn from the declared types. A body that kept state where tracing does not look would be refuted here, with the input and the first cell that differs.

scan's postcondition is decided by isl, from what the kernel's statements store. No cell of off is stored twice, or read before its store, which is the row above it, so what off[r + 1] = off[r] + cnt[r] stores at r is what off[r + 1] holds when scan returns, and the claim is that equation at every r. It is also evaluated at what every native run of scan left in its arguments, on the file's example inputs and on drawn ones, and held after each, which its provenance keeps. The split(j, 2) row is a cast of the schedule the file holds, rows_split, decided as the step was taken: lanky check lists the casts of a file's schedules beside the facts of the kernel they schedule. The last row belongs to solve, the @program that runs scan and then spmv. It restates scan's postcondition in the program's scope, decided by the call, and rests on scan's own postcondition fact, postcondition:spmv.scan@69, so it is worth what that fact is worth. The id is the kind, then scan's definition, the module the file's path gives it, its name and its line, so a scan imported from another file, or defined twice, is never taken for this one. loopty run compiles solve too, as one kernel: its term is the two kernels' statements in call order, and the compiled program is compared with the program run natively.

A postcondition is also what a later call may assume. A kernel's requirements on its inputs are its argument types, and where an earlier call of a program wrote the array, the requirement is decided under what held when the call was made, or checked by the compiled program between the two calls; examples/travel.py shows both.

And a transformation is a cast, checked before it is applied:

$ python examples/stencil_skew.py
...
Schedule(jacobi).tile(jacobi.t, jacobi.i, 8, 8) ->
  IllegalCast: tile(jacobi.t, jacobi.i, 8, 8) illegal: u[t + 1, i] = (u[t, i - 1] + u[t, i + 1]) / 2 (stencil_skew.py:61) at t=0, i=8 writes u[1, 8], read by the same statement at t=1, i=7, scheduled earlier (at nt=16, nx=16, as hinted)
...
accepted: Schedule(jacobi, target='c').skew(jacobi.i, by=jacobi.t).tile(jacobi.t, jacobi.i, 8, 8)
  loop nest: jacobi.t.outer jacobi.i.outer jacobi.t.inner jacobi.i.inner
  decided  isl  skew(jacobi.i, by=jacobi.t) renames the instances of jacobi one for one
  decided  isl  the order after skew(jacobi.i, by=jacobi.t) runs every dependence of jacobi forward
  decided  isl  tile(jacobi.t, jacobi.i, 8, 8) renames the instances of jacobi one for one
  decided  isl  the order after tile(jacobi.t, jacobi.i, 8, 8) runs every dependence of jacobi forward
...

The rejection names two real instances of the kernel, not an empty set or a failed pattern match: the statement as the file writes it, its line, and each instance's loop variables. It is the pair a person would find by hand, at sizes the message states, because which violating pair isl picks depends on them.

And an array's index set need not be a box. One value per pair of particles is an array over the lower triangle, and the triangle is its type:

@kernel
def pairs(
    x: Arr[Fin[n], Real],
    y: Arr[Fin[n], Real],
    q: Arr[Fin[n], Real],
    f: Arr[Where[i: Fin[n], j: Fin[n], j < i], Real],
    e: Arr[Fin[n], Real],
):
    for i, j in f.dom:
        dx = x[i] - x[j]
        dy = y[i] - y[j]
        f[i, j] = q[i] * q[j] / (1.0 + dx * dx + dy * dy)
    for p in e.dom:
        e[p] = reduce_sum(f[p, j] for j in f.dom[p]) + reduce_sum(
            f[k, p] for k in e.dom if k > p
        )

Where takes binders, written as slices, and then the constraints that cut their box; Sigma[i: Fin[n], Fin[i + 1]] is a sum with affine fibers, and Fin[n] + Fin[m] a union of pieces. Every in-bounds obligation is decided over the exact set:

$ lanky check examples/pairs.py
STATUS   BY           WHERE        OWNER  STATEMENT
-------  -----------  -----------  -----  ----------------------------------------------------------------------------------------------------
decided  isl          pairs.py:77  pairs  f[i, j] is in bounds for every instance of f[i, j] = q[i] * q[j] /
                                            (1.0 + dx * dx + dy * dy) (pairs.py:77)
...
decided  isl          pairs.py:79  pairs  f[p, j] is in bounds for every instance of e[p] = reduce_sum(...) + reduce_sum(...) (pairs.py:79)
decided  isl          pairs.py:79  pairs  f[k, p] is in bounds for every instance of e[p] = reduce_sum(...) + reduce_sum(...) (pairs.py:79)
...
17 facts: 16 decided, 1 tested

f[k, p] is read under k > p, which makes (k, p) a point of the triangle, while f[p, p] would be refused although the n x n box around the triangle has the cell. How f is stored is a separate choice that changes no fact: Schedule(pairs) keeps it in that box, and Schedule(pairs).pack(pairs.f) keeps its cells and no others, row after row, read as f[off_f[i] + j] through a table of row starts. Both compiled runs agree with the native one.

What nothing else does

• The ragged shape is a type, and it is the same type isl reasons about. A CSR matrix is a dependent sum, Arr[Fin[n], Fin[cnt], Real]. In-bounds and disjointness come from the shape, so a fiber loop is checked rather than trusted. Systems that model sparsity as offset arithmetic have to prove what this states.

• An indirection can be in bounds by type. col: Arr[..., Fin[m]] makes x[col[r, j]] sound with no proof obligation at all. This is where the polyhedral model normally gives up and inserts a runtime check or a "trust me". The type is what discharges it, so the executor is what enforces the type: every argument is checked against its declared element sort on the way in, and a col entry of -1 or of m is a ValueError naming the cell rather than an address outside x.

• Transformations are casts with witnesses. Every split, tile, interchange, skew and realize states its reindexing as an isl map, and affine takes the map from you, any injective affine one, the diamond (t, i) -> (t + i, t - i) included, or a map per statement. Each is checked for bijectivity on statement instances and for monotonicity on the dependence relation. A failure prints two instances and the array cell between them, before any code is generated. Fusing two kernels a program calls is one of them, a map per statement that runs both in one loop, and so the storage of what passes between them is a question of the schedule and not of the kernels: computed where it is read, kept in full, or kept in the cells live at once, with a checked fact that no read sees its slot overwritten.

• Reassociation is visible in the type. Splitting an accumulation and summing the pieces is not free on floating point. A trace reads the accumulation's class off what it sums, realize(spmv.y, tree=True) is what lowers it to reassoc, and the differential comparison against the native run widens its tolerance because of that fact and not because someone chose a number. The tolerance is per element and local: exact is bitwise, and reassoc and approx ask that |got - want| <= eps_class * (|want| + 1) at every cell, so it is the accuracy claimed for that cell and does not grow with the size of the output. Over an exact accumulation the same cast is refused, and a kernel with an exact output is compiled with floating-point contraction off, so a fused multiply-add cannot change its last bit.

• An assembly's race is named by the mesh. PyOP2 adds each cell's values into the cell's nodes, y[cell[c, k]] += yl[c, k], and two cells that share a node add into one cell of y. That is an update: either order gives the sum, so the kernel owes the class of the sum and not distinct writes, and a schedule may run two updates of one cell in another order but not at the same time. A split of the cells onto work items is refused with the two updates that collide on the example inputs, the node PyOP2's two triangles share, not with a pair isl picked out of every cell the map could reach, and accepted once the updates are atomic:

  cells = mass_action.c
  Schedule(mass_action, target="opencl").atomic(mass_action.y).split(
      cells, 32
  ).tag({cells.outer: "g.0", cells.inner: "l.0"})
  

• Legal and buildable are different questions, and both are answered. A transformation can preserve the meaning of a program and still be one the backend cannot generate. Every accepted step is asked whether the target can build it, and a failure is a refuted fact of kind buildable decided by loopy-target, with the limit in words as its reason, rather than a LoopyError thrown from inside code generation several steps later.

• The reference implementation is the kernel. The same body runs on numpy under plain python and traces to the term loopy compiles, so the differential test compares a program with itself rather than with a second implementation. That the trace is the body is checked too, not assumed: every kernel and every program has a trace-faithful fact, the traced term interpreted and compared with the native run, bit for bit when the output is exact.

• The index set of an argument is a type, and not a box. An array over the lower triangle, a band, a sum with affine fibers or a union of pieces is written as that set, Arr[Where[i: Fin[n], j: Fin[n], j < i], Real], and its in-bounds obligations are decided over the set itself, so a cell of the bounding box outside it is refused. Where the cells are kept is a layout, boxed or packed, chosen by a schedule step, and changes no fact.

• A layout is a checked map. Where an array's cells are stored is a Python function from a cell's indices to its address, and that it stores every cell inside its buffer, no two at one address and a cell at every address is a fact isl decides and Lean proves, as is that two layouts store every cell alike: sumpy's compressed-Taylor storage index, and grudge's dim_tags N3,N0,N1,N2 against modepy's reshape of the DOFArray (examples/storage_layouts.py), at every order in examples/hexahedra.py, where isl tests it at a few orders and Lean proves it for every one. loopy's shaped arrays (loopy #1053) take a layout's injectivity on trust.

• What isl decides, Lean proves, one dependence class or one place at a time, from loopty's description of the kernel and never from a set isl computed, so the two answers are independent (the top of this page); docs/lean-checks.md has the details and the coverage.

• A kernel is a rule with a domain of its own. A kernel body calls a kernel as Python calls a function, and the call is inlined with binders of its own, which is loopy's substitution rule with a loop of its own (loopy #1044) expanded where it is used. A rule returns its value:

  @kernel
  def dot(a: Arr[Fin[n], Real], b: Arr[Fin[n], Real]) -> Real:
      return reduce_sum(a[k] * b[k] for k in a.dom)
  
  @kernel
  def matvec(A: Arr[Fin[m], Fin[n], Real], x: Arr[Fin[n], Real],
             y: Arr[Fin[m], Real]):
      for i in y.dom:
          y[i] = dot(A[i], x)
  

  What dot's contract checks at each call natively, the row A[i] in bounds, an index argument in its sort, two rows apart where it writes one, matvec owes at every point it calls dot, as facts isl decides and Lean proves, and a callee's postcondition holds after each call, restated of what it was passed. An argument is evaluated once, when the call is made, as in Python, and a scalar is converted into the sort the callee declares, as its contract converts it natively: same(x[i]) of a rule taking s: np.float32 is np.float32(x[i]) in matvec's term, and a value the contract refuses (a real for Int) is refused where it is passed.

• It plugs into a proof host. loopty registers a theory, an isl oracle, an executor and a run verb with lanky, so a residual obligation an oracle cannot decide is an ordinary theorem a person can prove, in the same ledger.

Status

This is 0.4.0. The core works: kernels and programs run natively, trace, state their facts, lower and run compiled on the C target, and every schedule step is checked; the edges are sharp. New in 0.4.0, a schedule names what it moves by where the file defines it, l2p.unpack.box, and takes no string any more; a kernel body branches with Python's own if and else, and a loop nest is one line, for box, a, b in R.dom:; a statement is named by its text and line in the ledger and in every refusal; and the code a schedule generates computes a sum's operands that read none of its binders once, before the sum. The CHANGELOG has the migration from 0.3, step by step. The five demos at the top of this page came with 0.3.0, as did Lean's proof of each fact isl decides. Below is what works, what works in part and what does not work yet. After it, Known limits lists the open issues you are most likely to meet, each with what to do about it. None of them makes the ledger call a false fact decided, but in several the compiled run computes something other than the native run, and the ledger catches that only on the inputs it tries: read that section before you rely on a compiled run. CI runs the suite on every change, on Python 3.12 and 3.13 with lanky's main and on 3.13 with the released lanky, and reruns every command whose output this README, the quickstart and the examples' README show, holding them to what it prints. It runs no OpenCL device. Until 1.0 an interface can change in a minor release; the CHANGELOG says when one does.

Works.

• @kernel and @program: inert, registering, running natively on numpy.

• A program's term, and one kernel for it. Program.term is what the body does run once against placeholders: every kernel call is recorded instead of run, and the callees' terms follow one another in call order, in the program's names. Sizes are unified through the arrays passed (scale's n is scan's n + 1 when both are handed off), each call's loops get names of their own, and its statements' ids are named after it (scan.S1, step@2.S0), which a map per statement names them by; a message names a statement by its text and line. An array the body makes with Arr.zeros_like(u) is a temporary of the term, zeroed where it was made, and a loopy temporary of the lowered kernel, so the intermediate between two kernels is nobody's argument. The edges between the calls are not declared: an array one call writes and a later one reads is one array of the term, so the dependence is in the footprints, and it orders the instructions of the one loopy kernel the program lowers to. LoopyExecutor().run(program, ...) runs it, Schedule(program) schedules it, and loopty run compares every program in a file with its native run, as it does a kernel. A program may also make an array of a declared shape, whose sizes it reads off its arguments as it does natively, Arr.zeros(TRIANGLE, n=x.dom.size) over a domain or Arr.zeros((x.dom.size, 3)), or Arr.zeros((x.dom[0].size,)) where x has a row 0 (the compiled run is refused where it has none, as natively); a call may be given one array for two parameters it only reads, outer_product(u, u, z); and a callee may take an array over a Where, Sigma or union domain, whose sizes are unified across the calls and whose points two calls given it have to agree on, compared as isl sets (#74); an array the program makes over a domain is stored boxed, and Schedule.pack packs an argument over one. A body that does anything to an argument but pass it to a kernel, make an array like it, or read its sizes, is refused with a TraceError naming the fix.

• Facts that travel between a program's calls (loopty.hypotheses). A kernel's requirements on its inputs are its argument types, and three of them are about what an array's cells hold: an element of a Fin[m] sort is a point of it, one of the Nat sort is not negative, and the offsets a ragged family is read through are the ones its counts give. Where an earlier call wrote the array, the requirement is a requirement fact of the program, decided by isl under what held at the call: the earlier callees' postconditions that nothing has written over since, the zeros an Arr.zeros_like starts an array with, the types the program's contract checks of what nothing has written yet, and the theorems the program cites (@program(uses=[scan_monotone])), instantiated at the cells the requirement reads. The fact rests on the facts it used, so gather's requirement after number reads decided and is worth tested, as number's postcondition is. Where the hypotheses do not decide it, the compiled program checks the cells between the two calls and stops there, with the message of the native refusal; the requirement stays assumed, and its fact says why. The same hypotheses decide a callee's in-bounds fact its own term leaves assumed, a flat val[off[r] + j] after the scan, as a fact of the program's. examples/travel.py shows all three.

• Fusion, and the storage of what passes between calls. An array a program makes and one call writes and a later call reads is an internal edge, and the program has a definedness fact for it, decided by isl: every cell the reader reads is one a call before it stored, or one of the zeros the array was made with, which the fact's reason then lists (y[i] = f[i] + 1.0 (padding.py:19) reads the zeros at f[0] and f[n - 1]), so zero padding at a boundary is read as written and lanky check passes. Where no read sees the zeros, every schedule of the program runs without them: isl decides that when the schedule is made (every cell of burgers_rhs.f that burgers_rhs reads, burgers_rhs.flux has stored by the time it is read, so burgers_rhs runs without the zeros burgers_rhs.f is made with), the fact is the schedule's first, and a window or a substitution rests on it. Schedule(burgers_rhs).fuse(burgers_rhs.flux.j, burgers_rhs.divergence.i, shift=1) runs the reader's loop inside the writer's, one step behind, and is affine with the map per statement it builds, the divergence's i at the flux's step j = i + 1, checked on the dependences between the calls; without the shift it is refused with the pair of instances, and the message names the fusion that is accepted. fuse takes one loop of each call, as many calls as are fused, each loop at the same depth of its own call, and fuses those loops and every loop outside them, leaving the loops inside one after the other: two passes over elements, the second of which reads every cell of the element's row the first stores, fuse on the element and the row (fuse(weak_derivative.along_r.a, weak_derivative.along_s.a)) and no deeper, and the refusal of the deeper fusion names the shallower one (examples/sum_factorization.py). .substitute(burgers_rhs.f) then stores none of f: the one statement that writes it becomes a substitution rule (loopy's assignment_to_subst) computed again at every read, legal when every read is of a cell it stored before the read, since no zeros are left to read, and nothing writes what it read in between, and every later step is checked against the dependences of the program as it then runs. .contract(burgers_rhs.f) keeps f in the cells live at once instead: fused one step behind, the divergence reads f[j] and f[j - 2] at step j, so a window of three cells is found, each flux computed once:

  flux, divergence = burgers_rhs.flux, burgers_rhs.divergence
  fused = Schedule(burgers_rhs).fuse(flux.j, divergence.i, shift=1)
  windowed = fused.contract(burgers_rhs.f)          # lambda j: j % 3
  fused.contract(burgers_rhs.f, lambda j: j % 2)    # refused: a read sees its slot overwritten
  

  A placement is a Python function from a cell's indices to its slot, given or found axis by axis along the order as it stands (one row of the sum-factorized passes' v, lambda e, a, b: b); an array over a domain is placed by the domain's binders (one box's triangle of sumpy's L2P, lambda box, a, b: (a, b), in examples/l2p.py). The array may be stored by several statements, and read by one of them, a recurrence. The step is legal when every read is of a cell stored before it, the definedness fact, and no read sees its slot overwritten between that store and the read, the no-clobber fact, decided by isl over the schedule's order and proved by Lean, and refuted with the read, the cell and the instance that overwrote its slot. From then on every pair of accesses to one slot is a dependence of the schedule, so a later step that breaks the window is refused with its pair. A window of a number of cells is each work item's own, so the boxes of the L2P on g.0 each keep their triangle, and a read on a work item its cell's store is not on is refused; a window of a size, f_window[n], is one for every work item (in global memory on OpenCL). It lowers to double f_window[3] with modulo addressing, or, with registers=True (the default on OpenCL), to a shift register whose every access has a constant index, f_window[0] = f_window[1] at each step. examples/fusion.py is the Burgers flux and divergence of examples/composition.py, fused, then substituted, or kept in a window, each compiled run compared with the native one.

• Declared layouts, checked (loopty.layout). A layout is a function from a cell's indices to its address, its cells and its buffer in its annotations, and a layout of pieces is layout(first, second):

  def centre(a: Fin[1], b: Fin[1]) -> Fin[2 * P + 1]:
      return 0
  
  def derivatives(
      cell: Where[a: Fin[2], b: Fin[Q], (a + b >= 1) & (a + b < Q)],
  ) -> Fin[2 * P + 1]:
      a, b = cell
      return 2 * (a + b) - 1 + a
  
  storage_index = layout(centre, derivatives)        # sumpy's get_storage_index
  as_dim_tags = dof_array.same_as("N3,N0,N1,N2")     # grudge's output layout
  

  lanky check asks three facts of each, every cell at an address of the buffer, no two at one address, a cell at every address (whose witness for Lean is isl's inverse, with // and %), and same_as one more, and a kernel reads through the layout's own function. A wrong nesting, N3,N2,N1,N0, is refuted with the cell it moves (examples/storage_layouts.py). A layout's body may compute with a size, n**3 * e + (i + n * j + n**2 * k) over Fin[n] for every order: isl cannot state n * j, so it decides each fact at n from 1 to 6, which reads tested, and Lean proves it for every n, given the steps omega lacks (the products that bound each digit, the digits of an address by division); two such layouts the same is a polynomial identity Lean multiplies out (examples/hexahedra.py). A call gives the size by name, dof_array(e, i, j, k, n=4).

• Tracing a body to a typed term: accesses, statements, reductions, ragged fibers, branches, source locations, and a TraceError that names the fix when a while loops on a computed value, when a body has more paths through its branches than loopty.trace.MAX_RUNS, when a Python name, a global, or a list, dict or set carries state from one loop iteration to the next, when an array is used whole (y[:] = ..., x * 2, a numpy function of it), when a reduction's if clause is not a bound isl can state, when the trace changes Python state outside the arrays (a global, a closure variable, an object's attribute or slot, a list, dict, set or numpy array they hold, a ragged array's offsets, or the same in a helper the body calls), and when the body prints, reads input, opens a file or draws a random number. The kernel's own module and package count as its code wherever they are installed, site-packages included.

• Branches: if, elif and else, and, or and not, and conditional expressions a if c else b, on loop variables, sizes, scalars and data. The trace runs the body once per path through them (#266 for the select), each condition answered from the path, and joins the runs into one term: a statement after an if/else is one statement, a value two paths choose between on data, vn[i] * (a[i] if vn[i] > 0 else b[i]), is one statement computing a select (?: in C), and a path isl shows to be empty records nothing. The native run is Python's own: the branch not taken is not run. A path's conditions narrow a statement's isl domain where isl can state them: an affine comparison of loop variables, sizes and scalars of an integral sort. A condition that reads an array, compares with !=, or compares with a Real scalar is evaluated at run time only, the statement lists it in Stmt.unnarrowed, and the facts stated over its domain say so in their provenance. loopy's own bounds check reads a guard that names no array into isl only when it is affine in integers; under any other, a product of loop variables, a float literal or a Real scalar (if i * 0.5 < 2:, if i * 2.0 < n:), it checks an access at every point of the loop, and refuses one it cannot show in bounds there, one the guard keeps in bounds included. The refusal (loopty.refusal.GuardNotRead, a LoopyIndexError) names the statement, the guard, and a guard over integers that keeps the access in bounds, as the body writes it, if i <= y.dom.size - 5: for x[i + 4] (#148; note 23 in docs/loopy-notes.md). A condition has to be a truth value: natively ~(i > 0) is -2 or -1 (~ on a Python bool is bitwise), which Python's if takes as true, while the trace reads not (i > 0), so the trace refuses it, naming i <= 0 and not. What is stored into an array of Bool has to be a truth value too, a comparison, a connective of truth values or a Bool read, since natively the array is a bool and compiled a byte, into which C converts 0.5 as 0: b[i] = u[i] is a TraceError naming b[i] = u[i] != 0, and an integer arriving at a bool array natively, ~(i > 0) again, is refused there. Every operand of &, | and ~ has to be a truth value, wherever the connective is, the index of the cell written included, since the trace and the compiled kernel read them as and, or and not and natively they are bitwise on an integer: (k[i] & 1) * x[i] is a TraceError naming k[i] != 0, and k % 2 for k & 1. And ~(i > 0) used as a number (x[i] * ~(i > 0), in a sum's body or an index) is a TraceError naming i <= 0, since natively it is the -2 or -1 the native run refuses only where it is stored into a Bool array. More operations are refused, naming what to write: b[i] + c[i] of two truth values, which is or natively and 2 compiled (b[i] | c[i] for or, 1 * b[i] + c[i] for a count), and b[i] - c[i], which numpy refuses; an integer element or scalar to a negative integer power, which numpy refuses and loopy computed as an integer (1 / k[i] ** 2, a real); ^, << or >> of a real, and // or % of a complex value, which numpy refuses and C cannot compile; a loop variable divided by a literal zero or shifted by a negative literal, which Python refuses natively, or shifted left by a literal of 64 or more, which Python computes exactly and the compiled run to 0 (i * 2.0 ** 64); -b[i] of a truth value, which numpy refuses (~b[i] for not, -(1 * b[i]) for the integer); an integer literal beside a numpy integer whose type does not hold it, u[i] // -1 of a uint64 or (b[i] // c[i]) * 200 of two truth values (an int8 natively), which numpy refuses (np.int64(-1)); ^, << or >> of a uint64 and a signed integer, which numpy refuses; and an integer literal past 64 bits where numpy does not compute it in a real or a uint64 that holds it, 1.0 * (i + 2**64), which the compiled kernel has no integer type for (2.0 ** 64).

• The faithfulness fact. For each kernel and each program, the traced term is run by an interpreter (loopty.interpret: statement by statement in source order over each statement's isl domain, each loop enumerated when the run reaches it, so that a ragged loop runs to the length its row has when the loop starts, expressions evaluated with numpy's arithmetic, reductions summed in the order reduce_sum sums natively) and compared with the native run, on the file's example_inputs() and on three inputs drawn from the declared types. It is a trace-faithful fact, tested on agreement and refuted with the input and the first differing cell, which is where state hidden past every check above shows up, and, for a program, a call its body makes or skips by looking at an argument as no placeholder can be looked at (isinstance(x, Arr)).

• Typing rules and the ledger: in-bounds by isl or by type, write disjointness, ordering, reduction exactness, postconditions. A write into an array whose element sort is Fin[m] owes the fact that the value written is a point of Fin[m] (element-sort), decided by isl for a quasi-affine value, by type for a value read from an array of that sort, and assumed otherwise, and an index read from such an array is in bounds by type resting on those facts, so perm[i] = i + 1 followed by x[perm[j]] is refuted where it was decided. The reads a ragged layout makes are accesses like any other, in-bounds obligations and dependences every cast is checked against: the start of the row a ragged access is flattened through (off[r], when the kernel declares the offsets), and the length of the row a loop over a ragged fiber runs to (cnt[r], which the lowered code reads once per row). Each is listed where the lowered code reads it and nowhere else, so a kernel that writes its counts or its offsets is ordered against exactly the rows that use them.

• One meaning for a kernel that writes its own layout. The native run and the term interpreter read a ragged array through the counts and offsets the kernel declares, as the kernel has left them, which is what the lowered code does with the arguments it is handed; the contract checks that they are the array's own layout when the run starts.

• IslOracle: Empty, Subset, Bijective, Monotone, each refutation with a witness.

• isl decides, Lean checks (loopty.description, loopty.obligations; docs/lean-checks.md). Every set loopty builds keeps the text it was written from, on the set, with each call's renaming, and every schedule step a record of its maps as written; from those, and from no set isl computed, every fact isl decides is stated for lanky's Lean oracle in parts (one per place an access is made over, one per dependence class of a cast), which proves them with omega where Lean is installed. isl's work comes in as hints Lean checks: its functional form of a step's map as the witness that the step is total, its lexmax pieces as the writer of each cell a read reads, the instances it took of a program's universal hypotheses, and its counterexample to a refuted fact, which Lean proves false. loopty run asks Lean about the casts and layout facts it reports, and --show-lean prints the Lean. On the examples and the test kernels Lean proves every fact isl decides but those listed under Known limits; scripts/lean_coverage.py reports it.

• Definedness inside one kernel. An array a kernel reads only at cells an earlier instance of its own wrote (scan reads off[r] where off[0] = 0 or the row before stored it) gets a definedness fact where isl shows it, stated for Lean as a ∀∃ with isl's last writer as the witness; a kernel that reads what an array held when it was called owes nothing and gets none.

• Postconditions decided from the term (loopty.dataflow). Where no cell of an array a kernel writes is stored twice, or read before its store, a single-assignment fact says so, decided by isl and proved by Lean, and what each statement stores is what the kernel leaves: its value equation, off[r + 1] == off[r] + cnt[r] at every r of its statement's loop. The postcondition is decided under those equations, each taken at the cells the claim reads, and rests on the single-assignment facts. A claim about reals is decided the same way, each equation isl cannot read an atom, so R[a + 2, b] == -R[a, b + 2] is what the recurrence stores at (a + 2, b); Lean states it with Mathlib's reals and proves it by linarith. The native runs still test it, and a run that refutes it refutes it.

• A kernel cites theorems, @kernel(uses=[on_the_triangle]). A theorem whose hypotheses are conjuncts of the postcondition, up to the names of their binders (a family R(a)(b) stands for the array R[a, b]), gives its conclusion at the kernel's arrays, a by-theorem fact resting on the theorem and the postcondition: what is left unbound stays universal, and the hypotheses the postcondition does not give are its hypotheses. Lean checks the application, the theorem as it was proved applied with the arrays loopty matched: examples/taylor_reconstruction.py,

  @kernel(uses=[on_the_triangle])
  def reconstruct(
      c: Arr[Where[a : Fin[2], b : Fin[p + 1], a + b <= p], Real],
      R: Arr[Where[a : Fin[p + 1], b : Fin[p + 1], a + b <= p], Real],
  ) -> all(
      R[a, b] == c[a, b] for a in Fin[p + 1] for b in Fin[p + 1] if a < 2 if a + b <= p
  ) & all(
      R[a + 2, b] == -R[a, b + 2] for a in Fin[p + 1] for b in Fin[p + 1] if a + b + 2 <= p
  ):
      ...
  

  reads by on_the_triangle, for every D with pde and stored: ... R[a, b] == D(a)(b): whatever table D satisfies the PDE (pde) and has its stored entries in c (stored), R is D on the triangle. A theorem whose conclusion is a conjunct of the postcondition decides that conjunct instead, as a program's cited theorem decides a requirement. The application says something of what the kernel computes, or it is not made: every hypothesis about an array the kernel writes has to be a conjunct of the postcondition, so a kernel whose recurrence is off by one gets nothing from the theorem; and a conjunct about reals is a hypothesis only where nothing in it rounds (copies and sign changes), since the theorem is about the reals, and s[i] == a[i] + b[i] holds of the floats and not of them read as reals. The theorem is taken at the kernel's sizes in every sort as well: a table D: Fn[Fin[n], Real] that nothing binds, with n matched to the kernel's p, ranges over Fin[p].

• A cited theorem about a series, applied to the terms a kernel adds (loopty.summands). A theorem that concludes all(c(n) == E for n in Nat if ...) is about the n-th term of a series; where the kernel's term adds E at whole numbers n, with the theorem's variables bound to the kernel's expressions, a by-theorem row says so, resting on the theorem. An identity the kernel cites, a theorem that concludes L == R, rewrites R in E into L first, the kernel's way of computing it. examples/layer_potential.py's single_layer cites line_taylor and by_the_conjugate, ratio(alpha, beta) == -beta / alpha, and reads single_layer adds line_taylor's c(n) for n in 1..5, for every c with away, first and ode, read over the complex numbers: about the term it evaluates, read over the complex numbers, and not about the floats it produces. With Mathlib, Lean checks it a term at a time, the identity rewritten and the theorem applied.

• Equations of reals at their class (loopty.reals). An equation of a postcondition that rounds, a closed form the kernel computes another way, is compared after each native run as the compiled run is with the native one, at the class of what it reads: |a - b| <= 1e-6 (|b| + 1) for approx. One that rounds nothing (a copy, a sign change), one of integers and one whose class is exact are compared bit for bit, and the statement says which: ..., to the approx tolerance. Where the cells it reads are stored once, by statements whose sums run over concrete bounds, Lean proves it with Mathlib over the reals, the tables at the values they hold; that is a proof of the term read over the reals, and the fact rests on the kernel's trace-faithful fact. examples/assembly.py's mass_ffc promises the element's mass, |J| (1 + delta_jk) / 24, from FFC's 15-digit tables.

• lanky check lists a file's schedules (#191). The casts of every schedule a file holds at module level are rows of the ledger beside the facts of the kernel or program it schedules, decided by isl as each step was taken and proved by Lean where it is installed, each id once.

• Mixed subscripts. X[cell[c, k], d] is in bounds by type on its first axis, cell[c, k] : Fin(nv), and by isl on the second, where it was assumed.

• Schedule.assume(k.n % 4 == 0), loopy's assume, checked: loopy is told, every later cast is decided and stated where the condition holds and rests on an assumption fact, and a call whose sizes break it is refused.

• Schedule: tag, split, interchange, prioritize, tile, skew, affine, fuse, each checked as a cast, each emitting its fact; realize and atomic, which move no instance and record what they make of a sum; and substitute and contract, checked as storage decisions; retarget, which replays every step against another loopy target and re-checks it. prioritize and interchange also take the loop of a sum, named before loops of its statement, which it then runs outside, as sumpy's CPU tile of a layer potential runs 256 sources at a time outside 16 targets (.split(single_layer.s, 256).prioritize(single_layer.s.outer, single_layer.t.inner)): each sum still adds its terms in order, the loops it runs outside of order none of their instances, and loopy keeps the order where it keeps an accumulator per instance, on a hardware axis or an ilp loop, which the buildable fact asks (examples/layer_potential.py runs sumpy's tile on C with an ilp loop for its work items). A split sum adds its terms in order, outer loop first. A tag belongs to a loop, so split, tile and affine refuse a loop that carries one, and the loops they make are tagged after them. affine(map) takes an isl map from loops to the loops that replace them, refuses one that misses or merges an instance or runs a dependence backwards, and rewrites the kernel over the map's image; skew is that method with a particular map. A union map whose tuples name statements, { S0[t, i] -> [a, b] : ...; S1[t, i] -> [a, b] : ... }, moves each statement by its own map, checked on the dependences between the statements as well as within each, which is the time offset a diamond tiling of two statements that feed each other needs. A program's statement, flux.S0 or step@2.S0, which isl cannot read as a tuple name, is named with every other character spelled _, flux_S0 or step_2_S0, the id of its instruction in the lowered kernel. A loop over an image with holes, such as the diamond's b (acoustic.loops.b), counts its steps (b = 2*b_step - a, Schedule.strides) instead of testing a parity at every b. A loop on a hardware axis (g.*, l.*) is the launch grid, outside every other loop, and nothing in a kernel orders two of its work items through global memory, so the monotone cast also refuses any dependence between instances on two work items, with the dependence and the two work items as the witness (the stencil's jacobi with i on g.0, the acoustic pair likewise). A statement with no loop on the axis runs on every work item of it. A schedule runs without the zeros a program makes an array with where no read sees them, a definedness fact isl decides when the schedule is made, so examples/l2p.py's fused.on_groups(l2p.unpack.box), sumpy's schedule for a GPU, names only the file's loops and zeroes nothing. Zeros a read sees, padding at a boundary, run where the array is stored: their loop along an axis of the array goes on the axis of the loop that every statement storing the array stores it along. A sum on a local axis stays allowed, since loopy synchronizes its partial sums; its body reads on every work item, and its statement stores the result from one, so a sum whose body reads what another sum's statement stored is refused, as loopy refuses it for want of a global barrier.

• The target-capability check: a concurrent tag (a hardware axis, ilp or vec) inside a data-dependent (ragged) loop bound or its domain, a hardware axis on a reduction nested in another, a reduction loopy will not realize (partly in parallel and partly in sequence, across two local axes, on a group axis, or on a local axis whose extent has no numeric maximum), a hardware axis loopy will not assign (numbered past an unused one, shared by two loops of one statement, missing from an instruction the kernel runs beside it, or l.auto), an unr, ilp or vec loop whose length is not a number, a temporary loopy misreads once an ilp or vec loop has a copy of it per iteration (a ragged row's length), a loop ordered outside a loop loopy nests it inside, or a hardware axis on the C target, which has none, is reported as a refuted buildable fact and raises UnbuildableSchedule when something asks for code. It is asked of the schedule as it stands after every step, so an interchange can make a tiled ragged loop buildable again.

• Lowering to loopy, including a ragged axis as a flat buffer plus offsets, and running on lp.ExecutableCTarget. Every argument of LoopyExecutor.run is an argument of the kernel; the target is chosen by the schedule (Schedule(kernel, target="opencl")) or by the executor (LoopyExecutor(target="opencl")). Each operation is computed in the type numpy computes it in natively (loopty.promotion, by NEP 50), where C's conversions would pick another: a quotient of integers is a double, a Python float beside a float32 is single precision (0.1f), a float32 beside an integer array is a double, and a floating power calls pow, as numpy does. Any power compiles on the C target, of a complex base too. Integer arithmetic is 64 bits wide in both runs: Nat and Int are stored in 64 bits, and a product of a Fin[m] array's entries or of loop variables is computed in 64, a subscript's too unless loopy reads the subscript as affine with no division (note 23). //, %, << and >> are computed as numpy computes them, by functions loopty's targets define (loopty.operations): of reals too, by zero (0 for integers) and past the width of a shift. An operation numpy computes in a numpy integer type is computed in it compiled too: one narrower than int (np.int8, np.uint16), which C computes as an int, is converted back into it, and an unsigned integer beside a signed one is computed in the signed type numpy computes in (int64 for a uint32 and an int32), where C's would be unsigned; integers compare exactly, as numpy's do, a negative one with an unsigned one included. A sum is computed as the body wrote it: c[i] - a[i] as the difference, and c[i] + -a[i] as the sum of a[i] negated in its own type, which differ for an int8 a at -128 and an unsigned one at every value but zero; u[i] + -1 of an unsigned u and k[i] - 2**31 of an int32 k, which numpy refuses at every point, are refused by the trace (#151). A real less an integer keeps the sign of zero, x[i] - 0 and x[i] - i at x[i] = -0.0, and a truth value subtracted, k[i] - True, which pymbolic takes no operand of, is refused naming the integer. Loop arithmetic beside a numpy integer whose type does not hold it at some point, u[i] + (i - 3) of a uint32 u at i = 0, is refused by numpy natively there, and the compiled run is refused before it starts where the call's sizes reach such a point, naming it (loopty.refusal.RefusedNatively, #153). abs of an integer is numpy's, the smallest value of its type included, and abs of a real in a guard is written without a call, since loopy resolves no function only a guard calls (#154); another function only a guard calls is refused, naming the fix. abs of a term compares on either side, if abs(k[i]) > 2: and 127 == abs(c[i]), since lanky 0.2.0 (xywei/lanky#94). An integer literal past 64 bits is written as a real where numpy computes in one (x[i] * 2**70), or as a uint64 beside one (u[i] + 2**63), and refused elsewhere. A kernel named like a function, a macro or a type loopy, the C headers or OpenCL C know (floor, pow, cpow, get_global_id, clamp, NULL, size_t) is renamed in the generated code, and a parameter, size or loop variable named like a macro the headers define (I with complex values, NAN, INT32_MAX, M_PI, NULL) or a function the kernel calls (pow with a power) is refused, naming what the code means by it. So is a size, an integral scalar or a reduction binder named like a keyword of isl's (max, floor, mod), and a loop variable of such a name is renamed. See notes 19, 20 and 23 in docs/loopy-notes.md.

• Array arguments over polyhedral domains (loopty.domain): Where[...], binders written as slices and then the comparisons that cut their box, joined by &; Sigma[...], binders and an unnamed last fiber affine in them; and a union of pieces, Fin[n] + Fin[m] (lanky's SumType). .dom runs binder by binder, L.dom[i] over the points the constraints allow at i (a fiber at a point outside the domain is empty), and a union's pieces by number, which a trace runs as the Python loop it is. The in-bounds obligations are decided over the exact set. An argument over other points is refused at every entry point, and so is a plain ndarray. Both layouts lower and run on the C target: the box of the binders, addressed by loopy, and packed rows through a table of row starts (Schedule.pack), which the executor computes from the domain and passes. Arr.zeros(domain, n=..., storage=...) and Arr.from_cells build such an array, and Arr.cells() reads it in one order whatever its layout. examples/pairs.py is the demo.

• Elementary functions and complex parts, loopty.math: sqrt, exp, log, sin, cos, real, imag and conj, numpy's natively and the target library's compiled (csqrt, clog, creal on the C target, pyopencl's cdouble_log on OpenCL), each computed in the type numpy computes it in: sqrt(k[i]) of an integer is a double, and real, imag and conj of a real are the number, 0 and the number. One numpy computes in half precision, sqrt of an int8, is a TraceError naming sqrt(1.0 * k[i]). A value of one compares as any value does, if sqrt(x[i]) > 1.5:, and numpy's own np.sqrt, Python's math.sqrt and z[i].real are refused when traced, naming the function of loopty.math to write. examples/layer_potential.py is pytential's QBX single and double layer over complex points. ipow(x, n) is x to an integer power multiplied out as loopy multiplies one, by squaring, in both runs, so they agree bit for bit and the compiled kernel calls loopy's own loopy_pow_float64_int32; x ** n of a real x is numpy's power, the C library's pow, in both. In examples/l2p.py's Taylor sum that call was the whole cost on the C target.

• Agreement with the code a kernel mirrors, @upstream(kernel) (loopty.upstream): a function that computes, on the file's example inputs, what the upstream code computes there, and a ledger row under the kernel that compares it with the native run at the accuracy the kernel's types claim, tested by upstream, refuted with the cell, or assumed where the function cannot run (its package or its recording is not there). examples/l2p.py compares its program with sumpy's own E2PFromSingleBox, recorded on its inputs. @upstream(kernel, within=1e-12) compares each output computed in floating point within 1e-12 of the largest value the function gave for it, and the row says so (... agrees with sumpy_single_layer in out within 1e-12 of the largest value, on its example inputs): round-off between two orders of summation is far inside the class's allowance, and the row states how far.

• Constant tables, loopty.table: W = table([1 / 6, 1 / 6, 1 / 6]) next to the kernel, read by the body as an array is (W[q], W.dom). A table is a read-only copy of its values; its reads are in bounds over its shape, decided by isl, and the lowering writes it into the kernel as a read-only loopy temporary with the values as its initializer, static double W[3] = {...} on the C target and __constant on OpenCL, with every digit. A program holds one table for every call that reads it. Natively a table is read as an array is, and it is no argument: an array of its values is passed instead, Arr.from_numpy(W.values.copy()). examples/quadrature.py is firedrake's P1 mass kernel.

• Claims settled by running them, loopty.tested: a function that takes nothing and returns whether its claim held, alone or with what it measured, decorated @tested, is a row of the ledger, tested by run with its docstring's first line for statement, or refuted with the reason. It is for what no oracle reads, an accuracy that converges at its order say, and it never reads more than tested; the agreement with the code a kernel mirrors is @upstream's row. examples/layer_potential.py measures pytential's Green identity and the jump relations this way, and compares both its kernels with sumpy's own through @upstream.

• Updates and atomics. A statement that adds to the cell it writes a value that reads nothing else of the array, y[cell[c, k]] += yl[c, k] or y[r] += v[r, j], is an update. Its sum has the class of what it adds, as a reduction's has (approx for Real). An update whose instances can add to one cell owes that class, the accumulation into y[cell[c, k]] over c, k is approx, and no distinct writes; one that updates a cell of its own at every instance, y[i] += a * x[i], owes its distinct writes. Two updates of one cell are no dependence, so a schedule may reorder them, which it records as a reassociation. A step that runs two of them at the same time, on two work items or in no order (ilp, vec), is refused with the pair and the fix until Schedule.atomic(mass_action.y). Then each update is loopy's AtomicUpdate on OpenCL, a compare-and-swap, and the C target, which runs in one thread, keeps its code. A refused cast through a cell read from an array names two instances that touch one cell on the example inputs, with what the subscripts read, where cell[0, 2] = 3 and cell[1, 1] = 3. A write of that kind, y[perm[i]] = x[i], has its distinct writes decided on the native runs. examples/assembly.py is PyOP2's P1 mass action, assembled.

• The argument contract, enforced at every entry point that runs a kernel (compiled, differential and native): two distinct array parameters may not share storage where the kernel writes either, a ragged argument has to agree with the counts array its type names and with any offsets passed alongside it, an argument over a domain has to have the declared domain's points at the sizes of the call, and a value of a refined sort such as Fin[m] has to be one — an array element and a scalar argument alike, and being one means being a finite whole number in range, not merely passing two comparisons. A value of Nat or Int is also inside the 64 bits both runs store it in, and one of Fin[m] inside the 32 the compiled run stores it in, since 2**32 + 5 would run natively as it was and be 5 compiled. An array the kernel writes has to be stored as its element sort is natively (float64 for Real, bool for Bool, int64 for Nat and Int, a signed integer of 32 bits or more for Fin[m], a numpy sort as itself), since an integer x for a Real parameter truncates every write the compiled run keeps; an array it only reads is read by the native run in that dtype, as the compiled run converts it, so an integer x no longer overflows natively where the compiled double does not. A complex entry of a sort that is not complex has no imaginary part, and an entry of Bool stored as a number is 0 or 1. A scalar is asked the same, and is converted into the dtype of its sort in both runs, since it is passed by value: np.int64(2**32) for a Real is a double natively too, a Bool a numpy bool, on which ~ is not, and a Python number a numpy scalar, so a=0.7 and a=np.float64(0.7) compute alike beside a float32. These are the assumptions the typing rules and the two runs make about a call rather than about the term, and a violation is a ValueError naming the argument. Distinct parameters being disjoint storage is the load-bearing one: dependences are computed per array name, so a kernel reading x[i - 1] and writing y[i] may legally run i in parallel, and the same kernel called with x is y is a race that the differential test cannot see, because it copies each argument separately.

• loopty run FILE [--target c|opencl] [--emit-code] [--json OUT], loopty check FILE, and lanky run FILE through the entry point. --target retargets every schedule in the file, re-checking its casts, and says so by name when one cannot be retargeted; without it each schedule keeps the target it was written for. A refuted fact is repeated under the ledger with what explains it, the block lanky check prints (a compiled run that disagrees names the outputs and by how much), and the command exits 1. So it does when two kernels claim one fact id, as two kernels one definition makes (a factory) do: a DUPLICATE block names the kernel and its ids, as lanky check names them, and the --json ledger keeps each claim not in the table under refused_claims, with its status and what explains it. A schedule of a schedule, Schedule(Schedule(k).split(k.i, 2)), is that schedule built again from k: its steps replayed and checked again, on its target and with its sizes unless others are given. Its run is compared with k's body, and k is not run again through the identity.

Partial.

• Ragged bounds are reflected into isl as one parameter per distinct bound term (allocated once per kernel, so the same cnt[r] is one parameter everywhere and two different bounds are never given one name), so cnt[r] and cnt[r + 1] are unrelated to isl. Within one kernel nothing knows that counts sum to the offsets, or that off is monotone, so an access against flat storage, val[off[r] + j], is reported assumed in the kernel's own ledger, with the reason in its provenance. The ragged spelling val[r, j] over 0 <= j < cnt[r], which is what the tracer and the demos produce, is decided. In a program, the same access is decided where the kernel is called after the call that wrote off, under that call's postcondition: the cells off[r], off[r + 1] and cnt[r] become isl parameters, and the scan's recurrence, instantiated at them, is an affine constraint between them (loopty.hypotheses, examples/travel.py). The fact is the program's, resting on what it used.
• A kernel's postcondition is decided from what its statements store where every cell of the arrays it claims about is stored once at most, and read only after its store (loopty.dataflow, below); elsewhere it is tested against its native runs, on the file's example inputs and on drawn ones, and its provenance says why it is not decided (undecided). A cell stored twice, a running value a loop overwrites, needs the last store before each read, which is not done; nor is an update, y[i] += v or a scatter-add y[cell[c, k]] += v, whose cell holds what it held on entry plus every update. @program restates a callee's postcondition as a fact in scope, decided by the call and resting on the callee's own fact (lanky's rests_on), and offers it as a hypothesis to the calls after it, until a call writes an array it mentions. A callee imported from another file has its fact in that file's ledger, so lanky check counts the id as an assumption and names it in an UNRESOLVED line; the id is the one that ledger holds, since every fact of a kernel is keyed by the kernel's definition (lanky.ledger.fact_id over the module the file's path gives it), and a kernel of the program's own file with the callee's name has an id of its own. A theorem the program cites (@program(uses=[scan_monotone])) is used at the arrays its hypotheses match, and only where each array's cells are points of the family's sort when the call is made: scan_monotone is about offsets of Nat, and nothing says the cells a scan wrote are, so it takes a theorem over Int to say the offsets a scan computed are monotone. Only affine facts reach isl; a hypothesis isl cannot state is dropped, which assumes less, never more. A program's own return annotation is its postcondition, a claim about what it leaves in its parameters. It is tested on the program's native runs, an equation of an output compared at that output's class (approx for reals), and decided by isl from what each call stores where that decides it; otherwise Lean is asked to prove it from the same equations and the theorems the program cites, unfolding a call at a time (loopty.unfolding). examples/hexahedra.py's weak_grad claims the dense Kronecker operators grudge's escape hatch applies, and Lean proves it with Mathlib from what the three passes store and sum_factorization.
• A program lowers sequentially: its calls' loops run one after another, as the program runs them, until a schedule fuses them. A fusion moves the outer loops of each side, as many as both have or as loops says, and a nest the lowering wrote as one domain is cut after them; loops a map takes that are not the outer loops of their domain leave the casts decided and the kernel unbuildable. Only a pointwise producer is substituted (one cell per instance, at its loop variables, no sum, no guard isl cannot state). A window is taken of an array its statements store a value in, not update, which they read in values, before any of their loops is tagged (#213). Every read gets the value as the store converted it, a read in a subscript included (y[(int16_t) (t[i])]), and a conversion inside the value stays ((t[j] + 1) % n of a 32-bit Fin[m] entry, computed in 64 bits); a store's widening keeps the value, and a read in a subscript reads the value inside it (x[t[i]], and x[n - 1 - i], which loopy then reads as affine and checks). The compiled program is one call, so the contract checks its arguments when it starts and not at every call. What a callee's contract checks of the cells of an array an earlier call wrote, or the program made (an element sort Fin[m] or Nat, the offsets a ragged family is read through), is a requirement of the program: decided by isl under the hypotheses that held at the call, or, where they do not decide it, checked by the compiled program between the two calls, which then stops with the requirement's message where the native callee is refused (examples/travel.py). An integral array such a check reads is stored in 64 bits, as the native run holds it, a Fin[m] one included, so that a value written outside 32 bits is read as written and not as the store narrowed it; the compiled run refuses one passed in fewer bits, which the native run would narrow. Offsets a scan computes and a later call declares Nat without reading rows through them are checked: that they are naturals follows from the scan's recurrence only by induction. The counts of a ragged family an earlier call wrote are still refused: its rows are laid out in a buffer the program is given, which no hypothesis about the program's arrays can speak of. The kernels an array is passed to have to declare the same element sort for it, and every call has to read a ragged family's rows through the same offsets; a loop in the body whose trip count is an argument (a host loop) is refused, and so are an array made like a ragged one, an array made over a union of pieces, which is stored as one box (#208), and an array over a domain of a Fin[m] or Nat element sort that an earlier call wrote and a later one is given, whose cells no checked point reads yet (#209). A temporary made like a parameter, Arr.zeros_like(u), has u's dtype natively, so the compiled program refuses a u whose dtype does not hold what the compiled temporary holds: float64 for Real, the dtype itself for a numpy one such as np.complex128, bool for Bool, int64 for Nat and Int, and a signed integer of 32 bits or more for Fin[m] (int64 for an array a checked point reads). On the C target the executor allocates a program's array on the heap and passes it after the program's own arguments, double *__restrict__ f, where loopy would declare it on the stack of the call, whose size bounded it (#62; note 16 in docs/loopy-notes.md); on OpenCL it is a global temporary, which the PyOpenCL host code allocates.
• Only a two-axis (row, fiber) ragged array, its rows counted by an array of one axis, is built and lowered; tracing refuses any other ragged type, a fiber after two dense axes, a dense axis after the fiber, or counts of two axes, before any fact is stated about it.
• Integers are 64 bits wide in both runs, and a result outside 64 bits wraps round compiled where numpy wraps or refuses, which is numpy's limit too. Arithmetic of loop variables, sizes and literals alone is Python's natively, which is exact: 1.0 * i ** 5 past i = 6208 is a wider number natively and wraps compiled, which the differential fact reports (#139). A loop bound, and a subscript loopy reads as affine with no division, are loopy's index arithmetic, 32 bits wide; under -fwrapv such a subscript's sums and products give the cell it names. Any other subscript is computed in 64 bits, as the rest is, its sums of loop variables and sizes included: x[(i * i) % n] (#129), x[(i * 499999) // 1000000], which loopy wrote with a division of the wrapped 499999 * i, and x[(i + n) // 2], whose sum wrapped round past 2**30 (#149). A sum of loop variables, sizes and small literals in a value or a guard on the loops is still 32 bits wide, y[i] = i + n past 2**31 (#162). A subscript over a narrow numpy integer type is computed in it, as numpy computes it: x[s + i] of a uint8 s wraps round at 256 in both runs (#157). One of a signed type that wraps round below zero is refused natively, as any negative index is, and read before the array compiled; its in-bounds fact is not decided, since isl reads a narrow scalar or element as any integer, and is refuted or assumed. A Fin[m] array the kernel writes may be an int32 one natively, whose entries the native run reads as int64, as the compiled run computes with them, unless a program's checked point reads it, which takes an int64 one. An integer to a negative power whose exponent is not a literal is left to numpy's refusal natively, and computed as an integer compiled. Loop arithmetic beside a numpy integer that does not hold it is refused compiled where isl reads it and the call's sizes fix its loops (#153); a product of loop variables, a ragged row's loop, a store under a guard isl does not read, and a sum's body under both a guard and the generator's if clause where the guard is false, are left to the native refusal and computed compiled (#166). So is a value stored into a signed integer array that does not hold it, a uint64 past 2**63 into Int, an Int into an np.int8 array, or a real past 2**63, or a NaN, into Int, whose value only the compiled kernel knows (#159). The written -1 * u[i] of an unsigned u, which numpy refuses, is read as -u[i], which it computes (#168). A numpy scalar on the left of a term is traced as the Python number it equals: np.float64(0.1) * x[i] as 0.1 * x[i], single precision beside a float32, and x[np.uint8(200) + i] as x[200 + i], not taken round (#167). Two statements of one kernel whose lowered arithmetic differs only by the dtype of an equal literal are generated with the first one's literal, since loopy maps a subexpression to the result of an equal one, and pymbolic finds a numpy scalar equal to any number of its value (#170). See notes 19, 20, 23 and 25 in docs/loopy-notes.md.
• A polyhedral domain is an array's whole index set, so it cannot sit beside a dense axis (Arr[Fin[k], Where[...], Real] is refused; write the axis as a binder of the domain). The pieces of a union have the same number of axes, and a piece is chosen by a Python integer, never by a loop variable. A constraint is a conjunction of comparisons: != and | are refused rather than widened (write a union instead), and the packed layout refuses a domain whose rows skip columns (a remainder in a constraint). An array over a domain is indexed, and its fibers taken, at quasi-affine expressions of loop variables and sizes; an indirect index such as L[p[k], j] is refused when traced. A size a binder's bound runs up to is never negative, so a scalar that makes one negative is refused. Natively, an array over a domain enumerates its points with isl when it is built and checks each access against them in Python. A run, native or compiled, is over the declared domain, its loops, its sizes and its layout, so an argument over the same points written otherwise, or stored the other way, is copied into it and back.
• A reduction nested in another one cannot take its bound from the outer binder when that bound is not affine: reduce_sum(reduce_sum(val[q, j] for j in val.dom[q]) for q in val.dom) is decided by the analysis and refused by the lowering with a LoweringError, because a row length is computed inside the loop over its row and a reduction binder has no such loop. Writing the outer reduction as a loop that accumulates into the output, or keeping each row's sum in a cell indexed by the row, lowers and runs. An affine inner bound (Fin[i + 1]) lowers as it is.
• An index expression isl cannot express widens the footprint to the whole array. That is sound, but it can reject a legal schedule.
• A constant table's values are read in its dtype and have no sort beyond it: an index read from a table of integers, x[perm[i]] with perm = table([2, 0, 1]), is not in bounds by type, as one read from an array of Fin[m] is, and its in-bounds fact is assumed (#182). A table is a kernel's or a program's own; a table the body reads through a container, which no frame names, is named table. A table holds finite numbers (#193), and a loop bounded by a value read from one is refused when traced (#194).
• A kernel called inside another's loops (loopty.calls) is given whole arrays, rows of dense, ragged and table arrays, the cells an index row picks, scalars and functions; a row of an array over a domain is not passed (a whole one is), and the cells a gather picks are a copy natively, so a called kernel does not write them (#205). Two arguments over one array the callee writes have to be two rows of it, decided at every call; a whole array and a part of it are refused, and so is an argument that reads an array the callee writes, since natively it is read once, when the call is made, a value a function given captured included. A rule's claim about its value (-> Real & (...)) is refused: state it as a theorem about the rule (#206). A call made inside a sum over Fin[...] of a size, rather than over an array's .dom, has points the trace does not see, and its facts are assumed (#207). A rule that takes a function is compiled where it is called, not on its own. A real given for a numpy integer sort (np.int8) is truncated in both runs, and one outside the type converts as C leaves undefined and numpy warns; a uint64 value given for Fin[m], Nat or Int is refused, since the contract refuses one past int64 (#232). A call that reads one cell of a row, first(A[i]), leaves A's second size in no loop, which fails to compile (#220). A domain declared outside an annotation, for Arr.zeros(TRIANGLE, ...) in a program, names its binders and sizes as lanky variables, i, j, n = Var("i"), Var("j"), Var("n") (#210). A program's array over a domain kept in a window is kept in the box its placement's slots span, 11 x 11 for a triangle of 66 cells, as plain loopy keeps one in a temporary.
• loopty.math has eight functions (#184), and numpy's own functions on a traced value are refused rather than read as them (#195). A function only a guard calls is refused when lowered: loopy resolves a function by the values that call it (note 20 in docs/loopy-notes.md).
• The in-bounds fact of a ragged access val[r, j] is decided against the row's length, and the dependences and disjoint writes of val are decided over [r, j]. Both assume the offsets lay the rows out inside the flat buffer and apart from each other, which the contract checks when a run starts. A kernel that writes its counts or its offsets can break that during the run, so such a kernel has a layout fact for each counts family it rewrites, and the in-bounds and disjoint-writes facts of the family's ragged arrays rest on it, as do a fact decided by type through an index read from one of them (x[col[r, j]]) and the monotone casts of its schedules: the ledger shows them decided under layout:.... isl decides the layout fact when the kernel writes only the offsets, each either as the counts lay it out (off[r + 1] = off[r] + cnt[r], which writes back what the contract checked) or as a value of the loop variables and the sizes, and refutes it with the instance that writes a start the counts can contradict (off[r] = 0 at r = 1). Any other write, a count, a start read from another array, one under a guard isl cannot state, is decided by induction over the run when isl shows it keeps the rows in order (each starting no earlier than the row before it ends, inside the buffer) wherever they were before it: cnt[r] = 0 does, and so does a start moved within the room its row leaves. One isl cannot show does leaves the fact assumed, and the facts on it worth an assumption, unless a native run of the kernel refutes it: one that reads a row off the buffer through the layout it has written (off[r + 1] = off[r] + cnt[r] + 1, which moves the last row past the end), or that leaves two rows on one cell (off[r] = s[r]). The order is more than the rows' staying apart, so a kernel that permutes its rows inside the buffer is left assumed; and it is less than some facts on the layout need, so the induction is not asked of a kernel that moves rows and writes the family's arrays, whose disjoint writes and dependences tell the cells apart as [r, j], one cell each only while no row moves, nor of one that writes a row's length inside a loop over the row, which runs to the length it read when it started. A read of the row inside such a loop is checked natively against the length the write left, which isl cannot read, so its in-bounds fact is assumed itself, with the reason (cnt[r] = 0 inside for j in val.dom[r] makes the native run refuse val[r, 1]); a sum or an inner loop over the row reads the length again where it starts, and its reads stay decided. A program's layout facts are not tried on runs.
• Schedule.affine and maps whose image has holes. The diamond (t, i) -> (t + i, t - i) reaches only the points of equal parity, and loopy's own map_domain refuses it, so loopty rewrites the kernel over the image itself, and the loop left with the holes counts its steps. The code is correct, bit for bit against the native run for the stencil and the acoustic pair, untiled and tiled in diamond coordinates. Statements moved by maps of their own keep sharing their loops, because loopy gives the statements of a loop one domain: the loops run over the union of the images, and each statement tests that a point is its own. For the acoustic pair every point is one statement's; for maps whose images leave holes between them, such as S0 at 2t and S1 at 2t + 1 along the diamond, the loops run over the hull of the union and visit the holes. A map the rewrite cannot write for loopy, such as one over a row and the ragged fiber inside it, or maps that move two statements of one ragged fiber, or of two fibers of one row, differently, is a refuted buildable fact. See note 13 in docs/loopy-notes.md.
• realize(var, tree=True) checks and marks the reassociation; the reduction tree itself comes from splitting and tagging the reduction iname, which is checked separately and not verified on the C target.
• An update's cell read from an array is any cell to isl, so a step that runs two updates through it at the same time is refused even where no input can make them collide (y[perm[i]] += x[i] of a permutation perm), with a pair from the example inputs where they have one and isl's otherwise, which the message says: make the updates atomic (#202). An exact update, m[cell[c, k]] += 1 of integers, keeps the order of its updates and is refused atomic, as an exact sum is refused a tree, although integer addition gives the same bits in any order (#200).
• The target-capability check knows the limits of loopy 2025.2 listed in notes 6, 11 and 14 of docs/loopy-notes.md and no others, so it is a list rather than a model of what the backend can do. A schedule it passes can still fail in code generation for a reason nobody has met yet.
• A statement with no loop on a hardware axis that a loop or a sum of another statement is on runs on every work item of it, so every dependence to or from it is refused, even where each work item reads back only what it wrote itself (every work item storing one value into one cell, say; loopy refuses such a statement for the axis it lacks in any case). The body of a sum on such an axis is taken to read on every work item of it, though each read happens on the work item of the sum's loop, and the work item a sum's statement stores its result from is taken to be one unknown work item, the same for every such statement, though loopy always uses the first. Both are coarser than the check could be, and refuse more. Without a kernel to read a loop's start off (after an affine step the kernel rewrite could not write), only two instances of one loop are taken to share a work item.
• Zeros a read sees follow a loop that stores the array, and not a split of one (#289): after split(k.i, 4).tag({k.i.outer: "g.0", k.i.inner: "l.0"}), a statement stores f[4*k.i.outer + k.i.inner], along no loop of its own, so zeros a later read sees (padding at a boundary) stay on every work item and the tag is refused, saying which read sees them. Zeros no read sees are not run at all (#290), so the split is legal for every array a program stores before it reads.
• A kernel called with a zero-length shape-bearing argument (a matrix with no rows at all) cannot run on the C target: loopy cannot pass an empty array, and the workaround that rescues the empty flat buffer of a ragged axis cannot be applied to an argument loopy reads a size from. A matrix whose rows are all empty does run. See docs/loopy-notes.md.
• Array shapes are not checked at the executor boundary, though element types, ragged layouts and aliasing now are. Lowering has to declare some arrays with shape=None (see docs/loopy-notes.md), so a wrongly sized array is undefined behaviour rather than an error.
• A reduction binder keeps its name in the generated kernel only where nothing else has it. loopy gives an iname one domain, a reduction cannot carry a predicate, and two statements cannot share a reduction's loop, so a reduction whose binder another statement already uses (as a loop variable or a binder), or that its own statement binds over another domain, is lowered under a fresh iname (j_0), and the binders nested in it follow. The name in the term is unchanged, and so is every message, and a Schedule step names such a reduction by its handle, split(k.j, 2), which reaches the iname Lowering.reduction_inames lists; the generated code has the fresh one.
• The trace-time refusals of hidden state look one level below a name: into the containers and the objects it holds (their __dict__ and their slots), the buffers of the arrays it holds (both of a ragged one), and the containers those objects hold. acc[0][0] += 1, holder.inner.s = ..., a deque, a loop over a generator that wraps a domain, and a dir() probe trace without an error, and the trace-faithful fact is what refutes them. That fact is a test, not a proof: it compares the runs on the inputs it tries, so hidden state no such input exercises goes unseen. It stays assumed, with the reason, when no input runs natively, when the term calls a function the interpreter has no numpy counterpart for, or when one ragged bound of a statement bounds two of its loops (a fiber loop inside another over the same row), or a loop and a sum inside it over the same row, and reads an array the kernel writes, which the body reads where each loop or sum starts and the interpreter reads once.

Not yet.

• The OpenCL target is never executed on a development machine: nothing here imports pyopencl, and the test suite asserts it. Device runs happen on real hardware elsewhere and are reported under docs/device-runs/.
• CUDA, and targets beyond C and OpenCL.
• The two-level FMM point-to-point shape, which is a dependent sum of dependent sums and needs more than the two ragged axes lowering handles.

Known limits

The open issues you are most likely to meet, and what to do about each. Most show in lanky check or loopty run as a refuted or assumed fact or as an error, but loopty run compares the compiled run with the native one only on the inputs it runs.

• A compiled run does not check indices (#171). loopty run and LoopyExecutor.run compile a kernel whatever its in-bounds facts say. Where an index leaves its array, the native run raises IndexError and the compiled run reads or writes outside the buffer, with no error. lanky check shows where this can happen: the access's row, x[i + 4] is in bounds for every instance of y[i] = x[i + 4] (shift.py:9), reads refuted or assumed. Put the access under a condition isl can read, which decides the fact: if i + 4 < x.dom.size: before y[i] = x[i + 4].
• A numpy scalar on the left of an operation loses its type (#167). np.float64(0.1) * x[i] is traced as 0.1 * x[i], so beside a float32 x the compiled run multiplies in single precision where numpy uses double. The difference is within the tolerance loopty run allows a float32 input, so it reads tested. Put the scalar on the right: x[i] * np.float64(0.1).
• Loop arithmetic outside a subscript is 32 bits wide compiled (#162). A sum of loop variables and sizes used as a value, y[i] = i + n, (i + n) // 2 or 1.0 * (i + n), wraps round compiled once it passes 2**31, which takes a size past about 2**30; natively it is exact. loopty run reports it only on an input that large. Make the first operand a 64-bit integer: y[i] = i * np.int64(1) + n. A loop bound and a condition on the loops, if i + n > 5:, stay 32 bits wide: keep their sums below 2**31.
• Integers past 64 bits (#139). Arithmetic of loop variables, sizes and literals alone is exact natively, as Python's integers are, and 64 bits wide compiled: 1.0 * i ** 5 wraps round compiled past i = 6208, which loopty run reports on an input that reaches it. Where you want a real, compute in one from the start: (1.0 * i) ** 5.
• Values too large for a narrow numpy integer (#166, #159, #177). numpy raises when a value does not fit an integer type such as np.int16, and loopty refuses a compiled run before it starts wherever isl finds such a value at the call's sizes. The check misses a product of loop variables, a[i] + i * i of an np.int16 a, loop arithmetic in a ragged row's loop, and a uint64 past 2**63 stored into an Int array: the compiled run wraps these round, and loopty run raises the native run's error only on an input that reaches one. It also refuses one case numpy completes, a sum whose if clause repeats the statement's guard. In each case, widen the narrow value first, a[i] * np.int64(1) + i * i, which both runs compute alike, and store a uint64 into an np.uint64 array.
• A few spellings fail without naming the fix (#172, #173, #175). A sum in a guard, if reduce_sum(x[j] for j in x.dom) > 1:, fails inside loopy with name already used as iname: compute the sum into an array in a statement of its own and read that array in the condition. A truth value first in arithmetic, True - k[i], fails while tracing with Python's TypeError: write the integer, 1 - k[i]. k[i] - -(2**63) of an Int is refused as a literal past 64 bits, though numpy computes it: write k[i] + -(2**63), the same value in int64.
• A refused program has written its outputs (#176). When a compiled program raises RefusedNatively for a call that follows one of its checks between calls, the whole program has already run: that call and the calls after it have written their outputs, where natively the program stopped inside that call. After the error, do not trust the arrays, and pass copies if you need them as they were.
• A window keeps values stored, not updated (#213, #214, #229, #133). Schedule.contract refuses an array a statement updates (f[j] += v adds to the zeros a window drops: store each cell once, or keep the array whole), and a loop tagged with no order before it; an update of another array may read a window, and atomic does not lift the order of its slots. A cell stored twice is asked about at both stores, not only at the last before the read, so a window that is safe can be refused there. A window of a size, f_window[n], is one for every work item (in global memory on OpenCL), so a loop it drops cannot go on a hardware axis after it; a window of a number of cells is each work item's own. An array kept in a window, and another computed from it where it is read, are refused together, in either order: substitute moves the reads of the window to where the other is read. A substitution through two intermediates, one computed from the other, is refused, unless the first is computed from arrays nothing in the program writes (grudge's v.n from the normals, in examples/advection.py): substitute one of them.
• Some facts are isl's alone (docs/lean-checks.md). Where Lean is installed, a cast through the part of a tile isl computes stays decided isl: loopty states it for isl only, and lanky check --verbose or loopty run --verbose says so. A cast refused through a substitution's carried read, or a window refused at a read in a sum, is refuted by isl with no counterexample stated for Lean (#215).
• A declared layout stands on its own (#216). Its buffer is one axis, -> Fin[k], and an argument of a kernel does not name the layout it is stored through: a kernel reads through the layout's function, and its facts do not rest on the layout's.
• The OpenCL target is not run in CI (#174). --target opencl, with the opencl extra, generates OpenCL and runs it through pyopencl, and nothing in CI runs a device. A schedule for "opencl" is checked anywhere: where pyopencl is not installed it is built on loopy's plain OpenCL target, its facts are in the ledger, loopty run --emit-code prints the OpenCL C a device runs, and loopty run says it ran none (examples/l2p.py's on_a_gpu). The demos were run on two device classes by hand, and docs/device-runs.md has the commands and the transcripts. Signed integer overflow is undefined in OpenCL C, which has no -fwrapv, so a signed integer sum that overflows, which numpy wraps round, can come out otherwise there: run such a kernel on the C target.
• A theorem's or a postcondition's statement is lanky's reading of it. loopty traces a kernel's body itself and checks that trace, the trace-faithful fact. A theorem, and the claim after a kernel's ->, are annotations lanky reads. lanky 0.2.0 checks a theorem's reading at drawn values, the row the term computes what the annotations compute under it, and its Known limits say what that check can miss. A postcondition has no such row: read its STATEMENT column in the ledger.

The Partial and Not yet lists under Status have the rest, and the open issues have every one.

Install

pip install loopty

That is the 0.4.0 release, whose spelling this page uses, with the five demos at the top of it (Run it). loopty runs on Python 3.12 and 3.13: the islpy it is held to has no wheels for 3.14. CI, and every run this README reports, is on Linux. In a uv project, uv add loopty, and the same for the extra below. The install pulls in:

• lanky 0.4.0 or later, the proof host: lanky check prints the ledger, and loopty plugs into it;
• loopy 2025.2, which generates the code. loopty builds on loopy's internals, so a release of loopty is held to the loopy release it is tested with;
• islpy from 2025.2.5 and below 2026, since islpy 2026 removed a method loopy 2025.2 calls (note 4 in docs/loopy-notes.md);
• numpy 2 or later, whose promotion rules (NEP 50) the compiled run follows, and pymbolic.

A compiled run, loopty run or LoopyExecutor, also needs a C compiler, whatever gcc is on the PATH: GCC on Linux, and on macOS the gcc of Apple's command-line tools, which is clang. lanky check needs none, and neither does calling a kernel natively; the demos under python compile a schedule as well, so they do.

Device execution is an extra, which installs pyopencl, for loopty run FILE --target opencl and LoopyExecutor(target="opencl"). The device and its OpenCL driver are yours to provide, and CI never runs this path:

pip install "loopty[opencl]"

lanky's Lean oracle is an extra of lanky, so a ledger that says proved lean needs pip install "lanky[lean]" and a Lean toolchain in the same environment, and a theorem about real numbers needs Mathlib too (lanky's README has the details). The induction steps of a theorem for every order are found by sympy, lanky's cas extra. Without them a fact isl decided reads decided isl, a theorem tested or assumed, and a WITHOUT LEAN, WITHOUT MATHLIB or WITHOUT SYMPY line under the table names the rows each piece would prove and the command that installs it.

The examples and the documentation are in the repository and not in the package; clone it to follow them.

For work on loopty itself, lanky is resolved from a sibling checkout:

git clone https://github.com/xywei/lanky.git     # next to this one
uv sync --group dev
uv run pytest -q
uv run ruff check .

The dev group has lanky's cas extra, as "Run it" installs it, and not its lean extra, as CI runs, so the transcripts in the documentation are what a reader's install prints before Lean is set up. uv run needs the sibling checkout: in a clone with none beside it, it stops with Distribution not found at: file:///.../lanky.

CI clones lanky's branch of the same name as the one under test when there is one, and lanky's main otherwise, so a change that needs both repositories is tested as one.

A kernel file needs from __future__ import annotations and a ruff F821 per-file ignore, because a size such as n in Arr[Fin[n], Real] is a symbolic variable lanky invents while evaluating the annotation. The same scope invents float and int, so a sort is written Real, Nat, Int, Fin[n] or a numpy type such as np.float64; a sort that is a free name is refused when the kernel is traced.

Architecture

Types are isl objects. A statement's type is its iteration domain, an isl set; its effects are its read, write and accumulation footprints, isl maps. Dependences are derived from the footprints, not declared. Every statement instance is a point of one padded, statement-tagged space, so schedules and dependences are plain maps and the checker is a handful of isl calls.

Index types carry shape. Fin[n] is an index type, Fin[a*b] normalizes to Fin[a] x Fin[b], and a ragged axis is a dependent sum whose bound is another array's entry. An array's index set may also be a polyhedral domain, an isl set it is compared with exactly. Layout and RaggedLayout are the maps from index space to storage, and a domain's box and its packed rows are two more.

Evaluate annotations, trace bodies. No Python parser and no AST pass. Annotations are evaluated with lanky's scope, and the body is run once against symbolic arrays. That is why the native run and the compiled run are the same program. Tracing is the only thing that gives a body its meaning, and whether the trace kept that meaning is a checked fact rather than an assumption: the term is interpreted on its own and compared with the native run.

Transformations are casts. Each states a reindexing map, which isl checks for bijectivity, and a new execution order, which isl checks for monotonicity on the dependence relation. Failure is an IllegalCast carrying the witness (by the statements' ids), its instances as the message names them (by the statements' text and line), and the refuted fact, whose reason is the exception's message. Casts that change floating-point semantics mark the result's exactness class instead of being refused.

loopty is a lanky plugin. It registers a theory (KernelTheory, which turns a kernel into facts), an oracle (IslOracle, trust class decision-procedure), an executor (LoopyExecutor, trust class test) and a verb (run), through the lanky.theories, lanky.oracles, lanky.executors and lanky.verbs entry-point groups. lanky never imports loopty; it finds these and asks each what it can do. Obligations loopty cannot decide stay ASSUMED in the ledger, or become theorems a person proves.

Name

loopty is loop + ty, for types: loops, typed. It follows loopy, sumpy, and pytato in the naming tradition of the loopy ecosystem.

Documentation

• docs/quickstart.md: the two demos end to end, with the output the commands actually print.
• examples/README.md: all nineteen demos, with every console block regenerated by scripts/refresh_example_outputs.py.
• docs/device-runs.md and docs/device-runs/: the demos run on real OpenCL devices, with the commands and the transcripts.
• docs/lean-checks.md: how a fact isl decides is stated for Lean and proved, what isl hands Lean as hints, what stays trusted, and the coverage over the examples.
• docs/loopy-notes.md: the loopy and islpy interactions that cost debugging time, each with its local workaround and why it is local.
• CHANGELOG.md.
• lanky: the proof host loopty plugs into, and where the ledger, the statuses and the oracle protocol are defined.

AI disclosure

This project is developed with substantial assistance from AI coding agents (Anthropic's Claude, via Claude Code). Design, direction, and review are by Xiaoyu Wei. Generated text and code are reviewed before release, but readers should assume AI involvement throughout.

License

MIT. See LICENSE.