We return the sqrt of the covariance matrix, which is nonsense...

Check the code at core.py. We may have to return the diagonal and reshape, to give the same return as for 'dif'. Maybe this has consequences elsewhere, so better check thoroughly...