`First function call produced NaNs` when using NUTS but not MH

[devmotion]

@devmotion

As yet another alternative, probably it would be sufficient to add a check for i == j || (xj = SVector(u[1, j], u[2, j])) == xi in the inner for loop in cuckersmale! since the contribution from these cases where xj == xi (which might occur even if i != j) to totx and toty is zero anyway, it seems.

On a second thought, that would only work if xi == xj would imply vi == vj (which is not guaranteed in general, I assume?). Generally, if xi == xj, the contribution to totx and toty would be K * (vj[1] - vi[1]) and K * (vj[2] - vi[2]), which might be nonzero if vi != vj (which, of course, is only possible if i != j).

So probably the safest and easiest option is to enable NaNsafe-mode. Skipping the cases i != j would still be reasonable, IMO, but it seems it would not ensure that you don’t run into these gradient issues.