Weighted LSQfit fails, Message: DomainError [...] sqrt will only return a complex result if called with a complex argument.

[StefanPofahl]

Here is my example script that fails :-(
Do you have an idea, what causes this problem?

using PlotlyJS, RobustModels
using LsqFit
using StatsBase # priovides rmsd (Return the root mean squared deviation between two optionally normalized arrays.)

## --- compose signals: clean and disturbed --------------------------------------------------------------------------------
sampling_rate   = 10000                    # Sampling frequency
n_signal        = 3000                     # Length of signal
off_set_sign    = 0.4                      # 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 _LSQfit_four_param(_SineFourPar::_MySineFourParStruct, _data_pts::Vector{<:Number}, _sampl_rate::Number, _weighing_vec::Vector{Float64}=Vector{Float64}(undef, 0))
    # ---
    n_data = length(_data_pts)
    _t_vec = collect(range(0, step = 1/_sampl_rate, length = n_data))

    # --- constants
    idx_frequ = 1;   idx_ampl = 2;     idx_phase = 3;     idx_offset = 4
    
    # --- fitting function
    # @. model(x, _ampl, _phi, _offset) = _ampl * cos(2 * pi * _frequ * x + _phi) + _offset
    @. model(x, p)                      = p[2]  * cos(2 * pi * p[1] * x + p[3]) + p[4]

    # --- 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 = [max(0.05, 0.8 * _p0[idx_frequ]),   max(0.001, 0.8 * _p0[idx_ampl]),    -2 * pi,    _p0[idx_offset] - 0.1]
    _upper_bounds = [1.2 * _p0[idx_frequ],              1.2 * _p0[idx_ampl],                +2 * pi,    _p0[idx_offset] + 0.1]
        
    # --- call optimizer:
    if isempty(_weighing_vec)
        _fit = LsqFit.curve_fit(model, _t_vec, _data_pts, _p0, lower = _lower_bounds, upper = _upper_bounds)
    else
        if size(_weighing_vec)[1] == n_data
            _fit = LsqFit.curve_fit(model, _t_vec, _data_pts, _weighing_vec, _p0, lower = _lower_bounds, upper = _upper_bounds)
        else
            error("Weight Array has not the right dimension!")
        end
    end
    # ---
    _frequ          = _fit.param[idx_frequ]
    _ampl           = _fit.param[idx_ampl]
    _phase          = _fit.param[idx_phase]
    _offset         = _fit.param[idx_offset]
    # ---
    _results = _MySineFourParStruct(_frequ, _ampl, _phase, _offset)
    return _results, _fit 
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     = PlotlyJS.scatter(; x = _time_vec[_range], y = _signal_clean[_range],     name = "clean")
    line_disturbed = PlotlyJS.scatter(; x = _time_vec[_range], y = _signal_disturbed[_range], name = "disturbed")
    line_fitted    = PlotlyJS.scatter(; x = _time_vec[_range], y = _signal_fitted[_range],    name = "fitted")
    _data = [line_clean, line_disturbed, line_fitted]
    _layout = PlotlyJS.Layout(
        title_text       = "Clean and Disturbed Signal",
        xaxis_title_text = "Time / s",
        yaxis_title_text = "Time Domaine Signal / -"
    )
    return PlotlyJS.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_ULS = _LSQfit_four_param(_SineFourPar, signal_disturbed, sampling_rate)

fitted_signal       = _build_fited_signal(fit_result, t_vec)

PlotlyJS.update!(plot_clean_disturbed_and_fited_signal(signal_clean, signal_disturbed, fitted_signal, sampling_rate, frequ_A))

quality_normed      = StatsBase.rmsd(signal_disturbed, fitted_signal;    normalize = true)
quality             = StatsBase.rmsd(signal_disturbed, fitted_signal;    normalize = false)

# --- calculate good weighing vector:
weight_vec        = 1 ./ fit_ULS.resid

fit_result_WLS, fit_WLS = _LSQfit_four_param(_SineFourPar, signal_disturbed, sampling_rate, weight_vec)
fitted_signal_WLS       = _build_fited_signal(fit_result_WLS, t_vec)
quality_normed_WLS      = StatsBase.rmsd(signal_disturbed, fitted_signal_WLS;    normalize = true)
quality_WLS             = StatsBase.rmsd(signal_disturbed, fitted_signal_WLS;    normalize = false)

@info(string("Q:     \t", quality, ",\t Q_normed: \t", quality_normed))
@info(string("Q_WLS: \t", quality_WLS, ",\t Q_WLSnmd: \t", quality_normed_WLS))

[KronosTheLate]

This is a relativly long and complicated script - honestly, it is longer than I can bother debugging in my free time. Please try to cut away as much as possible, and provide a minimal working example (MWE) demonstrating the error. Also please include the error.

[StefanPofahl]

Thanks for the fast response!
To come to the point, the following thing does not work:

fit_OLS = curve_fit(m, tdata, ydata, p0)
wt = 1 ./ fit_OLS.resid
fit_WLS = curve_fit(m, tdata, ydata, wt, p0)

In my script this are the code lines:

fit_result, fit_ULS = _LSQfit_four_param(_SineFourPar, signal_disturbed, sampling_rate)
[...]
weight_vec        = 1 ./ fit_ULS.resid
[...]
_fit = LsqFit.curve_fit(model, _t_vec, _data_pts, _weighing_vec, _p0, lower = _lower_bounds, upper = _upper_bounds)

The full error message is:

ERROR: DomainError with -14.755570448982017:
sqrt will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).
Stacktrace:
  [1] throw_complex_domainerror(f::Symbol, x::Float64)
    @ Base.Math ./math.jl:33
  [2] sqrt
    @ ./math.jl:582 [inlined]
  [3] _broadcast_getindex_evalf
    @ ./broadcast.jl:648 [inlined]
  [4] _broadcast_getindex
    @ ./broadcast.jl:621 [inlined]
  [5] getindex
    @ ./broadcast.jl:575 [inlined]
  [6] macro expansion
    @ ./broadcast.jl:984 [inlined]
  [7] macro expansion
    @ ./simdloop.jl:77 [inlined]
  [8] copyto!
    @ ./broadcast.jl:983 [inlined]
  [9] copyto!
    @ ./broadcast.jl:936 [inlined]
 [10] copy
    @ ./broadcast.jl:908 [inlined]
 [11] materialize
    @ ./broadcast.jl:883 [inlined]
 [12] curve_fit(model::var"#model#5", xdata::Vector{Float64}, ydata::Vector{Float64}, wt::Vector{Float64}, p0::Vector{Float64}; inplace::Bool, kwargs::Base.Iterators.Pairs{Symbol, Vector{Float64}, Tuple{Symbol, Symbol}, NamedTuple{(:lower, :upper), Tuple{Vector{Float64}, Vector{Float64}}}})
    @ LsqFit ~/.julia/packages/LsqFit/BBrNp/src/curve_fit.jl:143
 [13] _LSQfit_four_param(_SineFourPar::_MySineFourParStruct, _data_pts::Vector{Float64}, _sampl_rate::Int64, _weighing_vec::Vector{Float64})
    @ Main /media/stefan/DATA/data/julia/tmp/test.jl:67
 [14] top-level scope
    @ /media/stefan/DATA/data/julia/tmp/test.jl:130

[KronosTheLate]

I am sorry, I am quite pressed with school at the moment and can not find time to sit down with this. I would suggest making a discourse post, following the guidelines in https://discourse.julialang.org/t/please-read-make-it-easier-to-help-you/14757. Good luck!

[StefanPofahl]

Hi,
thank you for taking the time to replay. Here is my proposal

wt = collect(Diagonal(1 ./ fit_OLS.resid.^2))

instead of a vector:

wt = 1 ./ fit_OLS.resid

Unfortunately, this suggestion to derive a suitable weighing matrix change the result not significantly
and it does not improve the result, this is the situation at least for a simple sinusoidal 4-parameter fit
(as can be found in my example script).

Regards,

Stefan