loopty

loopy, with types.

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.

Three commands, all of them working today:

$ 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   EFFECTIVE  BY             WHERE        OWNER          STATEMENT
-------  ---------  -------------  -----------  -------------  ------------------------------------------------------------------------
decided  decided    isl            spmv.py:79   scan           off[0] is in bounds for every instance of S0
decided  decided    isl            spmv.py:81   scan           off[r + 1] is in bounds for every instance of S1
decided  decided    isl            spmv.py:79   scan           distinct instances of S0 write distinct cells of off
tested   tested     native         spmv.py:69   scan           off[0] == 0 and (forall r in Fin(n). off[r + 1] == off[r] + cnt[r])
tested   tested     property-test  spmv.py:84   scan_monotone  n : Nat, cnt : Fn[Fin(n), Nat], off : Fn[Fin(n + 1), Nat] | off(0) ==...
decided  decided    isl            spmv.py:112  spmv           y[r] is in bounds for every instance of S0
decided  decided    type           spmv.py:112  spmv           x[col[r, j]] is in bounds by type (col[r, j] : Fin(m))
decided  decided    isl            spmv.py:112  spmv           distinct instances of S0 write distinct cells of y
decided  decided    type           spmv.py:112  spmv           the accumulation into y[r] over j is approx
tested   tested     interpreter    spmv.py:102  spmv           the traced term computes what the body computes
decided  tested     call           spmv.py:115  solve          after scan(...) in solve: off[0] == 0 and (forall r in Fin(n). off[r ...
...
21 facts: 16 decided, 5 tested

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 tested by native: it is evaluated at what every native run of scan left in its arguments, on the file's example inputs and on the drawn ones, and held after each. 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 EFFECTIVE column reads tested. 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('t', 'i', 8, 8) ->
  IllegalCast: tile(t,i,8,8) illegal: instance S0[t=0, i=8] writes u[1, 8] read by S0[t=1, i=7] scheduled earlier (at nt=16, nx=16, as hinted)
...
accepted: Schedule(jacobi, target='c').skew(i, by='t').tile(t,i,8,8)
  loop nest: t_outer i_outer t_inner i_inner
  decided  isl  skew(i, by='t') renames the instances of jacobi one for one
  decided  isl  the order after skew(i, by='t') runs every dependence of jacobi forward
  decided  isl  tile(t,i,8,8) renames the instances of jacobi one for one
  decided  isl  the order after tile(t,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. 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 in f.dom:
        for j in f.dom[i]:
            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:78  pairs  f[i, j] is in bounds for every instance of S0
...
decided  isl          pairs.py:80  pairs  f[p, j] is in bounds for every instance of S1
decided  isl          pairs.py:80  pairs  f[k, p] is in bounds for every instance of S1
...
16 facts: 15 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("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, or kept.
• 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("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.
• 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.
• 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.1.0, the first release. 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. Below is what works, what works in part and what does not work yet, and after it, under Known limits, the open issues you are most likely to meet, one of them a fact the ledger calls decided that is false: read that before you rely on a ledger. CI runs the suite on every change, on Python 3.12 and 3.13 and against lanky's main, 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 are named after it (scan.S1, step@2.S0). 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 body that does anything to an argument but pass it to a kernel, or make an array like it, 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 (shifted.S0 reads the zeros at f[0] and f[n - 1]), so zero padding at a boundary is read as written and lanky check passes. Schedule(program).fuse("flux", "divergence", shift=1) runs the reader's loop inside the writer's, one step behind, and is affine with the map per statement it builds, { flux_S0[j] -> [j]; divergence_S0[i] -> [j] : 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 shift that is accepted. .substitute("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. examples/fusion.py is the Burgers flux and divergence of examples/composition.py, fused, then substituted, each compiled run compared with the native one.
• Tracing a body to a typed term: accesses, statements, reductions, ragged fibers, when guards, source locations, and a TraceError that names the fix when a Python if is used on a computed value, 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.
• when guards 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 guard 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 (when(i * 0.5 < 2), when(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 (#148; note 23 in docs/loopy-notes.md). A guard has to be a truth value, and one whose value is an integer is a TraceError, on a native run and under tracing alike. Natively that is what ~(i > 0) is (~ on a Python bool is bitwise), while the trace records not (i > 0), so such a kernel traces and its trace-faithful fact is refuted by the native refusal, which names the fix. 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 as it is by when. 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 a truth value is asked for. 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.
• Schedule: tag, split, interchange, prioritize, tile, skew, affine, fuse, realize, each checked as a cast, each emitting its fact, and substitute, checked as a storage decision; retarget, which replays every step against another loopy target and re-checks it. 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, 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 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. abs of an integer is numpy's, the smallest value of its type included; a comparison of one does not trace yet (abs(k[i]) > 2, 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.
• The argument contract, enforced at every entry point that runs a kernel (compiled, differential and native): two distinct array parameters may not share storage, 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("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 tested against its native runs, on the file's example inputs and on drawn ones, never decided: deciding it from the term, by which statement last writes each cell, is not done yet. @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 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, 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), and the array is either stored in full or not at all: contracting it to the cells live at once is not done. 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 and a callee with an array over a Where, Sigma or union domain. 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 a temporary is a variable-length array on the stack of the call, which bounds its size (note 16 in docs/loopy-notes.md); on OpenCL it is a global temporary, which is generated but, like every device path, not run from a development machine.
• 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, and so is a sum of loop variables, sizes and small literals elsewhere; 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: x[(i * i) % n] (#129), and x[(i * 499999) // 1000000], which loopy wrote with a division of the wrapped 499999 * i. A sum the plan leaves in 32 bits stays so in a division too, x[(i + n) // 2] past 2**30 (#149). A subscript over a narrow numpy integer type is not converted back into it (#157). 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. pymbolic builds u[i] - 1 and u[i] + -1 as one term, read as the difference: of an unsigned u numpy refuses the second at every point, and the compiled run computes it as the first; a - b and a + -b of an int8 b likewise differ natively at b = -128 alone, and of an Int a and an unsigned b at every b but zero (#151). A loop variable's arithmetic beside a numpy integer whose type does not hold its value at some point, u[i] + (i - 3) of a uint32 u, is refused natively there and computed compiled (#153), and abs of a real in a guard fails inside loopy (#154). See notes 19, 20 and 23 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.
• 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.
• The accumulation convention: a traced y[r] += ... under a parallel iname is reported as a disjointness refutation, which is the conservative reading. The reassoc fact is what should license it and nothing consumes that yet.
• 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.
• 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, but a Schedule step names such a reduction by its iname in the kernel (split("j_0", 2)), which Lowering.reduction_inames lists.
• 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 a user is most likely to meet. Each says what the ledger and the runs make of it, since that is what decides whether it can go unnoticed.

• Some index arithmetic is 32 bits wide compiled (#149, #157). A loop bound, and a subscript loopy reads as affine with no division, are loopy's index arithmetic, 32 bits wide. Every other subscript is computed in 64 bits, as natively. A sum of loop variables and sizes inside a dividing subscript stays in 32 bits, so x[(i + n) // 2] wraps round past 2**30 (#149). A subscript over a narrow numpy integer type, such as int8 or uint16, is computed in C's int compiled, and wraps round natively (#157).
• A guard loopy's bounds check cannot read (#148). An access that only a floating-point guard keeps in bounds, x[2 * i] under when(i * 2.0 < n), is refused by loopy's bounds check, with a message that names no fix. Write the guard over integers.
• Integers past 64 bits (#139, #159). Arithmetic of loop variables, sizes and literals alone is exact natively and 64 bits wide compiled: 1.0 * i ** 5 wraps round compiled past i = 6208, and the differential fact is refuted. A uint64 value past 2**63 stored into an Int array is refused natively and wraps round compiled.
• What passes between calls is stored in full or not at all (#132, #133). Schedule.substitute computes an intermediate again at every read; keeping only the cells live at once, a window of three cells of f after the Burgers fusion, is not done. Nor is a substitution through two intermediates, one computed from the other.
• A program's intermediates live on the C stack (#62). An array a program makes with Arr.zeros_like is a variable-length array on the stack of the compiled call, so one larger than the stack crashes the compiled run, and nothing checks the size first.
• The OpenCL target is not run in CI. --target opencl, with the opencl extra, generates OpenCL and runs it through pyopencl, and nothing in CI or in the test suite runs a device. The demos were run on two device classes by hand, and docs/device-runs.md has the commands and the transcripts. A comparison of complex values does not build on it (#158).
• A theorem's or a postcondition's statement is lanky's reading of it. loopty traces a kernel's body itself, and the trace-faithful fact checks that trace. A theorem, and the claim after a kernel's ->, are annotations lanky evaluates, and lanky's Known limits apply to them: an operation that answers from a term object rather than through its overloads can make the statement say something other than what was written. Read the ledger's STATEMENT column.

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

Install

pip install loopty

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.1.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 (lanky's README has the details). Without it the same theorem reads tested property-test, and every other row is unchanged.

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 .

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 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 nine 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/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.