Fix residuals space of modified King fits and add reinit!() checks

[JamesWrigley]

Some more things found by auditing with Claude:

• d683795: modified King fits are solved using a slightly different function (I assume for numerical stability) and the residuals were calculated in that function space. Now we return the residuals in velocity space, which is consistent with the rest of the codebase such that res == y - ŷ holds.
• c8fd4e7: removes some unused code and cleans up some tests.
• e14ecfd: adds checks to reinit!() for nonlinear fits so that it will explicitly error out if the input lengths are wrong or unsupported by the cache.

Checklist

• Appropriate tests were added
• Any code changes were done in a way that does not break public API
• All documentation related to code changes were updated
• The new code follows the
  contributor guidelines, in particular the SciML Style Guide and
  COLPRAC.
• Any new documentation only uses public API

[JamesWrigley]

Sneaked in one more fix in 49bbb66, which fixes #92.

[JamesWrigley]

Two more fixes found by Fable:

• 14b6ac6: expsum fits were modifying the inputs during the solve 🫠 It did invert the transform but it still causes floating point error. Now we copy and modify the inputs instead.
• 7976513: when passing y=nothing to nonlinear fits, the objective function is treated as returning the residual directly. But in such cases sigma would be silently dropped from but kept in the problem, which meant that the statistics methods would wrongly assume that sigma was used and calculate the wrong results. Now we always use sigma. Also made some other minor fixes for in-place/y=nothing solutions.