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)