Scipy.optimize.curve_fit(): "ValueError("`x0` is infeasible.")"

Yesterday I struggled with the error message:
"ValueError("x0 is infeasible.")"
of the function scipy.optimize.curve_fit(), this is my contribution for others that start to work with this well tested method. If you work with measured data and you would like to utilize the feature of boundary (which often make sense) you should spent effort to ensure three things:

  • The initial guess should make sense (best is a set of start values close to the expected optimum)

  • select sound boundaries (a wide as necessary as narrow as suitable)

  • Make sure that your initial guess lies inside your boundaries

To avoid the error message and a search in the net, I would like to share my example
which includes a check of boundaries and initial guess. I hope it helps others to get started:

# --- example for multy variable optimization via SciPy.optimize.curve_fit: ------------------------------------------------
using PlotlyJS, RobustModels
using SciPy

## --- compose signals: clean and disturbed --------------------------------------------------------------------------------
sampling_rate   = 10000                    # Sampling frequency
n_signal        = 3000                     # Length of signal
off_set_sign    = 1.0                      # vertical offset of signal
frequ_A         = 50.0                     # unit = Hz
ampl_A          = 0.7                      # Amplitude @f1
phase_A         = 0.1 * pi                 # phase offset
noise_ampl      = 0.1                      # max noise amplitude
b_symmetry_on   = false                    # substract arithmetic mean to make it symetric around zero

# --- create time vector:
delta_t         = 1.0 / sampling_rate      # Time step in sec
t_vec           = collect(range(0, step=delta_t, length=n_signal))  # Time vector

# create signal with sinus at frequ_A:
signal_clean     = off_set_sign .+ ampl_A .* sin.((2 * pi * frequ_A) .* t_vec .+ phase_A)  
signal_disturbed = signal_clean + noise_ampl * randn(size(t_vec))
if b_symmetry_on
    signal_disturbed = signal_disturbed .- mean(signal_disturbed)
end

## --- local structs: ------------------------------------------------------------------------------------------------------
mutable struct _MySineFourParStruct
    frequency :: Number
    amplitude :: Number
    phase     :: Number
    offset    :: Number
end
## --- local functions -----------------------------------------------------------------------------------------------------
function _PtsPerXperiods(_frequency::Real, _sampl_rate::Real, _num_periods::Int=1)
    if _sampl_rate < 10 * _frequency
        error("sampling rate too low, it needs to be at least ten times the investigated frequency!")
    end
    _pts_of_n_periods = round(Int, _num_periods * _sampl_rate / _frequency)
    return _pts_of_n_periods 
end


function _fit_cosinus(_SineFourPar::_MySineFourParStruct, _data_pts::Vector{<:Number}, _sampl_rate::Number)
    # ---
    n_data = length(_data_pts)
    t_vec = collect(range(0, step = 1/_sampl_rate, length = n_data))

    # --- constants
    i_frequ = 1;     i_ampl = 2;     i_phase = 3;     i_offset = 4
    
    # --- fitting function
    function _model_fun(_t, _frequ, _ampl, _phi, _offset)
        return _ampl .* cos.(2 * pi * _frequ .* _t .+ _phi) .+ _offset
    end
    # --- initial values
    _p0 = [_SineFourPar.frequency, _SineFourPar.amplitude, _SineFourPar.phase, _SineFourPar.offset] # initial values for fit
    # bounds        p1= frequ,               p2= ampl,                              p3= phase, p4= offset 
    _lower_bounds = [0.9999 * _SineFourPar.frequency, max(0.01, 0.9 * _SineFourPar.amplitude), -2 * pi, 0.0]
    _upper_bounds = [1.0001 * _SineFourPar.frequency, max(0.03, 1.1 * _SineFourPar.amplitude), +2 * pi, max(1.0, 2 * _SineFourPar.offset)]
    _pbounds = [_lower_bounds, _upper_bounds]

    # --- check boundaries
    any_equal_or_lower = _upper_bounds .<= _lower_bounds
    indx_upper_equal_lower = findfirst(any_equal_or_lower)
    if ~(indx_upper_equal_lower==nothing)
        @info("Upper is equal or lower as lower boundary at index: ",  indx_upper_equal_lower, ",\t _lower_bounds[]: ", _lower_bounds[indx_upper_equal_lower], ", \t _upper_bounds[]: ", _upper_bounds[indx_upper_equal_lower])
    end
    
    any_equal_or_lower = _p0 .<= _lower_bounds
    indx_p0_equal_lower = findfirst(any_equal_or_lower)
    if ~(indx_p0_equal_lower==nothing)
        @info("_p0 is equal or lower as lower boundry at index: ",  indx_p0_equal_lower, ",\t _p0[]: ", _p0[indx_p0_equal_lower], ", \t _lower_bounds[]: ", _lower_bounds[indx_p0_equal_lower])
    end
    
    any_equal_or_higher = _p0 .>= _upper_bounds
    indx_p0_equal_higher = findfirst(any_equal_or_higher)
    if ~(indx_p0_equal_higher==nothing)
        @info("_p0 is equal or higher as upper boundry at index: ",  indx_p0_equal_higher, ",\t _p0[]: ", _p0[indx_p0_equal_higher], ", \t _upper_bounds[]: ", _upper_bounds[indx_p0_equal_higher])
    end
    
    # --- call optimizer:
    _fit_python_param, _ = SciPy.optimize.curve_fit(_model_fun, t_vec, signal_disturbed, p0 = _p0, bounds = _pbounds)
    
    # ---
    _frequ_opt      = _fit_python_param[i_frequ]
    _ampl           = _fit_python_param[i_ampl]
    _phase          = _fit_python_param[i_phase]
    _offset         = _fit_python_param[i_offset]
    # ---
    _results = _MySineFourParStruct(_frequ_opt, _ampl, _phase, _offset)
    return _results
end

## --- local plot functions ---------------------------------------------------------------------
function plot_clean_disturbed_and_fited_signal(_signal_clean::Vector{<:Number}, _signal_disturbed::Vector{<:Number},
    _signal_fitted::Vector{<:Number}, _sampl_rate::Real, _frequency::Real, _num_periods::Int=1)
    # --- all need to have the same length:
    if ~(length(_signal_clean) ==  length(_signal_disturbed) )
        error("Missmatch of vector length!")
    end
    # --- select part of the signal:
    _num_pts = _PtsPerXperiods(_frequency, _sampl_rate, _num_periods)
    # --- build time vector:
    _time_vec = collect(range(0, step= 1 / _sampl_rate , length= _num_pts))
    # --- build index vector / index range:
    _range = range(1, step= 1, length= _num_pts)
    # --- construct plot objects:
    line_clean     = scatter(; x = _time_vec[_range], y = _signal_clean[_range],     name = "clean")
    line_disturbed = scatter(; x = _time_vec[_range], y = _signal_disturbed[_range], name = "disturbed")
    line_fitted    = scatter(; x = _time_vec[_range], y = _signal_disturbed[_range], name = "fitted")
    _data = [line_clean, line_disturbed, line_fitted]
    _layout = Layout(
        title_text       = "Clean and Disturbed Signal",
        xaxis_title_text = "Time / s",
        yaxis_title_text = "Time Domaine Signal / -"
    )
    return plot(_data, _layout)
end

function _build_fited_signal(_four_param::_MySineFourParStruct, _time_vec::Vector{<:Number})
    _signal_fited = _four_param.offset .+ _four_param.amplitude .* cos.((2 * pi * _four_param.frequency) .* _time_vec .+ _four_param.phase)
    return _signal_fited
end

## --- main ----------------------------------------------------------------------------------------------------------------

_SineFourPar = _MySineFourParStruct(frequ_A, ampl_A, phase_A, off_set_sign)

fit_result = _fit_cosinus(_SineFourPar, signal_disturbed, sampling_rate)

fitted_signal = _build_fited_signal(fit_result, t_vec)

plot_clean_disturbed_and_fited_signal(signal_clean, signal_disturbed, fitted_signal, sampling_rate, frequ_A)

1 Like