Find the best fit with some conditions

Question

I have this code to find the best values for $l$, $s$ and $j$:

Clear[j,l,s,norm,maxx,maxy];
data=Import["https://pastebin.com/raw/2DG5Xes6","Table"];
g=3/2+(s(s+1)-l(l+1))/(2j*(j+1));
\[Mu]=9.274*10^-24;k=1.380*10^-23;
y=\[Mu]*g*j*x/k;
maxy=Max[data[[All,2]]];maxx=Max[data[[All,1]]];minx=Min[data[[All,1]]];
conds={Mod[l,1]==0&&Mod[j,1/2]==0,j-s==0||j-(l+s)==0||j-Abs[l-s]==0};
b[x_]=maxy*(((2j+1)/(2j))Coth[(y(2j+1))/(2j)]-(1/(2j))Coth[y/(2j)]);
fit=FindFit[data,{b[x],conds},{l,j,s},x]

$l$ is an integer and $j$ and $s$ are half-integers. One of these conditions must hold:

j-s==0||j-(l+s)==0||j-Abs[l-s]==0

I tried to fit the data using all those conditions but the result was $l=s=j=1$, which is not the best fit. I happen to know the correct parameters for this case $(l=0,s=j=3/2)$ and if I use those as the initial guesses, I do find the correct fit. Is it possible to rewrite the conditions so that Mathematica gives the best fit automatically?

Answer 1

I'm not sure why FindFit can't manage it, but if I construct the objective function myself to minimize square residuals and I use Method -> "DifferentialEvolution" then I get the answer:

objective = Total[(#[[2]] - b[#[[1]]])^2 & /@ data];
fit = Last[NMinimize[{objective, conds}, {l, j, s}, 
  Method -> "DifferentialEvolution"]]//Chop

(* result: {l -> 0, j -> 1.5, s -> 1.5} *)

However, using FindFit's Method->NMinimize I was unable to achieve the above result.

Answer 2

I am not certain why there is a hiccup with the fitting, generally both FindFit and NonlinearModelFit are quite reliable if you set your conditions right. I found that the issue resides with the "or" conditionals. If you request only one of the 'or' conditions (j - Abs[l - s] == 0) then you can get your answer. The other conditions don't seem to satisfy the model even though, weirdly, both j-s == 0 and j - (l + s) == 0 give perfectly valid quantum number conditions for your established solution.

ifs = {j - Abs[l - s] == 0, j - s == 0, j - (l + s) == 0};
fit = NonlinearModelFit[
    data, {b[x], Mod[j, 1/2] == 0 && Mod[l, 1] == 0, ifs[[1]]}, {l, j,
      s}, x, Method -> "NMinimize"] // Chop ;
fit2 = FindFit[
    data, {b[x], Mod[j, 1/2] == 0 && Mod[l, 1] == 0, ifs[[1]]}, {l, j,
      s}, x, Method -> "NMinimize"] // Chop ;
fitvals = fit["BestFitParameters"];
Print["Using NonlinearModelFit: ", fitvals]
Print["Using FindFit: ", fit2]

Show[{Plot[b[x] /. fitvals, {x, minx, maxx}], data // ListPlot}]

Output:

You could in principle evaluate each "or" condition separately and only catch the output that doesn't return an error. I hope this helps.