
-- Is the global state the most elegant way to achieve the desired flexibility?
   Apparently julia 1.6 comes with scratch space, which might be better.

-- Restriced nll/REML (presumably rnll.jl file)?
  -- Numerically thoughtful way for obtaining contrast basis
  -- Numerically thoughtful way to minimize matrix-matrix arithmetic

-- Include a parametric mean. Obviously very easy for the linear algebra, but
   the question is really about the best interface for communicating what
   indices are for what part of the distribution.

-- Add tests for the many parameters and varying workspace count setting.

-- Add a "presolve!" function to compute [K\K_j for K_j in Kdv] in the exact
   derivative cases.

