Major rearrangement of code in advance of dramatic improvements to
computation time. Modular computation of density.basis which can be
computed outside of sum.log.density and sum.log.density.gradient then
be passed to optim() so that all optim does is i) simple
multiplication of P%*%beta and gradient sums and integrates. By
avoiding repeated computation of P this _must_ be orders of magnitude
faster. Perhaps we will approach the impressive computation time of
logspline yet?
By removing non-essential computation outside of optim reduced runtime
from 89 seconds to 6 seconds for n=1000 and ought to scale better to boot!