From 8098366c53cf0ba8a86c9c12f8c1e21bf4dcba29 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Tue, 17 Jan 2023 10:49:14 +0100 Subject: [PATCH 01/51] bug correction --- src/special_matrices.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/special_matrices.jl b/src/special_matrices.jl index c4a610be..4f2a36b6 100644 --- a/src/special_matrices.jl +++ b/src/special_matrices.jl @@ -30,7 +30,7 @@ function Base.:*(C::CovarIS, v::TV)::TV where {TV<:AbstractVector{Float64}} @debug "checksum $(sum(C.IS)) $(sum(v))" @debug "size $(size(C.IS)) $(size(v))" - log = false + log = true @debug begin log = true end @@ -47,8 +47,8 @@ function Base.:*(C::CovarIS, v::TV)::TV where {TV<:AbstractVector{Float64}} @debug begin @show norm(C.IS * x - v) end - #@show convergence_history - return x + @show convergence_history + return x elseif C.factors != nothing return C.factors \ v else From 728f012cdfd566626c4f431c2464105d21c1f5e9 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Wed, 18 Jan 2023 17:09:48 +0100 Subject: [PATCH 02/51] Update DIVAndjog.jl https://github.com/gher-uliege/DIVAnd.jl/issues/117 --- src/DIVAndjog.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/DIVAndjog.jl b/src/DIVAndjog.jl index 54cb95cb..e7ac66ab 100644 --- a/src/DIVAndjog.jl +++ b/src/DIVAndjog.jl @@ -1249,9 +1249,9 @@ function DIVAndjog( PC1 = 1 sc = 0 - if size(lmask)[2] > 2 + #if size(lmask)[2] > 2 fc, sc = DIVAndrun(maskc, pmnc, xic, x, f, Labsccut, 10000.0; otherargsc...) - end + #end if !(sc == 0) # TEST makes sure there are values in the coarse resolution solution From 4ab1c8216416c0ad6c7ee846110423c2b980fb82 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Wed, 25 Jan 2023 06:41:42 +0100 Subject: [PATCH 03/51] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index a42a6d02..fa9b9aa8 100644 --- a/README.md +++ b/README.md @@ -163,7 +163,7 @@ Tools to help you are included in ([DIVAnd_cv.jl](https://github.com/gher-ulieg ## Note about the error fields -`DIVAnd` allows the calculation of the analysis error variance, scaled by the background error variance. Though it can be calculated "exactly" using the diagonal of the error covariance matrix s.P, it is too costly and approximations are provided. Two version are recommended, `DIVAnd_cpme` for a quick estimate and `DIVAnd_aexerr` for a version closer the theoretical estimate (see [Beckers et al 2014](https://doi.org/10.1175/JTECH-D-13-00130.1) ) +`DIVAnd` allows the calculation of the analysis error variance, scaled by the background error variance. Though it can be calculated "exactly" using the diagonal of the error covariance matrix s.P, it is generally too costly and approximations are provided. All of them are accessible as options via `DIVAnd_errormap` or you can let `DIVAnd` decide which version to use (possibly by specifying if you just need a quick estimate or a version closer the theoretical estimate) (see [Beckers et al 2014](https://doi.org/10.1175/JTECH-D-13-00130.1) ) ## Advanced usage From 305f59503d778b6a83360b5dccea648ac3056494 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Wed, 25 Jan 2023 11:21:48 +0100 Subject: [PATCH 04/51] Update README.md --- README.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/README.md b/README.md index fa9b9aa8..8a0729a7 100644 --- a/README.md +++ b/README.md @@ -214,3 +214,8 @@ Note that only [official julia builds](https://julialang.org/downloads/) are sup # Fun An [educational web application](http://data-assimilation.net/Tools/divand_demo/html/) has been developed to reconstruct a field based on point "observations". The user must choose in an optimal way the location of 10 observations such that the analysed field obtained by `DIVAnd` based on these observations is as close as possible to the original field. + +# You do not want to use Julia + + +https://github.com/gher-uliege/DIVAnd.py From b750052d068f59c1ad5a7cb5eb85258f5fabd3c2 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Wed, 25 Jan 2023 11:25:53 +0100 Subject: [PATCH 05/51] Update README.md --- README.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/README.md b/README.md index 8a0729a7..f14dd9ff 100644 --- a/README.md +++ b/README.md @@ -217,5 +217,8 @@ An [educational web application](http://data-assimilation.net/Tools/divand_demo/ # You do not want to use Julia +You should really reconsider and try out Julia. It is easy to use and provides the native interface to `DIVAnd`. + +If you have a stable workflow using python, into which you want to integrate `DIVAnd`, you might try https://github.com/gher-uliege/DIVAnd.py From a4c9813f26be87fccd30de67e4e1da633444f192 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Sat, 28 Jan 2023 17:08:22 +0100 Subject: [PATCH 06/51] First step to fixing alphabc problem --- src/DIVAnd_bc_stretch.jl | 20 +++++++++++++------- src/DIVAndrun.jl | 2 +- 2 files changed, 14 insertions(+), 8 deletions(-) diff --git a/src/DIVAnd_bc_stretch.jl b/src/DIVAnd_bc_stretch.jl index 37e7fc71..c0aafedf 100644 --- a/src/DIVAnd_bc_stretch.jl +++ b/src/DIVAnd_bc_stretch.jl @@ -1,7 +1,13 @@ +function DIVAnd_bc_stretch(mask, pmnin, xiin, Lin, moddim, alphabc::Number = 1.0) + + return DIVAnd_bc_stretch(mask, pmnin, xiin, Lin, moddim, alphabc*ones(ndims(mask))) + +end + """ """ -function DIVAnd_bc_stretch(mask, pmnin, xiin, Lin, moddim, alphabc = 1) +function DIVAnd_bc_stretch(mask, pmnin, xiin, Lin, moddim, alphabc::Vector{Float64} = ones(ndims(mask))) # number of dimensions n = ndims(mask) @@ -18,12 +24,12 @@ function DIVAnd_bc_stretch(mask, pmnin, xiin, Lin, moddim, alphabc = 1) # Just used to fill the Labs tuple (so background will not fill it again) # - if alphabc == 0 + if sum(alphabc) == 0 #@warn "DIVAnd_bc_stretch was just used to fill in Labs" return pmnin, xiin, Labs end - if alphabc > 0 + if sum(alphabc) > 0 pmn = deepcopy(pmnin) xi = deepcopy(xiin) @@ -38,7 +44,7 @@ function DIVAnd_bc_stretch(mask, pmnin, xiin, Lin, moddim, alphabc = 1) max.( ( ( - 2.0 .* alphabc .* Labs[i][ind1...] .* pmnin[i][ind2...] .- + 2.0 .* alphabc[i] .* Labs[i][ind1...] .* pmnin[i][ind2...] .- 1.0 ) .* pmnin[i][ind1...] .- pmnin[i][ind2...] ) ./ (pmnin[i][ind1...] + pmnin[i][ind2...]), @@ -55,7 +61,7 @@ function DIVAnd_bc_stretch(mask, pmnin, xiin, Lin, moddim, alphabc = 1) max.( ( ( - 2.0 .* alphabc .* Labs[i][ind1...] .* pmnin[i][ind2...] .- + 2.0 .* alphabc[i] .* Labs[i][ind1...] .* pmnin[i][ind2...] .- 1.0 ) .* pmnin[i][ind1...] - pmnin[i][ind2...] ) ./ (pmnin[i][ind1...] + pmnin[i][ind2...]), @@ -78,7 +84,7 @@ function DIVAnd_bc_stretch(mask, pmnin, xiin, Lin, moddim, alphabc = 1) wjmb[ind1...] = 1.0 ./ max.( - (2 * alphabc .* Labs[i][ind1...] .- 1.0 ./ wjmb[ind2...]), + (2 * alphabc[i] .* Labs[i][ind1...] .- 1.0 ./ wjmb[ind2...]), 1.0 ./ wjmb[ind2...], ) @@ -88,7 +94,7 @@ function DIVAnd_bc_stretch(mask, pmnin, xiin, Lin, moddim, alphabc = 1) wjmb[ind1...] = 1.0 ./ max.( - (2 * alphabc .* Labs[i][ind1...] .- 1.0 ./ wjmb[ind2...]), + (2 * alphabc[i] .* Labs[i][ind1...] .- 1.0 ./ wjmb[ind2...]), 1.0 ./ wjmb[ind2...], ) end diff --git a/src/DIVAndrun.jl b/src/DIVAndrun.jl index 9dd575f9..f658e689 100644 --- a/src/DIVAndrun.jl +++ b/src/DIVAndrun.jl @@ -116,7 +116,7 @@ function DIVAndrun( # add all additional constrains for i = 1:length(constraints) #@info "Adding additional constrain - $(i)" - s = DIVAnd_addc(s, constraints[i]) + s = DIVAnd_addc(s, constraints[i]) end # factorize a posteriori error covariance matrix From db546292758e6d10feb68fb6b4b4e29357928db6 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Sat, 28 Jan 2023 17:42:16 +0100 Subject: [PATCH 07/51] Update DIVAndjog.jl --- src/DIVAndjog.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/DIVAndjog.jl b/src/DIVAndjog.jl index e7ac66ab..ab52ab29 100644 --- a/src/DIVAndjog.jl +++ b/src/DIVAndjog.jl @@ -1248,10 +1248,10 @@ function DIVAndjog( Labsccut = ([Labsc[i] * lmask1[i] for i = 1:n]...,) PC1 = 1 sc = 0 - - #if size(lmask)[2] > 2 + #@show size(lmask),size(lmaskl) + if size(lmask)[1] > 2 fc, sc = DIVAndrun(maskc, pmnc, xic, x, f, Labsccut, 10000.0; otherargsc...) - #end + end if !(sc == 0) # TEST makes sure there are values in the coarse resolution solution From 84b49b839be277227d623bfa1cbe0404b04543b7 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Mon, 30 Jan 2023 10:07:11 +0100 Subject: [PATCH 08/51] Update README.md --- README.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/README.md b/README.md index f14dd9ff..023e0fc4 100644 --- a/README.md +++ b/README.md @@ -139,6 +139,8 @@ More examples are available in the notebooks from the [Diva Workshop](https://gi `DIVAndgo` is only needed for very large problems when a call to `DIVAndrun` leads to memory or CPU time problems. This function tries to decide which solver (direct or iterative) to use and how to make an automatic domain decomposition. Not all options from `DIVAndrun` are available. +If you want to try out multivariate approaches, you can look at `DIVAnd_multivarEOF` and `DIVAnd_multivarJAC` + ## Note about the background field If zero is not a valid first guess for your variable (as it is the case for e.g. ocean temperature), you have to subtract the first guess from the observations before calling `DIVAnd` and then add the first guess back in. From 00f558881cfca9b9c15e9896a47f64307e138e3a Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Mon, 30 Jan 2023 10:10:02 +0100 Subject: [PATCH 09/51] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 023e0fc4..b4558b1d 100644 --- a/README.md +++ b/README.md @@ -130,7 +130,7 @@ More examples are available in the notebooks from the [Diva Workshop](https://gi `DIVAndrun` is the core analysis function in n dimensions. It does not know anything about the physical parameters or units you work with. Coordinates can also be very general. The only constraint is that the metrics `(pm,pn,po,...)` when multiplied by the corresponding length scales `len` lead to non-dimensional parameters. Furthermore the coordinates of the output grid `(xi,yi,zi,...)` need to have the same units as the observation coordinates `(x,y,z,...)`. -`DIVAndfun` is a version with a minimal set of parameters (the coordinates and values of observations) `(x,f)` and provides and interpolation function rather than an already gridded field. +`DIVAndfun` is a version with a minimal set of parameters (the coordinates and values of observations, i.e. `(x,f)`, the remaining parameters being optional) and provides an interpolation *function* rather than an already gridded field. `diva3D` is a higher-level function specifically designed for climatological analysis of data on Earth, using longitude/latitude/depth/time coordinates and correlations length in meters. It makes the necessary preparation of metrics, parameter optimizations etc you normally would program yourself before calling the analysis function `DIVAndrun`. From f98545dc887661007ab1de2edaa3a5df412cb438 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Wed, 1 Feb 2023 11:27:37 +0100 Subject: [PATCH 10/51] Update README.md --- README.md | 1 + 1 file changed, 1 insertion(+) diff --git a/README.md b/README.md index b4558b1d..d265905a 100644 --- a/README.md +++ b/README.md @@ -212,6 +212,7 @@ Please include the following information when reporting an issue: * The command and their arguments which produced the error Note that only [official julia builds](https://julialang.org/downloads/) are supported. +In all cases, if we provide a tentative solution, please provide a feedback in all cases (whether it solved your issue or not). # Fun From df4b25bbf92988f61333e943ebf322a313456c2d Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Wed, 1 Feb 2023 11:28:18 +0100 Subject: [PATCH 11/51] Update README.md --- README.md | 1 + 1 file changed, 1 insertion(+) diff --git a/README.md b/README.md index d265905a..74f4f060 100644 --- a/README.md +++ b/README.md @@ -212,6 +212,7 @@ Please include the following information when reporting an issue: * The command and their arguments which produced the error Note that only [official julia builds](https://julialang.org/downloads/) are supported. + In all cases, if we provide a tentative solution, please provide a feedback in all cases (whether it solved your issue or not). # Fun From 50af9054164546943b3394f62dccee0b75cf481f Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Wed, 1 Feb 2023 13:20:44 +0100 Subject: [PATCH 12/51] Update fit.jl Try to fix issues https://github.com/gher-uliege/DIVAnd.jl/issues/121#top https://github.com/gher-uliege/DIVAnd.jl/issues/123#issuecomment-1411889797 --- src/fit.jl | 22 ++++++++++++++++++++-- 1 file changed, 20 insertions(+), 2 deletions(-) diff --git a/src/fit.jl b/src/fit.jl index e607ed15..7c5c5f08 100644 --- a/src/fit.jl +++ b/src/fit.jl @@ -426,7 +426,18 @@ function fitlen( ) if length(d) == 0 @warn "no data is provided to fitlen" - return NaN, NaN, Dict{Symbol,Any}() + dbinfo = Dict{Symbol,Any}( + :covar => [], + :fitcovar => [], + :distx => [], + :rqual => 0.0, + :range => 1:0, + :covarweight => [], + :distx => [], + :sn => 0.0, + :meandist => 0.0, + ) + return NaN, NaN, dbinfo end # number of dimensions @@ -505,7 +516,14 @@ function fitlen( @debug "Number of probable active bins: $rnbins" ddist = meandist / rnbins - nbmax = floor(Int, maxdist / ddist + 1) + nbmax=1 + worktmp=maxdist / ddist + if !isnan(worktmp)&&isfinite(worktmp) + nbmax = floor(Int, worktmp + 1) + end + + + #nbmax = floor(Int, maxdist / ddist + 1) @debug "distance for binning: $ddist" @debug "maximum number of bins: $nbmax" From 2023e4e9e061dbfb2779d9ceb8c54fed64ab5844 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Fri, 3 Feb 2023 14:08:11 +0100 Subject: [PATCH 13/51] Update DIVAnd_errormap.jl --- src/DIVAnd_errormap.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/DIVAnd_errormap.jl b/src/DIVAnd_errormap.jl index 96cb2309..01b402ca 100644 --- a/src/DIVAnd_errormap.jl +++ b/src/DIVAnd_errormap.jl @@ -210,7 +210,7 @@ function DIVAnd_errormap( ScalebyB = false end - if errmethod == "exact" && Bscale + if errmethod == :exact && Bscale # Or maybe if all info is there run locally ? Yes probably possible as aexerr also needs all infos ? warn("You need to do that scaling by yourself, running diva again with a very high R matrix and divide by this second map") ScalebyB = false @@ -221,7 +221,7 @@ function DIVAnd_errormap( errmethod = cpme end - if errmethod == "exact" && noP + if errmethod == :exact && noP warn("Sorry, that method needs s.P to be available. Will use aexerr instead") errmethod = aexerr end @@ -269,7 +269,7 @@ function DIVAnd_errormap( return scpme, errmethod end - if errmethod == "exact" + if errmethod == :exact errormap, =statevector_unpack(s.sv,diag(s.P)) From 9518d49c6a52d007e7c75355559e8e26bd2a0382 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Fri, 24 Mar 2023 13:41:14 +0100 Subject: [PATCH 14/51] norm components --- src/DIVAnd.jl | 4 +++ src/DIVAnd_diagnostics.jl | 57 +++++++++++++++++++++++++++++++++++++++ test/test_2dvar_adv.jl | 10 +++++++ 3 files changed, 71 insertions(+) create mode 100644 src/DIVAnd_diagnostics.jl diff --git a/src/DIVAnd.jl b/src/DIVAnd.jl index 8eb4baa5..04b9f0b2 100644 --- a/src/DIVAnd.jl +++ b/src/DIVAnd.jl @@ -647,6 +647,10 @@ export DIVAnd_diagapp export DIVAnd_errormap +include("DIVAnd_diagnostics.jl") + +export DIVAnd_norms + # old function name (to be depreciated) const sparse_gradient = DIVAnd_gradient export sparse_gradient diff --git a/src/DIVAnd_diagnostics.jl b/src/DIVAnd_diagnostics.jl new file mode 100644 index 00000000..faf64544 --- /dev/null +++ b/src/DIVAnd_diagnostics.jl @@ -0,0 +1,57 @@ +""" + norm1,norm2,norm3,epsilon2=DIVAnd_norms(fi,s) + +### Input: + +* `fi` and `s` from the corresponding `DIVANDrun` exectution `fi,s = DIVAndrun(...` + +### Output: + +the different terms of the norm which is minimised + +* `norm1`: the regularity terms + +* `norm2`: the data-misfit terms + +* `norm3`: the remaining constraints + +* `epsilon2` : the average diagonal of `R`; that allows to look at the L-shape curve by looking at `Log(norm1)` vs `Log(norm2*epsilon2)` + + + +""" +function DIVAnd_norms(fi,s) + + R = s.R + iB = s.iB + H = s.H + fstate = statevector_pack(s.sv, (fi,)) + residuals=s.yo - (s.H) * fstate + residualsobs=s.obsconstrain.yo - (s.obsconstrain.H) *fstate + residualsobs[s.obsout] .= 0.0 + + norm2= residualsobs'*(s.obsconstrain.R\residualsobs) + norm3= residuals'*(R \ residuals)-norm2 + norm1=fstate'*(iB*fstate) + dR=diag(s.obsconstrain.R) + epsilon2=mean(dR[isfinite.(dR)]) + + return norm1,norm2,norm3,epsilon2 + +end + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/test/test_2dvar_adv.jl b/test/test_2dvar_adv.jl index 212c19f7..15e2f131 100644 --- a/test/test_2dvar_adv.jl +++ b/test/test_2dvar_adv.jl @@ -29,7 +29,17 @@ fi, s = DIVAndrun( @test abs(fi[18, 24] - 0.8993529043140029) < 1e-2 +norm1,norm2,norm3,epsilon2=DIVAnd_norms(fi,s) +#@show norm1,norm2,norm3,epsilon2 + +@test norm1 ≈ 4.864863745730252 + +@test norm2 ≈ 0.1276096541011819 + +@test norm3 ≈ 0.059450077426666054 + +@test epsilon2 ≈ 1 / 200 # Copyright (C) 2014, 2017 Alexander Barth # From b950448f2f83dc470ba9b22af9d8df0957ae18e5 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Fri, 24 Mar 2023 16:42:22 +0100 Subject: [PATCH 15/51] With inequalities --- src/DIVAnd.jl | 8 ++++ src/DIVAndrun.jl | 104 ++++++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 111 insertions(+), 1 deletion(-) diff --git a/src/DIVAnd.jl b/src/DIVAnd.jl index 04b9f0b2..87a4e9fa 100644 --- a/src/DIVAnd.jl +++ b/src/DIVAnd.jl @@ -71,6 +71,14 @@ mutable struct DIVAnd_constrain{ H::TH end +mutable struct DIVAnd_ineqconstrain{ + T<:AbstractFloat, + TH<:AbstractMatrix{<:Number}, +} + yo::Vector{T} + H::TH +end + # T is the type of floats and # Ti: the type of integers # N: the number of dimensions diff --git a/src/DIVAndrun.jl b/src/DIVAndrun.jl index f658e689..f6b03566 100644 --- a/src/DIVAndrun.jl +++ b/src/DIVAndrun.jl @@ -15,6 +15,7 @@ function DIVAndrun( maxit = 100, minit = 0, constraints = (), + ineqconstraints = (), inversion = :chol, moddim = [], fracindex = Matrix{T}(undef, 0, 0), @@ -39,6 +40,105 @@ function DIVAndrun( mean_Labs = nothing, ) where {N,T} +# Inequality constraints via loop around classical analysis (recursive call) + + + + if !isempty(ineqconstraints) + + + +function Ineqtoeq(yo,H,xo;tol1=1e-7,tol2=1e-8) + # Translates H x - yo > 0 into a quadratic weak constraint depending on a first guess of the statevector xo + vv=ones(size(yo)[1]).*Inf + Hxo=H*xo + @show sum(Hxo.0 + ineqok=false + @show violated,ntries + end + allconstraints=(allconstraints...,onemoreconstraint)# Add to the list of constraints + end + + + + end + return fi,s + end + # End of special treatment for inequality constraints + # check pmn .* len > 4 #checkresolution(mask,pmnin,lin) @@ -215,7 +315,9 @@ For oceanographic application, this is the land-sea mask where sea is true and l m=2 1 3 3 1 (n=3,4) ... -* `constraints`: a structure with user specified constraints (see `DIVAnd_addc`). +* `constraints`: a structure with user specified quandratic constraints (see `DIVAnd_addc`). + +* `ineqconstraints`: a structure with user specified inequality constraints `Hx >= y0`. There is no check if the inequality constraints make sense are compatible with each other or with the data. * `moddim`: modulo for cyclic dimension (vector with n elements). Zero is used for non-cyclic dimensions. One should not include a boundary From 3f79103579967a7a73d34a0865e30d0012552b98 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Fri, 24 Mar 2023 17:44:15 +0100 Subject: [PATCH 16/51] Update test_2dvar_adv.jl --- test/test_2dvar_adv.jl | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/test/test_2dvar_adv.jl b/test/test_2dvar_adv.jl index 15e2f131..02b5febd 100644 --- a/test/test_2dvar_adv.jl +++ b/test/test_2dvar_adv.jl @@ -41,6 +41,26 @@ norm1,norm2,norm3,epsilon2=DIVAnd_norms(fi,s) @test epsilon2 ≈ 1 / 200 +m=sum(mask) +m2c=DIVAnd.DIVAnd_ineqconstrain(0*ones(Float64,m),Diagonal(ones(Float64,m))) +x = [0.4,0.5] +y = [0.4,0.5] +f = [1.0,-0.1] +fi, s = DIVAndrun( + mask, + (pm, pn), + (xi, yi), + (x, y), + f, + len, + epsilon2; + velocity = (u, v), + alphabc = 0, + ineqconstraints=(m2c,) +); + +@test fi[18,24] ≈ 0.5381625233419283 + # Copyright (C) 2014, 2017 Alexander Barth # # This program is free software; you can redistribute it and/or modify it under From afab248748c8c4bd7e025143282f276efea55828 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Fri, 24 Mar 2023 18:16:52 +0100 Subject: [PATCH 17/51] Update DIVAndrun.jl --- src/DIVAndrun.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/DIVAndrun.jl b/src/DIVAndrun.jl index f6b03566..f959e7e5 100644 --- a/src/DIVAndrun.jl +++ b/src/DIVAndrun.jl @@ -16,6 +16,7 @@ function DIVAndrun( minit = 0, constraints = (), ineqconstraints = (), + ntriesmax = 10, inversion = :chol, moddim = [], fracindex = Matrix{T}(undef, 0, 0), @@ -67,7 +68,7 @@ end ntries=0 fi=() s=() - ntriesmax=10 + while !ineqok && ntries= y0`. There is no check if the inequality constraints make sense are compatible with each other or with the data. +* `ineqconstraints`: a structure with user specified inequality constraints `Hx >= y0`. There is no check if the inequality constraints make sense are compatible with each other or with the data. Inequalities will not be satisfied exactly everywhere unless they are already satisfied with a normal analysis. You can increase the number of iterations by increasing `ntriesmax`. * `moddim`: modulo for cyclic dimension (vector with n elements). Zero is used for non-cyclic dimensions. One should not include a boundary From caafbc223c6b62eec461957ffb7c31740ce73fe0 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Sat, 25 Mar 2023 07:50:30 +0100 Subject: [PATCH 18/51] Update DIVAndrun.jl --- src/DIVAndrun.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/DIVAndrun.jl b/src/DIVAndrun.jl index f959e7e5..a30d1158 100644 --- a/src/DIVAndrun.jl +++ b/src/DIVAndrun.jl @@ -318,7 +318,7 @@ For oceanographic application, this is the land-sea mask where sea is true and l * `constraints`: a structure with user specified quandratic constraints (see `DIVAnd_addc`). -* `ineqconstraints`: a structure with user specified inequality constraints `Hx >= y0`. There is no check if the inequality constraints make sense are compatible with each other or with the data. Inequalities will not be satisfied exactly everywhere unless they are already satisfied with a normal analysis. You can increase the number of iterations by increasing `ntriesmax`. +* `ineqconstraints`: a structure with user specified inequality constraints such that the analysis `x` satisfies`Hx >= y0`. There is no check if the inequality constraints make sense are compatible with each other or with the data. Inequalities will not be satisfied exactly everywhere unless they are already satisfied with a normal analysis. You can increase the number of iterations by increasing `ntriesmax`. * `moddim`: modulo for cyclic dimension (vector with n elements). Zero is used for non-cyclic dimensions. One should not include a boundary From 56eb29faff67081306955643d5a88f805ba662bb Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Tue, 28 Mar 2023 08:48:17 +0200 Subject: [PATCH 19/51] more tests (formal only) --- test/runtests.jl | 7 ++++ test/test_multivar.jl | 74 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 81 insertions(+) create mode 100644 test/test_multivar.jl diff --git a/test/runtests.jl b/test/runtests.jl index f7ff9688..c8a0ba62 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -29,6 +29,13 @@ using Downloads: download include("test_1dvar.jl") include("test_1D_seminormed.jl") + + + # multivar + + include("test_multivar.jl") + + # dynamical constraints include("test_2dvar_adv.jl") diff --git a/test/test_multivar.jl b/test/test_multivar.jl new file mode 100644 index 00000000..e3646481 --- /dev/null +++ b/test/test_multivar.jl @@ -0,0 +1,74 @@ +# A simple example of DIVAnd in 2 dimensions +# with observations from an analytical function. + +using Test +using StableRNGs +using Random +using DIVAnd + +# true error variance of observation +epsilon2_true = 1.0 + +rng = StableRNG(1234) +#rng = MersenneTwister(1234) + +# observations +nobs = 99 +x = rand(rng,nobs); +y = rand(rng,nobs); +v = rand(rng,[1,2],nobs); + +# true length-scale +len_true = 0.5 +f = (1.5 .-v ).*sin.(π * x / len_true) .* cos.(π * y / len_true); + +# add noise +f = f + sqrt(epsilon2_true) * randn(rng,nobs); + +# final grid +mask, (pm, pn,pv), (xi, yi,vi) = + DIVAnd_rectdom(range(0, stop = 1, length = 14), range(0, stop = 1, length = 13),1:2) + +# correlation length (first guess) +len = 0.1; + +# obs. error variance normalized by the background error variance (first guess) +epsilon2 = 2.0 + +# loop over all methods + +fi,s=DIVAnd_multivarEOF(mask,(pm,pn,pv),(xi,yi,vi),(x,y,v),f,len,epsilon2;velocity=(ones(size(mask)),ones(size(mask)),zeros(size(mask)))) + + +@test fi[10,10,2] ≈ 0.15748780548566835 + +fi,s=DIVAnd_multivarEOF(mask,(pm,pn,pv),(xi,yi,vi),(x,y,v),f,(len,len,0.),epsilon2;velocity=(ones(size(mask)),ones(size(mask)),zeros(size(mask))),eof=[-1.,1.]) + +@test fi[10,10,2] ≈ 0.1177009988865455 + +fi,s=DIVAnd_multivarJAC(mask,(pm,pn,pv),(xi,yi,vi),(x,y,v),f,len,epsilon2;velocity=(ones(size(mask)),ones(size(mask)),zeros(size(mask)))) + + +@test fi[10,10,2] ≈ 0.051194689308468086 + +fi,s=DIVAnd_multivarJAC(mask,(pm,pn,pv),(xi,yi,vi),(x,y,v),f,(len,len,0.),epsilon2;velocity=(ones(size(mask)),ones(size(mask)),zeros(size(mask)))) + +@test fi[10,10,2] ≈ 0.051194689308468086 + + +# Copyright (C) 2014, 2017, 2018 +# Alexander Barth +# Jean-Marie Beckers +# +# This program is free software; you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the Free Software +# Foundation; either version 2 of the License, or (at your option) any later +# version. +# +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. +# +# You should have received a copy of the GNU General Public License along with +# this program; if not, see . From c17c5baa1964db090983873b87c3ebe521fe2e5b Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Tue, 28 Mar 2023 10:23:24 +0200 Subject: [PATCH 20/51] more tests --- src/DIVAnd_aexerr.jl | 5 ++- src/DIVAnd_errormap.jl | 23 ++++++++---- test/runtests.jl | 2 +- test/test_errormaps.jl | 79 ++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 99 insertions(+), 10 deletions(-) create mode 100644 test/test_errormaps.jl diff --git a/src/DIVAnd_aexerr.jl b/src/DIVAnd_aexerr.jl index 9231674d..7e6f2941 100644 --- a/src/DIVAnd_aexerr.jl +++ b/src/DIVAnd_aexerr.jl @@ -200,11 +200,10 @@ function DIVAnd_aexerr(mask, pmn, xi, x, f, len, epsilon2; otherargs...) # The factor 1.70677 is the best one in 2D but should be slightly different for other dimensions # Could be a small improvement. Also used in DIVAnd_cpme - f1, s1 = - DIVAndrun(mask, pmn, xi, xfake, ffake, len ./ 1.70766, epsilonforB; otherargs...) + f1, s1 = DIVAndrun(mask, pmn, xi, xfake, ffake, len ./ 1.70766, epsilonforB; otherargs...) # Calculate final error - aexerr = max.(Bjmb - f1, 0) + aexerr = max.(Bjmb - f1, 0.0) #@show mean(aexerr[.!isnan.(aexerr)]) diff --git a/src/DIVAnd_errormap.jl b/src/DIVAnd_errormap.jl index 01b402ca..27d463ba 100644 --- a/src/DIVAnd_errormap.jl +++ b/src/DIVAnd_errormap.jl @@ -201,33 +201,33 @@ function DIVAnd_errormap( if errmethod == :cpme && Bscale - warn("Sorry, that method does not allow rescaling by spatial dependance of B ") + @warn "Sorry, that method does not allow rescaling by spatial dependance of B " ScalebyB = false end if errmethod == :scpme && Bscale - warn("Sorry, that method does not allow rescaling by spatial dependance of B ") + @warn "Sorry, that method does not allow rescaling by spatial dependance of B " ScalebyB = false end if errmethod == :exact && Bscale # Or maybe if all info is there run locally ? Yes probably possible as aexerr also needs all infos ? - warn("You need to do that scaling by yourself, running diva again with a very high R matrix and divide by this second map") + @warn "You need to do that scaling by yourself, running diva again with a very high R matrix and divide by this second map" ScalebyB = false end if errmethod == :scpme && noP - warn("Sorry, that method needs s.P to be available. Will use cpme instead") + @warn "Sorry, that method needs s.P to be available. Will use cpme instead" errmethod = cpme end if errmethod == :exact && noP - warn("Sorry, that method needs s.P to be available. Will use aexerr instead") + @warn "Sorry, that method needs s.P to be available. Will use aexerr instead" errmethod = aexerr end if errmethod == :diagapp && noP - warn("Sorry, that method needs s.P to be available. Will use aexerr instead") + @warn "Sorry, that method needs s.P to be available. Will use aexerr instead" errmethod = aexerr end @@ -297,7 +297,18 @@ function DIVAnd_errormap( epsilon2; otherargs... ) + if errormap==0 + @warn "too fine resolution for aexerr, using exact" + errormap, =statevector_unpack(s.sv,diag(s.P)) + + return errormap, errmethod + + end + if ScalebyB + return errormap./bi, errmethod + else return errormap, errmethod + end end @show "You should not be here" end diff --git a/test/runtests.jl b/test/runtests.jl index c8a0ba62..45207559 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -35,7 +35,7 @@ using Downloads: download include("test_multivar.jl") - + include("test_errormaps.jl") # dynamical constraints include("test_2dvar_adv.jl") diff --git a/test/test_errormaps.jl b/test/test_errormaps.jl new file mode 100644 index 00000000..6829b2a6 --- /dev/null +++ b/test/test_errormaps.jl @@ -0,0 +1,79 @@ +# Testing DIVAnd in 2 dimensions with advection. +using DIVAnd +using Test +using StableRNGs +using Random + +# grid of background field +mask, (pm, pn), (xi, yi) = DIVAnd_squaredom(2, range(-1, stop = 1, length = 30)) +epsilon2 = 1 / 200 +errmean=0 + +a = 5; +u = a * yi; +v = -a * xi; + + +for ND in (10, 1000, 10000) +for len in (0.05, 0.2 ,1.0) +for mymeth in (:auto, :cheap, :precise, :cpme, :scpme, :exact, :aexerr, :diagapp) +for myscale in (true, false) + +xdd = rand(ND) +ydd = rand(ND) +fdd = randn(ND) + + + + + +fil, sl = DIVAndrun( + mask, + (pm, pn), + (xi, yi), + (xdd, ydd), + fdd, + len, + epsilon2; + velocity = (u, v), + alphabc = 0, +); + +errorm,methodc=DIVAnd_errormap( + mask, + (pm, pn), + (xi, yi), + (xdd, ydd), + fdd, + len, + epsilon2,sl; + velocity = (u, v), + alphabc = 0, + method=mymeth, + Bscale=myscale + ) + +global errmean=errmean+errorm[10,10] + +end +end +end +end + +@test errmean ≈ 6.951700275663623 + + +# Copyright (C) 2014, 2017 Alexander Barth +# +# This program is free software; you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the Free Software +# Foundation; either version 2 of the License, or (at your option) any later +# version. +# +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. +# +# You should have received a copy of the GNU General Public License along with +# this program; if not, see . From e8c186826b78a39ce5135cd7c2c75b2d1420b358 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Tue, 28 Mar 2023 16:50:54 +0200 Subject: [PATCH 21/51] Update DIVAndfun.jl --- src/DIVAndfun.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/DIVAndfun.jl b/src/DIVAndfun.jl index e361bd7d..fcc75b67 100644 --- a/src/DIVAndfun.jl +++ b/src/DIVAndfun.jl @@ -117,5 +117,6 @@ function DIVAndfun(x,f;mask=nothing,pmn=nothing,xi=nothing,len=nothing,epsilon2= # Now initialize interpolations function and return that - return LinearInterpolation(tuple(collect.(coords)...), fi.+backg) + #return LinearInterpolation(tuple(collect.(coords)...), fi.+backg) + return linear_interpolation(tuple(collect.(coords)...), fi.+backg) end \ No newline at end of file From c322399ce70c76e9e542ec0dc31e8001e00e9132 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Tue, 28 Mar 2023 20:20:51 +0200 Subject: [PATCH 22/51] Update test_errormaps.jl --- test/test_errormaps.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/test_errormaps.jl b/test/test_errormaps.jl index 6829b2a6..a659c236 100644 --- a/test/test_errormaps.jl +++ b/test/test_errormaps.jl @@ -4,6 +4,8 @@ using Test using StableRNGs using Random +rng = StableRNG(1234) + # grid of background field mask, (pm, pn), (xi, yi) = DIVAnd_squaredom(2, range(-1, stop = 1, length = 30)) epsilon2 = 1 / 200 @@ -59,7 +61,7 @@ end end end end - +@show errmean @test errmean ≈ 6.951700275663623 From d057bc9d3393b43fee3ec332648b0b811828ac00 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Wed, 29 Mar 2023 09:51:02 +0200 Subject: [PATCH 23/51] Update test_errormaps.jl --- test/test_errormaps.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/test_errormaps.jl b/test/test_errormaps.jl index a659c236..9680dd98 100644 --- a/test/test_errormaps.jl +++ b/test/test_errormaps.jl @@ -56,6 +56,7 @@ errorm,methodc=DIVAnd_errormap( ) global errmean=errmean+errorm[10,10] +@show errmean,errorm[10,10] end end From a888e847ea7c8d46dbee4043841c9146e31dc8cc Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Wed, 29 Mar 2023 10:47:17 +0200 Subject: [PATCH 24/51] Update test_errormaps.jl --- test/test_errormaps.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_errormaps.jl b/test/test_errormaps.jl index 9680dd98..d4fd29ff 100644 --- a/test/test_errormaps.jl +++ b/test/test_errormaps.jl @@ -26,7 +26,7 @@ ydd = rand(ND) fdd = randn(ND) - +@show xdd[1],ND,len,mymeth,myscale fil, sl = DIVAndrun( @@ -56,7 +56,7 @@ errorm,methodc=DIVAnd_errormap( ) global errmean=errmean+errorm[10,10] -@show errmean,errorm[10,10] +@show errmean,errorm[10,10],fi1[10,10],xdd[1],ND,len,mymeth,myscale,methodc end end From 6ff7793435be6e9fe3816defcb51f567acc0b65f Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Wed, 29 Mar 2023 15:08:39 +0200 Subject: [PATCH 25/51] Update test_errormaps.jl --- test/test_errormaps.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/test_errormaps.jl b/test/test_errormaps.jl index d4fd29ff..8f3bf849 100644 --- a/test/test_errormaps.jl +++ b/test/test_errormaps.jl @@ -26,7 +26,7 @@ ydd = rand(ND) fdd = randn(ND) -@show xdd[1],ND,len,mymeth,myscale +@show xdd[1],ND,len,mymeth,myscale,size(mask),size(pm) fil, sl = DIVAndrun( @@ -55,6 +55,8 @@ errorm,methodc=DIVAnd_errormap( Bscale=myscale ) + @show size(errorm),size(fi1) + global errmean=errmean+errorm[10,10] @show errmean,errorm[10,10],fi1[10,10],xdd[1],ND,len,mymeth,myscale,methodc From d9455258139bd4e06e7f8389d179882746bc2d95 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Wed, 29 Mar 2023 16:02:15 +0200 Subject: [PATCH 26/51] Update test_errormaps.jl --- test/test_errormaps.jl | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/test/test_errormaps.jl b/test/test_errormaps.jl index 8f3bf849..1a8b7666 100644 --- a/test/test_errormaps.jl +++ b/test/test_errormaps.jl @@ -26,10 +26,10 @@ ydd = rand(ND) fdd = randn(ND) -@show xdd[1],ND,len,mymeth,myscale,size(mask),size(pm) +@show xdd[1],ND,len,mymeth,myscale,size(mask),size(pm),size(xi),size(yi) -fil, sl = DIVAndrun( +fi11,sl = DIVAndrun( mask, (pm, pn), (xi, yi), @@ -38,8 +38,10 @@ fil, sl = DIVAndrun( len, epsilon2; velocity = (u, v), - alphabc = 0, -); + alphabc = 0 +) + +@show size(fi11) errorm,methodc=DIVAnd_errormap( mask, @@ -55,10 +57,10 @@ errorm,methodc=DIVAnd_errormap( Bscale=myscale ) - @show size(errorm),size(fi1) + @show size(errorm),size(fi11) global errmean=errmean+errorm[10,10] -@show errmean,errorm[10,10],fi1[10,10],xdd[1],ND,len,mymeth,myscale,methodc +@show errmean,errorm[10,10],fi11[10,10],xdd[1],ND,len,mymeth,myscale,methodc end end From 2a982b12e71edbbcd99a798b8f8ae343600eff10 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Thu, 30 Mar 2023 07:56:07 +0200 Subject: [PATCH 27/51] Update test_errormaps.jl --- test/test_errormaps.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_errormaps.jl b/test/test_errormaps.jl index 1a8b7666..d118f78c 100644 --- a/test/test_errormaps.jl +++ b/test/test_errormaps.jl @@ -11,7 +11,7 @@ mask, (pm, pn), (xi, yi) = DIVAnd_squaredom(2, range(-1, stop = 1, length = 30)) epsilon2 = 1 / 200 errmean=0 -a = 5; +a = 5.0 u = a * yi; v = -a * xi; From 1883c556d8332b451987a9b8f2bf3ac44bf9df88 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Thu, 30 Mar 2023 10:56:51 +0200 Subject: [PATCH 28/51] Update test_errormaps.jl --- test/test_errormaps.jl | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/test/test_errormaps.jl b/test/test_errormaps.jl index d118f78c..cd76c366 100644 --- a/test/test_errormaps.jl +++ b/test/test_errormaps.jl @@ -1,8 +1,9 @@ # Testing DIVAnd in 2 dimensions with advection. -using DIVAnd + using Test using StableRNGs using Random +using DIVAnd rng = StableRNG(1234) @@ -21,12 +22,12 @@ for len in (0.05, 0.2 ,1.0) for mymeth in (:auto, :cheap, :precise, :cpme, :scpme, :exact, :aexerr, :diagapp) for myscale in (true, false) -xdd = rand(ND) -ydd = rand(ND) -fdd = randn(ND) +xdd = rand(rng,ND) +ydd = rand(rng,ND) +fdd = randn(rng,ND) -@show xdd[1],ND,len,mymeth,myscale,size(mask),size(pm),size(xi),size(yi) +#show xdd[1],ND,len,mymeth,myscale,size(mask),size(pm),size(xi),size(yi) fi11,sl = DIVAndrun( @@ -41,7 +42,7 @@ fi11,sl = DIVAndrun( alphabc = 0 ) -@show size(fi11) + errorm,methodc=DIVAnd_errormap( mask, @@ -60,14 +61,14 @@ errorm,methodc=DIVAnd_errormap( @show size(errorm),size(fi11) global errmean=errmean+errorm[10,10] -@show errmean,errorm[10,10],fi11[10,10],xdd[1],ND,len,mymeth,myscale,methodc +#show errmean,errorm[10,10],fi11[10,10],xdd[1],ND,len,mymeth,myscale,methodc end end end end -@show errmean -@test errmean ≈ 6.951700275663623 +##how errmean +@test errmean ≈ 7.541790086372698 # Copyright (C) 2014, 2017 Alexander Barth From e9f12efedfffac318fb0941f473182336664454d Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Thu, 30 Mar 2023 11:35:48 +0200 Subject: [PATCH 29/51] Update test_errormaps.jl --- test/test_errormaps.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_errormaps.jl b/test/test_errormaps.jl index cd76c366..a4d98efb 100644 --- a/test/test_errormaps.jl +++ b/test/test_errormaps.jl @@ -27,7 +27,7 @@ ydd = rand(rng,ND) fdd = randn(rng,ND) -#show xdd[1],ND,len,mymeth,myscale,size(mask),size(pm),size(xi),size(yi) +@show xdd[1],ND,len,mymeth,myscale,size(mask),size(pm),size(xi),size(yi) fi11,sl = DIVAndrun( @@ -58,7 +58,7 @@ errorm,methodc=DIVAnd_errormap( Bscale=myscale ) - @show size(errorm),size(fi11) + #@show size(errorm),size(fi11) global errmean=errmean+errorm[10,10] #show errmean,errorm[10,10],fi11[10,10],xdd[1],ND,len,mymeth,myscale,methodc From 6c9bf962f52f609d5c8e9f3d63ce043a876af733 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Thu, 30 Mar 2023 12:03:17 +0200 Subject: [PATCH 30/51] Update test_errormaps.jl --- test/test_errormaps.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_errormaps.jl b/test/test_errormaps.jl index a4d98efb..426b0184 100644 --- a/test/test_errormaps.jl +++ b/test/test_errormaps.jl @@ -61,7 +61,7 @@ errorm,methodc=DIVAnd_errormap( #@show size(errorm),size(fi11) global errmean=errmean+errorm[10,10] -#show errmean,errorm[10,10],fi11[10,10],xdd[1],ND,len,mymeth,myscale,methodc +@show errmean,errorm[10,10],fi11[10,10],xdd[1],ND,len,mymeth,myscale,methodc end end From 9edcfad5cf610fa0e0d993997c221b429068cdd1 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Thu, 30 Mar 2023 12:33:17 +0200 Subject: [PATCH 31/51] Update DIVAnd_scalecpme!.jl --- src/DIVAnd_scalecpme!.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/DIVAnd_scalecpme!.jl b/src/DIVAnd_scalecpme!.jl index 92c25226..9201ef3e 100644 --- a/src/DIVAnd_scalecpme!.jl +++ b/src/DIVAnd_scalecpme!.jl @@ -5,7 +5,7 @@ function DIVAnd_scalecpme!(cpme, P::CovarIS, nsamples = 7) fractionshift=0.5 - z = randn((size(P)[1], nsamples)) + z = randn(rng,(size(P)[1], nsamples)) errscale = 1 if P.factors != nothing ZZ = P.factors.PtL \ z From f632d1591674295100ac8b865e1396b04cea35fc Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Thu, 30 Mar 2023 12:44:11 +0200 Subject: [PATCH 32/51] Update DIVAnd_scalecpme!.jl --- src/DIVAnd_scalecpme!.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/DIVAnd_scalecpme!.jl b/src/DIVAnd_scalecpme!.jl index 9201ef3e..c899649f 100644 --- a/src/DIVAnd_scalecpme!.jl +++ b/src/DIVAnd_scalecpme!.jl @@ -5,7 +5,7 @@ function DIVAnd_scalecpme!(cpme, P::CovarIS, nsamples = 7) fractionshift=0.5 - z = randn(rng,(size(P)[1], nsamples)) + z = randn(rng,((size(P)[1], nsamples))) errscale = 1 if P.factors != nothing ZZ = P.factors.PtL \ z From 67ca1700c26a30a02e34ccda6771d659ad7d08c8 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Thu, 30 Mar 2023 13:10:11 +0200 Subject: [PATCH 33/51] rng --- src/DIVAnd_scalecpme!.jl | 4 ++-- test/test_errormaps.jl | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/DIVAnd_scalecpme!.jl b/src/DIVAnd_scalecpme!.jl index c899649f..70139942 100644 --- a/src/DIVAnd_scalecpme!.jl +++ b/src/DIVAnd_scalecpme!.jl @@ -4,8 +4,8 @@ function DIVAnd_scalecpme!(cpme, P::CovarIS, nsamples = 7) # nsamples is the number of random arrays used to estimate the value fractionshift=0.5 - - z = randn(rng,((size(P)[1], nsamples))) + rng = Random.GLOBAL_RNG + z = randn(rng,(size(P)[1], nsamples)) errscale = 1 if P.factors != nothing ZZ = P.factors.PtL \ z diff --git a/test/test_errormaps.jl b/test/test_errormaps.jl index 426b0184..3b97505c 100644 --- a/test/test_errormaps.jl +++ b/test/test_errormaps.jl @@ -68,7 +68,7 @@ end end end ##how errmean -@test errmean ≈ 7.541790086372698 +@test errmean ≈ 7.525552547782063 # Copyright (C) 2014, 2017 Alexander Barth From 0d0065ad85ec17c7499280345357e73ec945b869 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Thu, 30 Mar 2023 14:18:26 +0200 Subject: [PATCH 34/51] rng --- src/DIVAnd_errormap.jl | 7 ++++--- src/DIVAnd_scalecpme!.jl | 4 ++-- test/test_errormaps.jl | 5 +++-- 3 files changed, 9 insertions(+), 7 deletions(-) diff --git a/src/DIVAnd_errormap.jl b/src/DIVAnd_errormap.jl index 27d463ba..e33cc7f7 100644 --- a/src/DIVAnd_errormap.jl +++ b/src/DIVAnd_errormap.jl @@ -62,7 +62,8 @@ function DIVAnd_errormap( s; method = :auto, Bscale = false, - otherargs..., + rng=Random.GLOBAL_RNG, + otherargs... ) # Criteria to define which fraction of the domain size L can be to be called small @@ -260,11 +261,11 @@ function DIVAnd_errormap( f, len, epsilon2; - otherargs... + otherargs... ) scpme=deepcopy(errormap) - DIVAnd_scalecpme!(scpme,s.P) + DIVAnd_scalecpme!(scpme,s.P;rng=rng) return scpme, errmethod end diff --git a/src/DIVAnd_scalecpme!.jl b/src/DIVAnd_scalecpme!.jl index 70139942..c6c8310f 100644 --- a/src/DIVAnd_scalecpme!.jl +++ b/src/DIVAnd_scalecpme!.jl @@ -1,10 +1,10 @@ -function DIVAnd_scalecpme!(cpme, P::CovarIS, nsamples = 7) +function DIVAnd_scalecpme!(cpme, P::CovarIS, nsamples = 7;rng=Random.GLOBAL_RNG) # IN PLACE rescaling of the clever poor mans estimate using a randomized estimate of the analysis error covariance # P returned in structure s (so s.P) from a previous run # nsamples is the number of random arrays used to estimate the value fractionshift=0.5 - rng = Random.GLOBAL_RNG + z = randn(rng,(size(P)[1], nsamples)) errscale = 1 if P.factors != nothing diff --git a/test/test_errormaps.jl b/test/test_errormaps.jl index 3b97505c..88745b2f 100644 --- a/test/test_errormaps.jl +++ b/test/test_errormaps.jl @@ -55,7 +55,8 @@ errorm,methodc=DIVAnd_errormap( velocity = (u, v), alphabc = 0, method=mymeth, - Bscale=myscale + Bscale=myscale, + rng=StableRNG(123) ) #@show size(errorm),size(fi11) @@ -68,7 +69,7 @@ end end end ##how errmean -@test errmean ≈ 7.525552547782063 +@test errmean ≈ 7.488115758433231 # Copyright (C) 2014, 2017 Alexander Barth From 694fff40d9fe45fc86fce6d7835b978c48997044 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Thu, 30 Mar 2023 15:02:39 +0200 Subject: [PATCH 35/51] rng again --- src/DIVAnd_aexerr.jl | 4 ++-- src/DIVAnd_errormap.jl | 1 + test/test_errormaps.jl | 2 +- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/DIVAnd_aexerr.jl b/src/DIVAnd_aexerr.jl index 7e6f2941..5a3408b7 100644 --- a/src/DIVAnd_aexerr.jl +++ b/src/DIVAnd_aexerr.jl @@ -40,7 +40,7 @@ Compute a variational analysis of arbitrarily located observations to calculate the almost exact error """ -function DIVAnd_aexerr(mask, pmn, xi, x, f, len, epsilon2; otherargs...) +function DIVAnd_aexerr(mask, pmn, xi, x, f, len, epsilon2; rng=Random.GLOBAL_RNG, otherargs...) @@ -114,7 +114,7 @@ function DIVAnd_aexerr(mask, pmn, xi, x, f, len, epsilon2; otherargs...) # not sure that covers nicely the domain?? Check random idea? # randindexes = collect(1:nsa:npgrid) - randindexes = shuffle(collect(1:npgrid))[1:npongrid] + randindexes = shuffle(rng,collect(1:npgrid))[1:npongrid] ncv = size(randindexes)[1] diff --git a/src/DIVAnd_errormap.jl b/src/DIVAnd_errormap.jl index e33cc7f7..506f8f1b 100644 --- a/src/DIVAnd_errormap.jl +++ b/src/DIVAnd_errormap.jl @@ -296,6 +296,7 @@ function DIVAnd_errormap( f, len, epsilon2; + rng=rng, otherargs... ) if errormap==0 diff --git a/test/test_errormaps.jl b/test/test_errormaps.jl index 88745b2f..e6dd55f8 100644 --- a/test/test_errormaps.jl +++ b/test/test_errormaps.jl @@ -69,7 +69,7 @@ end end end ##how errmean -@test errmean ≈ 7.488115758433231 +@test errmean ≈ 7.524924869150642 # Copyright (C) 2014, 2017 Alexander Barth From 000f58aa2c30c9a3d9a9a799cecf08a1a2a64f98 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Thu, 30 Mar 2023 15:44:23 +0200 Subject: [PATCH 36/51] rng --- src/DIVAnd_aexerr.jl | 2 +- test/test_errormaps.jl | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/DIVAnd_aexerr.jl b/src/DIVAnd_aexerr.jl index 5a3408b7..e1ac93c6 100644 --- a/src/DIVAnd_aexerr.jl +++ b/src/DIVAnd_aexerr.jl @@ -147,7 +147,7 @@ function DIVAnd_aexerr(mask, pmn, xi, x, f, len, epsilon2; rng=Random.GLOBAL_RNG restrictedlist = falses(size(ffake)[1]) restrictedlist[size(f)[1]+1:end] .= true # Limit the number of data points to the number of additional fake points - samples = shuffle(collect(1:size(f)[1]))[1:min(size(f)[1], npongrid)] + samples = shuffle(rng,collect(1:size(f)[1]))[1:min(size(f)[1], npongrid)] restrictedlist[samples] .= true #restrictedlist[:].=true diff --git a/test/test_errormaps.jl b/test/test_errormaps.jl index e6dd55f8..f1be2e7b 100644 --- a/test/test_errormaps.jl +++ b/test/test_errormaps.jl @@ -27,7 +27,7 @@ ydd = rand(rng,ND) fdd = randn(rng,ND) -@show xdd[1],ND,len,mymeth,myscale,size(mask),size(pm),size(xi),size(yi) +#@show xdd[1],ND,len,mymeth,myscale,size(mask),size(pm),size(xi),size(yi) fi11,sl = DIVAndrun( @@ -62,14 +62,14 @@ errorm,methodc=DIVAnd_errormap( #@show size(errorm),size(fi11) global errmean=errmean+errorm[10,10] -@show errmean,errorm[10,10],fi11[10,10],xdd[1],ND,len,mymeth,myscale,methodc +#@show errmean,errorm[10,10],fi11[10,10],xdd[1],ND,len,mymeth,myscale,methodc end end end end -##how errmean -@test errmean ≈ 7.524924869150642 +@show errmean +@test abs(errmean - 7.527185124068747) < 0.00000001 # Copyright (C) 2014, 2017 Alexander Barth From 6e66b86d3316708a073bcb51dfc5f2dccc648f19 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Fri, 31 Mar 2023 09:51:23 +0200 Subject: [PATCH 37/51] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index a42a6d02..c7fc4381 100644 --- a/README.md +++ b/README.md @@ -5,7 +5,7 @@ [![Project Status: Active – The project has reached a stable, usable state and is being actively developed.](https://www.repostatus.org/badges/latest/active.svg)](https://www.repostatus.org/#active) [![Build Status](https://github.com/gher-uliege/DIVAnd.jl/workflows/CI/badge.svg)](https://github.com/gher-uliege/DIVAnd.jl/actions) -[![codecov.io](http://codecov.io/github/gher-uliege/DIVAnd.jl/coverage.svg?branch=master)](http://codecov.io/github/gher-uliege/DIVAnd.jl?branch=master) +[![codecov.io](http://codecov.io/github/gher-uliege/DIVAnd.jl/coverage.svg?branch=master)](http://codecov.io/github/gher-uliege/DIVAnd.jl?branch=JMB) [![documentation stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://gher-uliege.github.io/DIVAnd.jl/stable/) [![documentation latest](https://img.shields.io/badge/docs-latest-blue.svg)](https://gher-uliege.github.io/DIVAnd.jl/latest/) [![DOI](https://zenodo.org/badge/79277337.svg)](https://zenodo.org/badge/latestdoi/79277337) From fb21aaa3377d831ff70e224ff46d172e0608b998 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Fri, 31 Mar 2023 12:29:24 +0200 Subject: [PATCH 38/51] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index c7fc4381..986e95f2 100644 --- a/README.md +++ b/README.md @@ -5,7 +5,7 @@ [![Project Status: Active – The project has reached a stable, usable state and is being actively developed.](https://www.repostatus.org/badges/latest/active.svg)](https://www.repostatus.org/#active) [![Build Status](https://github.com/gher-uliege/DIVAnd.jl/workflows/CI/badge.svg)](https://github.com/gher-uliege/DIVAnd.jl/actions) -[![codecov.io](http://codecov.io/github/gher-uliege/DIVAnd.jl/coverage.svg?branch=master)](http://codecov.io/github/gher-uliege/DIVAnd.jl?branch=JMB) +[![codecov.io](http://codecov.io/github/gher-uliege/DIVAnd.jl/coverage.svg?branch=JMB)](http://codecov.io/github/gher-uliege/DIVAnd.jl?branch=JMB) [![documentation stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://gher-uliege.github.io/DIVAnd.jl/stable/) [![documentation latest](https://img.shields.io/badge/docs-latest-blue.svg)](https://gher-uliege.github.io/DIVAnd.jl/latest/) [![DOI](https://zenodo.org/badge/79277337.svg)](https://zenodo.org/badge/latestdoi/79277337) From bf74114fd1e8bb51759f65ea799c946c6c628054 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Fri, 31 Mar 2023 12:31:39 +0200 Subject: [PATCH 39/51] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 986e95f2..e1d94ae7 100644 --- a/README.md +++ b/README.md @@ -5,7 +5,7 @@ [![Project Status: Active – The project has reached a stable, usable state and is being actively developed.](https://www.repostatus.org/badges/latest/active.svg)](https://www.repostatus.org/#active) [![Build Status](https://github.com/gher-uliege/DIVAnd.jl/workflows/CI/badge.svg)](https://github.com/gher-uliege/DIVAnd.jl/actions) -[![codecov.io](http://codecov.io/github/gher-uliege/DIVAnd.jl/coverage.svg?branch=JMB)](http://codecov.io/github/gher-uliege/DIVAnd.jl?branch=JMB) +[![codecov.io](http://codecov.io/github/gher-uliege/DIVAnd.jl/coverage.svg?branch=JMB)](http://app.codecov.io/github/gher-uliege/DIVAnd.jl?tree=JMB) [![documentation stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://gher-uliege.github.io/DIVAnd.jl/stable/) [![documentation latest](https://img.shields.io/badge/docs-latest-blue.svg)](https://gher-uliege.github.io/DIVAnd.jl/latest/) [![DOI](https://zenodo.org/badge/79277337.svg)](https://zenodo.org/badge/latestdoi/79277337) From 06205774c225a8e430b20a8e1d9073779c72d8ab Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Fri, 31 Mar 2023 12:32:28 +0200 Subject: [PATCH 40/51] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index e1d94ae7..e8534455 100644 --- a/README.md +++ b/README.md @@ -5,7 +5,7 @@ [![Project Status: Active – The project has reached a stable, usable state and is being actively developed.](https://www.repostatus.org/badges/latest/active.svg)](https://www.repostatus.org/#active) [![Build Status](https://github.com/gher-uliege/DIVAnd.jl/workflows/CI/badge.svg)](https://github.com/gher-uliege/DIVAnd.jl/actions) -[![codecov.io](http://codecov.io/github/gher-uliege/DIVAnd.jl/coverage.svg?branch=JMB)](http://app.codecov.io/github/gher-uliege/DIVAnd.jl?tree=JMB) +[![codecov.io](http://codecov.io/github/gher-uliege/DIVAnd.jl/coverage.svg?branch=JMB)](https://app.codecov.io/github/gher-uliege/DIVAnd.jl/tree/JMB) [![documentation stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://gher-uliege.github.io/DIVAnd.jl/stable/) [![documentation latest](https://img.shields.io/badge/docs-latest-blue.svg)](https://gher-uliege.github.io/DIVAnd.jl/latest/) [![DOI](https://zenodo.org/badge/79277337.svg)](https://zenodo.org/badge/latestdoi/79277337) From 2bb88c8194f47cea3c2d3d5a171e6b7e17179475 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Fri, 31 Mar 2023 12:35:34 +0200 Subject: [PATCH 41/51] Update README.md --- README.md | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index e8534455..736d3dcb 100644 --- a/README.md +++ b/README.md @@ -169,7 +169,7 @@ Tools to help you are included in ([DIVAnd_cv.jl](https://github.com/gher-ulieg ### Additional constraint -An arbitrary number of additional constraints can be included to the cost function which should have the following form: +An arbitrary number of additional quadratic constraints can be included to the cost function which should have the following form: *J*(**x**) = ∑*i* (**C***i* **x** - **z***i*)ᵀ **Q***i*-1 (**C***i* **x** - **z***i*) @@ -181,6 +181,19 @@ For every constrain, a structure with the following fields is passed to `DIVAnd` Internally the observations are also implemented as constraint defined in this way. +### Additional inequaliuty constraint + +An arbitrary number of additional quadratic constraints can be included to the cost function which should have the following form: + +(**H***i* **x** > **yo***i*) + +For every constrain, a structure with the following fields is passed to `DIVAnd`: + +* `yo`: a vector +* `H`: a matrix + + + ## Run notebooks on a server which has no graphical interface On the server, launch the notebook with: From 490a6e793dfc998e8820782602cd802e0a37d3bc Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Fri, 31 Mar 2023 12:37:15 +0200 Subject: [PATCH 42/51] Update README.md --- README.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index 736d3dcb..bae4210f 100644 --- a/README.md +++ b/README.md @@ -181,13 +181,13 @@ For every constrain, a structure with the following fields is passed to `DIVAnd` Internally the observations are also implemented as constraint defined in this way. -### Additional inequaliuty constraint +### Additional inequality constraint -An arbitrary number of additional quadratic constraints can be included to the cost function which should have the following form: +An arbitrary number of additional inequality constraints can be included and which should have the following form: (**H***i* **x** > **yo***i*) -For every constrain, a structure with the following fields is passed to `DIVAnd`: +For every constraint, a structure with the following fields is passed to `DIVAnd`: * `yo`: a vector * `H`: a matrix From 2e77c6ad2af202714c08db108031b2d40bffe6f3 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Sun, 2 Apr 2023 10:47:11 +0200 Subject: [PATCH 43/51] Update README.md --- README.md | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/README.md b/README.md index 74f4f060..156030f7 100644 --- a/README.md +++ b/README.md @@ -24,6 +24,19 @@ Barth, A., Beckers, J.-M., Troupin, C., Alvera-Azcárate, A., and Vandenbulcke, (click [here](./data/DIVAnd.bib) for the BibTeX entry). +# Summary of features + +* N-Dimensional analysis/interpolation +* Scattered data +* Noise allowed +* Physical constraints can be added +* Inequality constraints can be added +* Topological constraints are handled naturally (barriers, holes) +* Analysis error maps can be estimated +* Periodicity in selected directions can be enforced +* Multivariate data can be used (experimental) +* The output grid can be curvilinear + # Installing From 9a65aca8a0c403852187ff6ac9da1ecd410efc29 Mon Sep 17 00:00:00 2001 From: Alexander Barth Date: Mon, 3 Apr 2023 13:01:17 +0200 Subject: [PATCH 44/51] doc formatting --- src/DIVAnd_background.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/DIVAnd_background.jl b/src/DIVAnd_background.jl index f2a89392..b6262e38 100644 --- a/src/DIVAnd_background.jl +++ b/src/DIVAnd_background.jl @@ -1,6 +1,6 @@ """ -Form the inverse of the background error covariance matrix. -s = DIVAnd_background(mask,pmn,Labs,alpha,moddim) + s = DIVAnd_background(mask,pmn,Labs,alpha,moddim) + Form the inverse of the background error covariance matrix with finite-difference operators on a curvilinear grid # Input: @@ -65,6 +65,8 @@ function DIVAnd_background( end end + @debug "alpha" alpha + @debug "scaling" len_scale # mean correlation length in every dimension Ld = if mean_Labs == nothing From 34cc027f606b96da0dc37cea75517524ec23192f Mon Sep 17 00:00:00 2001 From: Alexander Barth Date: Mon, 3 Apr 2023 13:37:03 +0200 Subject: [PATCH 45/51] remove trailing white-space --- examples/DIVAnd_SST_pluto.jl | 4 +- src/DIVAnd_bc_stretch.jl | 2 +- src/DIVAnd_diagnostics.jl | 47 +++++++------- src/DIVAnd_errormap.jl | 14 ++-- src/DIVAnd_multivarJAC.jl | 122 +++++++++++++++++------------------ src/DIVAnd_scalecpme!.jl | 2 +- src/DIVAndfun.jl | 76 +++++++++++----------- src/DIVAndrun.jl | 30 ++++----- src/derived.jl | 2 +- src/fit.jl | 6 +- test/runtests.jl | 8 +-- test/test_errormaps.jl | 4 +- 12 files changed, 158 insertions(+), 159 deletions(-) diff --git a/examples/DIVAnd_SST_pluto.jl b/examples/DIVAnd_SST_pluto.jl index 739a115f..cd66f185 100644 --- a/examples/DIVAnd_SST_pluto.jl +++ b/examples/DIVAnd_SST_pluto.jl @@ -106,7 +106,7 @@ begin ) continents = ifelse.(isfinite.(v),NaN,0) - + function plotmap(v,title;clim = extrema(filter(isfinite,v)), kwargs...) heatmap(lon,lat,continents',c = palette([:grey, :grey], 2)) heatmap!(lon,lat,v'; title=title,clim=clim, kwargs...) @@ -120,7 +120,7 @@ begin edgecolor = :none) end clim = extrema(filter(isfinite,v)) - + plot( plotmap(v,"True field", clim = clim), plotmap(vi,"Analysis field", clim = clim), diff --git a/src/DIVAnd_bc_stretch.jl b/src/DIVAnd_bc_stretch.jl index c0aafedf..ce3b3564 100644 --- a/src/DIVAnd_bc_stretch.jl +++ b/src/DIVAnd_bc_stretch.jl @@ -1,5 +1,5 @@ function DIVAnd_bc_stretch(mask, pmnin, xiin, Lin, moddim, alphabc::Number = 1.0) - + return DIVAnd_bc_stretch(mask, pmnin, xiin, Lin, moddim, alphabc*ones(ndims(mask))) end diff --git a/src/DIVAnd_diagnostics.jl b/src/DIVAnd_diagnostics.jl index faf64544..35fb1a64 100644 --- a/src/DIVAnd_diagnostics.jl +++ b/src/DIVAnd_diagnostics.jl @@ -1,11 +1,11 @@ """ norm1,norm2,norm3,epsilon2=DIVAnd_norms(fi,s) -### Input: +### Input: * `fi` and `s` from the corresponding `DIVANDrun` exectution `fi,s = DIVAndrun(...` -### Output: +### Output: the different terms of the norm which is minimised @@ -21,37 +21,36 @@ the different terms of the norm which is minimised """ function DIVAnd_norms(fi,s) - + R = s.R iB = s.iB - H = s.H + H = s.H fstate = statevector_pack(s.sv, (fi,)) residuals=s.yo - (s.H) * fstate residualsobs=s.obsconstrain.yo - (s.obsconstrain.H) *fstate residualsobs[s.obsout] .= 0.0 - + norm2= residualsobs'*(s.obsconstrain.R\residualsobs) norm3= residuals'*(R \ residuals)-norm2 norm1=fstate'*(iB*fstate) dR=diag(s.obsconstrain.R) epsilon2=mean(dR[isfinite.(dR)]) - - return norm1,norm2,norm3,epsilon2 - + + return norm1,norm2,norm3,epsilon2 + end - - - - - - - - - - - - - - - - \ No newline at end of file + + + + + + + + + + + + + + + diff --git a/src/DIVAnd_errormap.jl b/src/DIVAnd_errormap.jl index 506f8f1b..c83fc185 100644 --- a/src/DIVAnd_errormap.jl +++ b/src/DIVAnd_errormap.jl @@ -28,13 +28,13 @@ * `epsilon2`: error variance of the observations (normalized by the error variance of the background field). `epsilon2` can be a scalar (all observations have the same error variance and their errors are decorrelated), a vector (all observations can have a difference error variance and their errors are decorrelated) or a matrix (all observations can have a difference error variance and their errors can be correlated). If `epsilon2` is a scalar, it is thus the *inverse of the signal-to-noise ratio*. -* `s`: this is the structure returned from the analysis itself. +* `s`: this is the structure returned from the analysis itself. # Optional input arguments specified as keyword arguments also as for DIVAnd *`method` : the method to be used, valid are `:auto`, `:cheap`, `:precise`, `:cpme`, `:scpme`, `:exact`, `:aexerr`, `:diagapp` - auto will select depenting on data coverage and length scale + auto will select depenting on data coverage and length scale cheap will also select but restrict to cheaper methods precise will also select but prefer better approximations the other choices are just the ones among which the automatic choices will choose from. You can force the choice by specifying the method. @@ -120,7 +120,7 @@ function DIVAnd_errormap( if method == :auto - + # try to guess # small L @@ -156,7 +156,7 @@ function DIVAnd_errormap( if method == :cheap - + if smallL if Lowdata errmethod = :cpme @@ -178,7 +178,7 @@ function DIVAnd_errormap( end if method == :precise - + if smallL @@ -271,7 +271,7 @@ function DIVAnd_errormap( end if errmethod == :exact - + errormap, =statevector_unpack(s.sv,diag(s.P)) return errormap, errmethod @@ -304,7 +304,7 @@ function DIVAnd_errormap( errormap, =statevector_unpack(s.sv,diag(s.P)) return errormap, errmethod - + end if ScalebyB return errormap./bi, errmethod diff --git a/src/DIVAnd_multivarJAC.jl b/src/DIVAnd_multivarJAC.jl index 099ba282..ffd1578c 100644 --- a/src/DIVAnd_multivarJAC.jl +++ b/src/DIVAnd_multivarJAC.jl @@ -13,11 +13,11 @@ DIVAnd_multivarJAC(mask,pmn,xi,x,f,lenin,epsilon2in;epsilon2jacobian=1.0,kwargs. if you have a habitat variable which you want to influence a variable for which you have few data only. Then you assign a large epsilon2forjacobian to the habitat variable and a lower value to the other one. -# Output: +# Output: * `fi`: analysis -* `s` : structure +* `s` : structure * `emap` : univariate error field @@ -25,9 +25,9 @@ DIVAnd_multivarJAC(mask,pmn,xi,x,f,lenin,epsilon2in;epsilon2jacobian=1.0,kwargs. * `pseudovelocity` : the "advection" field that guided the analysis - - -# + + +# """ @@ -36,37 +36,37 @@ DIVAnd_multivarJAC(mask,pmn,xi,x,f,lenin,epsilon2in;epsilon2jacobian=1.0,kwargs. if ndims(mask)<3 @error "Need two dimensional fields and several variables" end - + epsilon2=deepcopy(epsilon2in) sz = size(mask) lxl=zeros(Float64,size(mask)[end]) lyl=zeros(Float64,size(mask)[end]) - + eps2jac=epsilon2jacobian if isa(eps2jac, Number) eps2jac=tuple(repeat([eps2jac],size(mask)[end])...) end @show eps2jac - + # Some saveguards here; len must be zero on last dimension # At the same time create the lenlow correlation length for the analysis # in the lower dimension (without "variable" dimension) for the eof amplitudes # len=deepcopy(lenin) - + variableL=false if isa(lenin, Number) @warn "Expecting len array with zero in last direction; will be enforced" - + len=tuple(repeat([lenin],ndims(mask)-1)...,0.0) @debug "Len enforced $len" elseif isa(lenin, Tuple) if isa(lenin[1], Number) - - + + if lenin[end]!=0.0 - @warn "Expecting zero len in last direction; enforcing" + @warn "Expecting zero len in last direction; enforcing" #@show lenin[end] #@show lenin[1:end-1] len=tuple(lenin[1:end-1]...,0.0) @@ -75,72 +75,72 @@ DIVAnd_multivarJAC(mask,pmn,xi,x,f,lenin,epsilon2in;epsilon2jacobian=1.0,kwargs. else variableL=true if sum(len[end])!=0.0 - @warn "Expecting zero len in last direction , enforcing" + @warn "Expecting zero len in last direction , enforcing" len=(len[1:end-1]...,0.0 .* len[1]) end end end - - - - - - # univariate analysis - fi,s=DIVAndrun(mask,pmn,xi,x,f,len,epsilon2;kwargs...) - # - - - - - emap,meth=DIVAnd_errormap(mask,pmn,xi,x,f,len,epsilon2,s;method=:cheap,kwargs...) + + + + + + # univariate analysis + fi,s=DIVAndrun(mask,pmn,xi,x,f,len,epsilon2;kwargs...) + # + + + + + emap,meth=DIVAnd_errormap(mask,pmn,xi,x,f,len,epsilon2,s;method=:cheap,kwargs...) sori=deepcopy(s) @show "error method in multivar $meth" - + vconstrain=[] ML=5 pseudovelocity=[zeros(Float64,size(mask)) for idx in 1:ndims(mask)] for mloop=1:ML # # calculate pseudo-velocities. Maybe best to do it with already - # existing operators + # existing operators # dfdy - - - + + + S = sparse_stagger(sz, 2, s.iscyclic[2]) m = (S * mask[:]) .== 1 (coucouc,) = statevector_unpack(s.sv, sparse_pack(mask) * S' * sparse_pack(m)' *s.Dx[2]*statevector_pack(s.sv, (fi,))) pseudovelocity[1] .= coucouc - + S = sparse_stagger(sz, 1, s.iscyclic[1]) m = (S * mask[:]) .== 1 (coucouc,) = statevector_unpack(s.sv, sparse_pack(mask) * S' * sparse_pack(m)' *s.Dx[1]*statevector_pack(s.sv, (fi,))) pseudovelocity[2].= -coucouc - + for jjj=3:ndims(mask) pseudovelocity[jjj].=zeros(Float64,size(mask)) end # now scale each "2D layer" of the velocity () - + uu=reshape(pseudovelocity[1],(prod(size(pseudovelocity[1])[1:2]),prod(size(pseudovelocity[1])[3:end]))) vv=reshape(pseudovelocity[2],(prod(size(pseudovelocity[2])[1:2]),prod(size(pseudovelocity[2])[3:end]))) - + for kk=1:size(uu)[2] vn=sqrt(sum(uu[:,kk].^2 .+ vv[:,kk].^2)/size(uu)[1]) if vn>0 uu[:,kk].=uu[:,kk]./vn - vv[:,kk].=vv[:,kk]./vn + vv[:,kk].=vv[:,kk]./vn end end - - - + + + uu=reshape(uu,(prod(size(pseudovelocity[1])[1:end-1]),prod(size(pseudovelocity[1])[end]))) vv=reshape(vv,(prod(size(pseudovelocity[2])[1:end-1]),prod(size(pseudovelocity[2])[end]))) lu=zeros(size(uu)[2]) lv=zeros(size(uu)[2]) - wwww=ones(size(uu)[2]) + wwww=ones(size(uu)[2]) wtot=0 for lll=1:size(uu)[2] wwww[lll]=sum(abs.(x[end].-lll).<0.4) @@ -148,8 +148,8 @@ DIVAnd_multivarJAC(mask,pmn,xi,x,f,lenin,epsilon2in;epsilon2jacobian=1.0,kwargs. end wwww.=wwww./wtot #@show wwww - - + + #now for each physical point for ii=1:size(uu)[1] lu.=uu[ii,:] @@ -157,10 +157,10 @@ DIVAnd_multivarJAC(mask,pmn,xi,x,f,lenin,epsilon2in;epsilon2jacobian=1.0,kwargs. vvv=lu.^2 .+ lv.^2 wwww.=1.0 ./max.(0.000001,emap[ii:size(uu)[1]:end]) wwww.=wwww/sum(wwww) - - + + sorti=sortperm(vvv;rev=true) - + meanu=lu[sorti[1]]*wwww[sorti[1]] meanv=lv[sorti[1]]*wwww[sorti[1]] for kk=2:size(uu)[2] @@ -168,14 +168,14 @@ DIVAnd_multivarJAC(mask,pmn,xi,x,f,lenin,epsilon2in;epsilon2jacobian=1.0,kwargs. if meanu*lu[kks]+meanv*lv[kks]> 0 meanu=meanu+lu[kks]*wwww[kks] meanv=meanv+lv[kks]*wwww[kks] - + else meanu=meanu-lu[kks]*wwww[kks] meanv=meanv-lv[kks]*wwww[kks] - + end end - # apply possible weighting and scaling by L2 here + # apply possible weighting and scaling by L2 here if variableL lxl[:].= len[1][ii:size(uu)[1]:end] lyl[:].= len[2][ii:size(uu)[1]:end] @@ -183,12 +183,12 @@ DIVAnd_multivarJAC(mask,pmn,xi,x,f,lenin,epsilon2in;epsilon2jacobian=1.0,kwargs. lxl[:].= len[1] lyl[:].= len[2] end - - - + + + uu[ii,:].=meanu.*sqrt.(lxl[:].*lyl[:]./eps2jac[:]) vv[ii,:].=meanv.*sqrt.(lxl[:].*lyl[:]./eps2jac[:]) - + end uu[isnan.(uu)].=0.0 vv[isnan.(vv)].=0.0 @@ -196,17 +196,17 @@ DIVAnd_multivarJAC(mask,pmn,xi,x,f,lenin,epsilon2in;epsilon2jacobian=1.0,kwargs. pseudovelocity[2]=reshape(vv,size(pseudovelocity[2])) # vconstrain=DIVAnd_constr_advec(sori, tuple(pseudovelocity...)) - # pass the corresponding constraint to divandrun + # pass the corresponding constraint to divandrun fi,s=DIVAndrun(mask,pmn,xi,x,f,len,epsilon2;constraints=[vconstrain],kwargs...) - - - - - - end + + + + + + end # Error maps for the multivariate approach emapm,methm=DIVAnd_errormap(mask,pmn,xi,x,f,len,epsilon2,s;method=:cheap,constraints=[vconstrain],kwargs...) @show methm - # Banzaii, finished + # Banzaii, finished return fi,s,emap,emapm,pseudovelocity end \ No newline at end of file diff --git a/src/DIVAnd_scalecpme!.jl b/src/DIVAnd_scalecpme!.jl index c6c8310f..547f2477 100644 --- a/src/DIVAnd_scalecpme!.jl +++ b/src/DIVAnd_scalecpme!.jl @@ -4,7 +4,7 @@ function DIVAnd_scalecpme!(cpme, P::CovarIS, nsamples = 7;rng=Random.GLOBAL_RNG) # nsamples is the number of random arrays used to estimate the value fractionshift=0.5 - + z = randn(rng,(size(P)[1], nsamples)) errscale = 1 if P.factors != nothing diff --git a/src/DIVAndfun.jl b/src/DIVAndfun.jl index fcc75b67..5eace12c 100644 --- a/src/DIVAndfun.jl +++ b/src/DIVAndfun.jl @@ -1,12 +1,12 @@ - - - + + + """ myfunny=DIVAndfun(x,f;mask=nothing,pmn=nothing,xi=nothing,len=nothing,epsilon2=nothing,kwargs...) - + Provides a simplified interface to DIVAndrun and in return provides an interpolation FUNCTION rather than the gridded field of DIVAndrun You can use ALL paramters you can use with DIVAndrun, but some of them are made optional here by using keywords. @@ -16,7 +16,7 @@ The necessary input is the tuple of coordinates at which you have data points an The output is an interpolation function you can call at any coordinate in a hypercube defined by the bounding box of your input data or the domain explicitely defined by keyword xi If the user want more control on the grid he needs at least to provide xi, here with option to just provide vectors of coordinates in each direction (only works anyway for this case) so xi can be a tuple of vectors - + You can use all keyword parameters of divand # Input: @@ -27,43 +27,43 @@ You can use all keyword parameters of divand # Output: * `myfunny`: the interpolation function. If you had two dimensional input (i.e. x was a tuple of two coordinates), you can evaluate the interpolation as myfunny(0.1,0.2) for example -""" -function DIVAndfun(x,f;mask=nothing,pmn=nothing,xi=nothing,len=nothing,epsilon2=nothing,kwargs...) - - - - - - - +""" +function DIVAndfun(x,f;mask=nothing,pmn=nothing,xi=nothing,len=nothing,epsilon2=nothing,kwargs...) + + + + + + + maxsize=[1000,100,30,20] maxs=10 if length(x)>4 @warn "Are you sure ?" - + else maxs=maxsize[length(x)] end - + if size(f)[1]!=size(x[1])[1] @error "Need same size for coordinates than values" end - # Determine grid resolution + # Determine grid resolution nd=size(f)[1] - - - + + + NX=min(5*Int(round(10^(log10(nd)/length(x)))),maxs)*ones(Int32,length(x)) - + # Allocate array LX=ones(Float64,length(x)) - - - - + + + + coords=[] if xi==nothing - + for i=1:length(x) ri=[extrema(x[i])...] LX[i]=ri[2]-ri[1] @@ -75,48 +75,48 @@ function DIVAndfun(x,f;mask=nothing,pmn=nothing,xi=nothing,len=nothing,epsilon2= end mask,pmn,xi = DIVAnd_rectdom(coords...) else - + # No check done just taking the first point... for i=1:length(x) if isa(xi[i],Vector) # Get the coordinates along that direction - coorxi=deepcopy(xi[i]) + coorxi=deepcopy(xi[i]) else # Extract grid position in each direction assuming separable coordinates. coorxi=deepcopy(xi[i][1:stride(xi[i],i):stride(xi[i],i)*size(xi[i])[i]]) end coords=[coords...,coorxi] end - + if isa(xi[1],Vector) mask,pmn,xi = DIVAnd_rectdom(coords...) end - - + + if pmn==nothing # metric (inverse of the resolution) pmn = ndgrid([1 ./ localresolution(coords[i]) for i = 1:length(coords)]...) end - - + + if mask==nothing mask=trues(size(xi[1])) - end + end end - + if len==nothing len=tuple(LX./NX .* 4.0 ...) end if epsilon2==nothing epsilon2=0.1 end - + # Now make the DIVAndrun call backg=sum(f)/size(f)[1] fi,s=DIVAndrun(mask,pmn,xi,x,f.-backg,len,epsilon2; kwargs...) - + # Now initialize interpolations function and return that - + #return LinearInterpolation(tuple(collect.(coords)...), fi.+backg) return linear_interpolation(tuple(collect.(coords)...), fi.+backg) end \ No newline at end of file diff --git a/src/DIVAndrun.jl b/src/DIVAndrun.jl index a30d1158..f5ecd90d 100644 --- a/src/DIVAndrun.jl +++ b/src/DIVAndrun.jl @@ -43,7 +43,7 @@ function DIVAndrun( # Inequality constraints via loop around classical analysis (recursive call) - + if !isempty(ineqconstraints) @@ -55,25 +55,25 @@ function Ineqtoeq(yo,H,xo;tol1=1e-7,tol2=1e-8) Hxo=H*xo @show sum(Hxo. Date: Mon, 3 Apr 2023 13:38:41 +0200 Subject: [PATCH 46/51] expand tabs to space --- examples/DIVAnd_SST_pluto.jl | 14 +++--- examples/DIVAnd_simple_test_fluxes.jl | 24 +++++----- src/DIVAnd.jl | 2 +- src/DIVAnd_aexerr.jl | 2 +- src/DIVAnd_diagnostics.jl | 16 +++---- src/DIVAnd_errormap.jl | 42 +++++++++--------- src/DIVAnd_heatmap.jl | 4 +- src/DIVAnd_multivarEOF.jl | 4 +- src/DIVAnd_multivarJAC.jl | 6 +-- src/DIVAnd_scalecpme!.jl | 16 +++---- src/DIVAndfun.jl | 4 +- src/DIVAndrun.jl | 64 +++++++++++++-------------- src/fit.jl | 4 +- src/special_matrices.jl | 2 +- test/runtests.jl | 10 ++--- test/test_2dvar_adv.jl | 2 +- test/test_errormaps.jl | 10 ++--- 17 files changed, 113 insertions(+), 113 deletions(-) diff --git a/examples/DIVAnd_SST_pluto.jl b/examples/DIVAnd_SST_pluto.jl index cd66f185..8fdabca8 100644 --- a/examples/DIVAnd_SST_pluto.jl +++ b/examples/DIVAnd_SST_pluto.jl @@ -60,10 +60,10 @@ begin pm = ones(sz) ./ ((xi[2,1]-xi[1,1]) .* cosd.(yi)); pn = ones(sz) / (yi[1,2]-yi[1,1]); - md"""### Illustration of DIVAnd + md"""### Illustration of DIVAnd - We use the Reynolds et al. 2002 [OI SST](https://www.psl.noaa.gov/data/gridded/data.noaa.oisst.v2.html) for the month January and remove the zonal average. We extract pseudo-observations at random locations and aim to reconstruct the field from these data points by using DIVAnd. - """ + We use the Reynolds et al. 2002 [OI SST](https://www.psl.noaa.gov/data/gridded/data.noaa.oisst.v2.html) for the month January and remove the zonal average. We extract pseudo-observations at random locations and aim to reconstruct the field from these data points by using DIVAnd. + """ end @@ -118,15 +118,15 @@ begin title = "Observations",label = :none, markersize = 2, clim = clim, edgecolor = :none) - end - clim = extrema(filter(isfinite,v)) + end + clim = extrema(filter(isfinite,v)) plot( plotmap(v,"True field", clim = clim), plotmap(vi,"Analysis field", clim = clim), - scattermap(vobs, clim = clim), + scattermap(vobs, clim = clim), plotmap(vi-v,"Analysis - true field", clim = (-4,4), c = :bluesreds), - size=(800,400) + size=(800,400) ) end diff --git a/examples/DIVAnd_simple_test_fluxes.jl b/examples/DIVAnd_simple_test_fluxes.jl index 199254d5..d056d360 100644 --- a/examples/DIVAnd_simple_test_fluxes.jl +++ b/examples/DIVAnd_simple_test_fluxes.jl @@ -71,9 +71,9 @@ fluxesafter=zeros(size(h)[2]) for j=1:size(h)[2] for i=2:size(h)[1]-2 - if mask[i,j]&& mask[i+1,j] - fluxesafter[j]=fluxesafter[j]+h[i,j]*(fi[i+1,j]-fi[i,j]) - end + if mask[i,j]&& mask[i+1,j] + fluxesafter[j]=fluxesafter[j]+h[i,j]*(fi[i+1,j]-fi[i,j]) + end end end @show var(fluxes1+fluxesafter) @@ -97,9 +97,9 @@ fluxesafter=zeros(size(h)[1]) for i=1:size(h)[1] for j=2:size(h)[2]-2 - if mask[i,j]&& mask[i,j+1] - fluxesafter[i]=fluxesafter[i]+h[i,j]*(fi[i,j+1]-fi[i,j]) - end + if mask[i,j]&& mask[i,j+1] + fluxesafter[i]=fluxesafter[i]+h[i,j]*(fi[i,j+1]-fi[i,j]) + end end end @@ -129,9 +129,9 @@ fluxesafter=zeros(size(h)[2]) for j=1:size(h)[2] for i=2:size(h)[1]-2 - if mask[i,j]&& mask[i+1,j] - fluxesafter[j]=fluxesafter[j]+h[i,j]*(fi[i+1,j]-fi[i,j]) - end + if mask[i,j]&& mask[i+1,j] + fluxesafter[j]=fluxesafter[j]+h[i,j]*(fi[i+1,j]-fi[i,j]) + end end end @show var(fluxes1+fluxesafter) @@ -141,9 +141,9 @@ fluxesafter=zeros(size(h)[1]) for i=1:size(h)[1] for j=2:size(h)[2]-2 - if mask[i,j]&& mask[i,j+1] - fluxesafter[i]=fluxesafter[i]+h[i,j]*(fi[i,j+1]-fi[i,j]) - end + if mask[i,j]&& mask[i,j+1] + fluxesafter[i]=fluxesafter[i]+h[i,j]*(fi[i,j+1]-fi[i,j]) + end end end diff --git a/src/DIVAnd.jl b/src/DIVAnd.jl index 87a4e9fa..a9c947b5 100644 --- a/src/DIVAnd.jl +++ b/src/DIVAnd.jl @@ -598,7 +598,7 @@ export sparse_stagger, sparse_diff, localize_separable_grid, ndgrid, - localresolution, + localresolution, sparse_pack, sparse_interp, sparse_trim, diff --git a/src/DIVAnd_aexerr.jl b/src/DIVAnd_aexerr.jl index e1ac93c6..4397cefd 100644 --- a/src/DIVAnd_aexerr.jl +++ b/src/DIVAnd_aexerr.jl @@ -203,7 +203,7 @@ function DIVAnd_aexerr(mask, pmn, xi, x, f, len, epsilon2; rng=Random.GLOBAL_RNG f1, s1 = DIVAndrun(mask, pmn, xi, xfake, ffake, len ./ 1.70766, epsilonforB; otherargs...) # Calculate final error - aexerr = max.(Bjmb - f1, 0.0) + aexerr = max.(Bjmb - f1, 0.0) #@show mean(aexerr[.!isnan.(aexerr)]) diff --git a/src/DIVAnd_diagnostics.jl b/src/DIVAnd_diagnostics.jl index 35fb1a64..8869f007 100644 --- a/src/DIVAnd_diagnostics.jl +++ b/src/DIVAnd_diagnostics.jl @@ -1,5 +1,5 @@ """ - norm1,norm2,norm3,epsilon2=DIVAnd_norms(fi,s) + norm1,norm2,norm3,epsilon2=DIVAnd_norms(fi,s) ### Input: @@ -25,18 +25,18 @@ function DIVAnd_norms(fi,s) R = s.R iB = s.iB H = s.H - fstate = statevector_pack(s.sv, (fi,)) - residuals=s.yo - (s.H) * fstate - residualsobs=s.obsconstrain.yo - (s.obsconstrain.H) *fstate + fstate = statevector_pack(s.sv, (fi,)) + residuals=s.yo - (s.H) * fstate + residualsobs=s.obsconstrain.yo - (s.obsconstrain.H) *fstate residualsobs[s.obsout] .= 0.0 - norm2= residualsobs'*(s.obsconstrain.R\residualsobs) + norm2= residualsobs'*(s.obsconstrain.R\residualsobs) norm3= residuals'*(R \ residuals)-norm2 norm1=fstate'*(iB*fstate) - dR=diag(s.obsconstrain.R) - epsilon2=mean(dR[isfinite.(dR)]) + dR=diag(s.obsconstrain.R) + epsilon2=mean(dR[isfinite.(dR)]) - return norm1,norm2,norm3,epsilon2 + return norm1,norm2,norm3,epsilon2 end diff --git a/src/DIVAnd_errormap.jl b/src/DIVAnd_errormap.jl index c83fc185..1b3c5d4e 100644 --- a/src/DIVAnd_errormap.jl +++ b/src/DIVAnd_errormap.jl @@ -1,6 +1,6 @@ """ error,method = DIVAnd_errormap(mask,pmn,xi,x,f,len,epsilon2, - s; + s; method = :auto, Bscale = false, otherargs...,); @@ -35,9 +35,9 @@ *`method` : the method to be used, valid are `:auto`, `:cheap`, `:precise`, `:cpme`, `:scpme`, `:exact`, `:aexerr`, `:diagapp` auto will select depenting on data coverage and length scale - cheap will also select but restrict to cheaper methods - precise will also select but prefer better approximations - the other choices are just the ones among which the automatic choices will choose from. You can force the choice by specifying the method. + cheap will also select but restrict to cheaper methods + precise will also select but prefer better approximations + the other choices are just the ones among which the automatic choices will choose from. You can force the choice by specifying the method. *`Bscale` : it `true` will try to take out the the boundary effects in the background error variance. Not possible with all methods @@ -62,7 +62,7 @@ function DIVAnd_errormap( s; method = :auto, Bscale = false, - rng=Random.GLOBAL_RNG, + rng=Random.GLOBAL_RNG, otherargs... ) @@ -141,11 +141,11 @@ function DIVAnd_errormap( end end else - if Bigdata + if Bigdata errmethod = :scpme else errmethod = :aexerr - end + end end @@ -168,11 +168,11 @@ function DIVAnd_errormap( end end else - if Bigdata + if Bigdata errmethod = :scpme else errmethod = :cpme - end + end end end @@ -192,11 +192,11 @@ function DIVAnd_errormap( end end else - if Bigdata + if Bigdata errmethod = :diagapp else errmethod = :aexerr - end + end end end @@ -261,7 +261,7 @@ function DIVAnd_errormap( f, len, epsilon2; - otherargs... + otherargs... ) scpme=deepcopy(errormap) @@ -296,21 +296,21 @@ function DIVAnd_errormap( f, len, epsilon2; - rng=rng, + rng=rng, otherargs... ) - if errormap==0 - @warn "too fine resolution for aexerr, using exact" - errormap, =statevector_unpack(s.sv,diag(s.P)) + if errormap==0 + @warn "too fine resolution for aexerr, using exact" + errormap, =statevector_unpack(s.sv,diag(s.P)) return errormap, errmethod - end - if ScalebyB - return errormap./bi, errmethod - else + end + if ScalebyB + return errormap./bi, errmethod + else return errormap, errmethod - end + end end @show "You should not be here" end diff --git a/src/DIVAnd_heatmap.jl b/src/DIVAnd_heatmap.jl index f9c49c4c..09cb499a 100644 --- a/src/DIVAnd_heatmap.jl +++ b/src/DIVAnd_heatmap.jl @@ -108,8 +108,8 @@ function DIVAnd_heatmap( #varxb=zeros(Float64,DIMS) LF = zeros(Float64, DIMS) #fromgausstodivaL=[0.8099,0.62,0.62,0.62,0.62,0.62,0.62] - #{\left(\sum_i \eta_i \right) ^2 \over \sum_i \eta_i^2} - NPEFF=(sum(inflation))^2/sum(inflation .^ 2) + #{\left(\sum_i \eta_i \right) ^2 \over \sum_i \eta_i^2} + NPEFF=(sum(inflation))^2/sum(inflation .^ 2) for i = 1:DIMS meanxo = sum(inflation .* x[i]) / sum(inflation) varx[i] = sum(inflation .* (x[i] .- meanxo) .^ 2) / sum(inflation) diff --git a/src/DIVAnd_multivarEOF.jl b/src/DIVAnd_multivarEOF.jl index 23008588..9d24c666 100644 --- a/src/DIVAnd_multivarEOF.jl +++ b/src/DIVAnd_multivarEOF.jl @@ -42,7 +42,7 @@ For oceanographic application, this is the land-sea mask where sea is true and l coefficients from a linear fit between the variables (typically obtained by preliminary statistics or a run of this multivariate routine on a larger data set and which produced the eof coefficients). `eof`=[1.0,1.0] for example means the two variables are positively correlated with a slope 1. - `eof`=[1.0,-0.5] means the variables are negatively correlated and that for a variable 1 value of 1, variable 2 is expected to have a value of -0.5 + `eof`=[1.0,-0.5] means the variables are negatively correlated and that for a variable 1 value of 1, variable 2 is expected to have a value of -0.5 @@ -285,7 +285,7 @@ function DIVAnd_multivarEOF(mask,pmn,xi,x,f,lenin,epsilon2in;eof=(),velocity=(), for k=2:nlay coo=ndims(mask) # just take the points from one layer and copy them into other layers - # add 0.0001 to make sure rounding is not a problem + # add 0.0001 to make sure rounding is not a problem xm[coo][(k-1)*nd+1:k*nd].= mod.(xm[coo][(k-2)*nd+1:(k-1)*nd] .+0.0001,nlay).+1 end # Values of data are unimportant for the error field. So just repeated diff --git a/src/DIVAnd_multivarJAC.jl b/src/DIVAnd_multivarJAC.jl index ffd1578c..1318c84d 100644 --- a/src/DIVAnd_multivarJAC.jl +++ b/src/DIVAnd_multivarJAC.jl @@ -10,8 +10,8 @@ DIVAnd_multivarJAC(mask,pmn,xi,x,f,lenin,epsilon2in;epsilon2jacobian=1.0,kwargs. * `epsilon2jacobian` : epsilon2 array (or value) which defines the strength of the multivariate constraint for each variable. If you use a very large value for one variable, it means that variable will not be modified by the coupling, but the other will be. So typically - if you have a habitat variable which you want to influence a variable for which you have few data only. Then you assign a large epsilon2forjacobian to the habitat variable and - a lower value to the other one. + if you have a habitat variable which you want to influence a variable for which you have few data only. Then you assign a large epsilon2forjacobian to the habitat variable and + a lower value to the other one. # Output: @@ -94,7 +94,7 @@ DIVAnd_multivarJAC(mask,pmn,xi,x,f,lenin,epsilon2in;epsilon2jacobian=1.0,kwargs. emap,meth=DIVAnd_errormap(mask,pmn,xi,x,f,len,epsilon2,s;method=:cheap,kwargs...) - sori=deepcopy(s) + sori=deepcopy(s) @show "error method in multivar $meth" vconstrain=[] diff --git a/src/DIVAnd_scalecpme!.jl b/src/DIVAnd_scalecpme!.jl index 547f2477..440f8787 100644 --- a/src/DIVAnd_scalecpme!.jl +++ b/src/DIVAnd_scalecpme!.jl @@ -3,9 +3,9 @@ function DIVAnd_scalecpme!(cpme, P::CovarIS, nsamples = 7;rng=Random.GLOBAL_RNG) # P returned in structure s (so s.P) from a previous run # nsamples is the number of random arrays used to estimate the value - fractionshift=0.5 + fractionshift=0.5 - z = randn(rng,(size(P)[1], nsamples)) + z = randn(rng,(size(P)[1], nsamples)) errscale = 1 if P.factors != nothing ZZ = P.factors.PtL \ z @@ -14,16 +14,16 @@ function DIVAnd_scalecpme!(cpme, P::CovarIS, nsamples = 7;rng=Random.GLOBAL_RNG) ZZ = P * z errscale = mean(diag(z' * ZZ) ./ diag(z' * z)) end - # errscale here is the target + # errscale here is the target - # oldvalue is what we had - oldvalue=mean(cpme[.!isnan.(cpme)]) + # oldvalue is what we had + oldvalue=mean(cpme[.!isnan.(cpme)]) - #try a shift of fraction fractionshift of the mismatch + #try a shift of fraction fractionshift of the mismatch - cpme[:] .= cpme[:] .+ fractionshift*(errscale-oldvalue) + cpme[:] .= cpme[:] .+ fractionshift*(errscale-oldvalue) - # the rest is corrected by the multiplicative scaling + # the rest is corrected by the multiplicative scaling cpme[:] .= cpme[:] * errscale / mean(cpme[.!isnan.(cpme)]) diff --git a/src/DIVAndfun.jl b/src/DIVAndfun.jl index 5eace12c..1a4c047f 100644 --- a/src/DIVAndfun.jl +++ b/src/DIVAndfun.jl @@ -5,7 +5,7 @@ """ - myfunny=DIVAndfun(x,f;mask=nothing,pmn=nothing,xi=nothing,len=nothing,epsilon2=nothing,kwargs...) + myfunny=DIVAndfun(x,f;mask=nothing,pmn=nothing,xi=nothing,len=nothing,epsilon2=nothing,kwargs...) Provides a simplified interface to DIVAndrun and in return provides an interpolation FUNCTION rather than the gridded field of DIVAndrun @@ -118,5 +118,5 @@ function DIVAndfun(x,f;mask=nothing,pmn=nothing,xi=nothing,len=nothing,epsilon2= # Now initialize interpolations function and return that #return LinearInterpolation(tuple(collect.(coords)...), fi.+backg) - return linear_interpolation(tuple(collect.(coords)...), fi.+backg) + return linear_interpolation(tuple(collect.(coords)...), fi.+backg) end \ No newline at end of file diff --git a/src/DIVAndrun.jl b/src/DIVAndrun.jl index f5ecd90d..0f815edf 100644 --- a/src/DIVAndrun.jl +++ b/src/DIVAndrun.jl @@ -15,8 +15,8 @@ function DIVAndrun( maxit = 100, minit = 0, constraints = (), - ineqconstraints = (), - ntriesmax = 10, + ineqconstraints = (), + ntriesmax = 10, inversion = :chol, moddim = [], fracindex = Matrix{T}(undef, 0, 0), @@ -34,7 +34,7 @@ function DIVAndrun( fluxes = (), epsfluxes = 0, epsilon2forfractions = 0, - RTIMESONESCALES = (), + RTIMESONESCALES = (), QCMETHOD = (), coeff_laplacian::Vector{Float64} = ones(ndims(mask)), coeff_derivative2::Vector{Float64} = zeros(ndims(mask)), @@ -45,7 +45,7 @@ function DIVAndrun( - if !isempty(ineqconstraints) + if !isempty(ineqconstraints) @@ -63,18 +63,18 @@ end - allconstraints=constraints - ineqok=false - ntries=0 - fi=() - s=() + allconstraints=constraints + ineqok=false + ntries=0 + fi=() + s=() - while !ineqok && ntries0 - ineqok=false - @show violated,ntries - end - allconstraints=(allconstraints...,onemoreconstraint)# Add to the list of constraints + if violated>0 + ineqok=false + @show violated,ntries + end + allconstraints=(allconstraints...,onemoreconstraint)# Add to the list of constraints end - end - return fi,s - end + end + return fi,s + end # End of special treatment for inequality constraints # check pmn .* len > 4 @@ -217,7 +217,7 @@ end # add all additional constrains for i = 1:length(constraints) #@info "Adding additional constrain - $(i)" - s = DIVAnd_addc(s, constraints[i]) + s = DIVAnd_addc(s, constraints[i]) end # factorize a posteriori error covariance matrix diff --git a/src/fit.jl b/src/fit.jl index 412eb260..79c637df 100644 --- a/src/fit.jl +++ b/src/fit.jl @@ -426,7 +426,7 @@ function fitlen( ) if length(d) == 0 @warn "no data is provided to fitlen" - dbinfo = Dict{Symbol,Any}( + dbinfo = Dict{Symbol,Any}( :covar => [], :fitcovar => [], :distx => [], @@ -516,7 +516,7 @@ function fitlen( @debug "Number of probable active bins: $rnbins" ddist = meandist / rnbins - nbmax=1 + nbmax=1 worktmp=maxdist / ddist if !isnan(worktmp)&&isfinite(worktmp) nbmax = floor(Int, worktmp + 1) diff --git a/src/special_matrices.jl b/src/special_matrices.jl index 5bde71bf..48915bd6 100644 --- a/src/special_matrices.jl +++ b/src/special_matrices.jl @@ -53,7 +53,7 @@ function Base.:*(C::CovarIS, v::TV)::TV where {TV<:AbstractVector{Float64}} @show norm(C.IS * x - v) end @show convergence_history - return x + return x elseif C.factors != nothing return C.factors \ v else diff --git a/test/runtests.jl b/test/runtests.jl index b59b07e9..459c79cc 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -31,11 +31,11 @@ using Downloads: download include("test_1D_seminormed.jl") - # multivar + # multivar - include("test_multivar.jl") + include("test_multivar.jl") - include("test_errormaps.jl") + include("test_errormaps.jl") # dynamical constraints include("test_2dvar_adv.jl") @@ -134,8 +134,8 @@ using Downloads: download # Heatmap include("test_heatmap.jl") - # DIVAndfun - include("test_DIVAndfun.jl") + # DIVAndfun + include("test_DIVAndfun.jl") # test DIVAnd_filter3 A = zeros(5, 5, 5, 5, 5) diff --git a/test/test_2dvar_adv.jl b/test/test_2dvar_adv.jl index 02b5febd..900b2f81 100644 --- a/test/test_2dvar_adv.jl +++ b/test/test_2dvar_adv.jl @@ -56,7 +56,7 @@ fi, s = DIVAndrun( epsilon2; velocity = (u, v), alphabc = 0, - ineqconstraints=(m2c,) + ineqconstraints=(m2c,) ); @test fi[18,24] ≈ 0.5381625233419283 diff --git a/test/test_errormaps.jl b/test/test_errormaps.jl index 43e1e020..786b634f 100644 --- a/test/test_errormaps.jl +++ b/test/test_errormaps.jl @@ -54,12 +54,12 @@ errorm,methodc=DIVAnd_errormap( epsilon2,sl; velocity = (u, v), alphabc = 0, - method=mymeth, - Bscale=myscale, - rng=StableRNG(123) - ) + method=mymeth, + Bscale=myscale, + rng=StableRNG(123) + ) - #@show size(errorm),size(fi11) + #@show size(errorm),size(fi11) global errmean=errmean+errorm[10,10] #@show errmean,errorm[10,10],fi11[10,10],xdd[1],ND,len,mymeth,myscale,methodc From 6a7be8f2f6787b9eb0a433db34d832b62b3af2e3 Mon Sep 17 00:00:00 2001 From: jmbeckers Date: Mon, 3 Apr 2023 22:29:39 +0200 Subject: [PATCH 47/51] Update README.md --- README.md | 1 + 1 file changed, 1 insertion(+) diff --git a/README.md b/README.md index c1c3dcf5..beae5495 100644 --- a/README.md +++ b/README.md @@ -36,6 +36,7 @@ Barth, A., Beckers, J.-M., Troupin, C., Alvera-Azcárate, A., and Vandenbulcke, * Periodicity in selected directions can be enforced * Multivariate data can be used (experimental) * The output grid can be curvilinear +* Instead of interpolating scattered data you can also peform Kernel Density Estimations with the points. # Installing From d6cb1dca13f26ba01e2748594c5558e271fa50c8 Mon Sep 17 00:00:00 2001 From: Alexander Barth Date: Thu, 27 Apr 2023 12:02:30 +0200 Subject: [PATCH 48/51] use name from PROJECTS as project name --- src/SDNMetadata.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/SDNMetadata.jl b/src/SDNMetadata.jl index da3b3f44..2d59ab17 100644 --- a/src/SDNMetadata.jl +++ b/src/SDNMetadata.jl @@ -137,6 +137,8 @@ function SDNMetadata( ], "-", ) + elseif k == "project" + v = project["name"] end if typeof(v) <: Vector From 603c65d58d462de095a94227f7d884edbe53ed83 Mon Sep 17 00:00:00 2001 From: Alexander Barth Date: Thu, 27 Apr 2023 12:10:03 +0200 Subject: [PATCH 49/51] update attributes for data access and visualization --- src/SDNMetadata.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/SDNMetadata.jl b/src/SDNMetadata.jl index 2d59ab17..befda50a 100644 --- a/src/SDNMetadata.jl +++ b/src/SDNMetadata.jl @@ -5,9 +5,9 @@ const PROJECTS = Dict( "EMODNET-chemistry" => Dict( "name" => "EMODnet Chemistry", "URL" => "https://emodnet.ec.europa.eu/en/chemistry", - "baseurl_visualization" => "http://ec.oceanbrowser.net/emodnet/", + "baseurl_visualization" => "https://emodnet.ec.europa.eu/geoviewer", "baseurl_wms" => "http://ec.oceanbrowser.net/emodnet/Python/web/wms", - "baseurl_http" => "http://ec.oceanbrowser.net/data/emodnet-domains", + "baseurl_http" => "https://emodnet.ec.europa.eu/geoviewer", "baseurl_opendap" => "http://opendap.oceanbrowser.net/thredds/dodsC/data/emodnet-domains", "template" => joinpath(pathname, "templates", "emodnet-chemistry.mustache"), ), From e2863a511e73c76afe28535b30cf7b10bf0c7e42 Mon Sep 17 00:00:00 2001 From: Alexander Barth Date: Thu, 27 Apr 2023 12:18:42 +0200 Subject: [PATCH 50/51] test for EMODnet Chemistry --- test/test_product.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_product.jl b/test/test_product.jl index 2cb10067..386d4753 100644 --- a/test/test_product.jl +++ b/test/test_product.jl @@ -103,8 +103,8 @@ varname = "Salinity" filename = tempname() metadata = OrderedDict( - # Name of the project (SeaDataCloud, SeaDataNet, EMODNET-Chemistry, ...) - "project" => "SeaDataCloud", + # Name of the project (SeaDataCloud, SeaDataNet, EMODNET-chemistry, ...) + "project" => "EMODNET-chemistry", # URN code for the institution EDMO registry, # e.g. SDN:EDMO::1579 From 9af09ef7494d03df1e261f65841be707940b700e Mon Sep 17 00:00:00 2001 From: CompatHelper Julia Date: Thu, 1 Jun 2023 01:56:38 +0000 Subject: [PATCH 51/51] CompatHelper: add new compat entry for DelimitedFiles at version 1, (keep existing compat) --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index 12471f00..8e5bb9fe 100644 --- a/Project.toml +++ b/Project.toml @@ -33,6 +33,7 @@ ZipFile = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea" [compat] AlgebraicMultigrid = "0.2, 0.3, 0.4, 0.5" DataStructures = "0.17, 0.18" +DelimitedFiles = "1" EzXML = "0.8, 0.9, 1.1" HTTP = "1" Interpolations = "0.12, 0.13, 0.14"