text:
### A Pluto.jl notebook ###
# v0.16.1

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 60941eaa-1aea-11eb-1277-97b991548781
begin
    using PlutoUI,GridVisualize,SimplexGridFactory,Triangulate,ExtendableGrids,VoronoiFVM,PlutoVista,PlutoUI,PyPlot
end

# ╔═╡ de468cb9-b34d-4d2e-b911-9b93920caca1
md"""
# Flux reconstruction and visualization for the Laplace operator
"""

# ╔═╡ 928a70c5-4706-40a1-9387-abcb71c09443
"""
Define a "Swiss cheese domain" with punched-out holes, where each hole boundary corresponds to a different boundary condition.
"""
function swiss_cheese_2d()

    function circlehole!(builder, center, radius; n=20)
        points=[point!(builder, center[1]+radius*sin(t),center[2]+radius*cos(t)) for t in range(0,2π,length=n)]
        for i=1:n-1
            facet!(builder,points[i],points[i+1])
        end
        facet!(builder,points[end],points[1])
        holepoint!(builder,center)
    end

    builder=SimplexGridBuilder(Generator=Triangulate)
    cellregion!(builder,1)
    maxvolume!(builder,0.1)
    regionpoint!(builder,0.1,0.1)

    p1=point!(builder,0,0)
    p2=point!(builder,10,0)
    p3=point!(builder,10,10)
    p4=point!(builder,0,10)

    facetregion!(builder,1)
    facet!(builder,p1,p2)
    facet!(builder,p2,p3)
    facet!(builder,p3,p4)
    facet!(builder,p4,p1)

    holes=[1.0 2.0;
           8.0 9.0;
           2.0 8.0;
  		   8.0 4.0;
           9.0 1.0;
           3.0 4.0;
           4.0 6.0;
           7.0 9.0;
           4.0 7.0;
           7.0 5.0;
           2.0 1.0;
           4.0 1.0;
           4.0 8.0;
           3.0 6.0;
           4.0 9.0;
           6.0 9.0;
           3.0 5.0;
           1.0 4.0]'

    radii=[0.15, 0.15, 0.1, 0.35, 0.2, 0.3, 0.1, 0.4, 0.1, 0.4,  0.2, 0.2, 0.2, 0.35, 0.15, 0.25, 0.15, 0.25]

    for i=1:length(radii)
        facetregion!(builder,i+1)
        circlehole!(builder,holes[:,i], radii[i])
    end

    simplexgrid(builder)
end

# ╔═╡ 49f425a4-455c-40f5-b939-da704032cb88
md"""
Generate & plot grid
"""

# ╔═╡ bc304085-69c3-4974-beb4-6f2b981ac0f1
grid=swiss_cheese_2d()

# ╔═╡ ad0367e6-7308-46bc-a7ff-9c17dcfe470d
gridplot(grid,Plotter=PlutoVista)

# ╔═╡ 236d6cd5-c190-49c7-98fe-021fae224455
"""
Simple flux function for Laplace operator
"""
flux(y,u,edge)= y[1]=u[1,1]-u[1,2]

# ╔═╡ f4ebe6ad-4e04-4f33-9a66-6bec977adf4d
md"""
Define a finite volume system with Dirichlet boundary conditions at some of the holes
"""

# ╔═╡ 934a98f1-9734-45bc-848e-4d4d8ad80e4e
physics=VoronoiFVM.Physics(flux=flux)

# ╔═╡ f3af20fb-c3cc-41af-9782-d40adf2371f7
system=VoronoiFVM.System(grid,physics)

# ╔═╡ 50ea0a7e-97e0-428d-b4ec-6bd4d7aef31f
enable_species!(system,1,[1])

# ╔═╡ 49d3c96b-f6b4-4581-8600-8aaf50bb07fc
boundary_dirichlet!(system,1,2,10.0)

# ╔═╡ fa5200ff-713e-48fb-b452-c6ee96935cfc
boundary_dirichlet!(system,1,3,0.0)

# ╔═╡ 6b09f75a-521b-40d9-9636-155e62ab0cac
md"""
At hole #11, the value will be bound to a slider defined below
"""

# ╔═╡ d86c43f3-ec2f-4fae-88d2-1068603e7044
md"""
Solve, and trigger solution upon boundary value change
"""

# ╔═╡ 41bd1230-c87c-47b0-8e58-67ad55609fd3
md"""
Reconstruct the node flux. It is a ``d\times n_{spec}\times n_{nodes}`` tensor.
`nf[:,ispec,:]` then is a vector function representing the flux density of species `ispec` in each node of the domain. This readily can be fed into `GridVisualize.vectorplot`.
"""

# ╔═╡ 17be52fb-f55b-4b3d-85e5-33f36134046b
vis=GridVisualizer(Plotter=PlutoVista,dim=2,resolution=(400,400));vis

# ╔═╡ 03f582ec-4e95-4ca4-8482-9c797027810d
md"""
``v_{11}:`` $(@bind  val11 Slider(0:0.1:10,default=5,show_value=true))
"""

# ╔═╡ e8ae2b22-ac60-449b-952f-d2c55852677c
boundary_dirichlet!(system,1,11,val11)

# ╔═╡ 27efdcac-8ee4-47e4-905d-8d8c7313ddd1
val11; sol=solve(unknowns(system,inival=0),system)

# ╔═╡ 41f427c1-b6ad-46d4-9151-1d872b4efeb6
nf=nodeflux(system,sol)

# ╔═╡ c9b9fdb1-5734-430f-949a-6e6126d2f091
"""
Joint plot of solution and nodal flux of the solution
"""
function myplot(vis)
	scalarplot!(vis,grid,sol[1,:],levels=9,colormap=:summer,clear=true)
	vectorplot!(vis,grid,nf[:,1,:],clear=false,spacing=0.5,vscale=1.5)
	reveal(vis)
end

# ╔═╡ 531edb71-6d32-4231-b117-5e36416d2fb1
myplot(vis)

# ╔═╡ e0b92e00-dd0b-448b-baae-5cd4e0a2059a
md"""
Plot with PyPlot: $(@bind use_pyplot CheckBox(default=false))
"""

# ╔═╡ 1495eb67-12ec-4662-b0e7-049ffef569c0
if use_pyplot
	pyvis=GridVisualizer(Plotter=PyPlot,resolution=(400,400))
	myplot(pyvis)
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
ExtendableGrids = "cfc395e8-590f-11e8-1f13-43a2532b2fa8"
GridVisualize = "5eed8a63-0fb0-45eb-886d-8d5a387d12b8"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
PlutoVista = "646e1f28-b900-46d7-9d87-d554eb38a413"
PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
SimplexGridFactory = "57bfcd06-606e-45d6-baf4-4ba06da0efd5"
Triangulate = "f7e6ffb2-c36d-4f8f-a77e-16e897189344"
VoronoiFVM = "82b139dc-5afc-11e9-35da-9b9bdfd336f3"

[compat]
ExtendableGrids = "~0.8.3"
GridVisualize = "~0.3.5"
PlutoUI = "~0.7.16"
PlutoVista = "~0.8.2"
PyPlot = "~2.10.0"
SimplexGridFactory = "~0.5.8"
Triangulate = "~2.1.0"
VoronoiFVM = "~0.13.1"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.0-rc1"
manifest_format = "2.0"

[[deps.AbstractTrees]]
git-tree-sha1 = "03e0550477d86222521d254b741d470ba17ea0b5"
uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
version = "0.3.4"

[[deps.Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "84918055d15b3114ede17ac6a7182f68870c16f7"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.1"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.ArnoldiMethod]]
deps = ["LinearAlgebra", "Random", "StaticArrays"]
git-tree-sha1 = "f87e559f87a45bece9c9ed97458d3afe98b1ebb9"
uuid = "ec485272-7323-5ecc-a04f-4719b315124d"
version = "0.1.0"

[[deps.ArrayInterface]]
deps = ["Compat", "IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"]
git-tree-sha1 = "b8d49c34c3da35f220e7295659cd0bab8e739fed"
uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
version = "3.1.33"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.Cassette]]
git-tree-sha1 = "6ce3cd755d4130d43bab24ea5181e77b89b51839"
uuid = "7057c7e9-c182-5462-911a-8362d720325c"
version = "0.3.9"

[[deps.ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "2f294fae04aa5069a67964a3366e151e09ea7c09"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.9.0"

[[deps.ColorSchemes]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random"]
git-tree-sha1 = "a851fec56cb73cfdf43762999ec72eff5b86882a"
uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
version = "3.15.0"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[deps.CommonSubexpressions]]
deps = ["MacroTools", "Test"]
git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7"
uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
version = "0.3.0"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "31d0151f5716b655421d9d75b7fa74cc4e744df2"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.39.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.Conda]]
deps = ["JSON", "VersionParsing"]
git-tree-sha1 = "299304989a5e6473d985212c28928899c74e9421"
uuid = "8f4d0f93-b110-5947-807f-2305c1781a2d"
version = "1.5.2"

[[deps.DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

[[deps.DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.DiffResults]]
deps = ["StaticArrays"]
git-tree-sha1 = "c18e98cba888c6c25d1c3b048e4b3380ca956805"
uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
version = "1.0.3"

[[deps.DiffRules]]
deps = ["NaNMath", "Random", "SpecialFunctions"]
git-tree-sha1 = "7220bc21c33e990c14f4a9a319b1d242ebc5b269"
uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
version = "1.3.1"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "a32185f5428d3986f47c2ab78b1f216d5e6cc96f"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.5"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.EarCut_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "3f3a2501fa7236e9b911e0f7a588c657e822bb6d"
uuid = "5ae413db-bbd1-5e63-b57d-d24a61df00f5"
version = "2.2.3+0"

[[deps.ElasticArrays]]
deps = ["Adapt"]
git-tree-sha1 = "a0fcc1bb3c9ceaf07e1d0529c9806ce94be6adf9"
uuid = "fdbdab4c-e67f-52f5-8c3f-e7b388dad3d4"
version = "1.2.9"

[[deps.ExtendableGrids]]
deps = ["AbstractTrees", "Dates", "DocStringExtensions", "ElasticArrays", "InteractiveUtils", "LinearAlgebra", "Printf", "Random", "SparseArrays", "Test"]
git-tree-sha1 = "85dbe70afc7153ad510577f158214b45be22593b"
uuid = "cfc395e8-590f-11e8-1f13-43a2532b2fa8"
version = "0.8.3"

[[deps.ExtendableSparse]]
deps = ["DocStringExtensions", "LinearAlgebra", "Printf", "Requires", "SparseArrays", "SuiteSparse", "Test"]
git-tree-sha1 = "0341e41e45e6c5e46be89c33543082a0a867707c"
uuid = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
version = "0.6.5"

[[deps.FileIO]]
deps = ["Pkg", "Requires", "UUIDs"]
git-tree-sha1 = "3c041d2ac0a52a12a27af2782b34900d9c3ee68c"
uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
version = "1.11.1"

[[deps.FillArrays]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
git-tree-sha1 = "8756f9935b7ccc9064c6eef0bff0ad643df733a3"
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
version = "0.12.7"

[[deps.FiniteDiff]]
deps = ["ArrayInterface", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays"]
git-tree-sha1 = "8b3c09b56acaf3c0e581c66638b85c8650ee9dca"
uuid = "6a86dc24-6348-571c-b903-95158fe2bd41"
version = "2.8.1"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.ForwardDiff]]
deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions", "StaticArrays"]
git-tree-sha1 = "63777916efbcb0ab6173d09a658fb7f2783de485"
uuid = "f6369f11-7733-5829-9624-2563aa707210"
version = "0.10.21"

[[deps.GeometryBasics]]
deps = ["EarCut_jll", "IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"]
git-tree-sha1 = "58bcdf5ebc057b085e58d95c138725628dd7453c"
uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
version = "0.4.1"

[[deps.GridVisualize]]
deps = ["ColorSchemes", "Colors", "DocStringExtensions", "ElasticArrays", "ExtendableGrids", "GeometryBasics", "LinearAlgebra", "Observables", "OrderedCollections", "PkgVersion", "Printf", "Requires", "StaticArrays"]
git-tree-sha1 = "063aa184cec6aff22f7c2f6a498075a93b5fa304"
uuid = "5eed8a63-0fb0-45eb-886d-8d5a387d12b8"
version = "0.3.5"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "f6532909bf3d40b308a0f360b6a0e626c0e263a8"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.1"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.IfElse]]
git-tree-sha1 = "28e837ff3e7a6c3cdb252ce49fb412c8eb3caeef"
uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
version = "0.1.0"

[[deps.Inflate]]
git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c"
uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
version = "0.1.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

[[deps.IterTools]]
git-tree-sha1 = "05110a2ab1fc5f932622ffea2a003221f4782c18"
uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
version = "1.3.0"

[[deps.IterativeSolvers]]
deps = ["LinearAlgebra", "Printf", "Random", "RecipesBase", "SparseArrays"]
git-tree-sha1 = "1a8c6237e78b714e901e406c096fc8a65528af7d"
uuid = "42fd0dbc-a981-5370-80f2-aaf504508153"
version = "0.9.1"

[[deps.IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[deps.JLD2]]
deps = ["DataStructures", "FileIO", "MacroTools", "Mmap", "Pkg", "Printf", "Reexport", "TranscodingStreams", "UUIDs"]
git-tree-sha1 = "46b7834ec8165c541b0b5d1c8ba63ec940723ffb"
uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
version = "0.4.15"

[[deps.JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[deps.LaTeXStrings]]
git-tree-sha1 = "c7f1c695e06c01b95a67f0cd1d34994f3e7db104"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.2.1"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LightGraphs]]
deps = ["ArnoldiMethod", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"]
git-tree-sha1 = "432428df5f360964040ed60418dd5601ecd240b6"
uuid = "093fc24a-ae57-5d10-9952-331d41423f4d"
version = "1.3.5"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.LogExpFunctions]]
deps = ["ChainRulesCore", "DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "34dc30f868e368f8a17b728a1238f3fcda43931a"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.3"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "5a5bc6bf062f0f95e62d0fe0a2d99699fed82dd9"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.8"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NaNMath]]
git-tree-sha1 = "bfe47e760d60b82b66b61d2d44128b62e3a369fb"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "0.3.5"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.Observables]]
git-tree-sha1 = "fe29afdef3d0c4a8286128d4e45cc50621b1e43d"
uuid = "510215fc-4207-5dde-b226-833fc4488ee2"
version = "0.4.0"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"

[[deps.OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "98f59ff3639b3d9485a03a72f3ab35bab9465720"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.6"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PkgVersion]]
deps = ["Pkg"]
git-tree-sha1 = "a7a7e1a88853564e551e4eba8650f8c38df79b37"
uuid = "eebad327-c553-4316-9ea0-9fa01ccd7688"
version = "0.1.1"

[[deps.PlutoUI]]
deps = ["Base64", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "4c8a7d080daca18545c56f1cac28710c362478f3"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.16"

[[deps.PlutoVista]]
deps = ["ColorSchemes", "Colors", "DocStringExtensions", "GridVisualize", "UUIDs"]
git-tree-sha1 = "ce59550b517ab41b5e687174f93f57c76f872b6d"
uuid = "646e1f28-b900-46d7-9d87-d554eb38a413"
version = "0.8.2"

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.PyCall]]
deps = ["Conda", "Dates", "Libdl", "LinearAlgebra", "MacroTools", "Serialization", "VersionParsing"]
git-tree-sha1 = "169bb8ea6b1b143c5cf57df6d34d022a7b60c6db"
uuid = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
version = "1.92.3"

[[deps.PyPlot]]
deps = ["Colors", "LaTeXStrings", "PyCall", "Sockets", "Test", "VersionParsing"]
git-tree-sha1 = "14c1b795b9d764e1784713941e787e1384268103"
uuid = "d330b81b-6aea-500a-939a-2ce795aea3ee"
version = "2.10.0"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.RecipesBase]]
git-tree-sha1 = "44a75aa7a527910ee3d1751d1f0e4148698add9e"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.1.2"

[[deps.RecursiveArrayTools]]
deps = ["ArrayInterface", "ChainRulesCore", "DocStringExtensions", "FillArrays", "LinearAlgebra", "RecipesBase", "Requires", "StaticArrays", "Statistics", "ZygoteRules"]
git-tree-sha1 = "ff7495c78a192ff7d59531d9f14db300c847a4bc"
uuid = "731186ca-8d62-57ce-b412-fbd966d074cd"
version = "2.19.1"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.SimpleTraits]]
deps = ["InteractiveUtils", "MacroTools"]
git-tree-sha1 = "5d7e3f4e11935503d3ecaf7186eac40602e7d231"
uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
version = "0.9.4"

[[deps.SimplexGridFactory]]
deps = ["DocStringExtensions", "ElasticArrays", "ExtendableGrids", "GridVisualize", "LinearAlgebra", "Printf", "Test"]
git-tree-sha1 = "e55646af1d0f68763eaca28f666d381982eeb65f"
uuid = "57bfcd06-606e-45d6-baf4-4ba06da0efd5"
version = "0.5.8"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.SparseDiffTools]]
deps = ["Adapt", "ArrayInterface", "Compat", "DataStructures", "FiniteDiff", "ForwardDiff", "LightGraphs", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays", "VertexSafeGraphs"]
git-tree-sha1 = "36a4d27a02af48a1eafd2baff58b32deeeb68926"
uuid = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
version = "1.16.5"

[[deps.SparsityDetection]]
deps = ["Cassette", "DocStringExtensions", "LinearAlgebra", "SparseArrays", "SpecialFunctions"]
git-tree-sha1 = "9e182a311d169cb9fe0c6501aa252983215fe692"
uuid = "684fba80-ace3-11e9-3d08-3bc7ed6f96df"
version = "0.3.4"

[[deps.SpecialFunctions]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "793793f1df98e3d7d554b65a107e9c9a6399a6ed"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "1.7.0"

[[deps.Static]]
deps = ["IfElse"]
git-tree-sha1 = "a8f30abc7c64a39d389680b74e749cf33f872a70"
uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
version = "0.3.3"

[[deps.StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "3c76dde64d03699e074ac02eb2e8ba8254d428da"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.2.13"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.StructArrays]]
deps = ["Adapt", "DataAPI", "StaticArrays", "Tables"]
git-tree-sha1 = "2ce41e0d042c60ecd131e9fb7154a3bfadbf50d3"
uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
version = "0.6.3"

[[deps.SuiteSparse]]
deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[deps.Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "fed34d0e71b91734bf0a7e10eb1bb05296ddbcd0"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.6.0"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.6"

[[deps.Triangle_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "bfdd9ef1004eb9d407af935a6f36a4e0af711369"
uuid = "5639c1d2-226c-5e70-8d55-b3095415a16a"
version = "1.6.1+0"

[[deps.Triangulate]]
deps = ["DocStringExtensions", "Libdl", "Printf", "Test", "Triangle_jll"]
git-tree-sha1 = "2b4f716b192c0c615d96d541ee029e85666388cb"
uuid = "f7e6ffb2-c36d-4f8f-a77e-16e897189344"
version = "2.1.0"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.VersionParsing]]
git-tree-sha1 = "80229be1f670524750d905f8fc8148e5a8c4537f"
uuid = "81def892-9a0e-5fdd-b105-ffc91e053289"
version = "1.2.0"

[[deps.VertexSafeGraphs]]
deps = ["LightGraphs"]
git-tree-sha1 = "b9b450c99a3ca1cc1c6836f560d8d887bcbe356e"
uuid = "19fa3120-7c27-5ec5-8db8-b0b0aa330d6f"
version = "0.1.2"

[[deps.VoronoiFVM]]
deps = ["DiffResults", "DocStringExtensions", "ExtendableGrids", "ExtendableSparse", "ForwardDiff", "GridVisualize", "IterativeSolvers", "JLD2", "LinearAlgebra", "Printf", "RecursiveArrayTools", "SparseArrays", "SparseDiffTools", "SparsityDetection", "StaticArrays", "SuiteSparse", "Test"]
git-tree-sha1 = "db381645bec1d798c560fe415593663d4a2194ac"
uuid = "82b139dc-5afc-11e9-35da-9b9bdfd336f3"
version = "0.13.1"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.ZygoteRules]]
deps = ["MacroTools"]
git-tree-sha1 = "8c1a8e4dfacb1fd631745552c8db35d0deb09ea0"
uuid = "700de1a5-db45-46bc-99cf-38207098b444"
version = "0.2.2"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─de468cb9-b34d-4d2e-b911-9b93920caca1
# ╠═60941eaa-1aea-11eb-1277-97b991548781
# ╠═928a70c5-4706-40a1-9387-abcb71c09443
# ╟─49f425a4-455c-40f5-b939-da704032cb88
# ╠═bc304085-69c3-4974-beb4-6f2b981ac0f1
# ╠═ad0367e6-7308-46bc-a7ff-9c17dcfe470d
# ╠═236d6cd5-c190-49c7-98fe-021fae224455
# ╟─f4ebe6ad-4e04-4f33-9a66-6bec977adf4d
# ╠═934a98f1-9734-45bc-848e-4d4d8ad80e4e
# ╠═f3af20fb-c3cc-41af-9782-d40adf2371f7
# ╠═50ea0a7e-97e0-428d-b4ec-6bd4d7aef31f
# ╠═49d3c96b-f6b4-4581-8600-8aaf50bb07fc
# ╠═fa5200ff-713e-48fb-b452-c6ee96935cfc
# ╟─6b09f75a-521b-40d9-9636-155e62ab0cac
# ╠═e8ae2b22-ac60-449b-952f-d2c55852677c
# ╟─d86c43f3-ec2f-4fae-88d2-1068603e7044
# ╠═27efdcac-8ee4-47e4-905d-8d8c7313ddd1
# ╟─41bd1230-c87c-47b0-8e58-67ad55609fd3
# ╠═41f427c1-b6ad-46d4-9151-1d872b4efeb6
# ╠═c9b9fdb1-5734-430f-949a-6e6126d2f091
# ╠═17be52fb-f55b-4b3d-85e5-33f36134046b
# ╟─03f582ec-4e95-4ca4-8482-9c797027810d
# ╠═531edb71-6d32-4231-b117-5e36416d2fb1
# ╟─e0b92e00-dd0b-448b-baae-5cd4e0a2059a
# ╠═1495eb67-12ec-4662-b0e7-049ffef569c0
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using Clustering
using Distances: SqEuclidean, pairwise
using Test

@testset "seeding" begin

Random.seed!(34568)

@test RandSeedAlg <: SeedingAlgorithm
@test KmppAlg <: SeedingAlgorithm
@test KmCentralityAlg <: SeedingAlgorithm

alldistinct(x::Vector{Int}) = (length(Set(x)) == length(x))

function min_interdist(X::AbstractMatrix)
    dists = pairwise(SqEuclidean(), X, dims=2)
    n = size(X, 2)
    r = Inf
    for i = 1:n, j = 1:n
        if i != j && dists[i,j] < r
            r = dists[i,j]
        end
    end
    return r
end

d = 3
n = 100
k = 5
X = rand(d, n)
C = pairwise(SqEuclidean(), X, dims=2)

Xt = copy(transpose(X))
Ct = copy(transpose(C))

md0 = min_interdist(X)

@testset "Argument checks" begin
    @test_throws ArgumentError initseeds([1, 2], X, 3)
    @test initseeds([1, 2], X, 2) == [1, 2]
    @test_throws ArgumentError initseeds([-1, 2, 3], X, 3)
    @test_throws ArgumentError initseeds([1, n+2, 3], X, 3)

    @test_throws ArgumentError initseeds_by_costs([1, 2], C, 3)
    @test initseeds_by_costs([1, 2], C, 2) == [1, 2]

    @test_throws ArgumentError initseeds(:myseeding, X, 2)
    iseeds = initseeds(:kmpp, X, k)
    @test_throws DimensionMismatch copyseeds!(Matrix{Float64}(undef, 3, 6), X, iseeds)
    @test_throws DimensionMismatch copyseeds!(Matrix{Float64}(undef, 4, 5), X, iseeds)
    @test copyseeds!(Matrix{Float64}(undef, 3, 5), X, iseeds) isa Matrix{Float64}

    @testset "Seeds number check for $(typeof(alg))" for alg in
            (RandSeedAlg(), KmppAlg(), KmCentralityAlg())
        @test_throws ArgumentError initseeds(alg, X, 0)
        @test_throws ArgumentError initseeds(alg, X, n + 1)
        @test_throws ArgumentError initseeds_by_costs(alg, C, 0)
        @test_throws ArgumentError initseeds_by_costs(alg, C, n + 1)
        @test initseeds(alg, X, 4) isa Vector{Int}
        @test initseeds_by_costs(alg, C, 4) isa Vector{Int}
    end
end

@testset "RandSeed" begin
    Random.seed!(34568)
    iseeds = initseeds(RandSeedAlg(), X, k)
    @test length(iseeds) == k
    @test alldistinct(iseeds)
    Random.seed!(34568)
    iseeds_t = initseeds(RandSeedAlg(), Xt', k)
    @test iseeds == iseeds_t

    Random.seed!(34568)
    iseeds2 = initseeds(:rand, X, k)
    @test iseeds2 == iseeds

    Random.seed!(34568)
    iseeds = initseeds_by_costs(RandSeedAlg(), C, k)
    @test length(iseeds) == k
    @test alldistinct(iseeds)
    Random.seed!(34568)
    iseeds_t = initseeds_by_costs(RandSeedAlg(), Ct', k)
    @test iseeds == iseeds_t

    R = copyseeds!(Matrix{Float64}(undef, d, k), X, iseeds)
    @test isa(R, Matrix{Float64})
    @test R == X[:, iseeds]
    R_t = copyseeds!(Matrix{Float64}(undef, d, k), Xt', iseeds)
    @test R == R_t
end

@testset "Kmpp" begin
    Random.seed!(34568)
    iseeds = initseeds(KmppAlg(), X, k)
    @test length(iseeds) == k
    @test alldistinct(iseeds)
    Random.seed!(34568)
    iseeds_t = initseeds(KmppAlg(), Xt', k)
    @test iseeds == iseeds_t

    Random.seed!(34568)
    iseeds2 = initseeds(:kmpp, X, k)
    @test iseeds2 == iseeds
    Random.seed!(34568)
    iseeds_t2 = initseeds(:kmpp, Xt', k)
    @test iseeds_t2 == iseeds_t

    Random.seed!(34568)
    iseeds = initseeds_by_costs(KmppAlg(), C, k)
    @test length(iseeds) == k
    @test alldistinct(iseeds)
    Random.seed!(34568)
    iseeds_t = initseeds_by_costs(KmppAlg(), Ct', k)
    @test iseeds == iseeds_t

    @test min_interdist(X[:, iseeds]) > 20 * md0
    @test min_interdist((Xt')[:, iseeds]) > 20 * md0

    Random.seed!(34568)
    iseeds = initseeds_by_costs(:kmpp, C, k)
    @test length(iseeds) == k
    @test alldistinct(iseeds)
    Random.seed!(34568)
    iseeds_t = initseeds_by_costs(:kmpp, Ct', k)
    @test iseeds_t == iseeds
end

@testset "Kmcentrality" begin
    Random.seed!(34568)
    iseeds = initseeds(KmCentralityAlg(), X, k)
    @test length(iseeds) == k
    @test alldistinct(iseeds)
    Random.seed!(34568)
    iseeds_t = initseeds(KmCentralityAlg(), Xt', k)
    @test iseeds == iseeds_t

    Random.seed!(34568)
    iseeds2 = initseeds(:kmcen, X, k)
    @test iseeds2 == iseeds
    Random.seed!(34568)
    iseeds_t2 = initseeds(:kmcen, Xt', k)
    @test iseeds_t2 == iseeds_t

    Random.seed!(34568)
    iseeds = initseeds_by_costs(KmCentralityAlg(), C, k)
    @test length(iseeds) == k
    @test alldistinct(iseeds)
    Random.seed!(34568)
    iseeds_t = initseeds_by_costs(KmCentralityAlg(), Ct', k)
    @test iseeds == iseeds_t

    @test min_interdist(X[:, iseeds]) > 2 * md0
    @test min_interdist((Xt')[:, iseeds]) > 2 * md0
end

end

ext: jl


text:
function get_body_frame(obj::BodyShape)
    return jcall(obj, "getBodyFrame", Frame, ())
end

function get_ellipsoid(obj::EarthShape)
    return jcall(obj, "getEllipsoid", ReferenceEllipsoid, ())
end

function get_intersection_point(obj::BodyShape, arg0::FieldLine, arg1::FieldVector3D, arg2::Frame, arg3::FieldAbsoluteDate)
    return jcall(obj, "getIntersectionPoint", FieldGeodeticPoint, (FieldLine, FieldVector3D, Frame, FieldAbsoluteDate), arg0, arg1, arg2, arg3)
end

function get_intersection_point(obj::BodyShape, arg0::Line, arg1::Vector3D, arg2::Frame, arg3::AbsoluteDate)
    return jcall(obj, "getIntersectionPoint", GeodeticPoint, (Line, Vector3D, Frame, AbsoluteDate), arg0, arg1, arg2, arg3)
end

function project_to_ground(obj::BodyShape, arg0::TimeStampedPVCoordinates, arg1::Frame)
    return jcall(obj, "projectToGround", TimeStampedPVCoordinates, (TimeStampedPVCoordinates, Frame), arg0, arg1)
end

function project_to_ground(obj::BodyShape, arg0::Vector3D, arg1::AbsoluteDate, arg2::Frame)
    return jcall(obj, "projectToGround", Vector3D, (Vector3D, AbsoluteDate, Frame), arg0, arg1, arg2)
end

function transform(obj::BodyShape, arg0::FieldGeodeticPoint)
    return jcall(obj, "transform", FieldVector3D, (FieldGeodeticPoint,), arg0)
end

function transform(obj::BodyShape, arg0::FieldVector3D, arg1::Frame, arg2::FieldAbsoluteDate)
    return jcall(obj, "transform", FieldGeodeticPoint, (FieldVector3D, Frame, FieldAbsoluteDate), arg0, arg1, arg2)
end

function transform(obj::BodyShape, arg0::GeodeticPoint)
    return jcall(obj, "transform", Vector3D, (GeodeticPoint,), arg0)
end

function transform(obj::BodyShape, arg0::Vector3D, arg1::Frame, arg2::AbsoluteDate)
    return jcall(obj, "transform", GeodeticPoint, (Vector3D, Frame, AbsoluteDate), arg0, arg1, arg2)
end


text:
const plt = PyPlot

"Plots the contour of an airfoil given in x,y"
function plot_airfoil(x::Array{Float64,1}, y::Array{Float64,1};
                      label="Airfoil", style="-k", figfactor=1.0,
                      title_str="Airfoil geometry", legend_loc="best",
                      side_legend=(1.15, 1.1), zoom_factor=1.0, alpha=1.0)
  # Sizes the figure
  figsize = [7*1.5,5*0.5]*figfactor
  xmin, xmax = [-0.05, 1.05]*zoom_factor
  yrange = (xmax-xmin)/figsize[1] * figsize[2]
  ymin, ymax = -yrange/2, yrange/2
  fig1 = plt.figure(title_str, figsize=(figsize[1], figsize[2]))
  plt.xlim([xmin, xmax])
  plt.ylim([ymin, ymax])

  if label!="Airfoil"
    plt.plot(x,y, style, label=label, alpha=alpha)
    if side_legend!=nothing
      plt.legend(bbox_to_anchor=side_legend)
    else
      plt.legend(loc="best")
    end
  else
    plt.plot(x,y, style)
  end
  plt.xlabel("x")
  plt.ylabel("y")
  plt.grid(true, color="0.8", linestyle="--")
  plt.title(title_str)
end

"Plots loft splines"
function plot_loft_splines(b_pos, LE_x, LE_z, new_LE_x, new_LE_z, twists, new_twists,
  tilt_z, new_tilt_z, chords, new_chords)

     fig = plt.figure("spl_verif", figsize=(7*2,5*1))

     plt.subplot(121)
     plt.plot(LE_x[:,1], LE_x[:,2], "og", label="Org LE x", alpha=0.5)
     plt.plot(LE_z[:,1], LE_z[:,2], "ob", label="Org LE z", alpha=0.5)
     plt.plot(b_pos, new_LE_x, "--g", label="Spline LE x")
     plt.plot(b_pos, new_LE_z, "--b", label="Spline LE z")
     plt.xlabel(plt.L"y/b_{scale}")
     plt.ylabel(plt.L"x/b_{scale}, z/b_{scale}")
     plt.grid(true, color="0.8", linestyle="--")
     plt.legend(loc="best")

     plt.subplot(122)
     p1 = plt.plot(twists[:,1], twists[:,2], "og", label="Org Twist", alpha=0.5)
     p2 = plt.plot(b_pos, new_twists, "--g", label="Spline twist")
     pextra = []
     if tilt_z!=nothing
       pextra1 = plt.plot(tilt_z[:,1], tilt_z[:,2], "or", label="Org tilt z",
                                                                       alpha=0.5)
       pextra2 = plt.plot(b_pos, new_tlt_z, "--r", label="Spline tilt z")
       pextra = vcat(pextra, [pextra1[1], pextra2[1]])
     end
     plt.ylabel("Twist (deg)")

     plt.grid(true, color="0.8", linestyle="--")
     plt.xlabel(plt.L"y/b_{scale}")

     plt.twinx()
     p3 = plt.plot(chords[:,1], chords[:,2], "ob", label="Org Chord", alpha=0.5)
     p4 = plt.plot(b_pos, new_chords, "--b", label="Spline chord")
     plt.ylabel(plt.L"c/b_{scale}")

     ps = vcat([p1[1], p2[1], p3[1], p4[1]], pextra)
     plt.legend(ps, [p.get_label() for p in ps], loc="best")
     return nothing
end


text:
using MPSKit,TensorKit,Test

function bondmanage(state,ham,pars)
    final_bonddim = 12;

    upperbound = max_Ds(state);
    shouldincrease = reduce((a,i) -> a && dim(virtualspace(state,i))>=final_bonddim || upperbound[i+1] == dim(virtualspace(state,i)),1:length(state),init=true);

    if (shouldincrease)
        (state,pars) = changebonds(state, ham, OptimalExpand(),pars);
    end

    return (state,pars)
end

let
    #the operator used to evolve is the anticommutator
    th = nonsym_ising_ham()

    ham = anticommutator(th)

    ts = FiniteMPS(repeat(infinite_temperature(th),10))
    ca = params(ts,ham);

    betastep=0.1;endbeta=2;betas=collect(0:betastep:endbeta);

    for beta in betas
        (ts,ca) = bondmanage(ts,ham,ca)
        (ts,ca) = timestep(ts,ham,-betastep*0.25im,Tdvp(),ca) # find exp(-beta*H)
    end

end


text:
using CompressedSensing
using PyPlot
n=1000     #Signal Dimensions
m=400      #Measurement Dimension
k=100      #Approximate Sparsity
stdNoise=.3#standard deviation of noise

srand(0)
Signal = zeros(n)        #Create a random sparse signal
Signal[int(rand(k)*(n-1)+1)] = randn(k)

MM = randn(m,n)           #Create a random measurement matrix
Measurement = MM*Signal   #Measure the signal using the measurement matrix

Noise = randn(m)*stdNoise #Create measurement noise
Measurement += Noise      #add noise to measurement

#Reconstruct using IRLS
IRLS_Reconstruction = IRLS(MM,Measurement)

#reconstruct using UIRLS at multiple lambdas
DistToSig=[]
for lam = 0.0001:.25:5
    print("\r"*string("Calculating lambda=",lam))
    UIRLS_Reconstruction = UIRLS(MM,Measurement,lambda=lam)
    #Save the Euclidean Distance to original signal
    DistToSig =[DistToSig,sum(abs(UIRLS_Reconstruction-Signal).^2).^.5]
end

#plot the results
plot(0.0001:.25:5,DistToSig,label="UIRLS Distance to Signal")
DistIRLS=sum(abs(IRLS_Reconstruction-Signal).^2).^.5
plot([0,5],[DistIRLS,DistIRLS],label="IRLS Distance to Signal")
legend(loc=5)
xlabel("Lambda")
ylabel("Euclidean Distance to Original Signal")
savefig("Example_2_Fig_1.png")

max_stars_repo_path: examples/Example 2.jl

max_stars_repo_name: dahlend/CompressedSensing

max_issues_repo_path: examples/Example 2.jl

max_forks_repo_path: examples/Example 2.jl


text:
@testset "setequalconstant.jl" begin
    @testset "SetEqualConstant()" begin
        trailer = SeaPearl.Trailer()
        a = SeaPearl.IntSetVar(2, 5, "a", trailer)
        s = Set{Int}([5, 6, 7])

        constraint = SeaPearl.SetEqualConstant(a, s, trailer)

        @test constraint in a.onDomainChange
        @test constraint.active.value
    end

    @testset "propagate!(::SetEqualConstant)" begin
        trailer = SeaPearl.Trailer()
        a = SeaPearl.IntSetVar(2, 5, "a", trailer)
        s = Set{Int}([4, 5])
        constraint = SeaPearl.SetEqualConstant(a, s, trailer)

        toPropagate = Set{SeaPearl.Constraint}()
        prunedDomains = SeaPearl.CPModification()

        @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        @test SeaPearl.required_values(a.domain) == Set{Int}([4, 5])
        @test SeaPearl.possible_not_required_values(a.domain) == Set{Int}()
        @test !constraint.active.value

        # Infeasibility
        trailer = SeaPearl.Trailer()
        a = SeaPearl.IntSetVar(2, 5, "a", trailer)
        s = Set{Int}([4, 7])
        constraint = SeaPearl.SetEqualConstant(a, s, trailer)

        toPropagate = Set{SeaPearl.Constraint}()
        prunedDomains = SeaPearl.CPModification()

        @test !SeaPearl.propagate!(constraint, toPropagate, prunedDomains)
    end
end


text:
"""
Returns a set of all sub-expressions in an expression that look like \$some_expression
"""
function find_observable_expressions(obj::Expr)
    observable_expressions = Set()
    if is_interpolated_observable(obj)
        push!(observable_expressions, obj)
    else
        for a in obj.args
            observable_expressions = union(observable_expressions, find_observable_expressions(a))
        end
    end
    observable_expressions
end

# empty dict if x is not an Expr
find_observable_expressions(x) = Set()

is_interpolated_observable(x) = false
function is_interpolated_observable(e::Expr)
    e.head == Symbol(:$) && length(e.args) == 1
end

"""
Replaces every subexpression that looks like a observable expression with a substitute symbol stored in `exprdict`.
"""
function replace_observable_expressions!(exp::Expr, exprdict)
    if is_interpolated_observable(exp)
        error("You can't @lift an expression that only consists of a single observable.")
    else
        for (i, arg) in enumerate(exp.args)
            if is_interpolated_observable(arg)
                exp.args[i] = exprdict[arg]
            else
                replace_observable_expressions!(arg, exprdict)
            end
        end
    end
    return exp
end

replace_observable_expressions!(x, exprdict) = nothing

"""
Replaces an expression with lift(argtuple -> expression, args...), where args
are all expressions inside the main one that begin with \$.

# Example:

x = Observable(rand(100))
y = Observable(rand(100))

## before
z = lift((x, y) -> x .+ y, x, y)

## after
z = @lift(\$x .+ \$y)

You can also use parentheses around an expression if that expression evaluates to an observable.

```julia
nt = (x = Observable(1), y = Observable(2))
@lift(\$(nt.x) + \$(nt.y))
```
"""
macro lift(exp)

    observable_expr_set = find_observable_expressions(exp)

    if length(observable_expr_set) == 0
        error("Did not find any interpolated observables. Use '\$(observable)' to interpolate it into the macro.")
    end

    # store expressions with their substitute symbols, gensym them manually to be
    # able to escape the expression later
    observable_expr_arg_dict = Dict(expr => gensym("arg$i") for (i, expr) in enumerate(observable_expr_set))

    replace_observable_expressions!(exp, observable_expr_arg_dict)

    # keep an array for ordering
    observable_expressions_array = collect(keys(observable_expr_arg_dict))
    observable_substitutes_array = [observable_expr_arg_dict[expr] for expr in observable_expressions_array]
    observable_expressions_without_dollar = [n.args[1] for n in observable_expressions_array]

    # the arguments to the lifted function
    argtuple = Expr(Symbol(:tuple), observable_substitutes_array...)

    # the lifted function itself
    function_expression = Expr(Symbol(:->), argtuple, exp)

    # the full expression
    lift_expression = Expr(
        Symbol(:call),
        Symbol(:lift),
        esc(function_expression),
        esc.(observable_expressions_without_dollar)...
    )

    lift_expression
end


text:
using Test
using SequentialRaytrace

@testset "Dispersion" begin
    @test refractive_index(Air, 1.070) ≈ 1.0 atol = 1e-10
    @test refractive_index(Air, 0.94) ≈ 1.0 atol = 1e-10
    @test refractive_index(Silica, 1.070) ≈ 1.4495591322295867 atol = 1e-18
    @test refractive_index(Silica, 0.94) ≈ 1.4511992021095206 atol = 1e-18
end


text:
using RoadRunner
using Test

ant_str = """
    const Xo, X1
    Xo -> S1; k1*Xo - k2*S1
    S1 -> S2; k3*S2
    S2 -> X1; k4*S2

    Xo = 1;   X1 = 0
    S1 = 0;   S2 = 0
    k1 = 0.1; k2 = 0.56
    k3 = 1.2; k4 = 0.9
"""

rr = RoadRunner.loada(ant_str)

@testset "parameters" begin
    @test RoadRunner.getNumberOfGlobalParameters(rr) == 4
    @test RoadRunner.getGlobalParameterIds(rr) == ["k1", "k2", "k3", "k4"]
end

@testset "species" begin
    @test RoadRunner.getNumberOfFloatingSpecies(rr) == 2
    @test RoadRunner.getFloatingSpeciesIds(rr) == ["S1", "S2"]
    @test RoadRunner.getNumberOfBoundarySpecies(rr) == 2
    @test RoadRunner.getBoundarySpeciesIds(rr) == ["Xo", "X1"]
end

@testset "steadystates" begin
    @test RoadRunner.steadyState(rr) < 1e-6
end


text:
#=
    `CoeffBuffer` contains BigInt buffers used to minimize memory allocations
    during coefficient scaling, rational reconstruction and CRT reconstruction
=#
mutable struct CoeffBuffer
    # wow, 40 bytes per one BigInt

    # buffers that should be used
    # only for coefficient scaling
    scalebuf1::BigInt
    scalebuf2::BigInt

    # buffers that should be used
    # only for coefficient reduction
    reducebuf1::BigInt
    reducebuf2::BigInt
    reducebuf3::BigInt

    # buffers that should be used
    # only for coefficient reconstrction
    reconstructbuf1::BigInt
    reconstructbuf2::BigInt
    reconstructbuf3::BigInt
    reconstructbuf4::BigInt
    reconstructbuf5::BigInt
    reconstructbuf6::BigInt
    reconstructbuf7::BigInt
    reconstructbuf8::BigInt
    reconstructbuf9::BigInt
    reconstructbuf10::BigInt

    function CoeffBuffer()
        new(BigInt(), BigInt(), BigInt(), BigInt(), BigInt(),
            BigInt(), BigInt(), BigInt(), BigInt(), BigInt(), BigInt(),
            BigInt(), BigInt(), BigInt(), BigInt())
    end
end

#------------------------------------------------------------------------------

#=
    The structure stores accumulated integer coefficients and
    accumulated rational coefficients of a groebner basis
=#
mutable struct CoeffAccum
    gb_coeffs_zz::Vector{Vector{CoeffZZ}}
    prev_gb_coeffs_zz::Vector{Vector{CoeffZZ}}
    gb_coeffs_qq::Vector{Vector{CoeffQQ}}

    function CoeffAccum()
        new(Vector{Vector{CoeffZZ}}(undef, 0),
            Vector{Vector{CoeffZZ}}(undef, 0),
            Vector{Vector{CoeffQQ}}(undef, 0))
    end
end

#------------------------------------------------------------------------------

# compute gcd of denominators of `coeffs`
function common_denominator(coeffbuff::CoeffBuffer, coeffs::Vector{CoeffQQ})
    den = coeffbuff.scalebuf1
    Base.GMP.MPZ.set_si!(den, 1)
    for c in coeffs
        # TODO
        # use buffer?
        Base.GMP.MPZ.lcm!(den, denominator(c))
    end
    den
end

# compute gcd of denominators of `coeffs`
function common_denominator(coeffs::Vector{CoeffQQ})
    den = BigInt()
    Base.GMP.MPZ.set_si!(den, 1)
    for c in coeffs
        # TODO
        # use buffer?
        Base.GMP.MPZ.lcm!(den, denominator(c))
    end
    den
end

# scales denominators inplace
function scale_denominators!(
            coeffbuff::CoeffBuffer,
            coeffs_qq::Vector{Vector{CoeffQQ}},
            coeffs_zz::Vector{Vector{CoeffZZ}})

    @assert length(coeffs_zz) == length(coeffs_qq)

    buf = coeffbuff.scalebuf2
    for i in 1:length(coeffs_qq)
        @assert length(coeffs_zz[i]) == length(coeffs_qq[i])

        den = common_denominator(coeffbuff, coeffs_qq[i])
        sz  = Base.GMP.MPZ.sizeinbase(den, 2)
        @inbounds for j in 1:length(coeffs_qq[i])
            num = numerator(coeffs_qq[i][j])
            Base.GMP.MPZ.tdiv_q!(buf, den, denominator(coeffs_qq[i][j]))
            Base.GMP.MPZ.realloc2!(coeffs_zz[i][j], sz)
            Base.GMP.MPZ.mul!(coeffs_zz[i][j], num, buf)
        end
    end

    coeffs_zz
end

function scale_denominators(
            coeffbuff::CoeffBuffer,
            coeffs_qq::Vector{Vector{CoeffQQ}})
    coeffs_zz = [[CoeffZZ(0) for _ in 1:length(c)] for c in coeffs_qq]
    scale_denominators!(coeffbuff, coeffs_qq, coeffs_zz)
end

function scale_denominators(coeffs_qq::Vector{Vector{CoeffQQ}})
    coeffs_zz = [[CoeffZZ(0) for _ in 1:length(c)] for c in coeffs_qq]
    buf = BigInt()
    for i in 1:length(coeffs_qq)
        @assert length(coeffs_zz[i]) == length(coeffs_qq[i])

        den = common_denominator(coeffs_qq[i])
        @inbounds for j in 1:length(coeffs_qq[i])
            num = numerator(coeffs_qq[i][j])
            Base.GMP.MPZ.tdiv_q!(buf, den, denominator(coeffs_qq[i][j]))
            Base.GMP.MPZ.mul!(coeffs_zz[i][j], num, buf)
        end
    end

    coeffs_zz
end

function scale_denominators(coeffs_qq::Vector{CoeffQQ})
    coeffs_zz = [CoeffZZ(0) for _ in 1:length(coeffs_qq)]
    buf = BigInt()
    den = common_denominator(coeffs_qq)
    for i in 1:length(coeffs_qq)
        num = numerator(coeffs_qq[i])
        Base.GMP.MPZ.tdiv_q!(buf, den, denominator(coeffs_qq[i]))
        Base.GMP.MPZ.mul!(coeffs_zz[i], num, buf)
    end
    coeffs_zz
end

#------------------------------------------------------------------------------

function reduce_modulo!(
        coeffbuff::CoeffBuffer,
        coeffs_zz::Vector{Vector{CoeffZZ}},
        coeffs_ff::Vector{Vector{CoeffFF}},
        prime::CoeffFF)

    p   = coeffbuff.reducebuf1
    buf = coeffbuff.reducebuf2
    c   = coeffbuff.reducebuf3
    Base.GMP.MPZ.set_ui!(p, prime)

    for i in 1:length(coeffs_zz)
        cfs_zz_i = coeffs_zz[i]
        @inbounds for j in 1:length(cfs_zz_i)
            Base.GMP.MPZ.set!(c, cfs_zz_i[j])
            # TODO: rewrite this
            if Base.GMP.MPZ.cmp_ui(c, 0) < 0
                Base.GMP.MPZ.fdiv_q!(buf, c, p)
                Base.GMP.MPZ.neg!(buf)
                Base.GMP.MPZ.mul_ui!(buf, prime)
                Base.GMP.MPZ.add!(c, buf)
            end
            # @assert c >= 0 # TODO
            Base.GMP.MPZ.tdiv_r!(buf, c, p)
            coeffs_ff[i][j] = UInt64(buf)
        end
    end
end

function reduce_modulo(
        coeffbuff::CoeffBuffer,
        coeffs_zz::Vector{Vector{CoeffZZ}},
        prime::UInt64)
    coeffs_ff =  [Vector{CoeffFF}(undef, length(c)) for c in coeffs_zz]
    reduce_modulo!(coeffbuff, coeffs_zz, coeffs_ff, prime)
    coeffs_ff
end

#------------------------------------------------------------------------------

function resize_accum!(coeffaccum::CoeffAccum, gb_coeffs)
    resize!(coeffaccum.gb_coeffs_zz, length(gb_coeffs))
    resize!(coeffaccum.prev_gb_coeffs_zz, length(gb_coeffs))
    resize!(coeffaccum.gb_coeffs_qq, length(gb_coeffs))
    @inbounds for i in 1:length(gb_coeffs)
        coeffaccum.gb_coeffs_zz[i] = [CoeffZZ(0) for _ in 1:length(gb_coeffs[i])]
        coeffaccum.prev_gb_coeffs_zz[i] = [CoeffZZ(0) for _ in 1:length(gb_coeffs[i])]
        coeffaccum.gb_coeffs_qq[i] = [CoeffQQ(1) for _ in 1:length(gb_coeffs[i])]
    end
end

function reconstruct_trivial_crt!(coeffbuff::CoeffBuffer,
                                    coeffaccum::CoeffAccum,
                                    gb_coeffs_ff::Vector{Vector{CoeffFF}})
    gb_coeffs_zz = coeffaccum.gb_coeffs_zz
    for i in 1:length(gb_coeffs_ff)
        @inbounds for j in 1:length(gb_coeffs_ff[i])
            Base.GMP.MPZ.set_ui!(gb_coeffs_zz[i][j], gb_coeffs_ff[i][j])
        end
    end
end

function assure_structure(coeffaccum, gb_coeffs_ff)
    if length(coeffaccum.gb_coeffs_zz) != gb_coeffs_ff
        return false
    end
    @inbounds for i in 1:length(gb_coeffs_ff)
        if length(coeffaccum.gb_coeffs_zz[i]) != gb_coeffs_ff[i]
            return false
        end
    end
    return true
end

function reconstruct_crt!(
        coeffbuff::CoeffBuffer,
        coeffaccum::CoeffAccum,
        primetracker::PrimeTracker,
        gb_coeffs_ff::Vector{Vector{CoeffFF}},
        ch::UInt64)

    if isempty(coeffaccum.gb_coeffs_zz)
        # if first time reconstruction
        resize_accum!(coeffaccum, gb_coeffs_ff)
        reconstruct_trivial_crt!(coeffbuff, coeffaccum, gb_coeffs_ff)
    else
        gb_coeffs_zz = coeffaccum.gb_coeffs_zz
        prev_gb_coeffs_zz = coeffaccum.prev_gb_coeffs_zz

        # copy to previous gb coeffs
        for i in 1:length(gb_coeffs_zz)
            for j in 1:length(gb_coeffs_zz[i])
                Base.GMP.MPZ.set!(prev_gb_coeffs_zz[i][j], gb_coeffs_zz[i][j])
            end
        end

        buf = coeffbuff.reconstructbuf1
        n1, n2 = coeffbuff.reconstructbuf2, coeffbuff.reconstructbuf3
        M = coeffbuff.reconstructbuf4
        bigch = coeffbuff.reconstructbuf5
        invm1, invm2 = coeffbuff.reconstructbuf6, coeffbuff.reconstructbuf7

        Base.GMP.MPZ.set_ui!(bigch, ch)
        Base.GMP.MPZ.mul_ui!(M, primetracker.modulo, ch)
        Base.GMP.MPZ.gcdext!(buf, invm1, invm2, primetracker.modulo, bigch)

        for i in 1:length(gb_coeffs_ff)
            @inbounds for j in 1:length(gb_coeffs_ff[i])
                ca = gb_coeffs_zz[i][j]
                cf = gb_coeffs_ff[i][j]
                CRT!(M, buf, n1, n2, ca, invm1, cf, invm2, primetracker.modulo, bigch)
                # TODO: faster set!??
                Base.GMP.MPZ.set!(gb_coeffs_zz[i][j], buf)
            end
        end
    end
    updatemodulo!(primetracker)
end

#------------------------------------------------------------------------------

function reconstruct_modulo!(
        coeffbuff::CoeffBuffer,
        coeffaccum::CoeffAccum,
        primetracker::PrimeTracker)

    modulo = primetracker.modulo

    bnd = rational_reconstruction_bound(modulo)
    buf, buf1  = coeffbuff.reconstructbuf1, coeffbuff.reconstructbuf2
    buf2, buf3 = coeffbuff.reconstructbuf3, coeffbuff.reconstructbuf4
    u1, u2     = coeffbuff.reconstructbuf5, coeffbuff.reconstructbuf6
    u3, v1     = coeffbuff.reconstructbuf7, coeffbuff.reconstructbuf8
    v2, v3     = coeffbuff.reconstructbuf9, coeffbuff.reconstructbuf10

    gb_coeffs_zz = coeffaccum.gb_coeffs_zz
    gb_coeffs_qq = coeffaccum.gb_coeffs_qq

    for i in 1:length(gb_coeffs_zz)
        @inbounds for j in 2:length(gb_coeffs_zz[i])
            cz = gb_coeffs_zz[i][j]
            cq = gb_coeffs_qq[i][j]
            num, den = numerator(cq), denominator(cq)
            success = rational_reconstruction!(num, den, bnd, buf,
                                        buf1, buf2, buf3,
                                        u1, u2, u3, v1, v2, v3,
                                        cz, modulo)
            if !success
                # @error "not success"
                return false
            end
            # TODO
            # @assert gcd(numerator(cq), denominator(cq)) == 1
        end
    end

    return true
end


text:
# # State debug statistics
#
# This page shows how to use the `StateCheck` functions to get basic
# statistics for nodal values of fields held in ClimateMachine `MPIStateArray`
# data structures. The `StateCheck` functions can be used to
#
# 1. Generate statistics on `MPIStateArray` holding the state of a ClimateMachine experiment.
#
#    and to
#
# 2. Compare against saved reference statistics from ClimateMachine `MPIStateArray`
#    variables. This can enable simple automated regression test checks for
#    detecting unexpected changes introduced into numerical experiments
#    by code updates.
#
# These two cases are shown below:

# ## 1. Generating statistics for a set of MPIStateArrays
#
# Here we create a callback that can generate statistics for an arbitrary
# set of the MPIStateArray type variables of the sort that hold persistent state for
# ClimateMachine models. We then invoke the call back to show the statistics.
#
# In regular use the `MPIStateArray` variables will come from model configurations.
# Here we create a dummy set of `MPIStateArray` variables for use in stand alone
# examples.

# ### Create a dummy set of MPIStateArrays
#
# First we set up two `MPIStateArray` variables. This need a few packages to be in placeT,
# and utilizes some utility functions to create the array and add named
# persistent state variables.
# This is usually handled automatically as part of model definition in regular
# ClimateMachine activity.
# Calling `ClimateMachine.init()` includes initializing GPU CUDA and MPI parallel
# processing options that match the hardware/software system in use.

# Set up a basic environment
using MPI
using StaticArrays
using Random
using ClimateMachine
using ClimateMachine.VariableTemplates
using ClimateMachine.MPIStateArrays
using ClimateMachine.StateCheck
ClimateMachine.init()
FT = Float64

# Define some dummy vector and tensor abstract variables with associated types
# and dimensions
F1 = @vars begin
    ν∇u::SMatrix{3, 2, FT, 6}
    κ∇θ::SVector{3, FT}
end
F2 = @vars begin
    u::SVector{2, FT}
    θ::SVector{1, FT}
end
nothing # hide

# Create `MPIStateArray` variables with arrays to hold elements of the
# vectors and tensors
Q1 = MPIStateArray{Float32, F1}(
    MPI.COMM_WORLD,
    ClimateMachine.array_type(),
    4,
    9,
    8,
)
Q2 = MPIStateArray{Float64, F2}(
    MPI.COMM_WORLD,
    ClimateMachine.array_type(),
    4,
    6,
    8,
)
nothing # hide

# ### Create a call-back
#
# Now we can create a `StateCheck` call-back, _cb_, tied to the `MPIStateArray`
# variables _Q1_ and _Q2_. Each `MPIStateArray` in the array
# of `MPIStateArray` variables tracked is paired with a label
# to identify it. The call-back is also given a frequency (in time step numbers) and
# precision for printing summary tables.
cb = ClimateMachine.StateCheck.sccreate(
    [(Q1, "My gradients"), (Q2, "My fields")],
    1;
    prec = 15,
)
nothing # hide

# ### Invoke the call-back
#
# The call-back is of type `ClimateMachine.GenericCallbacks.EveryXSimulationSteps`
# and in regular use is designed to be passed to a ClimateMachine timestepping
# solver e.g.
typeof(cb)

# Here, for demonstration purposes, we can invoke
# the call-back after simply initializing the `MPIStateArray` fields to a random
# set of values e.g.
Q1.data .= rand(MersenneTwister(0), Float32, size(Q1.data))
Q2.data .= rand(MersenneTwister(0), Float64, size(Q2.data))
cb()

# ## 2. Comparing to reference values

# ### Generate arrays of reference values
#
# StateCheck functions can generate text that can be used to set the value of stored
# arrays that can be used in a reference test for subsequent regression testing. This
# involves 3 steps.
#
# **Step 1.** First a reference array setting program code is generated from the latest
# state of a given callback e.g.

ClimateMachine.StateCheck.scprintref(cb)

# **Step 2.** Next the array setting program code is executed (see below). At this stage the _parr[]_ array
# context may be hand edited. The parr[] array sets a target number of decimal places for
# matching against reference values in _varr[]_. For different experiments and different fields
# the degree of precision that constitutes failing a regression test may vary. Choosing the
# _parr[]_ values requires some sense as to the stability of the particular numerical
# and physical scenario an experiment represents. In the example below some precision
# settings have been hand edited from the default of 16 to illustrate the process.

#! format: off
varr = [
 [ "My gradients", "ν∇u[1]",  1.34348869323730468750e-04,  9.84732866287231445313e-01,  5.23545503616333007813e-01,  3.08209930764271777814e-01 ],
 [ "My gradients", "ν∇u[2]",  1.16317868232727050781e-01,  9.92088317871093750000e-01,  4.83800649642944335938e-01,  2.83350456014221541157e-01 ],
 [ "My gradients", "ν∇u[3]",  1.05845928192138671875e-03,  9.51775908470153808594e-01,  4.65474426746368408203e-01,  2.73615551085745090099e-01 ],
 [ "My gradients", "ν∇u[4]",  5.97668886184692382813e-02,  9.68048095703125000000e-01,  5.42618036270141601563e-01,  2.81570862027933854765e-01 ],
 [ "My gradients", "ν∇u[5]",  8.31030607223510742188e-02,  9.35931921005249023438e-01,  5.05405902862548828125e-01,  2.46073509972619536290e-01 ],
 [ "My gradients", "ν∇u[6]",  3.09681892395019531250e-02,  9.98341441154479980469e-01,  4.54375565052032470703e-01,  3.09461067853178561915e-01 ],
 [ "My gradients", "κ∇θ[1]",  8.47448110580444335938e-02,  9.94180679321289062500e-01,  5.27157366275787353516e-01,  2.92455951648181833313e-01 ],
 [ "My gradients", "κ∇θ[2]",  1.20514631271362304688e-02,  9.93527650833129882813e-01,  4.71063584089279174805e-01,  2.96449027197666359346e-01 ],
 [ "My gradients", "κ∇θ[3]",  8.14980268478393554688e-02,  9.55443382263183593750e-01,  5.05038917064666748047e-01,  2.77201022741208891187e-01 ],
 [    "My fields",   "u[1]",  4.31410233294131639781e-02,  9.97140933049696531754e-01,  4.62139750850942054861e-01,  3.23076684924287371725e-01 ],
 [    "My fields",   "u[2]",  1.01416659908237782872e-02,  9.14712023896926407218e-01,  4.76160523012988778913e-01,  2.71443440757963339038e-01 ],
 [    "My fields",   "θ[1]",  6.58965491052394547467e-02,  9.73216404386510802738e-01,  4.60007166313864512830e-01,  2.87310472114545079059e-01 ],
]
parr = [
 [ "My gradients", "ν∇u[1]",    16,     7,    16,     0 ],
 [ "My gradients", "ν∇u[2]",    16,     7,    16,     0 ],
 [ "My gradients", "ν∇u[3]",    16,     7,    16,     0 ],
 [ "My gradients", "ν∇u[4]",    16,     7,    16,     0 ],
 [ "My gradients", "ν∇u[5]",    16,     7,    16,     0 ],
 [ "My gradients", "ν∇u[6]",    16,     7,    16,     0 ],
 [ "My gradients", "κ∇θ[1]",    16,    16,    16,     0 ],
 [ "My gradients", "κ∇θ[2]",    16,    16,    16,     0 ],
 [ "My gradients", "κ∇θ[3]",    16,    16,    16,     0 ],
 [    "My fields",   "u[1]",    16,    16,    16,     0 ],
 [    "My fields",   "u[2]",    16,    16,    16,     0 ],
 [    "My fields",   "θ[1]",    16,    16,    16,     0 ],
]
#! format: on

# **Step 3.** Finally a call-back stored value can be compared for consistency to with _parr[]_ decimal places

ClimateMachine.StateCheck.scdocheck(cb, (varr, parr))
nothing # hide

# In this trivial case the match is guaranteed. The function will return _true_ to the calling
# routine and this can be passed to an `@test` block.
#
# However we can modify the reference test values to
# see the effect of a mismatch e.g.
varr[1][3] = varr[1][3] * 10.0
ClimateMachine.StateCheck.scdocheck(cb, (varr, parr))
nothing # hide

# Here the mis-matching field is highlighted with _N(0)_ indicating that the precision
# was not met and actual match length was (in this case) 0. If any field fails the test returns false
# for use in any regression testing control logic.

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
    LDLt <: Factorization

Matrix factorization type of the `LDLt` factorization of a real [`SymTridiagonal`](@ref)
matrix `S` such that `S = L*Diagonal(d)*L'`, where `L` is a [`UnitLowerTriangular`](@ref)
matrix and `d` is a vector. The main use of an `LDLt` factorization `F = ldlt(S)`
is to solve the linear system of equations `Sx = b` with `F\\b`. This is the
return type of [`ldlt`](@ref), the corresponding matrix factorization function.

The individual components of the factorization `F::LDLt` can be accessed via `getproperty`:

| Component | Description                                 |
|:---------:|:--------------------------------------------|
| `F.L`     | `L` (unit lower triangular) part of `LDLt`  |
| `F.D`     | `D` (diagonal) part of `LDLt`               |
| `F.Lt`    | `Lt` (unit upper triangular) part of `LDLt` |
| `F.d`     | diagonal values of `D` as a `Vector`        |

# Examples
```jldoctest
julia> S = SymTridiagonal([3., 4., 5.], [1., 2.])
3×3 SymTridiagonal{Float64, Vector{Float64}}:
 3.0  1.0   ⋅
 1.0  4.0  2.0
  ⋅   2.0  5.0

julia> F = ldlt(S)
LDLt{Float64, SymTridiagonal{Float64, Vector{Float64}}}
L factor:
3×3 UnitLowerTriangular{Float64, SymTridiagonal{Float64, Vector{Float64}}}:
 1.0        ⋅         ⋅
 0.333333  1.0        ⋅
 0.0       0.545455  1.0
D factor:
3×3 Diagonal{Float64, Vector{Float64}}:
 3.0   ⋅        ⋅
  ⋅   3.66667   ⋅
  ⋅    ⋅       3.90909
```
"""
struct LDLt{T,S<:AbstractMatrix{T}} <: Factorization{T}
    data::S

    function LDLt{T,S}(data) where {T,S<:AbstractMatrix{T}}
        require_one_based_indexing(data)
        new{T,S}(data)
    end
end
LDLt(data::AbstractMatrix{T}) where {T} = LDLt{T,typeof(data)}(data)
LDLt{T}(data::AbstractMatrix) where {T} = LDLt(convert(AbstractMatrix{T}, data)::AbstractMatrix{T})

size(S::LDLt) = size(S.data)
size(S::LDLt, i::Integer) = size(S.data, i)

LDLt{T,S}(F::LDLt{T,S}) where {T,S<:AbstractMatrix{T}} = F
LDLt{T,S}(F::LDLt) where {T,S<:AbstractMatrix{T}} = LDLt{T,S}(convert(S, F.data)::S)
LDLt{T}(F::LDLt{T}) where {T} = F
LDLt{T}(F::LDLt) where {T} = LDLt(convert(AbstractMatrix{T}, F.data)::AbstractMatrix{T})

Factorization{T}(F::LDLt{T}) where {T} = F
Factorization{T}(F::LDLt) where {T} = LDLt{T}(F)

function getproperty(F::LDLt, d::Symbol)
    Fdata = getfield(F, :data)
    if d === :d
        return Fdata.dv
    elseif d === :D
        return Diagonal(Fdata.dv)
    elseif d === :L
        return UnitLowerTriangular(Fdata)
    elseif d === :Lt
        return UnitUpperTriangular(Fdata)
    else
        return getfield(F, d)
    end
end

adjoint(F::LDLt{<:Real,<:SymTridiagonal}) = F
adjoint(F::LDLt) = LDLt(copy(adjoint(F.data)))

function show(io::IO, mime::MIME{Symbol("text/plain")}, F::LDLt)
    summary(io, F); println(io)
    println(io, "L factor:")
    show(io, mime, F.L)
    println(io, "\nD factor:")
    show(io, mime, F.D)
end

# SymTridiagonal
"""
    ldlt!(S::SymTridiagonal) -> LDLt

Same as [`ldlt`](@ref), but saves space by overwriting the input `S`, instead of creating a copy.

# Examples
```jldoctest
julia> S = SymTridiagonal([3., 4., 5.], [1., 2.])
3×3 SymTridiagonal{Float64, Vector{Float64}}:
 3.0  1.0   ⋅
 1.0  4.0  2.0
  ⋅   2.0  5.0

julia> ldltS = ldlt!(S);

julia> ldltS === S
false

julia> S
3×3 SymTridiagonal{Float64, Vector{Float64}}:
 3.0       0.333333   ⋅
 0.333333  3.66667   0.545455
  ⋅        0.545455  3.90909
```
"""
function ldlt!(S::SymTridiagonal{T,V}) where {T,V}
    n = size(S,1)
    d = S.dv
    e = S.ev
    @inbounds for i in 1:n-1
        iszero(d[i]) && throw(ZeroPivotException(i))
        e[i] /= d[i]
        d[i+1] -= e[i]^2*d[i]
    end
    return LDLt{T,SymTridiagonal{T,V}}(S)
end

"""
    ldlt(S::SymTridiagonal) -> LDLt

Compute an `LDLt` factorization of the real symmetric tridiagonal matrix `S` such that `S = L*Diagonal(d)*L'`
where `L` is a unit lower triangular matrix and `d` is a vector. The main use of an `LDLt`
factorization `F = ldlt(S)` is to solve the linear system of equations `Sx = b` with `F\\b`.

# Examples
```jldoctest
julia> S = SymTridiagonal([3., 4., 5.], [1., 2.])
3×3 SymTridiagonal{Float64, Vector{Float64}}:
 3.0  1.0   ⋅
 1.0  4.0  2.0
  ⋅   2.0  5.0

julia> ldltS = ldlt(S);

julia> b = [6., 7., 8.];

julia> ldltS \\ b
3-element Vector{Float64}:
 1.7906976744186047
 0.627906976744186
 1.3488372093023255

julia> S \\ b
3-element Vector{Float64}:
 1.7906976744186047
 0.627906976744186
 1.3488372093023255
```
"""
function ldlt(M::SymTridiagonal{T}; shift::Number=false) where T
    S = typeof((zero(T)+shift)/one(T))
    Mₛ = SymTridiagonal{S}(copy_oftype(M.dv, S), copy_oftype(M.ev, S))
    if !iszero(shift)
        Mₛ.dv .+= shift
    end
    return ldlt!(Mₛ)
end

factorize(S::SymTridiagonal) = ldlt(S)

function ldiv!(S::LDLt{<:Any,<:SymTridiagonal}, B::AbstractVecOrMat)
    require_one_based_indexing(B)
    n, nrhs = size(B, 1), size(B, 2)
    if size(S,1) != n
        throw(DimensionMismatch("Matrix has dimensions $(size(S)) but right hand side has first dimension $n"))
    end
    d = S.data.dv
    l = S.data.ev
    @inbounds begin
        for i = 2:n
            li1 = l[i-1]
            @simd for j = 1:nrhs
                B[i,j] -= li1*B[i-1,j]
            end
        end
        dn = d[n]
        @simd for j = 1:nrhs
            B[n,j] /= dn
        end
        for i = n-1:-1:1
            di = d[i]
            li = l[i]
            @simd for j = 1:nrhs
                B[i,j] /= di
                B[i,j] -= li*B[i+1,j]
            end
        end
    end
    return B
end

rdiv!(B::AbstractVecOrMat, S::LDLt{<:Any,<:SymTridiagonal}) =
    transpose(ldiv!(S, transpose(B)))

function logabsdet(F::LDLt{<:Any,<:SymTridiagonal})
    it = (F.data[i,i] for i in 1:size(F, 1))
    return sum(log∘abs, it), prod(sign, it)
end

# Conversion methods
function SymTridiagonal(F::LDLt)
    e = copy(F.data.ev)
    d = copy(F.data.dv)
    e .*= d[1:end-1]
    d[2:end] += e .* F.data.ev
    SymTridiagonal(d, e)
end
AbstractMatrix(F::LDLt) = SymTridiagonal(F)
AbstractArray(F::LDLt) = AbstractMatrix(F)
Matrix(F::LDLt) = Array(AbstractArray(F))
Array(F::LDLt) = Matrix(F)

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia