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

-- Some way to conditionally load the ForwardDiff code if that module is in
   scope? I'd really like to keep this package completely dependency-free.

