[🐞] DomainError when use background

[GeoTuxMan]

Describe the bug
I try to make a background analysis inspired by this notebook:
https://github.com/gher-uliege/Diva-Workshops/blob/master/notebooks/6-UnderConstruction/2do_background_notebook.ipynb
After that, I try to make a new analysis with previous background, but I got an error:
ERROR: LoadError: DomainError with -0.12812551856040955

Environment

• Julia version: 1.8.0
• Operating system: Windows 10

Installed modules

Type here the result of the command `import Pkg; Pkg.status()`
⌃ [336ed68f] CSV v0.10.4
⌃ [34da2185] Compat v4.2.0
  [efc8151c] DIVAnd v2.7.10 `https://github.com/gher-ulg/DIVAnd.jl.git#master`
⌃ [a93c6f00] DataFrames v1.3.6
  [864edb3b] DataStructures v0.18.13
  [c27321d9] Glob v1.3.0
⌃ [7073ff75] IJulia v1.23.3
⌃ [85f8d34a] NCDatasets v0.12.8
⌃ [3725be50] PhysOcean v0.6.6
  [d330b81b] PyPlot v2.11.0
⌃ [295af30f] Revise v3.4.0
  [ade2ca70] Dates
  [8bb1440f] DelimitedFiles
  [56ddb016] Logging
  [de0858da] Printf
  [9a3f8284] Random
  [10745b16] Statistics

To Reproduce
The command and their arguments which produced the error
include("test_background.jl")

Full screen output
Preferably obtained by setting ENV["JULIA_DEBUG"] = "DIVAnd"
julia> include("test_background.jl")
3.880550 seconds (11.86 M allocations: 1.397 GiB, 10.63% gc time, 80.00% compilation time)
extrema(b) = (NaN, NaN)
size(mask) = (388, 201, 13)
[ Info: Output dir already exists
[ Info: Starting computations
[ Info: grid size
[ Info: lon: 388
[ Info: lat: 201
[ Info: depth: 13
[ Info: N seasons: 1
[ Info: will compute 1 climatologies
yearlist = UnitRange{Int64}[1970:1986]
[ Info: year boundaries : [1970, 1986])
[ Info: Start reading NetCDF file
[ Info: 2023-02-03T09:57:16.401
23430 out of 63027 - 37.1745442429435 %
48510 out of 63027 - 76.96701413679851 %
[ Info: 2023-02-03T09:57:25.194
[ Info: End reading NetCDF file
31.637252 seconds (19.58 M allocations: 2.759 GiB, 1.37% gc time, 5.52% compilation time)
[ Info: Number of possible duplicates: 131
[ Info: Percentage of duplicates: 0.01%
[ Info: Checking ranges for dimensions and observations
minimum and maximum of obs. dimension 1: (26.382999420166016, 41.715999603271484)
minimum and maximum of obs. dimension 2: (40.16600036621094, 47.20000076293945)
minimum and maximum of obs. dimension 3: (0.0, 1682.0)
minimum and maximum of obs. dimension 4: (DateTime("1970-03-04T00:00:00"), DateTime("1986-12-27T05:00:00"))
minimum and maximum of data: (0.4466079771518707, 738.2000122070312)
[ Info: Computing weights using 1 CPU thread(s)
extrema(diagR) = (0.008203047821484443, 1.0)
lenz = [10, 10, 20, 20, 20, 20, 35, 50, 50, 50, 100, 100, 100]
lenxy = [33300.0, 44400.0, 55500.0, 66600.0, 77700.0, 88800.0, 99900.0, 111000.0, 122100.0, 133200.0, 144300.0, 149850.0, 155400.0]
[ Info: starting DIVA computations for Summer
[ Info: 2023-02-03T09:58:02.076
[ Info: Creating netCDF file Water_body_dissolved_oxygen_concentration_AllBS_background_3Feb2023_OneClim.nc
[ Info: Time step 1 / 1
[ Info: scaled correlation length (min,max) in dimension 1: (33300.0, 155400.0)
[ Info: scaled correlation length (min,max) in dimension 2: (33300.0, 155400.0)
[ Info: scaled correlation length (min,max) in dimension 3: (10.0, 100.0)
[ Info: number of windows: 1
220.924479 seconds (192.26 M allocations: 65.486 GiB, 7.06% gc time, 50.12% compilation time: 0% of which was recompilation)
[ Info: will compute 17 climatologies
[ Info: Will write results in C:/Users/Vostro/Documents/my_project_dir_2023/2023/Water_body_dissolved_oxygen_concentration_AllBS_Summer_withBackgr_3Feb2023.nc
[ Info: Creating netCDF file C:/Users/Vostro/Documents/my_project_dir_2023/2023/Water_body_dissolved_oxygen_concentration_AllBS_Summer_withBackgr_3Feb2023.nc
[ Info: Time step 1 / 17
[ Info: analysis time index 1 uses the background time index 1
ERROR: LoadError: DomainError with -0.12812551856040955:
log will only return a complex result if called with a complex argument. Try log(Complex(x)).
Stacktrace:
[1] throw_complex_domainerror(f::Symbol, x::Float64)
@ Base.Math .\math.jl:33
[2] _log(x::Float64, base::Val{:ℯ}, func::Symbol)
@ Base.Math .\special\log.jl:301
[3] log
@ .\special\log.jl:267 [inlined]
[4] trans
@ C:\Users\Vostro.julia\packages\DIVAnd\SijZ2\src\anamorphosis.jl:32 [inlined]
[5] trans
@ C:\Users\Vostro.julia\packages\DIVAnd\SijZ2\src\anamorphosis.jl:30 [inlined]
[6] _broadcast_getindex_evalf
@ .\broadcast.jl:670 [inlined]
[7] _broadcast_getindex
@ .\broadcast.jl:643 [inlined]
[8] getindex
@ .\broadcast.jl:597 [inlined]
[9] macro expansion
@ .\broadcast.jl:961 [inlined]
[10] macro expansion
@ .\simdloop.jl:77 [inlined]
[11] copyto!
@ .\broadcast.jl:960 [inlined]
[12] copyto!
@ .\broadcast.jl:913 [inlined]
[13] materialize!
@ .\broadcast.jl:871 [inlined]
[14] materialize!(dest::Array{Float64, 3}, bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{3}, Nothing, DIVAnd.Anam.var"#trans#10"{DIVAnd.Anam.var"#trans#8#11"{Float64, Float64}}, Tuple{Array{Float64, 3}}})
@ Base.Broadcast .\broadcast.jl:868
[15] (::DIVAnd.var"#615#618"{DIVAnd.var"#615#616#619"{TimeSelectorYearListMonthList{Vector{UnitRange{Int64}}, Vector{Vector{Int64}}}, Bool, NCDatasets.CFVariable{Union{Missing, Float32}, 4, NCDatasets.Variable{Float32, 4, NCDataset{Nothing}}, NCDatasets.Attributes{NCDataset{Nothing}}, NamedTuple{(:fillvalue, :missing_values, :scale_factor, :add_offset, :calendar, :time_origin, :time_factor), Tuple{Float32, Tuple{Float32}, Nothing, Nothing, Nothing, Nothing, Nothing}}}}})(xi::Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}, n::Int64, value::Vector{Float64}, trans::Function; selection::BitVector, obstime::Vector{DateTime})
@ DIVAnd C:\Users\Vostro.julia\packages\DIVAnd\SijZ2\src\utils.jl:775
[16] (::DIVAnd.var"#447#463"{Vector{Float64}, OrderedDict{String, String}, OrderedDict{String, String}, typeof(DIVAnd.distfun_m), DIVAnd.var"#615#618"{DIVAnd.var"#615#616#619"{TimeSelectorYearListMonthList{Vector{UnitRange{Int64}}, Vector{Vector{Int64}}}, Bool, NCDatasets.CFVariable{Union{Missing, Float32}, 4, NCDatasets.Variable{Float32, 4, NCDataset{Nothing}}, NCDatasets.Attributes{NCDataset{Nothing}}, NamedTuple{(:fillvalue, :missing_values, :scale_factor, :add_offset, :calendar, :time_origin, :time_factor), Tuple{Float32, Tuple{Float32}, Nothing, Nothing, Nothing, Nothing, Nothing}}}}}, Int64, Bool, Bool, Dict{Any, Any}, Dict{Any, Any}, Int64, Float64, Int64, Float64, Float64, Tuple{}, Bool, Vector{Tuple{String, Float64}}, typeof(DIVAndgo), Base.Pairs{Symbol, Symbol, Tuple{Symbol}, NamedTuple{(:solver,), Tuple{Symbol}}}, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{DateTime}}, Vector{Float64}, String, String, Bool, Matrix{DateTime}, Dict{Symbol, Any}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Vector{Float64}, Vector{Float64}, Int64, Tuple{Int64, Int64, Int64}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, DIVAnd.Anam.var"#invtrans#12"{DIVAnd.Anam.var"#invtrans#9#13"{Float64, Float64}}, DIVAnd.Anam.var"#trans#10"{DIVAnd.Anam.var"#trans#8#11"{Float64, Float64}}, Vector{DateTime}, TimeSelectorYearListMonthList{Vector{UnitRange{Int64}}, Vector{Vector{Int64}}}, Vector{Float64}, Int64})(ds::NCDataset{Nothing})
@ DIVAnd C:\Users\Vostro.julia\packages\DIVAnd\SijZ2\src\diva.jl:527
[17] NCDataset(::DIVAnd.var"#447#463"{Vector{Float64}, OrderedDict{String, String}, OrderedDict{String, String}, typeof(DIVAnd.distfun_m), DIVAnd.var"#615#618"{DIVAnd.var"#615#616#619"{TimeSelectorYearListMonthList{Vector{UnitRange{Int64}}, Vector{Vector{Int64}}}, Bool, NCDatasets.CFVariable{Union{Missing, Float32}, 4, NCDatasets.Variable{Float32, 4, NCDataset{Nothing}}, NCDatasets.Attributes{NCDataset{Nothing}}, NamedTuple{(:fillvalue, :missing_values, :scale_factor, :add_offset, :calendar, :time_origin, :time_factor), Tuple{Float32, Tuple{Float32}, Nothing, Nothing, Nothing, Nothing, Nothing}}}}}, Int64, Bool, Bool, Dict{Any, Any}, Dict{Any, Any}, Int64, Float64, Int64, Float64, Float64, Tuple{}, Bool, Vector{Tuple{String, Float64}}, typeof(DIVAndgo), Base.Pairs{Symbol, Symbol, Tuple{Symbol}, NamedTuple{(:solver,), Tuple{Symbol}}}, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{DateTime}}, Vector{Float64}, String, String, Bool, Matrix{DateTime}, Dict{Symbol, Any}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Vector{Float64}, Vector{Float64}, Int64, Tuple{Int64, Int64, Int64}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, DIVAnd.Anam.var"#invtrans#12"{DIVAnd.Anam.var"#invtrans#9#13"{Float64, Float64}}, DIVAnd.Anam.var"#trans#10"{DIVAnd.Anam.var"#trans#8#11"{Float64, Float64}}, Vector{DateTime}, TimeSelectorYearListMonthList{Vector{UnitRange{Int64}}, Vector{Vector{Int64}}}, Vector{Float64}, Int64}, ::String, ::Vararg{String}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ NCDatasets C:\Users\Vostro.julia\packages\NCDatasets\ipGBH\src\dataset.jl:241
[18] NCDataset(::Function, ::String, ::Vararg{String})
@ NCDatasets C:\Users\Vostro.julia\packages\NCDatasets\ipGBH\src\dataset.jl:238
[19] diva3d(xi::Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Vector{Float64}, TimeSelectorYearListMonthList{Vector{UnitRange{Int64}}, Vector{Vector{Int64}}}}, x::Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{DateTime}}, value::Vector{Float64}, len::Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Int64, 3}}, epsilon2::Vector{Float64}, filename::String, varname::String; datadir::String, bathname::String, bathisglobal::Bool, plotres::typeof(plotres2), timeorigin::DateTime, moddim::Vector{Float64}, zlevel::Symbol, ncvarattrib::OrderedDict{String, String}, ncglobalattrib::OrderedDict{String, String}, transform::Tuple{DIVAnd.Anam.var"#trans#10"{DIVAnd.Anam.var"#trans#8#11"{Float64, Float64}}, DIVAnd.Anam.var"#invtrans#12"{DIVAnd.Anam.var"#invtrans#9#13"{Float64, Float64}}}, distfun::typeof(DIVAnd.distfun_m), mask::BitArray{3}, background::DIVAnd.var"#615#618"{DIVAnd.var"#615#616#619"{TimeSelectorYearListMonthList{Vector{UnitRange{Int64}}, Vector{Vector{Int64}}}, Bool, NCDatasets.CFVariable{Union{Missing, Float32}, 4, NCDatasets.Variable{Float32, 4, NCDataset{Nothing}}, NCDatasets.Attributes{NCDataset{Nothing}}, NamedTuple{(:fillvalue, :missing_values, :scale_factor, :add_offset, :calendar, :time_origin, :time_factor), Tuple{Float32, Tuple{Float32}, Nothing, Nothing, Nothing, Nothing, Nothing}}}}}, background_epsilon2_factor::Nothing, background_lenz::Nothing, background_len::Nothing, background_lenz_factor::Int64, filterbackground::Int64, fitcorrlen::Bool, fithorzcorrlen::Bool, fitvertcorrlen::Bool, fithorz_param::Dict{Any, Any}, fitvert_param::Dict{Any, Any}, memtofit::Int64, overlapfactor::Float64, niter_e::Int64, minfield::Float64, maxfield::Float64, surfextend::Bool, velocity::Tuple{}, stat_per_timeslice::Bool, error_thresholds::Vector{Tuple{String, Float64}}, divamethod::typeof(DIVAndgo), kwargs::Base.Pairs{Symbol, Symbol, Tuple{Symbol}, NamedTuple{(:solver,), Tuple{Symbol}}})
@ DIVAnd C:\Users\Vostro.julia\packages\DIVAnd\SijZ2\src\diva.jl:375
[20] top-level scope
@ C:\Users\Vostro\Documents\my_project_dir_2023\test_background.jl:320
[21] include(fname::String)
@ Base.MainInclude .\client.jl:476
[22] top-level scope
@ REPL[6]:1
in expression starting at C:\Users\Vostro\Documents\my_project_dir_2023\test_background.jl:320
background.zip

** Full stack trace with error message**

[ctroupin]

Thanks for the full bug report, it really helps.

I think it comes again from negative values as argument of a logarithm functions.
Maybe the negative values are not in the original data but in the anomalies once the background is subtracted.

If that's the case we need to ensure that no negative anomalies are present...

[jmbeckers]

If you work with anomalies, I think you cannot apply the log transform because by definitions, anomalies can be negative.

If you still think the distribution is lognormal, what you should do is
work the whole process using log(data) as input and at the very end (once you combined your background and anomaly analysis) take the inverse transformation (but then you deactivate the transformation in DIVAnd during the whole process).

[GeoTuxMan]

Yes, without the line transform = Anam.loglin(slog), the code run fine with any errors
Thank you very much! 👍

[ctroupin]

Thanks @GeoTuxMan !
I've added a little explanation to the 2do_background_notebook.ipynb notebook.

Let's close the issue.