Merge branch 'master' of github.com:gher-ulg/DIVAnd.jl

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"

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=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)
@@ -24,6 +24,20 @@ 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
+* Instead of interpolating scattered data you can also peform Kernel Density Estimations with the points.
+
 
 # Installing
 
@@ -130,7 +144,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`.
 
@@ -139,6 +153,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.
@@ -163,13 +179,13 @@ 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
 
 ### 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**) = ∑<sub>*i*</sub> (**C**<sub>*i*</sub> **x**  - **z**<sub>*i*</sub>)ᵀ **Q**<sub>*i*</sub><sup>-1</sup> (**C**<sub>*i*</sub> **x** - **z**<sub>*i*</sub>)
 
@@ -181,6 +197,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 inequality constraint
+
+An arbitrary number of additional inequality constraints can be included and which should have the following form:
+
+(**H**<sub>*i*</sub> **x**  > **yo**<sub>*i*</sub>)
+
+For every constraint, 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:
@@ -211,6 +240,16 @@ Please include the following information when reporting an issue:
 
 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
 
 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
+
+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

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
 
 
@@ -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...)
@@ -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
 

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
 

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
@@ -590,7 +598,7 @@ export sparse_stagger,
     sparse_diff,
     localize_separable_grid,
     ndgrid,
-	localresolution,
+    localresolution,
     sparse_pack,
     sparse_interp,
     sparse_trim,
@@ -647,6 +655,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

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]
 
@@ -147,7 +147,7 @@ function DIVAnd_aexerr(mask, pmn, xi, x, f, len, epsilon2; otherargs...)
     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
@@ -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)])
 

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

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