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

MethodError when trying to run diva3d with fitcorrlen = true #109

Closed qcazenave closed 1 year ago

qcazenave commented 1 year ago

Hello, I am running DIVAnd v2.7.9 on JupyterLab with Julia 1.8.2. I did a first run of the function diva3d with predefined correlation length and it worked.

I then tried to run the function setting len=() and fitcorrlen = true and I get the following error message (see below). I checked in the diva.jl file, line 485 : "if all(background_len[3] .== 0)" So it seems "background_len" is not defined... What did I do wrong ?

Thank you

MethodError: no method matching getindex(::Nothing, ::Int64)

Stacktrace: [1] (::DIVAnd.var"#447#463"{Vector{Float64}, OrderedDict{String, String}, OrderedDict{Symbol, String}, typeof(DIVAnd.distfun_m), Nothing, Int64, Bool, Bool, Dict{Symbol, Function}, Dict{Symbol, Function}, Int64, Float64, Int64, Float64, Float64, Tuple{}, Bool, Vector{Tuple{String, Float64}}, typeof(DIVAndgo), Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, 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#5"{DIVAnd.Anam.var"#invtrans#2#6"}, DIVAnd.Anam.var"#trans#3"{DIVAnd.Anam.var"#trans#1#4"}, Vector{DateTime}, TimeSelectorYearListMonthList{Vector{UnitRange{Int64}}, Vector{Vector{Int64}}}, Vector{Float64}, Int64})(ds::NCDataset{Nothing}) @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\diva.jl:485 [2] NCDataset(::DIVAnd.var"#447#463"{Vector{Float64}, OrderedDict{String, String}, OrderedDict{Symbol, String}, typeof(DIVAnd.distfun_m), Nothing, Int64, Bool, Bool, Dict{Symbol, Function}, Dict{Symbol, Function}, Int64, Float64, Int64, Float64, Float64, Tuple{}, Bool, Vector{Tuple{String, Float64}}, typeof(DIVAndgo), Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, 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#5"{DIVAnd.Anam.var"#invtrans#2#6"}, DIVAnd.Anam.var"#trans#3"{DIVAnd.Anam.var"#trans#1#4"}, 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\qcazenav.julia\packages\NCDatasets\h1epE\src\dataset.jl:241 [3] NCDataset(::Function, ::String, ::Vararg{String}) @ NCDatasets C:\Users\qcazenav.julia\packages\NCDatasets\h1epE\src\dataset.jl:238 [4] 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{}, epsilon2::Vector{Float64}, filename::String, varname::String; datadir::String, bathname::String, bathisglobal::Bool, plotres::typeof(plotres), timeorigin::DateTime, moddim::Vector{Float64}, zlevel::Symbol, ncvarattrib::OrderedDict{String, String}, ncglobalattrib::OrderedDict{Symbol, String}, transform::Tuple{DIVAnd.Anam.var"#trans#3"{DIVAnd.Anam.var"#trans#1#4"}, DIVAnd.Anam.var"#invtrans#5"{DIVAnd.Anam.var"#invtrans#2#6"}}, distfun::typeof(DIVAnd.distfun_m), mask::BitArray{3}, background::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{Symbol, Function}, fitvert_param::Dict{Symbol, Function}, 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, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}) @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\diva.jl:375 [5] top-level scope @ .\timing.jl:262 [inlined] [6] top-level scope @ .\In[49]:0 [7] eval @ .\boot.jl:368 [inlined] [8] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String) @ Base .\loading.jl:1428

ctroupin commented 1 year ago

Hello, I cannot really tell for the moment, maybe if you could create a minimal example that fails (small dataset, region etc) so we can try the same code?

qcazenave commented 1 year ago

Test.zip

Here is a zip file containing input data and the Jupyter notebook. I did not add the gebco bathymetry file because it is too big and I assume you can easily get it.

ctroupin commented 1 year ago

thanks! Meanwhile we discussed with @jmbeckers about that issue and maybe you try to run diva3d not with len=() but with a real (i.e. not empty) tuple, with the expected dimension, to see if the fit is performed. For instance 

len = (5., 5., 200.)
qcazenave commented 1 year ago

So I tried with len = (5.0,5.0,200.0) and fitcorrlen=true and without changing any of the default fit parameters (fithorz_param and fitvert_param are not specified, contrary to what is written in the version I sent you, which I believe is invalid regarding the dictionaries I tried to define to specify fithorz_param and fitvert_param). It started working (it seems I have a horizontal correlation length for every depth level) but it still ended with another error message.

the way I call the function :

@time info = diva3d( (lonr,latr,depthr,TS), (obslon,obslat,obsdepth,obstime), log.(obsval), len, epsilon2, outputfile, varname, bathname = bathfile, plotres = plotres, mask = mask, fitcorrlen = fitcorrlen, ncvarattrib = ncvarattrib, ncglobalattrib = ncglobalattrib, surfextend = true );

the error message :

┌ Info: Creating netCDF file C:\Users\qcazenav\Documents\EMODnet_Chem_data\seasonal_maps\fitcorrlen_test_Water_body_dissolved_inorganicnitrogen(DIN)_maps_19712017_seasonal.4Danl.nc └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\diva.jl:374 ┌ Info: Time step 1 / 4 └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\diva.jl:427 ┌ Warning: Observations equal to NaN: 6 └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\DIVAnd_obs.jl:45 ┌ Info: Data points at z=-10.0: 7287, horz. correlation length: 70725.06209640756 (preliminary) └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1085 ┌ Info: Data points at z=0.0: 7585, horz. correlation length: 69350.53409283079 (preliminary) └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1085 ┌ Info: Data points at z=10.0: 7735, horz. correlation length: 83061.2197016761 (preliminary) └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1085 ┌ Info: Data points at z=20.0: 7791, horz. correlation length: 82228.08554316644 (preliminary) └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1085 ┌ Info: Data points at z=50.0: 7636, horz. correlation length: 69908.74813321588 (preliminary) └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1085 ┌ Info: Data points at z=100.0: 743, horz. correlation length: 183669.63952495012 (preliminary) └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1085 ┌ Info: Data points at z=150.0: 237, horz. correlation length: 89325.88839254994 (preliminary) └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1085 ┌ Info: Data points at z=200.0: 96, horz. correlation length: 22426.812735955686 (preliminary) └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1085 ┌ Warning: Too few data. Will use guesses (np = 9, RLz = 128333.53468397829, ) └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:900 ┌ Info: Data points at z=300.0: 46, horz. correlation length: 128333.53468397829 (preliminary) └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1085 ┌ Warning: Too few data. Will use guesses (np = 6, RLz = 76909.29767654871, ) └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:900 ┌ Info: Data points at z=400.0: 42, horz. correlation length: 76909.29767654871 (preliminary) └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1085 ┌ Warning: Too few data. Will use guesses (np = 4, RLz = 23570.275283835213, ) └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:900 ┌ Info: Data points at z=500.0: 95, horz. correlation length: 23570.275283835213 (preliminary) └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1085 ┌ Info: Data points at z=1000.0: 60, horz. correlation length: 217422.05576764955 (preliminary) └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1085 ┌ Warning: Too few data. Will use guesses (np = 5, RLz = 97211.20866295644, ) └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:900 ┌ Info: Data points at z=1500.0: 35, horz. correlation length: 97211.20866295644 (preliminary) └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1085 ┌ Warning: Too few data. Will use guesses (np = 5, RLz = 59459.245917581575, ) └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:900 ┌ Info: Data points at z=2000.0: 40, horz. correlation length: 59459.245917581575 (preliminary) └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1085 ┌ Warning: Too few data. Will use guesses (np = 4, RLz = 217614.64714484406, ) └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:900 ┌ Info: Data points at z=2500.0: 19, horz. correlation length: 217614.64714484406 (preliminary) └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1085 ┌ Warning: Too few data. Will use guesses (np = 1, RLz = 199574.50146485653, ) └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:900 ┌ Info: Data points at z=3000.0: 11, horz. correlation length: 199574.50146485653 (preliminary) └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1085 ┌ Info: Smoothed horz. correlation length at z=-10.0: 92855.11999059113 └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1113 ┌ Info: Smoothed horz. correlation length at z=0.0: 92858.4666022335 └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1113 ┌ Info: Smoothed horz. correlation length at z=10.0: 92858.18474983455 └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1113 ┌ Info: Smoothed horz. correlation length at z=20.0: 92840.75487077267 └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1113 ┌ Info: Smoothed horz. correlation length at z=50.0: 92567.81103659586 └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1113 ┌ Info: Smoothed horz. correlation length at z=100.0: 89596.57355187176 └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1113 ┌ Info: Smoothed horz. correlation length at z=150.0: 80170.55959796978 └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1113 ┌ Info: Smoothed horz. correlation length at z=200.0: 67305.52758974109 └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1113 ┌ Info: Smoothed horz. correlation length at z=300.0: 52109.20351135579 └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1113 ┌ Info: Smoothed horz. correlation length at z=400.0: 44282.90028713939 └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1113 ┌ Info: Smoothed horz. correlation length at z=500.0: 38921.652257863745 └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1113 ┌ Info: Smoothed horz. correlation length at z=1000.0: 35671.36953356311 └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1113 ┌ Info: Smoothed horz. correlation length at z=1500.0: 34470.49268077742 └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1113 ┌ Info: Smoothed horz. correlation length at z=2000.0: 57430.1790404207 └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1113 ┌ Info: Smoothed horz. correlation length at z=2500.0: 214061.43122342948 └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1113 ┌ Info: Smoothed horz. correlation length at z=3000.0: 199898.77654016748 └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\fit.jl:1113

MethodError: no method matching setindex!(::Float64, ::Float64, ::Int64, ::Int64, ::Int64)

Stacktrace: [1] (::DIVAnd.var"#447#463"{Vector{Float64}, OrderedDict{String, String}, OrderedDict{Symbol, String}, typeof(DIVAnd.distfun_m), 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, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{DateTime}}, Vector{Float64}, String, String, Bool, Matrix{DateTime}, Dict{Symbol, Any}, Tuple{Float64, Float64, Float64}, Tuple{Float64, Float64, Float64}, 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#5"{DIVAnd.Anam.var"#invtrans#2#6"}, DIVAnd.Anam.var"#trans#3"{DIVAnd.Anam.var"#trans#1#4"}, Vector{DateTime}, TimeSelectorYearListMonthList{Vector{UnitRange{Int64}}, Vector{Vector{Int64}}}, Vector{Float64}, Int64})(ds::NCDataset{Nothing}) @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\diva.jl:576 [2] NCDataset(::DIVAnd.var"#447#463"{Vector{Float64}, OrderedDict{String, String}, OrderedDict{Symbol, String}, typeof(DIVAnd.distfun_m), 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, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{DateTime}}, Vector{Float64}, String, String, Bool, Matrix{DateTime}, Dict{Symbol, Any}, Tuple{Float64, Float64, Float64}, Tuple{Float64, Float64, Float64}, 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#5"{DIVAnd.Anam.var"#invtrans#2#6"}, DIVAnd.Anam.var"#trans#3"{DIVAnd.Anam.var"#trans#1#4"}, 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\qcazenav.julia\packages\NCDatasets\h1epE\src\dataset.jl:241 [3] NCDataset(::Function, ::String, ::Vararg{String}) @ NCDatasets C:\Users\qcazenav.julia\packages\NCDatasets\h1epE\src\dataset.jl:238 [4] 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{Float64, Float64, Float64}, epsilon2::Vector{Float64}, filename::String, varname::String; datadir::String, bathname::String, bathisglobal::Bool, plotres::typeof(plotres), timeorigin::DateTime, moddim::Vector{Float64}, zlevel::Symbol, ncvarattrib::OrderedDict{String, String}, ncglobalattrib::OrderedDict{Symbol, String}, transform::Tuple{DIVAnd.Anam.var"#trans#3"{DIVAnd.Anam.var"#trans#1#4"}, DIVAnd.Anam.var"#invtrans#5"{DIVAnd.Anam.var"#invtrans#2#6"}}, distfun::typeof(DIVAnd.distfun_m), mask::BitArray{3}, background::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, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}) @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\diva.jl:375 [5] top-level scope @ .\timing.jl:262 [inlined] [6] top-level scope @ .\In[40]:0 [7] eval @ .\boot.jl:368 [inlined] [8] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String) @ Base .\loading.jl:1428

ctroupin commented 1 year ago

thanks, I'm trying now on my machine, a priori not easy to identify the issue, at least the first problem seems to be solved, we would need to change something in the code so that len=() doesn't generate that problem.

ctroupin commented 1 year ago

OK got the same error message as you -- always a good start --, now I've tested by setting 

fitcorrlen = false
len = (2., 2., 10.)

and I got it working without error. So it would be an issue with fitcorrlen.

Edit: Something that works is to take the correlation length you've created in the cell starting with 

sz = (length(lonr),length(latr),length(depthr));

lenh = max(10*dx*100_000.,300_000.)
...

then run the diva3d with 

fitcorrlen = true

and commenting the lines:

#fithorz_param = fithorz_param, 
#fitvert_param = fitvert_param,

Doing so, len is defined as a 3-element tuple, each of them having the same dimension as the grid (17 X 13 X 15) and it seems this allows the fitting to work properly.

Now I have to see why the fithorz_param's are creating an issue.

qcazenave commented 1 year ago

OK, thank you for your help ! I wonder if the fithorz_param is creating an issue because I defined limitfun as a function of (dx,dy)...

ctroupin commented 1 year ago

that's possible because I think the function doesn't know about dx and dy... Anyway it's easy to test/

qcazenave commented 1 year ago

Hi again ! So I tried the code again with: fitcorrlen = true len = (fill(300_000.0,sz),fill(300_000.0,sz),fill(200.0,sz)) and using diva3d as follows: `@time info = diva3d( (lonr,latr,depthr,TS), (obslon,obslat,obsdepth,obstime), obsval, len, epsilon2, outputfile, varname, bathname = bathfile, plotres = plotres, mask = mask, transform = Anam.loglin(200.0), fitcorrlen = fitcorrlen,

fithorz_param = fithorz_param,

#fitvert_param = fitvert_param,
niter_e = 1,
ncvarattrib = ncvarattrib,
ncglobalattrib = ncglobalattrib,
surfextend = true

);`

the fit works since I get the following message:

image

but then, a new error occurs:

┌ Info: scaled correlation length (min,max) in dimension 1: (3.673417765874348e10, 1.5490036994975006e11) └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\diva.jl:612 ┌ Info: scaled correlation length (min,max) in dimension 2: (3.673417765874348e10, 1.5490036994975006e11) └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\diva.jl:612 ┌ Info: scaled correlation length (min,max) in dimension 3: (37172.86865015531, 375733.049001446) └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\diva.jl:612 ┌ Info: number of windows: 1 └ @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\DIVAndgo.jl:110 Unhandled Task ERROR: PosDefException: matrix is not positive definite; Cholesky factorization failed. Stacktrace: [1] cholesky!(F::SuiteSparse.CHOLMOD.Factor{Float64}, A::SuiteSparse.CHOLMOD.Sparse{Float64}; shift::Float64, check::Bool) @ SuiteSparse.CHOLMOD C:\Users\qcazenav\Programmes\Julia-1.8.2\share\julia\stdlib\v1.8\SuiteSparse\src\cholmod.jl:1150 [2] #cholesky#8 @ C:\Users\qcazenav\Programmes\Julia-1.8.2\share\julia\stdlib\v1.8\SuiteSparse\src\cholmod.jl:1185 [inlined] [3] cholesky @ C:\Users\qcazenav\Programmes\Julia-1.8.2\share\julia\stdlib\v1.8\SuiteSparse\src\cholmod.jl:1178 [inlined] [4] #cholesky#9 @ C:\Users\qcazenav\Programmes\Julia-1.8.2\share\julia\stdlib\v1.8\SuiteSparse\src\cholmod.jl:1297 [inlined] [5] cholesky @ C:\Users\qcazenav\Programmes\Julia-1.8.2\share\julia\stdlib\v1.8\SuiteSparse\src\cholmod.jl:1297 [inlined] [6] factorize!(C::CovarIS{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}}) @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\special_matrices.jl:117 [7] DIVAnd_factorize!(s::DIVAnd.DIVAnd_struct{Float64, Int64, 3, SparseArrays.SparseMatrixCSC{Float64, Int64}}) @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\DIVAnd_factorize.jl:40 [8] DIVAndrun(operatortype::Type, mask::BitArray{3}, pmnin::Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, xiin::Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, x::Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}, f::Vector{Float64}, lin::Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, epsilon2::Vector{Float64}; velocity::Tuple{}, primal::Bool, factorize::Bool, tol::Float64, maxit::Int64, minit::Int64, constraints::Tuple{}, inversion::Symbol, moddim::Vector{Float64}, fracindex::Matrix{Float64}, alpha::Vector{Any}, keepLanczosVectors::Int64, compPC::typeof(DIVAnd_pc_none), progress::DIVAnd.var"#286#288", fi0::Array{Float64, 3}, f0::Vector{Float64}, alphabc::Float64, scale_len::Bool, btrunc::Vector{Any}, MEMTOFIT::Int64, topographyforfluxes::Tuple{}, fluxes::Tuple{}, epsfluxes::Int64, epsilon2forfractions::Int64, RTIMESONESCALES::Tuple{}, QCMETHOD::Tuple{}, coeff_laplacian::Vector{Float64}, coeff_derivative2::Vector{Float64}, mean_Labs::Vector{Float64}) @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\DIVAndrun.jl:125 [9] #DIVAndrun#290 @ C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\DIVAndrun.jl:298 [inlined] [10] macro expansion @ C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\DIVAndgo.jl:252 [inlined] [11] (::DIVAnd.var"#292#301"{3, Vector{Float64}, Int64, Tuple{}, Tuple{}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, BitArray{3}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}, Vector{Float64}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Vector{Float64}, Symbol, Vector{Float64}, SharedArrays.SharedVector{Float32}, SharedArrays.SharedVector{Float32}, SharedArrays.SharedVector{Float32}, SharedArrays.SharedArray{Float32, 3}, SharedArrays.SharedArray{Float32, 3}, Vector{Int64}, Vector{Float64}, Vector{Int64}, Vector{NTuple{6, Vector{Int64}}}, Bool, Bool})(R::UnitRange{Int64}, lo::Int64, hi::Int64) @ DIVAnd C:\Users\qcazenav\Programmes\Julia-1.8.2\share\julia\stdlib\v1.8\Distributed\src\macros.jl:303 [12] (::Distributed.var"#178#180"{UnitRange{Int64}, DIVAnd.var"#292#301"{3, Vector{Float64}, Int64, Tuple{}, Tuple{}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, BitArray{3}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}, Vector{Float64}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Vector{Float64}, Symbol, Vector{Float64}, SharedArrays.SharedVector{Float32}, SharedArrays.SharedVector{Float32}, SharedArrays.SharedVector{Float32}, SharedArrays.SharedArray{Float32, 3}, SharedArrays.SharedArray{Float32, 3}, Vector{Int64}, Vector{Float64}, Vector{Int64}, Vector{NTuple{6, Vector{Int64}}}, Bool, Bool}, UnitRange{Int64}})() @ Distributed C:\Users\qcazenav\Programmes\Julia-1.8.2\share\julia\stdlib\v1.8\Distributed\src\macros.jl:83 [13] #invokelatest#2 @ .\essentials.jl:729 [inlined] [14] invokelatest @ .\essentials.jl:726 [inlined] [15] #153 @ C:\Users\qcazenav\Programmes\Julia-1.8.2\share\julia\stdlib\v1.8\Distributed\src\remotecall.jl:425 [inlined] [16] run_work_thunk(thunk::Distributed.var"#153#154"{Distributed.var"#178#180"{UnitRange{Int64}, DIVAnd.var"#292#301"{3, Vector{Float64}, Int64, Tuple{}, Tuple{}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, BitArray{3}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}, Vector{Float64}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Vector{Float64}, Symbol, Vector{Float64}, SharedArrays.SharedVector{Float32}, SharedArrays.SharedVector{Float32}, SharedArrays.SharedVector{Float32}, SharedArrays.SharedArray{Float32, 3}, SharedArrays.SharedArray{Float32, 3}, Vector{Int64}, Vector{Float64}, Vector{Int64}, Vector{NTuple{6, Vector{Int64}}}, Bool, Bool}, UnitRange{Int64}}, Tuple{}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, print_error::Bool) @ Distributed C:\Users\qcazenav\Programmes\Julia-1.8.2\share\julia\stdlib\v1.8\Distributed\src\process_messages.jl:70 [17] run_work_thunk(rv::Distributed.RemoteValue, thunk::Function) @ Distributed C:\Users\qcazenav\Programmes\Julia-1.8.2\share\julia\stdlib\v1.8\Distributed\src\process_messages.jl:79 [18] (::Distributed.var"#100#102"{Distributed.RemoteValue, Distributed.var"#153#154"{Distributed.var"#178#180"{UnitRange{Int64}, DIVAnd.var"#292#301"{3, Vector{Float64}, Int64, Tuple{}, Tuple{}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, BitArray{3}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}, Vector{Float64}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Vector{Float64}, Symbol, Vector{Float64}, SharedArrays.SharedVector{Float32}, SharedArrays.SharedVector{Float32}, SharedArrays.SharedVector{Float32}, SharedArrays.SharedArray{Float32, 3}, SharedArrays.SharedArray{Float32, 3}, Vector{Int64}, Vector{Float64}, Vector{Int64}, Vector{NTuple{6, Vector{Int64}}}, Bool, Bool}, UnitRange{Int64}}, Tuple{}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}})() @ Distributed .\task.jl:484 Stacktrace: [1] sync_end(c::Channel{Any}) @ Base .\task.jl:436 [2] (::Distributed.var"#177#179"{DIVAnd.var"#292#301"{3, Vector{Float64}, Int64, Tuple{}, Tuple{}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, BitArray{3}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}, Vector{Float64}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Vector{Float64}, Symbol, Vector{Float64}, SharedArrays.SharedVector{Float32}, SharedArrays.SharedVector{Float32}, SharedArrays.SharedVector{Float32}, SharedArrays.SharedArray{Float32, 3}, SharedArrays.SharedArray{Float32, 3}, Vector{Int64}, Vector{Float64}, Vector{Int64}, Vector{NTuple{6, Vector{Int64}}}, Bool, Bool}, UnitRange{Int64}})() @ Distributed .\task.jl:455

TaskFailedException

nested task error: PosDefException: matrix is not positive definite; Cholesky factorization failed.
Stacktrace:
  [1] cholesky!(F::SuiteSparse.CHOLMOD.Factor{Float64}, A::SuiteSparse.CHOLMOD.Sparse{Float64}; shift::Float64, check::Bool)
    @ SuiteSparse.CHOLMOD C:\Users\qcazenav\Programmes\Julia-1.8.2\share\julia\stdlib\v1.8\SuiteSparse\src\cholmod.jl:1150
  [2] #cholesky#8
    @ C:\Users\qcazenav\Programmes\Julia-1.8.2\share\julia\stdlib\v1.8\SuiteSparse\src\cholmod.jl:1185 [inlined]
  [3] cholesky
    @ C:\Users\qcazenav\Programmes\Julia-1.8.2\share\julia\stdlib\v1.8\SuiteSparse\src\cholmod.jl:1178 [inlined]
  [4] #cholesky#9
    @ C:\Users\qcazenav\Programmes\Julia-1.8.2\share\julia\stdlib\v1.8\SuiteSparse\src\cholmod.jl:1297 [inlined]
  [5] cholesky
    @ C:\Users\qcazenav\Programmes\Julia-1.8.2\share\julia\stdlib\v1.8\SuiteSparse\src\cholmod.jl:1297 [inlined]
  [6] factorize!(C::CovarIS{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}})
    @ DIVAnd C:\Users\qcazenav\.julia\packages\DIVAnd\MV3j9\src\special_matrices.jl:117
  [7] DIVAnd_factorize!(s::DIVAnd.DIVAnd_struct{Float64, Int64, 3, SparseArrays.SparseMatrixCSC{Float64, Int64}})
    @ DIVAnd C:\Users\qcazenav\.julia\packages\DIVAnd\MV3j9\src\DIVAnd_factorize.jl:40
  [8] DIVAndrun(operatortype::Type, mask::BitArray{3}, pmnin::Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, xiin::Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, x::Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}, f::Vector{Float64}, lin::Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, epsilon2::Vector{Float64}; velocity::Tuple{}, primal::Bool, factorize::Bool, tol::Float64, maxit::Int64, minit::Int64, constraints::Tuple{}, inversion::Symbol, moddim::Vector{Float64}, fracindex::Matrix{Float64}, alpha::Vector{Any}, keepLanczosVectors::Int64, compPC::typeof(DIVAnd_pc_none), progress::DIVAnd.var"#286#288", fi0::Array{Float64, 3}, f0::Vector{Float64}, alphabc::Float64, scale_len::Bool, btrunc::Vector{Any}, MEMTOFIT::Int64, topographyforfluxes::Tuple{}, fluxes::Tuple{}, epsfluxes::Int64, epsilon2forfractions::Int64, RTIMESONESCALES::Tuple{}, QCMETHOD::Tuple{}, coeff_laplacian::Vector{Float64}, coeff_derivative2::Vector{Float64}, mean_Labs::Vector{Float64})
    @ DIVAnd C:\Users\qcazenav\.julia\packages\DIVAnd\MV3j9\src\DIVAndrun.jl:125
  [9] #DIVAndrun#290
    @ C:\Users\qcazenav\.julia\packages\DIVAnd\MV3j9\src\DIVAndrun.jl:298 [inlined]
 [10] macro expansion
    @ C:\Users\qcazenav\.julia\packages\DIVAnd\MV3j9\src\DIVAndgo.jl:252 [inlined]
 [11] (::DIVAnd.var"#292#301"{3, Vector{Float64}, Int64, Tuple{}, Tuple{}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, BitArray{3}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}, Vector{Float64}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Vector{Float64}, Symbol, Vector{Float64}, SharedArrays.SharedVector{Float32}, SharedArrays.SharedVector{Float32}, SharedArrays.SharedVector{Float32}, SharedArrays.SharedArray{Float32, 3}, SharedArrays.SharedArray{Float32, 3}, Vector{Int64}, Vector{Float64}, Vector{Int64}, Vector{NTuple{6, Vector{Int64}}}, Bool, Bool})(R::UnitRange{Int64}, lo::Int64, hi::Int64)
    @ DIVAnd C:\Users\qcazenav\Programmes\Julia-1.8.2\share\julia\stdlib\v1.8\Distributed\src\macros.jl:303
 [12] (::Distributed.var"#178#180"{UnitRange{Int64}, DIVAnd.var"#292#301"{3, Vector{Float64}, Int64, Tuple{}, Tuple{}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, BitArray{3}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}, Vector{Float64}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Vector{Float64}, Symbol, Vector{Float64}, SharedArrays.SharedVector{Float32}, SharedArrays.SharedVector{Float32}, SharedArrays.SharedVector{Float32}, SharedArrays.SharedArray{Float32, 3}, SharedArrays.SharedArray{Float32, 3}, Vector{Int64}, Vector{Float64}, Vector{Int64}, Vector{NTuple{6, Vector{Int64}}}, Bool, Bool}, UnitRange{Int64}})()
    @ Distributed C:\Users\qcazenav\Programmes\Julia-1.8.2\share\julia\stdlib\v1.8\Distributed\src\macros.jl:83
 [13] #invokelatest#2
    @ .\essentials.jl:729 [inlined]
 [14] invokelatest
    @ .\essentials.jl:726 [inlined]
 [15] #153
    @ C:\Users\qcazenav\Programmes\Julia-1.8.2\share\julia\stdlib\v1.8\Distributed\src\remotecall.jl:425 [inlined]
 [16] run_work_thunk(thunk::Distributed.var"#153#154"{Distributed.var"#178#180"{UnitRange{Int64}, DIVAnd.var"#292#301"{3, Vector{Float64}, Int64, Tuple{}, Tuple{}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, BitArray{3}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}, Vector{Float64}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Vector{Float64}, Symbol, Vector{Float64}, SharedArrays.SharedVector{Float32}, SharedArrays.SharedVector{Float32}, SharedArrays.SharedVector{Float32}, SharedArrays.SharedArray{Float32, 3}, SharedArrays.SharedArray{Float32, 3}, Vector{Int64}, Vector{Float64}, Vector{Int64}, Vector{NTuple{6, Vector{Int64}}}, Bool, Bool}, UnitRange{Int64}}, Tuple{}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, print_error::Bool)
    @ Distributed C:\Users\qcazenav\Programmes\Julia-1.8.2\share\julia\stdlib\v1.8\Distributed\src\process_messages.jl:70
 [17] run_work_thunk(rv::Distributed.RemoteValue, thunk::Function)
    @ Distributed C:\Users\qcazenav\Programmes\Julia-1.8.2\share\julia\stdlib\v1.8\Distributed\src\process_messages.jl:79
 [18] (::Distributed.var"#100#102"{Distributed.RemoteValue, Distributed.var"#153#154"{Distributed.var"#178#180"{UnitRange{Int64}, DIVAnd.var"#292#301"{3, Vector{Float64}, Int64, Tuple{}, Tuple{}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, BitArray{3}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}, Vector{Float64}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Vector{Float64}, Symbol, Vector{Float64}, SharedArrays.SharedVector{Float32}, SharedArrays.SharedVector{Float32}, SharedArrays.SharedVector{Float32}, SharedArrays.SharedArray{Float32, 3}, SharedArrays.SharedArray{Float32, 3}, Vector{Int64}, Vector{Float64}, Vector{Int64}, Vector{NTuple{6, Vector{Int64}}}, Bool, Bool}, UnitRange{Int64}}, Tuple{}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}})()
    @ Distributed .\task.jl:484
Stacktrace:
 [1] sync_end(c::Channel{Any})
   @ Base .\task.jl:436
 [2] (::Distributed.var"#177#179"{DIVAnd.var"#292#301"{3, Vector{Float64}, Int64, Tuple{}, Tuple{}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, BitArray{3}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}, Vector{Float64}, Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, Vector{Float64}, Symbol, Vector{Float64}, SharedArrays.SharedVector{Float32}, SharedArrays.SharedVector{Float32}, SharedArrays.SharedVector{Float32}, SharedArrays.SharedArray{Float32, 3}, SharedArrays.SharedArray{Float32, 3}, Vector{Int64}, Vector{Float64}, Vector{Int64}, Vector{NTuple{6, Vector{Int64}}}, Bool, Bool}, UnitRange{Int64}})()
   @ Distributed .\task.jl:455

Stacktrace: [1] sync_end(c::Channel{Any}) @ Base .\task.jl:436 [2] macro expansion @ .\task.jl:455 [inlined] [3] DIVAndgo(mask::BitArray{3}, pmn::Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, xi::Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, x::Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}, f::Vector{Float64}, Labs::Tuple{Array{Float64, 3}, Array{Float64, 3}, Array{Float64, 3}}, epsilon2::Vector{Float64}, errormethod::Symbol; moddim::Vector{Float64}, velocity::Tuple{}, MEMTOFIT::Int64, QCMETHOD::Tuple{}, RTIMESONESCALES::Tuple{}, solver::Symbol, overlapfactor::Float64, filteranom::Int64, filtererr::Int64, otherargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}) @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\DIVAndgo.jl:115 [4] (::DIVAnd.var"#447#463"{Vector{Float64}, OrderedDict{String, String}, OrderedDict{Symbol, String}, typeof(DIVAnd.distfun_m), 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, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, 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\qcazenav.julia\packages\DIVAnd\MV3j9\src\diva.jl:644 [5] NCDataset(::DIVAnd.var"#447#463"{Vector{Float64}, OrderedDict{String, String}, OrderedDict{Symbol, String}, typeof(DIVAnd.distfun_m), 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, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, 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\qcazenav.julia\packages\NCDatasets\h1epE\src\dataset.jl:241 [6] NCDataset(::Function, ::String, ::Vararg{String}) @ NCDatasets C:\Users\qcazenav.julia\packages\NCDatasets\h1epE\src\dataset.jl:238 [7] 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{Float64, 3}}, epsilon2::Vector{Float64}, filename::String, varname::String; datadir::String, bathname::String, bathisglobal::Bool, plotres::typeof(plotres), timeorigin::DateTime, moddim::Vector{Float64}, zlevel::Symbol, ncvarattrib::OrderedDict{String, String}, ncglobalattrib::OrderedDict{Symbol, 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::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, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}) @ DIVAnd C:\Users\qcazenav.julia\packages\DIVAnd\MV3j9\src\diva.jl:375 [8] top-level scope @ .\timing.jl:262 [inlined] [9] top-level scope @ .\In[36]:0 [10] eval @ .\boot.jl:368 [inlined] [11] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String) @ Base .\loading.jl:1428

qcazenave commented 1 year ago

Well, looking more closely at the result of the fit, I wonder if I can consider that it did indeed work... 1e10 for the horizontal scaled correlation length and 1e4 for the vertical corr. length does not look that good, after all...

ctroupin commented 1 year ago

The error with the factorisation can happen when the analysis parameters are not correct, i.e., too large values for example. https://gher-uliege.github.io/DIVAnd.jl/latest/index.html#Error-in-the-factorisation

During the estimation of the parameters there is a way to set bounds (lower and upper values). I need to check the option and let you know.

ctroupin commented 1 year ago

https://gher-uliege.github.io/DIVAnd.jl/latest/#DIVAnd.fithorzlen so you have the option limitfun which is a function that limits the values of the estimated correlation length.

qcazenave commented 1 year ago

Yes I tried with the following: image

And then uncommented the lines: fithorz_param = fithorz_param, fitvert_param = fitvert_param,

But I get the same error and the same orders of magnitude for the corrleation lengths... But I don't understand because the final information about the extrema of the correlation lengths: image do not correspond to the values mentioned before at each depth level: image image

qcazenave commented 1 year ago

I think I will run fithorzlen and fitvertlen independently beforehand to see if I get the same values.

ctroupin commented 1 year ago

That was the suggestion I was writing...

ctroupin commented 1 year ago

In diva3d the coordinates are in meters, while for the functions fithorzlen and fitvertlen called outside diva3d, the units for len should be the same as lon and lat.

qcazenave commented 1 year ago

Hi again again,

So I had it work with the following configuration:

image

(I used distfun=DIVAnd.distfun_m to force the correlation length to be in meters... )

And then:

image

And here is what I end up with:

image

GeoTuxMan commented 1 year ago

Hi @qcazenave !

Very nice plots; how can I do these plots ?

Thank you very much ! George

qcazenave commented 1 year ago

Hello @GeoTuxMan ,

Like this (more or less):

figure(figsize=[15,8])
subplot(121); title("Horizontal correlation lengths")
plot(info[:fithorzlen][:lenf][2:end,1]./1000.,-depthr,ls="--",color="c",marker="*",label="WINTER")
plot(info[:fithorzlen][:lenf][2:end,2]./1000.,-depthr,ls="--",color="g",marker="o",label="SPRING")
plot(info[:fithorzlen][:lenf][2:end,3]./1000.,-depthr,ls="--",color="gold",marker="d",label="SUMMER")
plot(info[:fithorzlen][:lenf][2:end,4]./1000.,-depthr,ls="--",color="brown",marker="v",label="AUTUMN")
legend()
xlabel("Corr. length [km]"),ylabel("Depth [m]")
grid()

subplot(122); title("Vertical correlation lengths")
plot(info[:fitvertlen][:lenf][2:end,1],-depthr,ls="--",color="c",marker="*",label="WINTER")
plot(info[:fitvertlen][:lenf][2:end,2],-depthr,ls="--",color="g",marker="o",label="SPRING")
plot(info[:fitvertlen][:lenf][2:end,3],-depthr,ls="--",color="gold",marker="d",label="SUMMER")
plot(info[:fitvertlen][:lenf][2:end,4],-depthr,ls="--",color="brown",marker="v",label="AUTUMN")
legend()
xlabel("Corr. length [m]"),ylabel("Depth [m]")
grid()

info is the output of the diva3d function (see previous comments) TS = DIVAnd.TimeSelectorYearListMonthList(yearlist,monthlist) withmonthlist = [1:3,4:6,7:9,10:12];

GeoTuxMan commented 1 year ago

Thank you very much ! Unfortunately, diva3d give me a lot of errors (e.g.: LoadError: TaskFailedException) I think the errors come from fithorz_param and fitvert_param.

@show fithorz_param fithorz_param = Dict{Symbol, Any}(:searchz => 50, :limitfun => var"#13#15"(), :distfun => DIVAnd.distfun_m)

My input file it's here: http://86.127.36.60/netcdf/index.html

My code:

using NCDatasets
using PhysOcean
using DataStructures
using DIVAnd
using PyPlot
using Dates
using Statistics
using Random
using Printf

basedir = "C:/Users/Vostro/Documents/my_project_dir_2022/"
bathname = joinpath(basedir,"gebco_30sec_4xxx.nc")
bathisglobal = true

# domain and grid definition
dx = 0.05 #  longitude resolution in degrees
dy = 0.05 # latitude resolution in degrees
# MSFD Black Sea (with Marmara)
lonr = 26.5:dx:42 # the range of longitudes (start:step:end)
latr = 40:dy:48 # the range of latitudes (start:step:end) ;

@time bx,by,b = load_bath(bathname,true,lonr,latr);
@show extrema(b)

# working directory where results will be written
wkdir=joinpath(basedir,"2022")
if isdir(wkdir)
    @info("Output dir already exists")
else
    mkdir(wkdir)
    @info("Creating output dir")
end

# variable name and ODV files to read
varname = "Water body dissolved oxygen concentration"
wkdir

# season definition
saisons=["Winter","Spring","Summer","Autumn"]
months=["(Dec-Feb)","(March-May)","(June-Aug)","(Sept-Nov)"]
monthlists = [[12,1,2],[3,4,5],[6,7,8],[9,10,11]]
year_start=1970:2016
year_window = 6 #  2018-1970=48 + 1

# depth levels
depthr = [0.,5.,10.,20.,30.,40.,50.,75.,100.,125.,150.,200.,250.];

srZ=50.
epsilon2=(1/1.6962)
slog=200.
#-----------------------------
@info("Starting computations")
@info("grid size")
@info("lon: $(length(lonr))")
@info("lat: $(length(latr))")
@info("depth: $(length(depthr))")
@info("N seasons:  $(length(monthlists))")

# Time origin for the NetCDF file
timeorigin = DateTime(1900,1,1,0,0,0) ;

# attributes for the NETCDF file 

    metadata = OrderedDict(
        "project" => "EMODNET-chemistry",    
    "institution_urn" => "SDN:EDMO::697",    
    "Author_e-mail" => ["Gheorghe Sarbu <ghesarbu@gmail.com>", "Luminita Buga <lbuga@web.de>"],   
    "source" => "Observational data from SeaDataNet/EMODNet Chemistry Data Network",    
    "comment" => "All Seasons",    
"parameter_keyword_urn" => "SDN:P35::EPC00002", # Water body dissolved oxygen concentration    
"search_keywords_urn" => ["SDN:P02::DOXY"], # Water body dissolved oxygen concentration    
    "area_keywords_urn" => ["SDN:C19::3_3"], # Black Sea
    "product_version" => "1.0",    
    "netcdf_standard_name" => "$(replace(varname," " => "_"))",
    "netcdf_long_name" => "Water body dissolved oxygen concentration",
    "netcdf_units" => "umol/l",
    # Abstract for the product
    "abstract" => "Dissolved oxygen concentration - Seasonal Climatology for NW Black Sea (region in front of Danube Delta months) for the period 1970-2018 on the domain: Lon E28.5-30.5°E, Lat N43.7-45.6°N.  (winter: December-February, - spring: March-May, - summer: June-August, - autumn: September-November)
Data Sources: observational data from SeaDataNet/EMODNet Chemistry Data Network.
Description of DIVA analysis: The computation was done with the DIVAnd (Data-Interpolating Variational Analysis in n dimensions), version 2.7.1, using GEBCO 30sec topography for the spatial connectivity of water masses. Horizontal correlation length of 47.730km and vertical correlation length of 45m were applied. Depth range: 0, 2.0, 4.0, 6.0, 8.0, 10.0, 15.0, 20.0, 25.0, 30.0, 40.0, 50.0 m. Units: umol/l. 
The horizontal resolution of the produced DIVAnd maps grids is 0.01 degrees.",    
    "doi" => "10.6092/7qst-j180")

# number of climatologies to compute
nc=year_start[end]-year_start[1]+1
@info("will compute $(nc) climatologies")
#yearlist=Array{Range}(nc)
#yearlist=UnitRange{Int64}(nc)
yearlist = Array{UnitRange{Int64},1}(undef, nc)

for i in 1:length(year_start)
    #@show year_start[i]
    yearlist[i]=year_start[i]:(year_start[i]+year_window-1)
end
@show yearlist
year_lim=[year_start[1],(year_start[end]+year_window-1)] ;
@info("year boundaries : $(year_lim))")

#-------------------------------
#reading NC files and removing duplicates
#------------------------------
filename1="data_from_Black_and_Marmara_DIVA_O2-test.nc"
@info("Start reading NetCDF file")
@info(Dates.now())
value,lon,lat,depth,obstime,ids = NCODV.load(Float64, filename1, "Water body dissolved oxygen concentration"; qv_flags = ["good_value", "probably_good_value", "changed_value", "value_below_detection", "interpolated_value", "value_below_limit_of_quantification"]);
@info(Dates.now())
@info("End reading NetCDF file")

######## FIND duplicates ######################
@time dupl = DIVAnd.Quadtrees.checkduplicates((lon,lat,depth,obstime), value, (0.01,0.01,0.01,1/(24*60)),0.01);
index = findall(.!isempty.(dupl));
ndupl = length(index);
pcdupl = round(ndupl / length(lon) * 100; digits=2);
@info("Number of possible duplicates: $ndupl")
@info("Percentage of duplicates: $pcdupl%")
######## REMOVE DUPLICATES from filename1 #####
@show(dupl);
index2keep=[dd[1] for dd in dupl]
lon=lon[index2keep]
lat=lat[index2keep]
depth=depth[index2keep]
obstime=obstime[index2keep]
value=value[index2keep]
ids=ids[index2keep]
###### END ###################################
# epsilon
eps1=fill(epsilon2[1],length(lon))

# select data over the requested time period
sel = (year_lim[1] .<= Dates.year.(obstime) .<= year_lim[end])
value = value[sel]
lon = lon[sel]
lat = lat[sel]
depth = depth[sel]
obstime = obstime[sel]
ids = ids[sel];
eps1=eps1[sel];
checkobs((lon,lat,depth,obstime),value,ids)
#--------------------------------------------------------------
# Correlation length
#--------------------------------------------------------------
# size for correlation length tables
sz = (length(lonr),length(latr),length(depthr)) ;
lenx = fill(166_500.,sz)
leny = fill(166_500.,sz)
#  Vertical Corr Length
ddepthr = depthr[2:end]-depthr[1:end-1] ;
ddepthr_below = [ddepthr[1], ddepthr...] ;
ddepthr_above = [ddepthr...,ddepthr[end]] ;
lend = [max(25,5*i,5*j) for (i,j) in zip(ddepthr_below,ddepthr_above)] ;
lenz = zeros(sz); 
for i in 1:sz[1]
    for j in 1:sz[2]
        lenz[i,j,:] = lend ;
    end;
end;
len = (lenx, leny, lenz);

# CL computed with diva3d; set boundaries
fitcorrlen = true;

if fitcorrlen
   corr_length_info = "Correlation lengths optimized by DIVAnd"
   len = (fill(1.,sz),fill(1.,sz),fill(1.,sz))
   fithorz_param = Dict(:limitfun => (z,len) -> max(min(len,166_500.),11_000.),
                        :searchz => 50,
                        :distfun => DIVAnd.distfun_m)
   fitvert_param = Dict(:limitfun => (z,len) -> max(min(len,200.),15.),
                        :searchz => 50,
                        :distfun => DIVAnd.distfun_m)
else
    fithorz_param = nothing
    fitvert_param = nothing
end;

@show fithorz_param

function plotres2(timeindex,sel,fit,erri)
    println("ok")
end

#--------------------------------------
# compute climatologies
#-------------------------------
    @info("starting DIVA computations for $(saisons[1])")
    @info(Dates.now())

    # Time selection
    TS = DIVAnd.TimeSelectorYearListMonthList(yearlist,monthlists) # toti anii
    #TS = DIVAnd.TimeSelectorYearListMonthList(yearlist[20:20],[monthlists[i]])  # un an,de test

    # NETCFD file
    # File name based on the variable (but all spaces are replaced by _)
    filename = joinpath(wkdir,"$(replace(varname," " => "_"))_ALLBS_testCL.4Danl.nc")
    if isfile(filename)
       rm(filename) # delete the previous analysis
    end
    @info("Will write results in $filename")
    # create attributes for the netcdf file (need an internet connexion)
    ncglobalattrib,ncvarattrib = SDNMetadata(metadata,filename,varname,lonr,latr)

    dbinfo = diva3d(
    (lonr,latr,depthr,TS),
    (lon,lat,depth,obstime), value,
    len, eps1,
    filename, varname,
    bathname = bathname,
    bathisglobal = bathisglobal,
    plotres = plotres2,
    timeorigin = timeorigin,
    fitcorrlen = fitcorrlen,
    fithorz_param = fithorz_param, 
    fitvert_param = fitvert_param,
    niter_e = 1,
    ncvarattrib = ncvarattrib,
    ncglobalattrib = ncglobalattrib,
    error_thresholds=[("L1", 0.3), ("L2", 0.5)]
);

    # Add used observations in the netcfd file
     DIVAnd.saveobs(filename,(lon,lat,depth,obstime),ids,used = dbinfo[:used])
     DIVAnd.derived(filename,varname,filename;error_thresholds=[("L1",0.3),("L2",0.5)],)
ctroupin commented 1 year ago

Hi George, could you put your code inside ``` symbols, for instance:

function myfun()
...

so it's properly formatted?

GeoTuxMan commented 1 year ago

I put the code here: http://86.127.36.60/netcdf/index.html

ctroupin commented 1 year ago

thanks, I've downloaded all and am now testing here. Let's close this issue and open another one if needed.

Alexander-Barth commented 1 year ago

For your information, I added some explications about the relative correlation length to the documentation:

https://github.com/gher-uliege/DIVAnd.jl/blob/3e67008985ae7ed2c6cf6ca904bcfa0bf0cf8134/src/diva.jl#L24-L31

ctroupin commented 11 months ago

@Alexander-Barth @jmbeckers here there is a PosDefException: matrix is not positive definite; Cholesky factorization failed. error with the data file available.

jmbeckers commented 11 months ago

https://sparsearrays.juliasparse.org/dev/solvers/#LinearAlgebra.cholesky

Maybe we can play with the shift parameter

Alexander-Barth commented 11 months ago

Do we have a reproducer with some realistic parameters?

I can only reproduce it with a vert. correlation of 40 with a vertical domain from [0,1]:

using DIVAnd
using Test
fun(x, y, z) = sin(6x) * cos(6y) * sin(6z)
mask, (pm, pn, po), (xi, yi, zi) = DIVAnd_squaredom(3, range(0, stop = 1, length = 50))

fi_ref = fun.(xi, yi, zi)

ϵ = eps()
x, y, z = ndgrid(
    range(ϵ, stop = 1 - ϵ, length = 10),
    range(ϵ, stop = 1 - ϵ, length = 10),
    range(ϵ, stop = 1 - ϵ, length = 10),
);
x = x[:];
y = y[:];
z = z[:];

f = fun.(x, y, z)
epsilon2 = 0.01;

# works with len[3] == 30
len = (0.1, 0.1, 40);

# fi is the interpolated field
fi, s =
    @time DIVAndrun(mask, (pm, pn, po), (xi, yi, zi), (x, y, z), f, len, epsilon2);

The diagonal elements of the matrix iB scale like (len[3]*po[1])^5 (here 1e16) while off-diagonal elements are of the order of 1. Any solver (in double precision) will likely to have a problem with such matrix.

jmbeckers commented 10 months ago

Maybe also add a test in function checkresolution to check when resolution is way too fine or length scale is much larger than domain size ? In the latter case, maybe suggest to drop that dimension in the analysis