I guess there is a bug in the code causing sporadic failures. And the failure is likely captured by the error handling either here:

I suggest commenting out those except xpsiError as e conditions so that we will know where in the code the failure actually happens.

It could be a similar issue as here:
#53

It seems to come from the default_background_marginalisation.pyx
https://github.com/xpsi-group/xpsi/blob/db1a05795d1e30ab4eb55ca3baa535b5a1266e16/xpsi/likelihoods/default_background_marginalisation.pyx#L397C1-L403C42
The function gsl_interp_eval_integ seems to sometimes be giving huge values that then get added to STAR, leading to a likelihood evaluation that is wrong.
One workaround would be to put a condition of len(phases)>2 to compute that function, and to just assign the component value if not ( STAR[i,j] = components[i][j] ), to deal with Everywhere and time_invariant=True where there is only one phase bin. But then using time_invariant = False leads to an error and needs to be solved