add inequality constraints

[pkienzle]

Bumps supports constraints such as parameter a < parameter b by having a penalty cost function evaluated for each parameter set. In this case, the function would be something like:

((b-a)^{2 + penalty) * (a>b)}

That way, if the fit finds itself outside the valid range a⇐b , it will be pushed back into the range with a quadratic function; the penalty value forces a step at the boundary so that a lower chisq out of bounds is not preferred to a value at the bounds. The penalty value should be large enough to remove any minima outside the bounds, otherwise algorithms like DREAM may be able to stumble onto it. Using a faster growing function, such as e^{(b-a) will reduce the problem. If the boundary can be a little less strict, then penalty can be zero.}

The constraints panel needs an extended syntax to handle inequality constraints. From this syntax, a constraint function could be automatically generated for use with bumps.

[pkienzle]

You can define a custom model for soft inequality constraints, returning "x - (C<0)*C". Use it with a data file containing the point (x,y,dy) = (1,1,1). For example, if you have M1.radius_core and M1.radius_shell, you could set

M2.C = M1.radius + M1.radius_core - 100.

If the total radius is less than 100 then C < 0 is false and x=1 is returned. If it is greater than 100, then (C<0) is true and the function returns the theory value 1 - C. The fitter then computes (yth - ydata)²/dy² which is ((1-C)-1)²/1² = C² and adds this as the residual to the fit. Note that we are using (x,y) = (1,1) so that it plots on a log-log graph.

You can extend this to allow multiple constraints by defining a file with N points (1,1,1), (1,1,1), … and fit it using the following plugin model:

parameters = [
   ["C1",  "",  0, [-inf, inf], "", "Constraint 1"],
   ["C2",  "",  0, [-inf, inf], "", "Constraint 2"],
   ["C3",  "",  0, [-inf, inf], "", "Constraint 3"],
   ["C4",  "",  0, [-inf, inf], "", "Constraint 4"],
   ["L1", "",   0, [-inf, inf], "", "Latent parameter 1"],
   ["L2", "",   0, [-inf, inf], "", "Latent parameter 2"],
   ["L3", "",   0, [-inf, inf], "", "Latent parameter 3"],
   ["L4", "",   0, [-inf, inf], "", "Latent parameter 4"],
   ["L5", "",   0, [-inf, inf], "", "Latent parameter 5"],
]

def Iq(q, C1, C2, C3, C4, L1, L2, L3, L4, L5):
     return np.array([1-(C1<0)*C1, 1-(C2<0)*C2, 1-(C3<0)*C3, 1-(C4<0)*C4])
Iq.vectorized = True

The data file would be:

   1 1 1
   1 1 1
   1 1 1
   1 1 1

The parameters L1, …, Lm are available as underlying parameters that can be fitted rather than fitting the model parameters directly. They don't actually contribute to the model output, so they can have an arbitrary value. For example, you can fit the total radius and the core radius by setting

M1.radius_core = M2.L1 - M1.radius_shell

and fitting M2.L1.

You can add hard constraints Hk as well by returning (inf if Hk<0 else 0). This won't work for the Levenberg-Marquardt fitter since it is not differentiable, but DREAM, DE and Nelder-Mead should completely avoid points in the infeasible region.