search

gher-uliege / DIVAnd.jl

DIVAnd performs an n-dimensional variational analysis of arbitrarily located observations
GNU General Public License v2.0
70 stars 11 forks source link

[🐞] DomainError when use background #128

Closed GeoTuxMan closed 1 year ago

GeoTuxMan commented 1 year ago

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

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) 10, 10, 20, 20, 20, 20, 35, 50, 50, 50, 100, 100, 100 = (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 commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

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

ctroupin commented 1 year ago

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

Let's close the issue.