Object function can return None with certain parameters, how to skip or avoid while converging in Curvefit？

Question

In my situation, the objective function is a numerical process contains a root finding process for an equation by bisection method. With certain set of parameters, the equation does not have root for a intermediate variable. I thought making the bisection root finding routine return None can solve such problem. As the object function with a set of date being regressed by scipy.optimize.curve_fit with p0 separate by this situation in between, error is then stop the process.

To study this case, a simplified case is shown.

import numpy as np

#Define object function:
def f(x,a1,a2):
    if a1 < 0:
        return None
    elif a2 < 0:
        return np.inf
    else:
        return a1 * x**2 + a2

#Making data:
x = np.linspace(-5,5,10)
i = 0
y = np.empty_like(x)
for xi in x:
    y[i] = f(xi,1,1)
    i += 1

import scipy.optimize as sp

para,pvoc = sp.curve_fit(f,x,y,p0=(-1,1))
#TypeError: unsupported operand type(s) for -: 'NoneType' and 'float'

para,pvoc = sp.curve_fit(f,x,y,p0=(1,-1))
#RuntimeError: Optimal parameters not found: Number of calls to function has reached maxfev = 600.

I also tried inf, and it is obviously not working. What should I return to continue the curve_fit process? Imagine it is trying to converge, what happen does the curve_fit do when it meets such situation.

Additional thinking: I tried the try...except... to ignore the error and also simulate a case that the p0 is in a solvable range, but will pass the unsolvable segment to the true fit.

import numpy as np

def f(x,a1,a2):
    if a1 < 0:
        return None
    elif a1 < 2:
        return a1 * x**2 + a2
    elif a2 < 0:
        return np.inf
    else:
        return a1 * x**2 + a2

def ff(x,a1,a2):
    output = f(x,a1,a2)
    if output == None:
        return 0
    else:
        return output

x = np.linspace(-5,5,10)
i = 0
y = np.empty_like(x)
for xi in x:
    y[i] = f(xi,1,1)
    i += 1


import scipy.optimize as sp

#para,pvoc = sp.curve_fit(f,x,y,p0=(-1,1))
#TypeError: unsupported operand type(s) for -: 'NoneType' and 'float':
#para,pvoc = sp.curve_fit(f,x,y,p0=(1,-1))

try:
    para,pvoc = sp.curve_fit(f,x,y,p0=(-3,1))
except TypeError:
    pass

Obviously error was met during converging and had been reported and was excepted. What should I do to continue the curve_fit with the original converging direction? Even I can make concession, how can I tell the curve_fit to return the last attempt to the a1?

On the other hand, I tried put this try... except... in the object function to return 0 when there is the error. The result is then as I expect:

para,pvoc = sp.curve_fit(ff,x,y,p0=(-3,1))

#OptimizeWarning: Covariance of the parameters could not be estimated
  category=OptimizeWarning)

Answer 1

I think you want to take a different approach. That is, you have written your objective function to return None or Inf to signal when the value for a1 or a2 is out of bounds: a1<0 and a2<0 are not acceptable values for the objective function.

If that is a correct interpretation of what you are trying to do, it would be better to place bounds on both a1 and a2 so that the objective function never gets those values at all. To do that with curve_fit you would need to create a tuple of arrays for lower and upper bounds, with an order matching your p0, so

pout, pcov = sp.curve_fit(f, x, y, p0=(1, 1), bounds=([0, 0], [np.inpf, np.inf])

As an aside: I have no idea why you're using a starting value for a1 that is < 0, and so out of bounds. That seems like you're asking for trouble.

For an even better experience setting bounds on fitting parameters, you might consider using the lmfit, which would allow you to write:

import numpy as np
from lmfit import Model

def f(x, a1, a2):
    return a1 * x**2 + a2

fmod = Model(f)

params = fmod.make_params(a1=1, a2=0.5)
params['a1'].min = 0
params['a2'].min = 0

x = np.linspace(-5, 5, 10)

np.random.seed(0)
y = f(x, 1, 1) + np.random.normal(size=len(x), scale=0.02)

result = fmod.fit(y, params, x=x)
print(result.fit_report())

which will print out

[[Model]]
    Model(f)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 13
    # data points      = 10
    # variables        = 2
    chi-square         = 0.00374066
    reduced chi-square = 4.6758e-04
    Akaike info crit   = -74.9107853
    Bayesian info crit = -74.3056151
[[Variables]]
    a1:  0.99998038 +/- 7.6225e-04 (0.08%) (init = 1)
    a2:  1.01496025 +/- 0.01034565 (1.02%) (init = 0.5)
[[Correlations]] (unreported correlations are < 0.100)
    C(a1, a2) = -0.750

Hope that helps.