diff --git a/.github/ISSUE_TEMPLATE/bug_report.md b/.github/ISSUE_TEMPLATE/bug_report.md index 674c8f27..0edd35ce 100644 --- a/.github/ISSUE_TEMPLATE/bug_report.md +++ b/.github/ISSUE_TEMPLATE/bug_report.md @@ -2,7 +2,7 @@ name: Bug report about: 'Tell us what''s going wrong ' title: "[\U0001F41E]" -labels: bug +labels: bug, enhancement assignees: '' --- diff --git a/.github/ISSUE_TEMPLATE/documentation-improvement.md b/.github/ISSUE_TEMPLATE/documentation-improvement.md new file mode 100644 index 00000000..25fb5593 --- /dev/null +++ b/.github/ISSUE_TEMPLATE/documentation-improvement.md @@ -0,0 +1,10 @@ +--- +name: Documentation improvement +about: Indicate the part of the code (functions, module) to improve +title: "[DOC]" +labels: documentation +assignees: '' + +--- + + 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" diff --git a/README.md b/README.md index a42a6d02..beae5495 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=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**) = ∑*i* (**C***i* **x** - **z***i*)ᵀ **Q***i*-1 (**C***i* **x** - **z***i*) @@ -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***i* **x** > **yo***i*) + +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 diff --git a/examples/DIVAnd_SST_pluto.jl b/examples/DIVAnd_SST_pluto.jl index beb268be..453aecc1 100644 --- a/examples/DIVAnd_SST_pluto.jl +++ b/examples/DIVAnd_SST_pluto.jl @@ -1,5 +1,5 @@ ### A Pluto.jl notebook ### -# v0.19.14 +# v0.19.26 using Markdown using InteractiveUtils @@ -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 @@ -154,18 +154,18 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" [compat] DIVAnd = "~2.7.9" -NCDatasets = "~0.12.8" -Plots = "~1.36.0" -PlutoUI = "~0.7.48" +NCDatasets = "~0.12.11" +Plots = "~1.38.0" +PlutoUI = "~0.7.49" """ # ╔═╡ 00000000-0000-0000-0000-000000000002 PLUTO_MANIFEST_TOML_CONTENTS = """ # This file is machine-generated - editing it directly is not advised -julia_version = "1.8.2" +julia_version = "1.9.0" manifest_format = "2.0" -project_hash = "50c7174e499f8549ff0ba5726f3dbef914d81bce" +project_hash = "3e4d1ab59f6e7abd4d49381315d59ab819200b77" [[deps.AbstractPlutoDingetjes]] deps = ["Pkg"] @@ -174,10 +174,14 @@ uuid = "6e696c72-6542-2067-7265-42206c756150" version = "1.1.4" [[deps.Adapt]] -deps = ["LinearAlgebra"] -git-tree-sha1 = "195c5505521008abea5aee4f96930717958eac6f" +deps = ["LinearAlgebra", "Requires"] +git-tree-sha1 = "76289dc51920fdc6e0013c872ba9551d54961c24" uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" -version = "3.4.0" +version = "3.6.2" +weakdeps = ["StaticArrays"] + + [deps.Adapt.extensions] + AdaptStaticArraysExt = "StaticArrays" [[deps.AlgebraicMultigrid]] deps = ["CommonSolve", "LinearAlgebra", "Printf", "Reexport", "SparseArrays"] @@ -202,9 +206,9 @@ version = "1.0.1" uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" [[deps.BitFlags]] -git-tree-sha1 = "84259bb6172806304b9101094a7cc4bc6f56dbc6" +git-tree-sha1 = "43b1a4a8f797c1cddadf60499a8a077d4af2cd2d" uuid = "d1d4a3ce-64b1-5f1a-9ba4-7e7e69966f35" -version = "0.1.5" +version = "0.1.7" [[deps.Bzip2_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] @@ -219,34 +223,28 @@ uuid = "179af706-886a-5703-950a-314cd64e0468" version = "0.1.2" [[deps.Cairo_jll]] -deps = ["Artifacts", "Bzip2_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"] +deps = ["Artifacts", "Bzip2_jll", "CompilerSupportLibraries_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"] git-tree-sha1 = "4b859a208b2397a7a623a03449e4636bdb17bcf2" uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a" version = "1.16.1+1" [[deps.ChainRulesCore]] deps = ["Compat", "LinearAlgebra", "SparseArrays"] -git-tree-sha1 = "e7ff6cadf743c098e08fca25c91103ee4303c9bb" +git-tree-sha1 = "e30f2f4e20f7f186dc36529910beaedc60cfa644" uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" -version = "1.15.6" - -[[deps.ChangesOfVariables]] -deps = ["ChainRulesCore", "LinearAlgebra", "Test"] -git-tree-sha1 = "38f7a08f19d8810338d4f5085211c7dfa5d5bdd8" -uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0" -version = "0.1.4" +version = "1.16.0" [[deps.CodecZlib]] deps = ["TranscodingStreams", "Zlib_jll"] -git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da" +git-tree-sha1 = "9c209fb7536406834aa938fb149964b985de6c83" uuid = "944b1d66-785c-5afd-91f1-9de20f533193" -version = "0.7.0" +version = "0.7.1" [[deps.ColorSchemes]] -deps = ["ColorTypes", "ColorVectorSpace", "Colors", "FixedPointNumbers", "Random"] -git-tree-sha1 = "1fd869cc3875b57347f7027521f561cf46d1fcd8" +deps = ["ColorTypes", "ColorVectorSpace", "Colors", "FixedPointNumbers", "PrecompileTools", "Random"] +git-tree-sha1 = "be6ab11021cd29f0344d5c4357b163af05a48cba" uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4" -version = "3.19.0" +version = "3.21.0" [[deps.ColorTypes]] deps = ["FixedPointNumbers", "Random"] @@ -256,31 +254,61 @@ version = "0.11.4" [[deps.ColorVectorSpace]] deps = ["ColorTypes", "FixedPointNumbers", "LinearAlgebra", "SpecialFunctions", "Statistics", "TensorCore"] -git-tree-sha1 = "d08c20eef1f2cbc6e60fd3612ac4340b89fea322" +git-tree-sha1 = "600cc5508d66b78aae350f7accdb58763ac18589" uuid = "c3611d14-8923-5661-9e6a-0046d554d3a4" -version = "0.9.9" +version = "0.9.10" [[deps.Colors]] deps = ["ColorTypes", "FixedPointNumbers", "Reexport"] -git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40" +git-tree-sha1 = "fc08e5930ee9a4e03f84bfb5211cb54e7769758a" uuid = "5ae59095-9a9b-59fe-a467-6f913c188581" -version = "0.12.8" +version = "0.12.10" + +[[deps.CommonDataModel]] +deps = ["CFTime", "DataStructures", "Dates", "Preferences", "Printf"] +git-tree-sha1 = "2678b3fc170d582655a14d22867b031b6e43c2d4" +uuid = "1fbeeb36-5f17-413c-809b-666fb144f157" +version = "0.2.4" [[deps.CommonSolve]] -git-tree-sha1 = "9441451ee712d1aec22edad62db1a9af3dc8d852" +git-tree-sha1 = "0eee5eb66b1cf62cd6ad1b460238e60e4b09400c" uuid = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" -version = "0.2.3" +version = "0.2.4" [[deps.Compat]] -deps = ["Dates", "LinearAlgebra", "UUIDs"] -git-tree-sha1 = "3ca828fe1b75fa84b021a7860bd039eaea84d2f2" +deps = ["UUIDs"] +git-tree-sha1 = "7a60c856b9fa189eb34f5f8a6f6b5529b7942957" uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "4.3.0" +version = "4.6.1" +weakdeps = ["Dates", "LinearAlgebra"] + + [deps.Compat.extensions] + CompatLinearAlgebraExt = "LinearAlgebra" [[deps.CompilerSupportLibraries_jll]] deps = ["Artifacts", "Libdl"] uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" -version = "0.5.2+0" +version = "1.0.2+0" + +[[deps.ConcurrentUtilities]] +deps = ["Serialization", "Sockets"] +git-tree-sha1 = "96d823b94ba8d187a6d8f0826e731195a74b90e9" +uuid = "f0e56b4a-5159-44fe-b623-3e5288b988bb" +version = "2.2.0" + +[[deps.ConstructionBase]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "738fec4d684a9a6ee9598a8bfee305b26831f28c" +uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9" +version = "1.5.2" + + [deps.ConstructionBase.extensions] + ConstructionBaseIntervalSetsExt = "IntervalSets" + ConstructionBaseStaticArraysExt = "StaticArrays" + + [deps.ConstructionBase.weakdeps] + IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953" + StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [[deps.Contour]] git-tree-sha1 = "d05d9e7b7aedff4e5b51a029dced05cfb6125781" @@ -294,9 +322,9 @@ uuid = "efc8151c-67de-5a8f-9a35-d8f54746ae9d" version = "2.7.9" [[deps.DataAPI]] -git-tree-sha1 = "46d2680e618f8abd007bce0c3026cb0c4a8f2032" +git-tree-sha1 = "8da84edb865b0b5b0100c0666a9bc9a0b71c553c" uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" -version = "1.12.0" +version = "1.15.0" [[deps.DataStructures]] deps = ["Compat", "InteractiveUtils", "OrderedCollections"] @@ -315,7 +343,9 @@ uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" [[deps.DelimitedFiles]] deps = ["Mmap"] +git-tree-sha1 = "9e2f36d3c96a820c678f2f1f1782582fcf685bae" uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab" +version = "1.9.1" [[deps.Distributed]] deps = ["Random", "Serialization", "Sockets"] @@ -323,20 +353,26 @@ uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" [[deps.DocStringExtensions]] deps = ["LibGit2"] -git-tree-sha1 = "c36550cb29cbe373e95b3f40486b9a4148f89ffd" +git-tree-sha1 = "2fb1e02f2b635d0845df5d7c167fec4dd739b00d" uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -version = "0.9.2" +version = "0.9.3" [[deps.Downloads]] deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"] uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" version = "1.6.0" +[[deps.ExceptionUnwrapping]] +deps = ["Test"] +git-tree-sha1 = "e90caa41f5a86296e014e148ee061bd6c3edec96" +uuid = "460bff9d-24e4-43bc-9d9f-a8973cb893f4" +version = "0.1.9" + [[deps.Expat_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "bad72f730e9e91c08d9427d5e8db95478a3c323d" +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "4558ab818dcceaab612d1bb8c19cee87eda2b83c" uuid = "2e619515-83b5-522b-bb60-26c02a35a201" -version = "2.4.8+0" +version = "2.5.0+0" [[deps.EzXML]] deps = ["Printf", "XML2_jll"] @@ -396,16 +432,16 @@ uuid = "0656b61e-2033-5cc2-a64a-77c0f6c09b89" version = "3.3.8+0" [[deps.GR]] -deps = ["Base64", "DelimitedFiles", "GR_jll", "HTTP", "JSON", "Libdl", "LinearAlgebra", "Pkg", "Preferences", "Printf", "Random", "Serialization", "Sockets", "Test", "UUIDs"] -git-tree-sha1 = "00a9d4abadc05b9476e937a5557fcce476b9e547" +deps = ["Artifacts", "Base64", "DelimitedFiles", "Downloads", "GR_jll", "HTTP", "JSON", "Libdl", "LinearAlgebra", "Pkg", "Preferences", "Printf", "Random", "Serialization", "Sockets", "TOML", "Tar", "Test", "UUIDs", "p7zip_jll"] +git-tree-sha1 = "8b8a2fd4536ece6e554168c21860b6820a8a83db" uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71" -version = "0.69.5" +version = "0.72.7" [[deps.GR_jll]] -deps = ["Artifacts", "Bzip2_jll", "Cairo_jll", "FFMPEG_jll", "Fontconfig_jll", "GLFW_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pixman_jll", "Pkg", "Qt5Base_jll", "Zlib_jll", "libpng_jll"] -git-tree-sha1 = "bc9f7725571ddb4ab2c4bc74fa397c1c5ad08943" +deps = ["Artifacts", "Bzip2_jll", "Cairo_jll", "FFMPEG_jll", "Fontconfig_jll", "GLFW_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pixman_jll", "Qt5Base_jll", "Zlib_jll", "libpng_jll"] +git-tree-sha1 = "19fad9cd9ae44847fe842558a744748084a722d1" uuid = "d2c73de3-f751-5644-a686-071e5b155ba9" -version = "0.69.1+0" +version = "0.72.7+0" [[deps.Gettext_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"] @@ -415,9 +451,9 @@ version = "0.21.0+0" [[deps.Glib_jll]] deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE2_jll", "Pkg", "Zlib_jll"] -git-tree-sha1 = "fb83fbe02fe57f2c068013aa94bcdf6760d3a7a7" +git-tree-sha1 = "d3b3624125c1474292d0d8ed0f65554ac37ddb23" uuid = "7746bdde-850d-59dc-9ae8-88ece973131d" -version = "2.74.0+1" +version = "2.74.0+2" [[deps.Graphite2_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] @@ -437,10 +473,10 @@ uuid = "0234f1f7-429e-5d53-9886-15a909be8d59" version = "1.12.2+2" [[deps.HTTP]] -deps = ["Base64", "CodecZlib", "Dates", "IniFile", "Logging", "LoggingExtras", "MbedTLS", "NetworkOptions", "OpenSSL", "Random", "SimpleBufferStream", "Sockets", "URIs", "UUIDs"] -git-tree-sha1 = "8556f4b387fcd1d9b3013d798eecbcfa0d985e66" +deps = ["Base64", "CodecZlib", "ConcurrentUtilities", "Dates", "ExceptionUnwrapping", "Logging", "LoggingExtras", "MbedTLS", "NetworkOptions", "OpenSSL", "Random", "SimpleBufferStream", "Sockets", "URIs", "UUIDs"] +git-tree-sha1 = "2613d054b0e18a3dea99ca1594e9a3960e025da4" uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3" -version = "1.5.2" +version = "1.9.7" [[deps.HarfBuzz_jll]] deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg"] @@ -462,14 +498,9 @@ version = "0.9.4" [[deps.IOCapture]] deps = ["Logging", "Random"] -git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a" +git-tree-sha1 = "d75853a0bdbfb1ac815478bacd89cd27b550ace6" uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89" -version = "0.2.2" - -[[deps.IniFile]] -git-tree-sha1 = "f550e6e32074c939295eb5ea6de31849ac2c9625" -uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f" -version = "0.5.1" +version = "0.2.3" [[deps.InteractiveUtils]] deps = ["Markdown"] @@ -477,20 +508,14 @@ uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" [[deps.Interpolations]] deps = ["Adapt", "AxisAlgorithms", "ChainRulesCore", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "Requires", "SharedArrays", "SparseArrays", "StaticArrays", "WoodburyMatrices"] -git-tree-sha1 = "842dd89a6cb75e02e85fdd75c760cdc43f5d6863" +git-tree-sha1 = "721ec2cf720536ad005cb38f50dbba7b02419a15" uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" -version = "0.14.6" - -[[deps.InverseFunctions]] -deps = ["Test"] -git-tree-sha1 = "49510dfcb407e572524ba94aeae2fced1f3feb0f" -uuid = "3587e190-3f89-42d0-90ee-14403ec27112" -version = "0.1.8" +version = "0.14.7" [[deps.IrrationalConstants]] -git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151" +git-tree-sha1 = "630b497eafcc20001bba38a4651b327dcfc491d2" uuid = "92d709cd-6900-40b7-9082-c6be49f344b6" -version = "0.1.1" +version = "0.2.2" [[deps.IterativeSolvers]] deps = ["LinearAlgebra", "Printf", "Random", "RecipesBase", "SparseArrays"] @@ -517,15 +542,15 @@ version = "1.4.1" [[deps.JSON]] deps = ["Dates", "Mmap", "Parsers", "Unicode"] -git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e" +git-tree-sha1 = "31e996f0a15c7b280ba9f76636b3ff9e2ae58c9a" uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" -version = "0.21.3" +version = "0.21.4" [[deps.JpegTurbo_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "b53380851c6e6664204efb2e62cd24fa5c47e4ba" +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "6f2675ef130a300a112286de91973805fcc5ffbc" uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8" -version = "2.1.2+0" +version = "2.1.91+0" [[deps.LAME_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] @@ -539,6 +564,12 @@ git-tree-sha1 = "bf36f528eec6634efc60d7ec062008f171071434" uuid = "88015f11-f218-50d7-93a8-a6af411a945d" version = "3.0.0+1" +[[deps.LLVMOpenMP_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "f689897ccbe049adb19a065c495e75f372ecd42b" +uuid = "1d63c593-3942-5779-bab2-d838dc0a180e" +version = "15.0.4+0" + [[deps.LZO_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] git-tree-sha1 = "e5b909bcf985c5e2605737d2ce278ed791b89be6" @@ -552,9 +583,17 @@ version = "1.3.0" [[deps.Latexify]] deps = ["Formatting", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "OrderedCollections", "Printf", "Requires"] -git-tree-sha1 = "ab9aa169d2160129beb241cb2750ca499b4e90e9" +git-tree-sha1 = "f428ae552340899a935973270b8d98e5a31c49fe" uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316" -version = "0.15.17" +version = "0.16.1" + + [deps.Latexify.extensions] + DataFramesExt = "DataFrames" + SymEngineExt = "SymEngine" + + [deps.Latexify.weakdeps] + DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" + SymEngine = "123dc426-2d89-5057-bbad-38513e3affd8" [[deps.LibCURL]] deps = ["LibCURL_jll", "MozillaCACerts_jll"] @@ -592,9 +631,9 @@ version = "1.8.7+0" [[deps.Libglvnd_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll", "Xorg_libXext_jll"] -git-tree-sha1 = "7739f837d6447403596a75d19ed01fd08d6f56bf" +git-tree-sha1 = "6f73d1dd803986947b2c750138528a999a6c7733" uuid = "7e76a0d4-f3c7-5321-8279-8d96eeed0f29" -version = "1.3.0+3" +version = "1.6.0+0" [[deps.Libgpg_error_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] @@ -604,9 +643,9 @@ version = "1.42.0+0" [[deps.Libiconv_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "42b62845d70a619f063a7da093d995ec8e15e778" +git-tree-sha1 = "c7cb1f5d892775ba13767a87c7ada0b980ea0a71" uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531" -version = "1.16.1+1" +version = "1.16.1+2" [[deps.Libmount_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] @@ -627,23 +666,33 @@ uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700" version = "2.36.0+0" [[deps.LinearAlgebra]] -deps = ["Libdl", "libblastrampoline_jll"] +deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"] uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" [[deps.LogExpFunctions]] -deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"] -git-tree-sha1 = "94d9c52ca447e23eac0c0f074effbcd38830deb5" +deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"] +git-tree-sha1 = "c3ce8e7420b3a6e071e0fe4745f5d4300e37b13f" uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" -version = "0.3.18" +version = "0.3.24" + + [deps.LogExpFunctions.extensions] + LogExpFunctionsChainRulesCoreExt = "ChainRulesCore" + LogExpFunctionsChangesOfVariablesExt = "ChangesOfVariables" + LogExpFunctionsInverseFunctionsExt = "InverseFunctions" + + [deps.LogExpFunctions.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + ChangesOfVariables = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0" + InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" [[deps.Logging]] uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" [[deps.LoggingExtras]] deps = ["Dates", "Logging"] -git-tree-sha1 = "5d4d2d9904227b8bd66386c1138cf4d5ffa826bf" +git-tree-sha1 = "cedb76b37bc5a6c702ade66be44f831fa23c681e" uuid = "e6f89c97-d47a-5376-807f-9c37f3926c36" -version = "0.4.9" +version = "1.0.0" [[deps.MIMEs]] git-tree-sha1 = "65f28ad4b594aebe22157d6fac869786a255b7eb" @@ -669,43 +718,43 @@ version = "1.1.7" [[deps.MbedTLS_jll]] deps = ["Artifacts", "Libdl"] uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" -version = "2.28.0+0" +version = "2.28.2+0" [[deps.Measures]] -git-tree-sha1 = "e498ddeee6f9fdb4551ce855a46f54dbd900245f" +git-tree-sha1 = "c13304c81eec1ed3af7fc20e75fb6b26092a1102" uuid = "442fdcdd-2543-5da2-b0f3-8c86c306513e" -version = "0.3.1" +version = "0.3.2" [[deps.Missings]] deps = ["DataAPI"] -git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f" +git-tree-sha1 = "f66bdc5de519e8f8ae43bdc598782d35a25b1272" uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28" -version = "1.0.2" +version = "1.1.0" [[deps.Mmap]] uuid = "a63ad114-7e13-5084-954f-fe012c677804" [[deps.MozillaCACerts_jll]] uuid = "14a3606d-f60d-562e-9121-12d972cd8159" -version = "2022.2.1" +version = "2022.10.11" [[deps.Mustache]] deps = ["Printf", "Tables"] -git-tree-sha1 = "1e566ae913a57d0062ff1af54d2697b9344b99cd" +git-tree-sha1 = "821e918c170ead5298ff84bffee41dd28929a681" uuid = "ffc61752-8dc7-55ee-8c37-f3e9cdd09e70" -version = "1.0.14" +version = "1.0.17" [[deps.NCDatasets]] -deps = ["CFTime", "DataStructures", "Dates", "NetCDF_jll", "NetworkOptions", "Printf"] -git-tree-sha1 = "1818d54a489492517841c2cb87475ba20918264c" +deps = ["CFTime", "CommonDataModel", "DataStructures", "Dates", "NetCDF_jll", "NetworkOptions", "Printf"] +git-tree-sha1 = "4263c4220f22e20729329838bf7e94a49d1ac32f" uuid = "85f8d34a-cbdd-5861-8df4-14fed0d494ab" -version = "0.12.8" +version = "0.12.17" [[deps.NaNMath]] deps = ["OpenLibm_jll"] -git-tree-sha1 = "a7c3d1da1189a1c2fe843a3bfa04d18d20eb3211" +git-tree-sha1 = "0877504529a3e5c3343c6f8b4c0381e57e4387e4" uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" -version = "1.0.1" +version = "1.0.2" [[deps.NetCDF_jll]] deps = ["Artifacts", "HDF5_jll", "JLLWrappers", "LibCURL_jll", "Libdl", "Pkg", "XML2_jll", "Zlib_jll"] @@ -719,9 +768,9 @@ version = "1.2.0" [[deps.OffsetArrays]] deps = ["Adapt"] -git-tree-sha1 = "f71d8950b724e9ff6110fc948dff5a329f901d64" +git-tree-sha1 = "82d7c9e310fe55aa54996e6f7f94674e2a38fcb4" uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" -version = "1.12.8" +version = "1.12.9" [[deps.Ogg_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] @@ -732,7 +781,7 @@ version = "1.3.5+1" [[deps.OpenBLAS_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" -version = "0.3.20+0" +version = "0.3.21+4" [[deps.OpenLibm_jll]] deps = ["Artifacts", "Libdl"] @@ -741,15 +790,15 @@ version = "0.8.1+0" [[deps.OpenSSL]] deps = ["BitFlags", "Dates", "MozillaCACerts_jll", "OpenSSL_jll", "Sockets"] -git-tree-sha1 = "5628f092c6186a80484bfefdf89ff64efdaec552" +git-tree-sha1 = "51901a49222b09e3743c65b8847687ae5fc78eb2" uuid = "4d8831e6-92b7-49fb-bdf8-b643e874388c" -version = "1.3.1" +version = "1.4.1" [[deps.OpenSSL_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "f6e9dba33f9f2c44e08a020b0caf6903be540004" +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "1aa4b74f80b01c6bc2b89992b861b5f210e665b5" uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95" -version = "1.1.19+0" +version = "1.1.21+0" [[deps.OpenSpecFun_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] @@ -764,20 +813,20 @@ uuid = "91d4177d-7536-5919-b921-800302f37372" version = "1.3.2+0" [[deps.OrderedCollections]] -git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c" +git-tree-sha1 = "d321bf2de576bf25ec4d3e4360faca399afca282" uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" -version = "1.4.1" +version = "1.6.0" [[deps.PCRE2_jll]] deps = ["Artifacts", "Libdl"] uuid = "efcefdf7-47ab-520b-bdef-62a2eaa19f15" -version = "10.40.0+0" +version = "10.42.0+0" [[deps.Parsers]] -deps = ["Dates", "SnoopPrecompile"] -git-tree-sha1 = "cceb0257b662528ecdf0b4b4302eb00e767b38e7" +deps = ["Dates", "PrecompileTools", "UUIDs"] +git-tree-sha1 = "4b2e829ee66d4218e0cef22c0a64ee37cf258c29" uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" -version = "2.5.0" +version = "2.7.1" [[deps.Pipe]] git-tree-sha1 = "6842804e7867b115ca9de748a0cf6b364523c16d" @@ -785,15 +834,15 @@ uuid = "b98c9c47-44ae-5843-9183-064241ee97a0" version = "1.3.0" [[deps.Pixman_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "b4f5d02549a10e20780a24fce72bea96b6329e29" +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LLVMOpenMP_jll", "Libdl"] +git-tree-sha1 = "64779bc4c9784fee475689a1752ef4d5747c5e87" uuid = "30392449-352a-5448-841d-b1acce4e97dc" -version = "0.40.1+0" +version = "0.42.2+0" [[deps.Pkg]] -deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] +deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" -version = "1.8.0" +version = "1.9.0" [[deps.PlotThemes]] deps = ["PlotUtils", "Statistics"] @@ -802,28 +851,48 @@ uuid = "ccf2f8ad-2431-5c83-bf29-c5338b663b6a" version = "3.1.0" [[deps.PlotUtils]] -deps = ["ColorSchemes", "Colors", "Dates", "Printf", "Random", "Reexport", "SnoopPrecompile", "Statistics"] -git-tree-sha1 = "21303256d239f6b484977314674aef4bb1fe4420" +deps = ["ColorSchemes", "Colors", "Dates", "PrecompileTools", "Printf", "Random", "Reexport", "Statistics"] +git-tree-sha1 = "f92e1315dadf8c46561fb9396e525f7200cdc227" uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043" -version = "1.3.1" +version = "1.3.5" [[deps.Plots]] -deps = ["Base64", "Contour", "Dates", "Downloads", "FFMPEG", "FixedPointNumbers", "GR", "JLFzf", "JSON", "LaTeXStrings", "Latexify", "LinearAlgebra", "Measures", "NaNMath", "Pkg", "PlotThemes", "PlotUtils", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "RelocatableFolders", "Requires", "Scratch", "Showoff", "SnoopPrecompile", "SparseArrays", "Statistics", "StatsBase", "UUIDs", "UnicodeFun", "Unzip"] -git-tree-sha1 = "ec23efe47c86da2c00dc5496e59cb3d36bbfce6d" +deps = ["Base64", "Contour", "Dates", "Downloads", "FFMPEG", "FixedPointNumbers", "GR", "JLFzf", "JSON", "LaTeXStrings", "Latexify", "LinearAlgebra", "Measures", "NaNMath", "Pkg", "PlotThemes", "PlotUtils", "PrecompileTools", "Preferences", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "RelocatableFolders", "Requires", "Scratch", "Showoff", "SparseArrays", "Statistics", "StatsBase", "UUIDs", "UnicodeFun", "UnitfulLatexify", "Unzip"] +git-tree-sha1 = "75ca67b2c6512ad2d0c767a7cfc55e75075f8bbc" uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" -version = "1.36.0" +version = "1.38.16" + + [deps.Plots.extensions] + FileIOExt = "FileIO" + GeometryBasicsExt = "GeometryBasics" + IJuliaExt = "IJulia" + ImageInTerminalExt = "ImageInTerminal" + UnitfulExt = "Unitful" + + [deps.Plots.weakdeps] + FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" + GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" + IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a" + ImageInTerminal = "d8c32880-2388-543b-8c61-d9f865259254" + Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [[deps.PlutoUI]] deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "FixedPointNumbers", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "MIMEs", "Markdown", "Random", "Reexport", "URIs", "UUIDs"] -git-tree-sha1 = "efc140104e6d0ae3e7e30d56c98c4a927154d684" +git-tree-sha1 = "b478a748be27bd2f2c73a7690da219d0844db305" uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8" -version = "0.7.48" +version = "0.7.51" + +[[deps.PrecompileTools]] +deps = ["Preferences"] +git-tree-sha1 = "9673d39decc5feece56ef3940e5dafba15ba0f81" +uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a" +version = "1.1.2" [[deps.Preferences]] deps = ["TOML"] -git-tree-sha1 = "47e5f437cc0e7ef2ce8406ce1e7e24d44915f88d" +git-tree-sha1 = "7eb1686b4f04b82f96ed7a4ea5890a4f0c7a09f1" uuid = "21216c6a-2e73-6563-6e65-726566657250" -version = "1.3.0" +version = "1.4.0" [[deps.Printf]] deps = ["Unicode"] @@ -845,21 +914,25 @@ uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" [[deps.Ratios]] deps = ["Requires"] -git-tree-sha1 = "dc84268fe0e3335a62e315a3a7cf2afa7178a734" +git-tree-sha1 = "1342a47bf3260ee108163042310d26f2be5ec90b" uuid = "c84ed2f1-dad5-54f0-aa8e-dbefe2724439" -version = "0.4.3" +version = "0.4.5" +weakdeps = ["FixedPointNumbers"] + + [deps.Ratios.extensions] + RatiosFixedPointNumbersExt = "FixedPointNumbers" [[deps.RecipesBase]] -deps = ["SnoopPrecompile"] -git-tree-sha1 = "d12e612bba40d189cead6ff857ddb67bd2e6a387" +deps = ["PrecompileTools"] +git-tree-sha1 = "5c3d09cc4f31f5fc6af001c250bf1278733100ff" uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" -version = "1.3.1" +version = "1.3.4" [[deps.RecipesPipeline]] -deps = ["Dates", "NaNMath", "PlotUtils", "RecipesBase", "SnoopPrecompile"] -git-tree-sha1 = "a030182cccc5c461386c6f055c36ab8449ef1340" +deps = ["Dates", "NaNMath", "PlotUtils", "PrecompileTools", "RecipesBase"] +git-tree-sha1 = "45cf9fd0ca5839d06ef333c8201714e888486342" uuid = "01d81517-befc-4cb6-b9ec-a95719d0359c" -version = "0.6.10" +version = "0.6.12" [[deps.Reexport]] git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b" @@ -884,9 +957,9 @@ version = "0.7.0" [[deps.Scratch]] deps = ["Dates"] -git-tree-sha1 = "f94f779c94e58bf9ea243e77a37e16d9de9126bd" +git-tree-sha1 = "30449ee12237627992a99d5e30ae63e4d78cd24a" uuid = "6c6a2e73-6563-6170-7368-637461726353" -version = "1.1.1" +version = "1.2.0" [[deps.Serialization]] uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" @@ -906,35 +979,34 @@ git-tree-sha1 = "874e8867b33a00e784c8a7e4b60afe9e037b74e1" uuid = "777ac1f9-54b0-4bf8-805c-2214025038e7" version = "1.1.0" -[[deps.SnoopPrecompile]] -git-tree-sha1 = "f604441450a3c0569830946e5b33b78c928e1a85" -uuid = "66db9d55-30c0-4569-8b51-7e840670fc0c" -version = "1.0.1" - [[deps.Sockets]] uuid = "6462fe0b-24de-5631-8697-dd941f90decc" [[deps.SortingAlgorithms]] deps = ["DataStructures"] -git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508" +git-tree-sha1 = "c60ec5c62180f27efea3ba2908480f8055e17cee" uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c" -version = "1.0.1" +version = "1.1.1" [[deps.SparseArrays]] -deps = ["LinearAlgebra", "Random"] +deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"] uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [[deps.SpecialFunctions]] -deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] -git-tree-sha1 = "d75bda01f8c31ebb72df80a46c88b25d1c79c56d" +deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] +git-tree-sha1 = "ef28127915f4229c971eb43f3fc075dd3fe91880" uuid = "276daf66-3868-5448-9aa4-cd146d93841b" -version = "2.1.7" +version = "2.2.0" +weakdeps = ["ChainRulesCore"] + + [deps.SpecialFunctions.extensions] + SpecialFunctionsChainRulesCoreExt = "ChainRulesCore" [[deps.StaticArrays]] deps = ["LinearAlgebra", "Random", "StaticArraysCore", "Statistics"] -git-tree-sha1 = "f86b3a049e5d05227b10e15dbb315c5b90f14988" +git-tree-sha1 = "832afbae2a45b4ae7e831f86965469a24d1d8a83" uuid = "90137ffa-7385-5640-81b9-e52037218182" -version = "1.5.9" +version = "1.5.26" [[deps.StaticArraysCore]] git-tree-sha1 = "6b7ba252635a5eff6a0b0664a41ee140a1c9e72a" @@ -944,27 +1016,33 @@ version = "1.4.0" [[deps.Statistics]] deps = ["LinearAlgebra", "SparseArrays"] uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +version = "1.9.0" [[deps.StatsAPI]] deps = ["LinearAlgebra"] -git-tree-sha1 = "f9af7f195fb13589dd2e2d57fdb401717d2eb1f6" +git-tree-sha1 = "45a7769a04a3cf80da1c1c7c60caf932e6f4c9f7" uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0" -version = "1.5.0" +version = "1.6.0" [[deps.StatsBase]] deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"] -git-tree-sha1 = "d1bf48bfcc554a3761a133fe3a9bb01488e06916" +git-tree-sha1 = "75ebe04c5bed70b91614d684259b661c9e6274a4" uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" -version = "0.33.21" +version = "0.34.0" [[deps.SuiteSparse]] deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"] uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9" +[[deps.SuiteSparse_jll]] +deps = ["Artifacts", "Libdl", "Pkg", "libblastrampoline_jll"] +uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c" +version = "5.10.1+6" + [[deps.TOML]] deps = ["Dates"] uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" -version = "1.0.0" +version = "1.0.3" [[deps.TableTraits]] deps = ["IteratorInterfaceExtensions"] @@ -974,14 +1052,14 @@ version = "1.0.1" [[deps.Tables]] deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "OrderedCollections", "TableTraits", "Test"] -git-tree-sha1 = "c79322d36826aa2f4fd8ecfa96ddb47b174ac78d" +git-tree-sha1 = "1544b926975372da01227b382066ab70e574a3ec" uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" -version = "1.10.0" +version = "1.10.1" [[deps.Tar]] deps = ["ArgTools", "SHA"] uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" -version = "1.10.1" +version = "1.10.0" [[deps.TensorCore]] deps = ["LinearAlgebra"] @@ -995,19 +1073,19 @@ uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [[deps.TranscodingStreams]] deps = ["Random", "Test"] -git-tree-sha1 = "8a75929dcd3c38611db2f8d08546decb514fcadf" +git-tree-sha1 = "9a6ae7ed916312b41236fcef7e0af564ef934769" uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa" -version = "0.9.9" +version = "0.9.13" [[deps.Tricks]] -git-tree-sha1 = "6bac775f2d42a611cdfcd1fb217ee719630c4175" +git-tree-sha1 = "aadb748be58b492045b4f56166b5188aa63ce549" uuid = "410a4b4d-49e4-4fbc-ab6d-cb71b17b3775" -version = "0.1.6" +version = "0.1.7" [[deps.URIs]] -git-tree-sha1 = "e59ecc5a41b000fa94423a578d29290c7266fc10" +git-tree-sha1 = "074f993b0ca030848b897beff716d93aca60f06a" uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4" -version = "1.4.0" +version = "1.4.2" [[deps.UUIDs]] deps = ["Random", "SHA"] @@ -1022,6 +1100,24 @@ git-tree-sha1 = "53915e50200959667e78a92a418594b428dffddf" uuid = "1cfade01-22cf-5700-b092-accc4b62d6e1" version = "0.4.1" +[[deps.Unitful]] +deps = ["ConstructionBase", "Dates", "LinearAlgebra", "Random"] +git-tree-sha1 = "ba4aa36b2d5c98d6ed1f149da916b3ba46527b2b" +uuid = "1986cc42-f94f-5a68-af5c-568840ba703d" +version = "1.14.0" + + [deps.Unitful.extensions] + InverseFunctionsUnitfulExt = "InverseFunctions" + + [deps.Unitful.weakdeps] + InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" + +[[deps.UnitfulLatexify]] +deps = ["LaTeXStrings", "Latexify", "Unitful"] +git-tree-sha1 = "e2d817cc500e960fdbafcf988ac8436ba3208bfd" +uuid = "45397f5d-5981-4c77-b2b3-fc36d6e9b728" +version = "1.6.3" + [[deps.Unzip]] git-tree-sha1 = "ca0969166a028236229f63514992fc073799bb78" uuid = "41fe7b60-77ed-43a1-b4f0-825fd5a5650d" @@ -1029,9 +1125,9 @@ version = "0.2.0" [[deps.Wayland_jll]] deps = ["Artifacts", "Expat_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg", "XML2_jll"] -git-tree-sha1 = "3e61f0b86f90dacb0bc0e73a0c5a83f6a8636e23" +git-tree-sha1 = "ed8d92d9774b077c53e1da50fd81a36af3744c1c" uuid = "a2964d1f-97da-50d4-b82a-358c7fce9d89" -version = "1.19.0+0" +version = "1.21.0+0" [[deps.Wayland_protocols_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] @@ -1047,9 +1143,9 @@ version = "0.5.5" [[deps.XML2_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"] -git-tree-sha1 = "58443b63fb7e465a8a7210828c91c08b92132dff" +git-tree-sha1 = "93c41695bc1c08c46c5899f4fe06d6ead504bb73" uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a" -version = "2.9.14+0" +version = "2.10.3+0" [[deps.XSLT_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "Pkg", "XML2_jll", "Zlib_jll"] @@ -1185,20 +1281,20 @@ version = "1.4.0+3" [[deps.ZipFile]] deps = ["Libdl", "Printf", "Zlib_jll"] -git-tree-sha1 = "ef4f23ffde3ee95114b461dc667ea4e6906874b2" +git-tree-sha1 = "f492b7fe1698e623024e873244f10d89c95c340a" uuid = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea" -version = "0.10.0" +version = "0.10.1" [[deps.Zlib_jll]] deps = ["Libdl"] uuid = "83775a58-1f1d-513f-b197-d71354ab007a" -version = "1.2.12+3" +version = "1.2.13+0" [[deps.Zstd_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "e45044cd873ded54b6a5bac0eb5c971392cf1927" +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "49ce682769cd5de6c72dcf1b94ed7790cd08974c" uuid = "3161d3a3-bdf6-5164-811a-617609db77b4" -version = "1.5.2+0" +version = "1.5.5+0" [[deps.fzf_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] @@ -1219,9 +1315,9 @@ uuid = "0ac62f75-1d6f-5e53-bd7c-93b484bb37c0" version = "0.15.1+0" [[deps.libblastrampoline_jll]] -deps = ["Artifacts", "Libdl", "OpenBLAS_jll"] +deps = ["Artifacts", "Libdl"] uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" -version = "5.1.1+0" +version = "5.7.0+0" [[deps.libfdk_aac_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] 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 8eb4baa5..a9c947b5 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 @@ -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 diff --git a/src/DIVAnd_aexerr.jl b/src/DIVAnd_aexerr.jl index 9231674d..4397cefd 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] @@ -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)]) 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 diff --git a/src/DIVAnd_bc_stretch.jl b/src/DIVAnd_bc_stretch.jl index 37e7fc71..ce3b3564 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/DIVAnd_diagnostics.jl b/src/DIVAnd_diagnostics.jl new file mode 100644 index 00000000..8869f007 --- /dev/null +++ b/src/DIVAnd_diagnostics.jl @@ -0,0 +1,56 @@ +""" + 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 + + + + + + + + + + + + + + + diff --git a/src/DIVAnd_errormap.jl b/src/DIVAnd_errormap.jl index 8513cbd7..1b3c5d4e 100644 --- a/src/DIVAnd_errormap.jl +++ b/src/DIVAnd_errormap.jl @@ -1,7 +1,7 @@ """ error,method = DIVAnd_errormap(mask,pmn,xi,x,f,len,epsilon2, - s; - method = "auto", + s; + method = :auto, Bscale = false, otherargs...,); @@ -28,16 +28,16 @@ * `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` +*`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. + 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. *`Bscale` : it `true` will try to take out the the boundary effects in the background error variance. Not possible with all methods @@ -60,9 +60,10 @@ function DIVAnd_errormap( len, epsilon2, s; - method = "auto", + 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 @@ -81,8 +82,6 @@ function DIVAnd_errormap( Lowdata = false - if method == "auto" - Lpmnrange = DIVAnd_Lpmnrange(pmn, len) # L compared to domain size @@ -118,6 +117,10 @@ function DIVAnd_errormap( Lowdata = true end + + if method == :auto + + # try to guess # small L @@ -129,20 +132,20 @@ function DIVAnd_errormap( # very high data coverage: scpme if smallL if Lowdata - errmethod = "cpme" + errmethod = :cpme else if Bigdata - errmethod = "scpme" + errmethod = :scpme else - errmethod = "diagapp" + errmethod = :diagapp end end else - if Bigdata - errmethod = "scpme" + if Bigdata + errmethod = :scpme else - errmethod = "aexerr" - end + errmethod = :aexerr + end end @@ -152,142 +155,80 @@ function DIVAnd_errormap( - if method == "cheap" - Lpmnrange = DIVAnd_Lpmnrange(pmn, len) - # L compared to domain size - - LoverLdomain = zeros(Float64, ndims(mask)) - - - for i = 1:ndims(mask) - LoverLdomain[i] = Lpmnrange[i][2] / size(mask)[i] - end - - if sum(LoverLdomain .< LoverLlimit) == ndims(mask) - smallL = true - end + if method == :cheap - - # Now look at lower values to check for data coverage - realdims = ndims(mask) - for i = 1:ndims(mask) - LoverLdomain[i] = Lpmnrange[i][1] / size(mask)[i] - if Lpmnrange[i][1] == 0 - LoverLdomain[i] = 1.0 / size(mask)[i] - realdims = realdims - 1 - end - end - #nbdonnee size of f a revoir en fonction depsilon2 - if prod(LoverLdomain) * size(f)[1] > pointsperbubblelimit^realdims - Bigdata = true - end - - if prod(LoverLdomain) * size(f)[1] < pointsperbubblelimitlow^realdims - Lowdata = true - end if smallL if Lowdata - errmethod = "cpme" + errmethod = :cpme else if Bigdata - errmethod = "scpme" + errmethod = :scpme else - errmethod = "cpme" + errmethod = :cpme end end else - if Bigdata - errmethod = "scpme" + if Bigdata + errmethod = :scpme else - errmethod = "cpme" - end + errmethod = :cpme + end end end - if method == "precise" - Lpmnrange = DIVAnd_Lpmnrange(pmn, len) - # L compared to domain size - - LoverLdomain = zeros(Float64, ndims(mask)) - - - for i = 1:ndims(mask) - LoverLdomain[i] = Lpmnrange[i][2] / size(mask)[i] - end - - if sum(LoverLdomain .< LoverLlimit) == ndims(mask) - smallL = true - end - - - # Now look at lower values to check for data coverage - realdims = ndims(mask) - for i = 1:ndims(mask) - LoverLdomain[i] = Lpmnrange[i][1] / size(mask)[i] - if Lpmnrange[i][1] == 0 - LoverLdomain[i] = 1.0 / size(mask)[i] - realdims = realdims - 1 - end - end - #nbdonnee size of f a revoir en fonction depsilon2 - if prod(LoverLdomain) * size(f)[1] > pointsperbubblelimit^realdims - Bigdata = true - end + if method == :precise - if prod(LoverLdomain) * size(f)[1] < pointsperbubblelimitlow^realdims - Lowdata = true - end if smallL if Lowdata - errmethod = "aexerr" + errmethod = :aexerr else if Bigdata - errmethod = "diagapp" + errmethod = :diagapp else - errmethod = "diagapp" + errmethod = :diagapp end end else - if Bigdata - errmethod = "diagapp" + if Bigdata + errmethod = :diagapp else - errmethod = "aexerr" - end + errmethod = :aexerr + end end end - if errmethod == "cpme" && Bscale - warn("Sorry, that method does not allow rescaling by spatial dependance of B ") + if errmethod == :cpme && Bscale + @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 ") + if errmethod == :scpme && Bscale + @warn "Sorry, that method does not allow rescaling by spatial dependance of B " 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") + @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") + if errmethod == :scpme && noP + @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") + if errmethod == :exact && noP + @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") + if errmethod == :diagapp && noP + @warn "Sorry, that method needs s.P to be available. Will use aexerr instead" errmethod = aexerr end @@ -296,7 +237,7 @@ function DIVAnd_errormap( # @show errmethod, ScalebyB, pointsperbubblelimit - if errmethod == "cpme" + if errmethod == :cpme errormap = DIVAnd_cpme( mask, pmn, @@ -311,7 +252,7 @@ function DIVAnd_errormap( return errormap, errmethod end - if errmethod == "scpme" + if errmethod == :scpme errormap = DIVAnd_cpme( mask, pmn, @@ -324,18 +265,19 @@ function DIVAnd_errormap( ) scpme=deepcopy(errormap) - DIVAnd_scalecpme!(scpme,s.P) + DIVAnd_scalecpme!(scpme,s.P;rng=rng) return scpme, errmethod end - if errmethod == "exact" - errormap =statevector_unpack(s.sv,diag(s.P)) + if errmethod == :exact + + errormap, =statevector_unpack(s.sv,diag(s.P)) return errormap, errmethod end - if errmethod == "diagapp" + if errmethod == :diagapp errormap = DIVAnd_diagapp( s.P, pmn, @@ -345,7 +287,7 @@ function DIVAnd_errormap( return errormap, errmethod end - if errmethod == "aexerr" + if errmethod == :aexerr errormap,bi,c,d =DIVAnd_aexerr( mask, pmn, @@ -354,9 +296,21 @@ function DIVAnd_errormap( f, len, epsilon2; + rng=rng, 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/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_metric.jl b/src/DIVAnd_metric.jl index ae04f733..4f064933 100644 --- a/src/DIVAnd_metric.jl +++ b/src/DIVAnd_metric.jl @@ -74,6 +74,12 @@ function DIVAnd_metric(lon::AbstractVector, lat::AbstractVector) return DIVAnd_metric(ndgrid(lon, lat)...) end +""" + dy = DIVAnd.deg2m(dlat) + +Convert an increment of latitude to a increment in meters +using the the mean Earth radius. +""" function deg2m(dlat) # Mean radius (http://en.wikipedia.org/wiki/Earth_radius) R = 6371.009e3 diff --git a/src/DIVAnd_multivarEOF.jl b/src/DIVAnd_multivarEOF.jl index df59bb9e..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 @@ -134,7 +134,7 @@ function DIVAnd_multivarEOF(mask,pmn,xi,x,f,lenin,epsilon2in;eof=(),velocity=(), # univariate analysis fi,s=DIVAndrun(mask,pmn,xi,x,f,len,epsilon2;velocity=velocity,kwargs...) - emap,meth=DIVAnd_errormap(mask,pmn,xi,x,f,len,epsilon2,s;method="cheap",velocity=velocity,kwargs...) + emap,meth=DIVAnd_errormap(mask,pmn,xi,x,f,len,epsilon2,s;method=:cheap,velocity=velocity,kwargs...) @debug "error method in multivar $meth" # keep points where error is small enough limite=0.5 @@ -206,7 +206,7 @@ function DIVAnd_multivarEOF(mask,pmn,xi,x,f,lenin,epsilon2in;eof=(),velocity=(), for k=1:size(eof)[1] if abs(eof[k])<0.1/size(eof)[1] @warn "Sorry too weak correlations layer $k" - return fi,s,eof,eofamplitudes,emap + return fi,s,eof,eofamplitudes,emap,emap end end @@ -285,14 +285,14 @@ 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 fm=repeat(f,nlay) # Error maps for the multivariate approach - emapm,methm=DIVAnd_errormap(mask,pmn,xi,tuple(xm...),fm,len,epsilon2m,s;method="cheap",velocity=velocity,kwargs...) + emapm,methm=DIVAnd_errormap(mask,pmn,xi,tuple(xm...),fm,len,epsilon2m,s;method=:cheap,velocity=velocity,kwargs...) @debug methm # Banzaii, finished return fi,s,eof,eofamplitudes,emap,emapm diff --git a/src/DIVAnd_multivarJAC.jl b/src/DIVAnd_multivarJAC.jl index f8fae4ed..1318c84d 100644 --- a/src/DIVAnd_multivarJAC.jl +++ b/src/DIVAnd_multivarJAC.jl @@ -10,14 +10,14 @@ 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: +# 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,71 +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) @@ -147,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,:] @@ -156,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] @@ -167,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] @@ -182,30 +183,30 @@ 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 pseudovelocity[1]=reshape(uu,size(pseudovelocity[1])) pseudovelocity[2]=reshape(vv,size(pseudovelocity[2])) # - vconstrain=DIVAnd_constr_advec(s, tuple(pseudovelocity...)) - # pass the corresponding constraint to divandrun + vconstrain=DIVAnd_constr_advec(sori, tuple(pseudovelocity...)) + # 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...) + 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 92c25226..440f8787 100644 --- a/src/DIVAnd_scalecpme!.jl +++ b/src/DIVAnd_scalecpme!.jl @@ -1,11 +1,11 @@ -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 + 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 @@ -14,16 +14,16 @@ function DIVAnd_scalecpme!(cpme, P::CovarIS, nsamples = 7) 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 e361bd7d..1a4c047f 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...) - + 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,47 +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 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/DIVAndjog.jl b/src/DIVAndjog.jl index 54cb95cb..ab52ab29 100644 --- a/src/DIVAndjog.jl +++ b/src/DIVAndjog.jl @@ -1248,8 +1248,8 @@ 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 diff --git a/src/DIVAndrun.jl b/src/DIVAndrun.jl index 9dd575f9..0f815edf 100644 --- a/src/DIVAndrun.jl +++ b/src/DIVAndrun.jl @@ -15,6 +15,8 @@ function DIVAndrun( maxit = 100, minit = 0, constraints = (), + ineqconstraints = (), + ntriesmax = 10, inversion = :chol, moddim = [], fracindex = Matrix{T}(undef, 0, 0), @@ -32,13 +34,112 @@ function DIVAndrun( fluxes = (), epsfluxes = 0, epsilon2forfractions = 0, - RTIMESONESCALES = (), + RTIMESONESCALES = (), QCMETHOD = (), coeff_laplacian::Vector{Float64} = ones(ndims(mask)), coeff_derivative2::Vector{Float64} = zeros(ndims(mask)), 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 +316,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 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 @@ -229,8 +332,8 @@ For oceanographic application, this is the land-sea mask where sea is true and l * `inversion`: direct solver (`:chol` for Cholesky factorization), an interative solver (`:pcg` for preconditioned conjugate gradient [1]) can be - used or `:cg_amg_sa` for a multigrid method with preconditioned conjugate - gradient. The two last methods are iterative methods who a controlled by + used or `:cg_amg_sa` for a multigrid method with preconditioned conjugate + gradient. The two last methods are iterative methods who a controlled by the number of iterations `maxit` and the tolerance `tol`. * `compPC`: function that returns a preconditioner for the primal formulation diff --git a/src/NCODV.jl b/src/NCODV.jl index 163a9b4a..7f2efb0b 100644 --- a/src/NCODV.jl +++ b/src/NCODV.jl @@ -77,7 +77,10 @@ function loadprof( fillval_time, accepted_status_flag_values_time; nchunk = 10 -) where {T,Tz} + ) where {T,Tz} + + @debug "accepted_status_flag_values_time: $accepted_status_flag_values_time, time name $(name(nctime))" + n_samples = size(ncvar,1) n_stations = size(ncvar, 2) # n_stations = 100000 @@ -295,8 +298,9 @@ end Load all profiles in the file `fname` corresponding to netCDF variable with the `long_name` attribute equal to the parameter `long_name`. `qv_flags` is a list of strings -with the quality flags to be kept. `obsids` is a vector of strings with the -EDMO code and local CDI id concatenated by a hyphen. +with the quality flags to be kept. The filtering of the quality flags is applied +to the data variables, time and depth coordinates. `obsids` is a vector of +strings with the EDMO code and local CDI id concatenated by a hyphen. `nchunk` is the number of profiles read at a time. Large values of `nchunk` can increase performance but requirer also more memory. diff --git a/src/SDNMetadata.jl b/src/SDNMetadata.jl index da3b3f44..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"), ), @@ -137,6 +137,8 @@ function SDNMetadata( ], "-", ) + elseif k == "project" + v = project["name"] end if typeof(v) <: Vector diff --git a/src/derived.jl b/src/derived.jl index 99a707c5..7170b539 100644 --- a/src/derived.jl +++ b/src/derived.jl @@ -57,7 +57,7 @@ end ) Compute derived quantities from a DIVAnd analyse (in particular the deepest -value of an analysis) using the NetCDF file `filename` with the variable +value of an analysis) using the NetCDF file `filename` with the variable `varname`. The result will be written in `new_filename`. See `DIVAnd.diva3d` for the optional parameter `error_thresholds`. diff --git a/src/diva.jl b/src/diva.jl index 10235a98..fc42344b 100644 --- a/src/diva.jl +++ b/src/diva.jl @@ -7,6 +7,7 @@ the regular grid defined by the vectors `xi` using the scaled observational erro variance `epsilon2` and the correlation length `len`. The result will be saved in the netCDF file `filename` under the variable `varname`. + ## Inputs * `xi`: tuple with n elements. Every element represents a coordinate @@ -128,6 +129,7 @@ end !!! note At all vertical levels, there should at least one sea point. + The function assumes a spherical Earth with a radius equal to the mean Earth radius. """ function diva3d( xi, diff --git a/src/domain.jl b/src/domain.jl index eb07433e..35008312 100644 --- a/src/domain.jl +++ b/src/domain.jl @@ -73,7 +73,9 @@ end mask,(pm,pn),(xi,yi) = DIVAnd.domain(bathname,bathisglobal,lonr,latr) Generate a 2D geospatial domain based on the topography from the netCDF file -`bathname`. +`bathname`. `lonr` and `latr` are expressed in degrees East and North +respectively (as are `xi` and `yi`). `pm` and `pn` are in m⁻¹ (using the the +mean Earth radius). """ function domain(bathname, bathisglobal, lonr, latr) mask, (pm, pn), (xi, yi) = DIVAnd.DIVAnd_rectdom(lonr, latr) @@ -86,12 +88,16 @@ function domain(bathname, bathisglobal, lonr, latr) end """ - mask,(pm,pn,po),(xi,yi,zi) = DIVAnd.domain(bathname,bathisglobal,lonr,latr,depthr) + mask,(pm,pn,po),(xi,yi,zi) = DIVAnd.domain(bathname,bathisglobal,lonr,latr,depthr; zlevel = :surface) Generate a 3D geospatial domain based on the topography from the netCDF file `bathname`. If `zlevel` is `:surface`, then `depthr` is zero for the sea surface and positive in water (positive is down). If `zlevel` is `:floor`, then `depthr` is -zero for the sea floor and positive in water (positive is up) +zero for the sea floor and positive in water (positive is up). + +`lonr` and `latr` are expressed in degrees East and North +respectively (as are `xi` and `yi`). `depthr` and `zi` are in meters and +`pm`, `pn` and `po` are in m⁻¹ (using the the mean Earth radius). """ function domain(bathname, bathisglobal, lonr, latr, depthr; zlevel = :surface) diff --git a/src/fit.jl b/src/fit.jl index e607ed15..79c637df 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" diff --git a/src/special_matrices.jl b/src/special_matrices.jl index c4a610be..48915bd6 100644 --- a/src/special_matrices.jl +++ b/src/special_matrices.jl @@ -30,11 +30,11 @@ 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 - x,convergence_history = cg( + output = cg( C.IS, v, Pl = C.factors, verbose = C.verbose, log = log, @@ -42,12 +42,17 @@ function Base.:*(C::CovarIS, v::TV)::TV where {TV<:AbstractVector{Float64}} reltol = C.reltol, maxiter = C.maxiter) + if log + x,convergence_history = output + else + x = output + end @debug "Number of iterations: $(convergence_history.iters)" @debug "Final norm of residue: $(convergence_history.data[:resnorm][end])" @debug begin @show norm(C.IS * x - v) end - #@show convergence_history + @show convergence_history return x elseif C.factors != nothing return C.factors \ v diff --git a/src/statevector.jl b/src/statevector.jl index d2dbea9a..cce3f9e8 100644 --- a/src/statevector.jl +++ b/src/statevector.jl @@ -22,7 +22,7 @@ end """ - sv = statevector_init((mask1, mask2, ...)) + sv = statevector((mask1, mask2, ...)) Initialize structure for packing and unpacking multiple variables given their corresponding land-sea mask. @@ -32,16 +32,12 @@ Input: Every mask can have a different shape. Output: - sv: structure to be used with statevector_pack and statevector_unpack. + sv: structure to be used with `pack` and `unpack`. Note: -see also statevector_pack, statevector_unpack - -Author: Alexander Barth, 2009,2017 -License: GPL 2 or later +see also `pack`, `unpack` """ function statevector(masks::NTuple{nvar_,BitArray{N}}) where {nvar_} where {N} - numels = [sum(mask) for mask in masks] ind = [0, cumsum(numels)...] diff --git a/src/varanalysis.jl b/src/varanalysis.jl index 192d3ecb..1e7748e7 100644 --- a/src/varanalysis.jl +++ b/src/varanalysis.jl @@ -18,8 +18,29 @@ function diffusion!(ivol, nus, α, nmax, x0, x, end +""" + x = DIVAnd.diffusion(mask,pmn,len,x0; boundary_condition! = nothing) + +Diffuse the field `x0` (n-dimensional array) defined in a domain with a mask `mask` and +metric `pmn` iteratively upto a length-scale of `len`. The units of `pmn` +must be the inverse of the units of `len`. The optional function +`boundary_condition!(x)` is applied after every interation to apply the +boundary condition. This function expected to modify its argument. + +## Example +```julia +using DIVAnd +mask, (pm, pn), (xi, yi) = DIVAnd_squaredom(2, range(-1, stop = 1, length = 30)) +f = zeros(size(mask)) +f[15,15] = 1 +len = 0.3 +f = DIVAnd.diffusion(mask,(pm,pn),len,f) +``` + +See also DIVAndrun for more information about `mask`, `pmn` and `len`. +""" function diffusion(mask,pmn,len_,x0; boundary_condition! = nothing) n = ndims(mask) len = DIVAnd.len_harmonize(len_, mask) diff --git a/test/runtests.jl b/test/runtests.jl index f7ff9688..459c79cc 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -30,6 +30,13 @@ using Downloads: download include("test_1dvar.jl") include("test_1D_seminormed.jl") + + # multivar + + include("test_multivar.jl") + + include("test_errormaps.jl") + # dynamical constraints include("test_2dvar_adv.jl") include("test_2dvar_constcoast.jl") @@ -127,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 212c19f7..900b2f81 100644 --- a/test/test_2dvar_adv.jl +++ b/test/test_2dvar_adv.jl @@ -29,7 +29,37 @@ 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 + +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 # diff --git a/test/test_errormaps.jl b/test/test_errormaps.jl new file mode 100644 index 00000000..786b634f --- /dev/null +++ b/test/test_errormaps.jl @@ -0,0 +1,88 @@ +# Testing DIVAnd in 2 dimensions with advection. + +using Test +using StableRNGs +using Random +using DIVAnd + +rng = StableRNG(1234) + +# 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.0 +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(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) + + +fi11,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, + rng=StableRNG(123) + ) + + #@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 + +end +end +end +end +@show errmean +@test abs(errmean - 7.527185124068747) < 0.00000001 + + +# 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 . 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 . diff --git a/test/test_product.jl b/test/test_product.jl index 62712f62..f4264efd 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