This is very tentative.

* Useful building blocks for models in [lib/nn_blocks.ml](lib/nn_blocks.ml).

* A language model example.

* Port (translate or bind) the Python files from [llm.c](https://github.com/karpathy/llm.c) to implement tokenization, data loading and saving etc.

* At the end of 0.6.x, we should have an apples-to-apples benchmark comparing OCANNL to [llm.c](https://github.com/karpathy/llm.c) for both CPU and GPU.

* First harvested from [Fast Multidimensional Matrix Multiplication on CPU from Scratch](https://siboehm.com/articles/22/Fast-MMM-on-CPU).

* Then harvested from [How to Optimize a CUDA Matmul Kernel for cuBLAS-like Performance: a Worklog](https://siboehm.com/articles/22/CUDA-MMM).

* Finally from [llm.c](https://github.com/karpathy/llm.c).

* These will require splitting a routine into multiple CUDA kernels.

* E.g. host-device transfers: copy from host if host update is later than the previous device update.

* Concise syntax for transfers into the merge buffer since we know which tensor node is transferred and where to.

* At the end of 0.8.x, OCANNL has a REPL.

* 0.10: Optimize performance: program search.

* Instead of dynamic scheduling as in tinygrad, we can schedule statically by program search.

* We should also reproduce the search that tinygrad is doing.

@@ Shape_error

( "Cannot compare Prod with unresolved variables in inequality",

[ Dim_mismatch [ cur; subr ] ] )

match (Map.find env.dim_env cur_v, Map.find env.dim_env subr_v) with

| Some (Bounds_dim { cur = cur1; _ }), _ when List.mem ~equal:equal_dim_var cur1 subr_v ->

([ Dim_eq { d1 = cur; d2 = subr } ], env)

Map.set env.dim_env ~key:subr_v

~data:(Bounds_dim { lub = Some cur; cur = cur2; subr = subr2; constr = constr2 });

} ))

| Dim _, Dim _ ->

raise

@@ Shape_error ("dimension comparison for axis: mismatch", [ Dim_mismatch [ cur; subr ] ])

val dim_to_int_exn : dim -> int

val dim_to_string : [> `Only_labels ] -> dim -> string

val dim_vars : dim -> dim_var list

During the solution process, the constraints are incorporated, or propagated, into the environment `constr` entry fields, and into further `constraint_` constraints, as needed. This provides sufficient scaffolding to implement the other complex constraints as the need arises.

The constraints are solved by: unification of the equation constraints, unification-like simplification of the inequality constraints, propagation of the complex constraints. Simplification of an inequality, and constraint propagation, can generate more constraints, so we need to be careful to keep it terminating. The solution proceeds in stages.

* Stage 1 is online as tensors are composed, and conservatively performs unification and constraint propagation. Stages 2, 3, 4 are only performed once necessary: when projections or dimensions are requested.

* Stage 2, when solving the constraints, substitutes dim variables in terminal shapes that do not have a LUB or other constraints, by dimension-1. (This is generalized at stage 6 to all variables.) It substitutes row variables in terminal shapes that do not have a LUB by one axis if that's required to satisfy the variable's constraint.

* Stage 4 sets yet-unknown dimensions with >1 lower bounds from direct accesses, to their LUBs if they have any, otherwise to the lower bound.

* Stage 5 addresses `Total_elems` constraints with yet-unknown row variables. If the constraint can be satisfied by assuming the row variable is no-further-axes, it sets the row variable to `Broadcastable`, otherwise it sets it to one axis of the required dimension.

* Stage 6 sets row variables in the remaining inequalities to no-further-axes values. This can unlock further between-axis inequalities because of row variables sandwiched between leftmost axes from their side of the inequality and rightmost axes from the other side of the inequality.

* `finish_inference` is called right before some projections or array dimensions are required (typically, because of jitting). It performs a second round of `propagate_shapes`, and then once again attempts to solve any remaining constraints that `propagate_shapes` didn't solve. Then it "closes the shapes": substitutes out remaining shape variables by their LUBs if any, or dimension-1 / `Broadcastable` (no-more-axes). Then it resets the environment state, since the shapes are now guaranteed to not have variables.

* `derive_projections` starts by freshening the `proj_id`s in the `update_step`. Then it generates and solves shape inequalities, and then generates and solves projection equations, and constructs the `projections` record.

* `of_spec` constructs a shape record from an einsum slot spec. If `deduced = Input_equals_output`, it adds the corresponding equation to the global environment.