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

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
	using ExtendableGrids
	using PlutoVista
	using GridVisualize
	using VoronoiFVM
	default_plotter!(PlutoVista)
end

# ╔═╡ 556480e0-94f1-4e47-be9a-3e1e0e99555c
TableOfContents(title="")

# ╔═╡ 5e13b3db-570c-4159-939a-7e2268f0a102
md"""
# Some problems with Voronoi FVM

Draft. J. Fuhrmann, Oct. 29. 2021

"""

# ╔═╡ fae47c55-eef8-4428-bb5f-45824978753d
md"""
## Transient problem

This problem was suggested by R. Eymard.
"""

# ╔═╡ 8ba2300c-17ff-44e1-b33a-c5bdf1ce12fe
md"""
Regard the following problem coupling Darcy's equation with Fick's law and transport:
"""

# ╔═╡ 51c9517c-8797-4406-b053-301694fb0484
md"""
```math
  \begin{aligned}
    \vec v &= k \nabla p \\
    \nabla \cdot \vec v &= 0\\
    \partial_t (\phi c) - \nabla \cdot (D\nabla c + c \vec v) &= 0
  \end{aligned}
```
"""

# ╔═╡ 99341e32-9c78-4e31-bec0-d1ffbc85ec32
md"""
The domain is described by the following discretization grid:
"""

# ╔═╡ cd013964-f329-4d2c-ae4b-305093f0ac56
md"""
### Results

In the calculations, we ramp up the inlet concentration and  measure  the amount  of solute flowing  through the outlet - the breaktrough curve.
"""

# ╔═╡ afccbadb-2ca8-4c3e-8c6d-c78df59d8d7e
nref=1

# ╔═╡ dd9f8d38-8812-40ba-88c8-f873ec7d6121
tend=100

# ╔═╡ 5f6ac608-b1a0-450e-910e-d7d8ea2ffae0
ε_fix=1.0e-4

# ╔═╡ 5b60c7d4-7bdb-4989-b055-6695b9fdeedc
md"""
Here, we plot the solutions for the `grid_n` case and the `grid_f` case.
"""

# ╔═╡ f6abea66-1e42-4201-8433-5d092989749d
begin
	vis_n=GridVisualizer(dim=2,resolution=(210,150))
	vis_f=GridVisualizer(dim=2,resolution=(210,150))
	(vis_n,vis_f)
end

# ╔═╡ 98ae56dd-d42d-4a93-bb0b-5956b6e981a3
md"""
Time: $(@bind t Slider(1:tend/100:tend,show_value=true,default=tend*0.1))
"""

# ╔═╡ 99c3b54b-d458-482e-8aa0-d2c2b51fdf25
md"""
## Reaction-Diffusion problem

Here we solve the following problem:

```math
    -\nabla D \nabla u + R u = 0
```
where D is large in the high permeability region and small otherwise. R is a constant.

"""

# ╔═╡ eef85cd7-eba4-4c10-9e1d-38411179314d
md"""
### Results
"""

# ╔═╡ fcd066f1-bcd8-4479-a4e4-7b8c235336c4
md"""
## Discussion

### Transient case
As there will be nearly no flow
in  y-direction, we should  get the  very same  results in  all four
cases for small permeability values in the low permeability region.

In the `grid_n` case,  the heterogeneous control volumina  ovrestimate the storage
capacity which shows itself  in a underestimation  of the
transferred solute.

With  the high  permeability contrast,  the results  for heterogeneous
domain should be essentially equal to those for 1D domain.
 However,   with  a   coarse  resolution   in
y-direction, we see large  differences in the transient behaviour of
the breaktrough curve compared to the 1D case.
The introduction of a thin  protection layer leads  to  reasonable   results.

Alternatively, the porosity of the low permeability region can be modified.
Arguably, this is the case in practice, see e.g.
[Ackerer et al, Transport in Porous Media35:345–373, 1999](https://link.springer.com/content/pdf/10.1023/A:1006564309167.pdf)
(eq. 7).

### Reaction diffusion case
In this case, we look at a homogeneous reaction in a domain divided
into a high and low diffusion region. With high contrast in the diffusion
coefficients, the reasonable assumption is that the reaction takes place only
in the high diffusion region, and the un-consumed share of species leaves at the outlet.

In this case we observe a similar related problem which can be fixed by adding a thin layer
of control volumes at the boundary. No problem occurs if the reaction rate at in the low diffusion
region is zero.

### Conclusion

Here, we indeed observe problem with the Voronoi approach: care must be taken to handle the case
of hetero interfaces in connection with transient processes and/or homogeneous reactions.
In these cases it should be analyzed if the problem occurs, and why, and it appears, that the discussion
should not be had without reference to the correct physical models. A remedy based on meshing
is available at least for straight interfaces.

### Opinion

With standard ways of using finite elements, the issue described here will occur in a similar way, so
the discussion is indeed with the alternative "cell centered" finite volume approach which places interfaces
at the boundaries of the control volumes rather than along the edges of a underlying triangulation.

#### Drawbacks of two point flux Voronoi methods based on simplicial meshes (as tested here):
- Anisotropic diffusion is only correct with aligned meshes
- Reliance on boundary conforming Delaunay property of the underyling mesh, thus narrowing the available meshing strategies
- The issue described  in the present notebook. However, in both cases discussed here, IMHO it might  "go  away"  depending on the correct physics.
  There should be more discussions with relevant application problems at hand.

#### Advantages (compared to the cell centered approach placing collocation points away from interfaces)
- Availability of P1 interpolant on simplices for visualization, interpolation, coupling etc.
- Mesh generators tend to place interfaces at triangle edges.
- Dirichlet BC can be applied exactly
- There is a straighforward way to link interface processes with bulk processes, e.g. an adsorption reaction is easily described by a reaction term at the boundary which involves interface and bulk value available at the same mesh node.

"""

# ╔═╡ c9d92201-813c-499b-b863-b138c30eb634
md"""
## Appendix
"""

# ╔═╡ a372ac90-c871-4dc0-a44b-a5bddef71823
md"""
### Domain data
"""

# ╔═╡ 124b2a0a-ef19-453e-9e3a-5b5ce7db5fac
md"""
Sizes:
"""

# ╔═╡ 1ad18670-e7cb-4f7a-be0f-3db98cdeb6a4
begin
L=10   # length of the high perm layer
W=0.5  # width of high perm layer
Wlow=2 # width of adjacent low perm layers
end;

# ╔═╡ cc325b2c-6174-4b8d-8e39-202ac68b5705
md"""
In the center of the domain, we assume a layer with high permeability.

As  boundary  conditions for the pressure ``p`` we choose  fixed pressure values at  the left
and right boundaries of the  domain, triggering a constant pressure gradient throughout the domain.

At the inlet of the high  permeability layer, we set ``c=1``, and at the
outlet, we set ``c=0``.

The high permeability layer has length `L`=$( L) and width `W`= $(W).

We solve the time dependent problem on three types of  rectangular grids with the same
resolution in   $x$ direction and different variants to to handle the  high permeability
layer.

- `grid_n` - a "naive" grid which just resolves the permeability layer and the surrounding material with equally spaced (in y direction) grids
- `grid_1` - a 1D grid  of the high permeability layer. With high permeability contrast, the solution of the 2D case at y=0 should conincide with the 1D solution
- `grid_f` - a "fixed" 2D grid which resolves the permeability layer and the surrounding material with equally spaced (in y direction) grids and "protection layers" of width `ε_fix`=$(ε_fix)  correcting the size of high permeability control volumes

"""

# ╔═╡ 47bc8e6a-e296-42c9-bfc5-967edfb0feb7
md"""
Boundary conditions:
"""

# ╔═╡ d1d5bad2-d282-4e7d-adb9-baf21f58155e
begin
Γ_top=3
Γ_bot=1
Γ_left=4
Γ_right=2
Γ_in=5
Γ_out=2
end;

# ╔═╡ 9d736062-6821-46d9-9e49-34b43b78e814
begin
    Ω_low=1
    Ω_high=2
end;

# ╔═╡ 83b9931f-9020-4400-8aeb-31ad391184db
function grid_2d(;nref=0,ε_fix=0.0)
    nx=10*2^nref
    ny=1*2^nref
    nylow=3*2^nref
    xc=linspace(0,L,nx+1)
    y0=linspace(-W/2,W/2,ny+1)
    if ε_fix>0.0
        yfix=[W/2,W/2+ε_fix]
	ytop=glue(yfix,linspace(yfix[end],Wlow,nylow+1))
    else
        ytop=linspace(W/2,Wlow,nylow+1)
    end
    yc=glue(-reverse(ytop),glue(y0,ytop))
    grid=simplexgrid(xc,yc)
    cellmask!(grid, [0,-W/2],[L,W/2],Ω_high)
    bfacemask!(grid, [0,-W/2],[0,W/2],Γ_in)
    bfacemask!(grid, [L,-W/2],[L,W/2],Γ_out)
end

# ╔═╡ 46a0f078-4165-4e37-9e69-e69af8584f6e
gridplot(grid_2d(),resolution=(400,300))

# ╔═╡ 3f693666-4026-4c01-a7aa-8c7dcbc32372
gridplot(grid_2d(;ε_fix=1.0e-1),resolution=(400,300))

# ╔═╡ c402f03c-746a-45b8-aaac-902a2f196094
function grid_1d(;nref=0)
    nx=10*2^nref
    xc=linspace(0,L,nx+1)
    grid=simplexgrid(xc)
    cellmask!(grid, [0],[L],Ω_high)
    bfacemask!(grid, [0],[0],Γ_in)
    bfacemask!(grid, [L],[L],Γ_out)
    grid
end

# ╔═╡ e33bf4b7-7b9a-4c5c-be4f-874af949c5ef
md"""
### Solvers
"""

# ╔═╡ a63a655c-e48b-4969-9409-31cd3db3bdaa
md"""
Pressure index in solution
"""

# ╔═╡ d7009231-4b43-44bf-96ba-9a203c0b5f5a
ip=1;

# ╔═╡ 26965e38-91cd-4022-bdff-4c503f724bfe
md"""
Concentration index in solution
"""

# ╔═╡ c904c921-fa10-43eb-bd46-b2869fa7f431
ic=2;

# ╔═╡ b143c846-2294-47f7-a2d1-8a6eabe942a3
md"""
Generate breaktrough courve from transient solution
"""

# ╔═╡ 92e4e4ab-3485-4cb9-9b41-e702a211a477
function breakthrough(sys,tf,sol)
	of=similar(sol.t)
	t=sol.t
	of[1]=0
	for i=2:length(sol.t)
	 of[i]=-integrate(sys,tf,sol[i],sol[i-1],t[i]-t[i-1])[ic]
	end
	of
end

# ╔═╡ 3df8bace-b4f1-4052-84f7-dff21d3a35f0
md"""
Transient solver:
"""

# ╔═╡ e866db69-9388-4691-99f7-879cf0658418
function trsolve(grid;
	κ=[1.0e-3,5],
	D=[1.0e-12,1.0e-12],
	Δp=1,
	ϕ=[1,1],
	tend=100)

    function flux(y,u,edge)
        y[ip]=κ[edge.region]*(u[ip,1]-u[ip,2])
	bp,bm=fbernoulli_pm(y[ip]/D[edge.region])
        y[ic]=D[edge.region]*(bm*u[ic,1]-bp*u[ic,2])
    end

    function stor(y,u,node)
        y[ip]=0
        y[ic]=ϕ[node.region]*u[ic]
    end
    physics=VoronoiFVM.Physics(flux=flux,storage=stor)
    sys=VoronoiFVM.System(grid,physics)
    enable_species!(sys,ip,[1,2])
    enable_species!(sys,ic,[1,2])
	dim=dim_space(grid)

	for ibc ∈ (dim == 1 ? [Γ_in] : [Γ_left,Γ_in])
		boundary_dirichlet!(sys,ip,ibc,Δp)
    end

	for ibc ∈ (dim == 1 ? [Γ_out] : [Γ_right,Γ_out])
      boundary_dirichlet!(sys,ip,ibc,0)
    end

    boundary_dirichlet!(sys,ic,Γ_in,0)
    boundary_dirichlet!(sys,ic,Γ_out,0)

    control=VoronoiFVM.NewtonControl()
    control.Δt=1.0e-4
    control.Δt_min=1.0e-6
    control.verbose=false

    inival=VoronoiFVM.solve(unknowns(sys,inival=0),sys;control=control)
    factory=VoronoiFVM.TestFunctionFactory(sys)
    tfc=testfunction(factory,[Γ_in,Γ_left,Γ_top,Γ_bot],[Γ_out])

    function pre(sol,t)
	cin=min(1.0,1.0e3*t)
	boundary_dirichlet!(sys,ic,Γ_in,cin)
    end

    sol=VoronoiFVM.solve(inival,sys,[0,tend]; pre=pre,control=control)

    bt=breakthrough(sys,tfc,sol)
    if dim==1
		bt=bt*W
	end

    grid,sol,bt
end

# ╔═╡ cd88123a-b042-43e2-99b9-ec925a8794ed
grid_n,sol_n,bt_n=trsolve(grid_2d(nref=nref),tend=tend);

# ╔═╡ 61ed054e-7116-4a51-9d13-e07279acb2a4
scalarplot!(vis_n,grid_n,sol_n(t)[ic,:],levels=0:0.2:1,resolution=(400,150),show=true)

# ╔═╡ b0ad0adf-6f6c-4fb3-b58e-e05cc8c0c796
grid_1,sol_1,bt_1=trsolve(grid_1d(nref=nref),tend=tend);

# ╔═╡ e36d2aef-1b5a-45a7-9289-8d1e544bcedd
scalarplot(grid_1,sol_1(t)[ic,:],levels=0:0.2:1,resolution=(500,150),
xlabel="x",ylabel="c",title="1D calculation, t=$t")

# ╔═╡ 76b77ec0-27b0-4a02-9ae4-43d756eb09dd
grid_f,sol_f,bt_f=trsolve(grid_2d(nref=nref,ε_fix=ε_fix),tend=tend);

# ╔═╡ d1ce7f69-3703-4cb2-9bf3-4c5efccb479e
scalarplot!(vis_f,grid_f,sol_f(t)[ic,:],levels=0:0.2:1,resolution=(400,150),show=true)

# ╔═╡ 904b36f0-10b4-4db6-9252-21668305de9c
grid_ϕ,sol_ϕ,bt_ϕ=trsolve(grid_2d(nref=nref), ϕ=[1.0e-3,1],tend=tend);

# ╔═╡ ce49bb25-b2d0-4d17-a8fe-d7b62e9b20be
let
    p=PlutoVistaPlot(resolution=(500,200),xlabel="t",ylabel="outflow",
                     legend=:rb,
                     title="Breakthrough Curves")
    plot!(p, sol_n.t,bt_n,label="naive grid")
    plot!(p, sol_1.t,bt_1,label="1D grid",markertype=:x)
    plot!(p, sol_f.t,bt_f,label="grid with fix",markertype=:circle)
    plot!(p, sol_ϕ.t,bt_ϕ,label="modified ϕ",markertype=:cross)
end

# ╔═╡ 78d92b4a-bdb1-4117-ab9c-b422eac403b1
md"""
Reaction-Diffusion solver
"""

# ╔═╡ bb3a50ed-32e7-4305-87d8-4093c054a4d2
function rdsolve(grid;D=[1.0e-12,1.0],R=[1,0.1])

    function flux(y,u,edge)
        y[1]=D[edge.region]*(u[1,1]-u[1,2])
    end

	function rea(y,u,node)
        y[1]=R[node.region]*u[1]
    end

    physics=VoronoiFVM.Physics(flux=flux,reaction=rea)
    sys=VoronoiFVM.System(grid,physics)
    enable_species!(sys,1,[1,2])
  	dim=dim_space(grid)

    boundary_dirichlet!(sys,1,Γ_in,1)
    boundary_dirichlet!(sys,1,Γ_out,0)

    sol=VoronoiFVM.solve(unknowns(sys,inival=0),sys)
    factory=VoronoiFVM.TestFunctionFactory(sys)
    tf=testfunction(factory,[Γ_in,Γ_left,Γ_top,Γ_bot],[Γ_out])
   	of=integrate(sys,tf,sol)
	    fac=1.0
	if dim==1
		fac=W
	end
    grid,sol[1,:],of[1]*fac

end

# ╔═╡ 2f560406-d169-4027-9cfe-7689494edf45
rdgrid_1,rdsol_1,of_1=rdsolve(grid_1d(nref=nref))

# ╔═╡ 34228382-4b1f-4897-afdd-19db7d5a7c59
scalarplot(rdgrid_1,rdsol_1,resolution=(300,200))

# ╔═╡ a6714eac-9e7e-4bdb-beb7-aca354664ad6
rdgrid_n,rdsol_n,of_n=rdsolve(grid_2d(nref=nref))

# ╔═╡ 20d7624b-f43c-4ac2-bad3-383a9e4e1b42
 rdgrid_f,rdsol_f,of_f=rdsolve(grid_2d(nref=nref,ε_fix=ε_fix))

# ╔═╡ 6a6d0e94-8f0d-4119-945c-dd48ec0798fd
scalarplot(rdgrid_n,rdsol_n,resolution=(210,200)),
scalarplot(rdgrid_f,rdsol_f,resolution=(210,200))

# ╔═╡ c0fc1f71-52ba-41a9-92d1-74e82ac7826c
 rdgrid_r,rdsol_r,of_r=rdsolve(grid_2d(nref=nref),R=[0,0.1])

# ╔═╡ c08e86f6-b5c2-4762-af23-382b1b153f45
md"""
We measure the outflow at the outlet. As a result, we obtain:
   - 1D case: $(of_1)
   - 2D case, naive grid: $(of_n)
   - 2D case, grid with "protective layer": $(of_f)
   - 2D case, naive grid, "modified" R: $(of_r)

"""

# ╔═╡ 0cc1c511-f351-421f-991a-a27f26a8db4f
  html"<hr>"

# ╔═╡ 99c8458a-a584-4825-a983-ae1a05e50000
md"""
"Un-markdown" the next cell for debugging/developing
"""

# ╔═╡ b6b826a1-b52f-41d3-8feb-b6464f76352e
md"""
  begin
    using Pkg
    Pkg.activate(".testenv")
    Pkg.add("Revise")
    using Revise
    Pkg.develop("VoronoiFVM")
    Pkg.add(["PlutoVista", "PlutoUI","GridVisualize","ExtendableGrids"])
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"
VoronoiFVM = "82b139dc-5afc-11e9-35da-9b9bdfd336f3"

[compat]
ExtendableGrids = "~0.8.7"
GridVisualize = "~0.3.9"
PlutoUI = "~0.7.16"
PlutoVista = "~0.8.7"
VoronoiFVM = "~0.13.2"
"""

# ╔═╡ 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-rc2"
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 = "d9352737cef8525944bf9ef34392d756321cbd54"
uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
version = "3.1.38"

[[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 = "3533f5a691e60601fe60c90d8bc47a27aa2907ec"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.11.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 = "dce3e3fea680869eaa0b774b2e8343e9ff442313"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.40.0"

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

[[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 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

[[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 = "1e8e50f054057f23e908fbd6935766dca6293cc2"
uuid = "cfc395e8-590f-11e8-1f13-43a2532b2fa8"
version = "0.8.7"

[[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 = "2db648b6712831ecb333eae76dbfd1c156ca13bb"
uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
version = "1.11.2"

[[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 = "925ba2f11df005d894b113292d32fca9afe3f8c8"
uuid = "5eed8a63-0fb0-45eb-886d-8d5a387d12b8"
version = "0.3.9"

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

[[deps.HypertextLiteral]]
git-tree-sha1 = "5efcf53d798efede8fee5b2c8b09284be359bf24"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.2"

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

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

[[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.InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "f0c6489b12d28fb4c2103073ec7452f3423bd308"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.1"

[[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.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", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "6193c3815f13ba1b78a51ce391db8be016ae9214"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.4"

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

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

[[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 = "d911b6a12ba974dabe2291c6d450094a7226b372"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.1.1"

[[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 = "b99d4e38e7dba4535cee937e0444aed5912245d0"
uuid = "646e1f28-b900-46d7-9d87-d554eb38a413"
version = "0.8.7"

[[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.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "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 = "c944fa4adbb47be43376359811c0a14757bdc8a8"
uuid = "731186ca-8d62-57ce-b412-fbd966d074cd"
version = "2.20.0"

[[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.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 = "999444c932d3bd32e7cb5a349e219bec9cfa2bd5"
uuid = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
version = "1.17.1"

[[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 = "f0bccf98e16759818ffc5d97ac3ebf87eb950150"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "1.8.1"

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

[[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.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

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

[[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 = "65f1d72aa575d3d6c348fcfc0d0d905fdbd7f5c3"
uuid = "82b139dc-5afc-11e9-35da-9b9bdfd336f3"
version = "0.13.2"

[[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:
# ╠═60941eaa-1aea-11eb-1277-97b991548781
# ╠═556480e0-94f1-4e47-be9a-3e1e0e99555c
# ╟─5e13b3db-570c-4159-939a-7e2268f0a102
# ╟─fae47c55-eef8-4428-bb5f-45824978753d
# ╟─8ba2300c-17ff-44e1-b33a-c5bdf1ce12fe
# ╟─51c9517c-8797-4406-b053-301694fb0484
# ╟─99341e32-9c78-4e31-bec0-d1ffbc85ec32
# ╟─46a0f078-4165-4e37-9e69-e69af8584f6e
# ╟─cc325b2c-6174-4b8d-8e39-202ac68b5705
# ╠═3f693666-4026-4c01-a7aa-8c7dcbc32372
# ╟─cd013964-f329-4d2c-ae4b-305093f0ac56
# ╠═afccbadb-2ca8-4c3e-8c6d-c78df59d8d7e
# ╠═dd9f8d38-8812-40ba-88c8-f873ec7d6121
# ╠═5f6ac608-b1a0-450e-910e-d7d8ea2ffae0
# ╠═cd88123a-b042-43e2-99b9-ec925a8794ed
# ╠═b0ad0adf-6f6c-4fb3-b58e-e05cc8c0c796
# ╠═76b77ec0-27b0-4a02-9ae4-43d756eb09dd
# ╠═904b36f0-10b4-4db6-9252-21668305de9c
# ╟─ce49bb25-b2d0-4d17-a8fe-d7b62e9b20be
# ╟─5b60c7d4-7bdb-4989-b055-6695b9fdeedc
# ╟─f6abea66-1e42-4201-8433-5d092989749d
# ╟─e36d2aef-1b5a-45a7-9289-8d1e544bcedd
# ╟─98ae56dd-d42d-4a93-bb0b-5956b6e981a3
# ╟─61ed054e-7116-4a51-9d13-e07279acb2a4
# ╟─d1ce7f69-3703-4cb2-9bf3-4c5efccb479e
# ╟─99c3b54b-d458-482e-8aa0-d2c2b51fdf25
# ╟─eef85cd7-eba4-4c10-9e1d-38411179314d
# ╠═2f560406-d169-4027-9cfe-7689494edf45
# ╠═a6714eac-9e7e-4bdb-beb7-aca354664ad6
# ╠═20d7624b-f43c-4ac2-bad3-383a9e4e1b42
# ╠═c0fc1f71-52ba-41a9-92d1-74e82ac7826c
# ╟─c08e86f6-b5c2-4762-af23-382b1b153f45
# ╠═34228382-4b1f-4897-afdd-19db7d5a7c59
# ╟─6a6d0e94-8f0d-4119-945c-dd48ec0798fd
# ╟─fcd066f1-bcd8-4479-a4e4-7b8c235336c4
# ╟─c9d92201-813c-499b-b863-b138c30eb634
# ╟─a372ac90-c871-4dc0-a44b-a5bddef71823
# ╟─124b2a0a-ef19-453e-9e3a-5b5ce7db5fac
# ╠═1ad18670-e7cb-4f7a-be0f-3db98cdeb6a4
# ╟─47bc8e6a-e296-42c9-bfc5-967edfb0feb7
# ╠═d1d5bad2-d282-4e7d-adb9-baf21f58155e
# ╠═9d736062-6821-46d9-9e49-34b43b78e814
# ╠═83b9931f-9020-4400-8aeb-31ad391184db
# ╠═c402f03c-746a-45b8-aaac-902a2f196094
# ╟─e33bf4b7-7b9a-4c5c-be4f-874af949c5ef
# ╟─a63a655c-e48b-4969-9409-31cd3db3bdaa
# ╠═d7009231-4b43-44bf-96ba-9a203c0b5f5a
# ╟─26965e38-91cd-4022-bdff-4c503f724bfe
# ╠═c904c921-fa10-43eb-bd46-b2869fa7f431
# ╟─b143c846-2294-47f7-a2d1-8a6eabe942a3
# ╠═92e4e4ab-3485-4cb9-9b41-e702a211a477
# ╟─3df8bace-b4f1-4052-84f7-dff21d3a35f0
# ╠═e866db69-9388-4691-99f7-879cf0658418
# ╟─78d92b4a-bdb1-4117-ab9c-b422eac403b1
# ╠═bb3a50ed-32e7-4305-87d8-4093c054a4d2
# ╟─0cc1c511-f351-421f-991a-a27f26a8db4f
# ╟─99c8458a-a584-4825-a983-ae1a05e50000
# ╠═b6b826a1-b52f-41d3-8feb-b6464f76352e
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
# This merger stream is always routed to a single synapse.
mutable struct StreamMerger <: AbstractBitStream
    base::BaseData{UInt8}

    # Inputs
    streams::Array{AbstractBitStream,1}

    function StreamMerger()
        o = new()

        o.streams = Array{AbstractBitStream,1}()
        o.base = BaseData{UInt8}()
        o
    end
end

# Inputs
function add_stream!(merger::StreamMerger, stream::AbstractBitStream)
    push!(merger.streams, stream)
end

function output(merger::StreamMerger)
    merger.base.output
end

function step!(merger::StreamMerger)
    # Combine each stream's output into a single value and send directly to output.
    # Generally the merger holds two or more streams: Poisson and PatternChannel
    merger.base.output = 0
    for stream in merger.streams
        # Merge streams.
        merger.base.output = merger.base.output | output(stream)
    end
end

ext: jl


text:
"""
    ESGPtrain(esn::AbstractLeakyESN, mean::GaussianProcesses.Mean, kernel::GaussianProcesses.Kernel
    [, lognoise::Float64, optimize::Bool, optimizer::Optim.AbstractOptimizer, y_target::AbstractArray{Float64})

Train the ESN using Gaussian Processes, as described in [1]

[1] Chatzis, Sotirios P., and Yiannis Demiris. "Echo state Gaussian process." IEEE Transactions on Neural Networks 22.9 (2011): 1435-1445.
"""
function ESGPtrain(esn::AbstractLeakyESN,
    mean::GaussianProcesses.Mean,
    kernel::GaussianProcesses.Kernel;
    lognoise::Float64 = -2.0,
    optimize::Bool = false,
    optimizer::Optim.AbstractOptimizer = Optim.LBFGS(),
    y_target::AbstractArray{Float64} = esn.train_data)

    states_new = nla(esn.nla_type, esn.states)
    if size(y_target, 1) == 1
        gp = GP(states_new, vec(y_target), mean, kernel, lognoise)
        if optimize == true
            optimize!(gp; method=optimizer)
        end
    else
        gp = []
        for i=1:size(y_target, 1)
            push!(gp, GP(states_new, y_target[i,:], mean, kernel, lognoise))
            if optimize == true
                optimize!(gp[i]; method=optimizer)
            end
        end
    end
    return gp
end

#one variable timestep prediction
function ESGPpredict(esn::AbstractLeakyESN,
    predict_len::Int,
    gp::GaussianProcesses.GPE)

    output = zeros(Float64, esn.in_size, predict_len)
    sigmas = zeros(Float64, esn.in_size, predict_len)
    x = esn.states[:, end]

    if esn.extended_states == false
        for i=1:predict_len
            x_new = hcat(nla(esn.nla_type, x)...)'
            out, sigma = GaussianProcesses.predict_y(gp, x_new)
            output[:, i] = out
            sigmas[:,i] = sigma
            x = leaky_fixed_rnn(esn.activation, esn.alpha, esn.W, esn.W_in, x, out)
        end
    else
        for i=1:predict_len
            x_new = hcat(nla(esn.nla_type, x)...)'
            out, sigma = GaussianProcesses.predict_y(gp, x_new)
            output[:, i] = out
            sigmas[:,i] = sigma
            x = vcat(leaky_fixed_rnn(esn.activation, esn.alpha, esn.W, esn.W_in, x[1:esn.res_size], out), out)
        end
    end
    return output, sigmas
end

#multidimensional timestep prediction

"""
    ESGPpredict(esn::AbstractLeakyESN, predict_len::Int, gp::AbstractArray{Any})

Return the prediction for a given lenght of the constructed ESN struct using GPs.

"""
function ESGPpredict(esn::AbstractLeakyESN,
    predict_len::Int,
    gp::AbstractArray{Any})

    output = zeros(Float64, esn.in_size, predict_len)
    sigmas = zeros(Float64, esn.in_size, predict_len)
    x = esn.states[:, end]

    if esn.extended_states == false
        for i=1:predict_len
            x_new = hcat(nla(esn.nla_type, x)...)'
            out = Array{Float64}(undef, esn.out_size)
            sigma = Array{Float64}(undef, esn.out_size)
            for i=1:size(gp, 1)
                out[i] = GaussianProcesses.predict_y(gp[i], x_new)[1][1]
                sigma[i] = GaussianProcesses.predict_y(gp[i], x_new)[2][1]
            end
            output[:, i] = out
            sigmas[:,i] = sigma
            x = leaky_fixed_rnn(esn.activation, esn.alpha, esn.W, esn.W_in, x, out)
        end
    else
        for i=1:predict_len
            x_new = hcat(nla(esn.nla_type, x)...)'
            out = Array{Float64}(undef, esn.out_size)
            sigma = Array{Float64}(undef, esn.out_size)
            for i=1:size(gp, 1)
                out[i] = GaussianProcesses.predict_y(gp[i], x_new)[1][1]
                sigma[i] = GaussianProcesses.predict_y(gp[i], x_new)[2][1]
            end
            output[:, i] = out
            sigmas[:,i] = sigma
            x = vcat(leaky_fixed_rnn(esn.activation, esn.alpha, esn.W, esn.W_in, x[1:esn.res_size], out), out)
        end
    end
    return output, sigmas
end

"""
    ESGPpredict_h_steps(esn::AbstractLeakyESN, predict_len::Int, h_steps::Int, test_data::AbstractArray{Float64}, gp::GaussianProcesses.GPE)

Return the prediction for h steps ahead of the constructed ESN struct using GPs.
"""
function ESGPpredict_h_steps(esn::AbstractLeakyESN,
    predict_len::Int,
    h_steps::Int,
    test_data::AbstractArray{Float64},
    gp::GaussianProcesses.GPE)

    output = zeros(Float64, esn.in_size, predict_len)
    sigmas = zeros(Float64, esn.in_size, predict_len)
    x = esn.states[:, end]

    if esn.extended_states == false
        for i=1:predict_len
            x_new = hcat(nla(esn.nla_type, x)...)'
            out, sigma = GaussianProcesses.predict_y(gp, x_new)
            output[:, i] = out
            sigmas[:,i] = sigma
            if mod(i, h_steps) == 0
                x = leaky_fixed_rnn(esn.activation, esn.alpha, esn.W, esn.W_in, x, test_data[:,i])
            else
                x = leaky_fixed_rnn(esn.activation, esn.alpha, esn.W, esn.W_in, x, out)
            end
        end
    else
        for i=1:predict_len
            x_new = hcat(nla(esn.nla_type, x)...)'
            out, sigma = GaussianProcesses.predict_y(gp, x_new)
            output[:, i] = out
            sigmas[:,i] = sigma
            if mod(i, h_steps) == 0
                x = vcat(leaky_fixed_rnn(esn.activation, esn.alpha, esn.W, esn.W_in, x[1:esn.res_size], test_data[:,i]), test_data[:,i])
            else
                x = vcat(leaky_fixed_rnn(esn.activation, esn.alpha, esn.W, esn.W_in, x[1:esn.res_size], out), out)
            end
        end
    end
    return output, sigmas
end


text:
function order_disjoint{T<:AbstractArray}(m::T, n::T)
    rlen = length(n)
    rdis = zeros(Int, rlen)
    for (i, e) in enumerate(n)
        j = findfirst(m, e)
        while j in rdis && j != 0
            j = findnext(m, e, j+1)
        end
        rdis[i] = j
    end
    if 0 in rdis
        throw(DomainError())
    end
    sort!(rdis)
    p = copy(m)
    p[rdis] = n
    return p
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@recipe function f(solution::AbstractRungeKuttaSolution; vars = nothing, is_complex = false)
    framestyle     --> :box
    gridalpha      --> 0.2
    legend         --> :none
    linewidth      --> 1.5
    minorgrid      --> 0.1
    minorgridstyle --> :dash
    seriestype     --> :path
    # tick_direction --> :out
    (L, N) = size(solution)
    if vars isa Tuple && length(vars) ≥ 1
        vars = vars
    elseif vars isa Integer
        vars = (vars, )
    elseif vars isa Nothing
        vars = 1:L
    else
        error("Got $(typeof(vars)) instead of `Tuple` with `length(vars)` ≥ 1 or `Integer`.")
    end
    @↓ u, t = solution
    if is_complex
        xwiden --> true
        for i in vars
            @series [u[n][i] for n = 1:N]
        end
    else
        xwiden --> false
        [(t, [u[n][i] for n = 1:N]) for i in vars]
    end
end

@userplot PHASEPLOT
@recipe function f(h::PHASEPLOT; vars = nothing)
    length(h.args) == 1 ? true : error("Got too many arguments: $(length(h.args)).")
    solution = h.args[1] isa AbstractRungeKuttaSolution ? h.args[1] : error("Got $(typeof(h.args)) instead of `AbstractRungeKuttaSolution`.")
    framestyle     --> :box
    gridalpha      --> 0.2
    legend         --> :none
    linewidth      --> 1.5
    minorgrid      --> 0.1
    minorgridstyle --> :dash
    seriestype     --> :path
    # tick_direction --> :out
    (L, N) = size(solution)
    if vars isa Tuple && length(vars) ≥ 1
        vars = vars
    elseif vars isa Integer
        vars = (vars, )
    elseif vars isa Nothing
        vars = 1:L
    else
        error("Got $(typeof(vars)) instead of `Tuple` with `length(vars)` ≥ 1 or `Integer`.")
    end
    @↓ u, t = solution
    tuple([[u[n][i] for n = 1:N] for i in vars]...)
end

@userplot STABILITY
@recipe function f(h::STABILITY; lims = (-4, 4), xlims = lims, ylims = lims)
    length(h.args) == 1 ? true : error("Got too many arguments: $(length(h.args)).")
    R = h.args[1] isa ButcherTableau   ? z -> ℛ(z, h.args[1])         :
        h.args[1] isa AbstractRungeKuttaSolver ? z -> ℛ(z, h.args[1].tableau) :
        h.args[1] isa Function         ? h.args[1]                    :
        error("Got $(typeof(h.args)) instead of `ButcherTableau`, `AbstractRungeKuttaSolver` or `Function`.")
    framestyle     --> :box
    gridalpha      --> 0.2
    gridstyle      --> :dot
    legend         --> :none
    levels         --> [1.0]
    linewidth      --> 1.5
    minorgrid      --> 0.1
    minorgridstyle --> :dash
    seriestype     --> :contour
    # tick_direction --> :out
    Δx = max(1, abs(xlims[2] - xlims[1]))
    Δy = max(1, abs(ylims[2] - ylims[1]))
    x = length(xlims) == 2 ? LinRange(xlims..., Int(Δx * 101)) : LinRange(xlims...)
    y = length(ylims) == 2 ? LinRange(ylims..., Int(Δy * 101)) : LinRange(ylims...)
    function f(x, y)
        z = x + 1im * y
        p = abs(R(z))
    end
    x, y, f
end

@userplot STABILITYF
@recipe function f(h::STABILITYF; lims = (-4, 4), xlims = lims, ylims = lims)
    length(h.args) == 1 ? true : error("Got too many arguments: $(length(h.args)).")
    R = h.args[1] isa ButcherTableau   ? z -> ℛ(z, h.args[1])         :
        h.args[1] isa AbstractRungeKuttaSolver ? z -> ℛ(z, h.args[1].tableau) :
        h.args[1] isa Function         ? h.args[1]                    :
        error("Got $(typeof(h.args)) instead of `ButcherTableau`, `AbstractRungeKuttaSolver` or `Function`.")
    clims          --> (0, 1)
    colorbar       --> true
    framestyle     --> :box
    gridalpha      --> 0.2
    legend         --> :none
    minorgrid      --> 0.1
    minorgridstyle --> :dash
    seriestype     --> :heatmap
    # tick_direction --> :out
    Δx = max(1, abs(xlims[2] - xlims[1]))
    Δy = max(1, abs(ylims[2] - ylims[1]))
    x = length(xlims) == 2 ? LinRange(xlims..., Int(Δx * 101)) : LinRange(xlims...)
    y = length(ylims) == 2 ? LinRange(ylims..., Int(Δy * 101)) : LinRange(ylims...)
    function f(x, y)
        z = x + 1im * y
        p = abs(R(z))
        (p > 1 ? -Inf : p)
    end
    x, y, f
end

@userplot ORDERSTAR
@recipe function f(h::ORDERSTAR; lims = (-4, 4), xlims = lims, ylims = lims)
    length(h.args) == 1 ? true : error("Got too many arguments: $(length(h.args)).")
    R = h.args[1] isa ButcherTableau   ? z -> ℛ(z, h.args[1])         :
        h.args[1] isa AbstractRungeKuttaSolver ? z -> ℛ(z, h.args[1].tableau) :
        h.args[1] isa Function         ? h.args[1]                    :
        error("Got $(typeof(h.args)) instead of `ButcherTableau`, `AbstractRungeKuttaSolver` or `Function`.")
    framestyle     --> :box
    gridalpha      --> 0.2
    legend         --> :none
    levels         --> [1.0]
    linewidth      --> 1.5
    minorgrid      --> 0.1
    minorgridstyle --> :dash
    seriestype     --> :contour
    # tick_direction --> :out
    Δx = max(1, abs(xlims[2] - xlims[1]))
    Δy = max(1, abs(ylims[2] - ylims[1]))
    x = length(xlims) == 2 ? LinRange(xlims..., Int(Δx * 101)) : LinRange(xlims...)
    y = length(ylims) == 2 ? LinRange(ylims..., Int(Δy * 101)) : LinRange(ylims...)
    function f(x, y)
        z = x + 1im * y
        abs(R(z) * exp(-z))
    end
    x, y, f
end

@userplot ORDERSTARF
@recipe function f(h::ORDERSTARF; lims = (-4, 4), xlims = lims, ylims = lims)
    length(h.args) == 1 ? true : error("Got too many arguments: $(length(h.args)).")
    R = h.args[1] isa ButcherTableau   ? z -> ℛ(z, h.args[1])         :
        h.args[1] isa AbstractRungeKuttaSolver ? z -> ℛ(z, h.args[1].tableau) :
        h.args[1] isa Function         ? h.args[1]                    :
        error("Got $(typeof(h.args)) instead of `ButcherTableau`, `AbstractRungeKuttaSolver` or `Function`.")
    clims          --> (0, 1)
    colorbar       --> true
    framestyle     --> :box
    gridalpha      --> 0.2
    legend         --> :none
    minorgrid      --> 0.1
    minorgridstyle --> :dash
    seriestype     --> :heatmap
    # tick_direction --> :out
    Δx = max(1, abs(xlims[2] - xlims[1]))
    Δy = max(1, abs(ylims[2] - ylims[1]))
    x = length(xlims) == 2 ? LinRange(xlims..., Int(Δx * 101)) : LinRange(xlims...)
    y = length(ylims) == 2 ? LinRange(ylims..., Int(Δy * 101)) : LinRange(ylims...)
    function f(x, y)
        z = x + 1im * y
        p = abs(R(z) * exp(-z))
        (p > 1 ? NaN : p)
    end
    x, y, f
end


text:
using Pkg
Pkg.activate(".")
using CSV
include("./parser/julia103/powermodels_dict.jl")
include("./parser/julia103/hydro_dict.jl")

#Params
case = "brasil_4"
case_dir = dirname(@__FILE__)

# Load case data
hydro_pre = CSV.read(joinpath(case_dir,"hydro.csv"))
lines = CSV.read(joinpath(case_dir,"trans.csv"))
thermal_gen = CSV.read(joinpath(case_dir,"thermal.csv"))
demand = CSV.read(joinpath(case_dir,"demand.csv"))

#------------------------------------------------------------
# Parse Data to PowerModels JSON
#------------------------------------------------------------
# prepocess
busName = unique(vcat(lines[:f],lines[:t]))
gen2bus = [findall(x-> x== i,busName)[1] for i in thermal_gen[:bus]]
lines[:f] = [findall(x-> x== i,busName)[1] for i in lines[:f]]
lines[:t] = [findall(x-> x== i,busName)[1] for i in lines[:t]]
demand[:Tr] = fill(0, size(demand,1))
bus2load = [demand[1,Symbol(bus)] for bus in busName]
# Create Dicts
gen_array= [gen_dict( gen_bus= gen2bus[idx],
    pmax= thermal_gen[:max][idx]/100,
    index= idx,
    cost= [
        thermal_gen[:cost][idx]*100,
        0
    ])
for idx = 1:size(thermal_gen,1)]

bus_array = [bus_dict( name = busName[idx],
        bus_i= idx,
        bus_type= idx==1 ? 3 : idx in unique(gen2bus) ? 2 : 1,
        index= idx)
for idx = 1:size(busName,1)]

branch_array = [branch_dict(
            br_r= lines[:r][idx],
            rate_a= lines[Symbol(" max")][idx]/100,
            f_bus= lines[:f][idx],
            t_bus= lines[:t][idx],
            index= idx)
for idx = 1:size(lines,1)]

lod_array = [load_dict(load_bus= idx,
                    pd= bus2load[idx]/100,
                    index= idx)
    for idx=1:size(bus2load,1)]

#------------------------------------------------------------
# Parse Data to Hidro JSON
#------------------------------------------------------------
hydroarray = [hydrogenerator_dict(
        index = idx,
        index_grid= idx+size(thermal_gen,1),
        max_volume= hydro_pre[:Max_store][idx]*0.0036,
        min_volume= 0,
        max_turn = hydro_pre[:Max_discharge][idx],
        min_turn = 0,
        initial_volume= hydro_pre[:Init_store][idx]*0.0036,
        production_factor= 1,
        spill_cost= 0,
        downstream_turn  = Int64[],
        downstream_spill = Int64[])
        for idx=1:size(hydro_pre,1)]

# Create Dicts gen
hydrogen_array= [gen_dict(gen_bus= [findall(x-> x== i,busName)[1] for i in hydro_pre[:bus]][idx],
    pmax= hydro_pre[:Max_discharge][idx]/100,
    index= idx+size(thermal_gen,1),
    cost= [0])
for idx=1:size(hydro_pre,1)]

#------------------------------------------------------------
# Final dicts
#------------------------------------------------------------

# Final Dict PowerModels
pm= powermodels_dict(name=case,cost_deficit=2300,gen=[gen_array;hydrogen_array],bus=bus_array,branch=branch_array,load=lod_array)

# Final Dict Hydro
hm = hydro_dict(hydrogenarray=hydroarray)

#------------------------------------------------------------
# Save Files
#------------------------------------------------------------

using JSON
pathcase_out = case_dir

open(joinpath(pathcase_out,"PowerModels.json"),"w") do f
    JSON.print(f, pm)
end

open(joinpath(pathcase_out,"hydro.json"),"w") do f
    JSON.print(f, hm)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Raylib"
version = v"4.0.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/raysan5/raylib.git",
              "0851960397f02a477d80eda2239f90fae14dec64"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/raylib/src/
atomic_patch -p1 ../../patches/add-missing-header.patch
atomic_patch -p1 ../../patches/make-install-everywhere.patch
export CFLAGS="-D_POSIX_C_SOURCE=200112L -I${includedir}"
if [[ "${target}" == *-freebsd* ]]; then
    # Allow definition of `u_char`, `u_short`, `u_int`, and `u_long` in sys/types.h
    CFLAGS="${CFLAGS} -D__BSD_VISIBLE"
fi
CFLAGS="${CFLAGS} -DSUPPORT_EVENTS_AUTOMATION -DSUPPORT_FILEFORMAT_BMP -DSUPPORT_FILEFORMAT_JPG"
make -j${nproc} USE_EXTERNAL_GLFW=TRUE PLATFORM=PLATFORM_DESKTOP RAYLIB_LIBTYPE=SHARED RAYLIB_MODULE_RAYGUI=TRUE RAYLIB_MODULE_PHYSAC=TRUE
make install RAYLIB_LIBTYPE=SHARED DESTDIR="${prefix}" RAYLIB_INSTALL_PATH="${libdir}"
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true, exclude=p->arch(p)=="armv6l" || Sys.iswindows(p))

# The products that we will ensure are always built
products = [
    LibraryProduct("libraylib", :libraylib)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="alsa_jll", uuid="45378030-f8ea-5b20-a7c7-1a9d95efb90e"))
    Dependency(PackageSpec(name="Mesa_jll", uuid="78dcde23-ec64-5e07-a917-6fe22bbc0f45"))
    Dependency(PackageSpec(name="Xorg_libX11_jll", uuid="4f6342f7-b3d2-589e-9d20-edeb45f2b2bc"))
    Dependency(PackageSpec(name="Xorg_libXrandr_jll", uuid="ec84b674-ba8e-5d96-8ba1-2a689ba10484"))
    Dependency(PackageSpec(name="Xorg_libXi_jll", uuid="a51aa0fd-4e3c-5386-b890-e753decda492"))
    Dependency(PackageSpec(name="Xorg_libXcursor_jll", uuid="935fb764-8cf2-53bf-bb30-45bb1f8bf724"))
    Dependency(PackageSpec(name="Xorg_libXinerama_jll", uuid="d1454406-59df-5ea1-beac-c340f2130bc3"))
    BuildDependency(PackageSpec(name="Xorg_xorgproto_jll", uuid="c4d99508-4286-5418-9131-c86396af500b"))
    Dependency(PackageSpec(name="GLFW_jll", uuid="0656b61e-2033-5cc2-a64a-77c0f6c09b89"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# Sample a specific integral point in the grid
function GlobalApproximationValueIteration.sample_state(mdp::SimpleGridWorld, rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}

    x = rand(rng, 1:mdp.size[1])
    y = rand(rng, 1:mdp.size[2])

    return GWPos(x, y)
end

function test_against_full_grid()

    rng = MersenneTwister(2378)

    # Attempt to approximate globally with N samples and M iterations. As N increases, the average error should decrease
    MAX_ITERS = 100
    NUM_SAMPLES_LOW = 50
    NUM_SAMPLES_HI = 500

    # Grid probabilities
    SIZE_X = 10
    SIZE_Y = 10
    REWARD_COV_PROB = 0.3

    # Generate a block of reward states from 40,40 to 60,60
    rewards = Dict{GWPos, Float64}()
    for x = 1:SIZE_X
        for y = 1:SIZE_Y
            if rand(rng) < REWARD_COV_PROB
                rewards[GWPos(x, y)] = 10
            end
        end
    end

    # Create MDP
    mdp = SimpleGridWorld(size=(SIZE_X, SIZE_Y), rewards=rewards)

    # Define learning model
    model1 = Chain(
                Dense(2, 10, relu),
                Dense(10, 5, relu),
                Dense(5, 1))
    model2 = Chain(
                Dense(2, 10, relu),
                Dense(10, 5, relu),
                Dense(5, 1))
    opt = ADAM(0.001)

    nonlin_gfa_1 = NonlinearGlobalFunctionApproximator(model1, opt, Flux.mse)
    nonlin_gfa_2 = NonlinearGlobalFunctionApproximator(model2, opt, Flux.mse)

    solver_low = GlobalApproximationValueIterationSolver(nonlin_gfa_1; num_samples=NUM_SAMPLES_LOW, max_iterations=MAX_ITERS, verbose=true)
    solver_hi = GlobalApproximationValueIterationSolver(nonlin_gfa_2; num_samples=NUM_SAMPLES_HI, max_iterations=MAX_ITERS, verbose=true)

    policy_low = solve(solver_low, mdp)
    policy_hi = solve(solver_hi, mdp)

    # Now solve with dicrete VI
    solver = ValueIterationSolver(max_iterations=1000, verbose=true)
    policy = solve(solver, mdp)

    total_err_low = 0.0
    total_err_hi = 0.0

    for state in states(mdp)

        full_val = value(policy, state)

        approx_val_low = value(policy_low, state)
        approx_val_hi = value(policy_hi, state)

        total_err_low += abs(full_val-approx_val_low)
        total_err_hi += abs(full_val-approx_val_hi)
    end

    avg_err_low = total_err_low / length(states(mdp))
    avg_err_hi = total_err_hi / length(states(mdp))

    @show avg_err_low
    @show avg_err_hi

    return (avg_err_low > avg_err_hi)
end

@test_broken test_against_full_grid() == true


text:
using Images
using StaticArrays

@testset "Visualization" begin
    @info "Running flow visualization test."
    flow = Array{SVector{2, Float64}, 2}(undef,(3, 2))

    flow[1] = SVector{2, Float64}(0.5, 0.5)
    flow[2] = SVector{2, Float64}(-0.4, -0.8)
    flow[3] = SVector{2, Float64}(-0.5, 0.3)
    flow[4] = SVector{2, Float64}(-0.4, -0.8)
    flow[5] = SVector{2, Float64}(-0.5, 0.3)
    flow[6] = SVector{2, Float64}(0.5, 0.5)

    @time hsv = visualize_flow(flow, ColorBased(), RasterConvention())

    @test hsv[1] == hsv[6]
    @test hsv[2] == hsv[4]
    @test hsv[3] == hsv[5]

end


text:
module AutoPages

export gather_pages

function replace_reverse(x, p...; count::Union{Integer,Nothing}=nothing)
  from = getproperty.(p, :first)
  to = getproperty.(p, :second)
  reversed_pairs = Pair.(reverse.(from), reverse.(to))
  x_reversed = reverse(x)
  if count==nothing
    return String(collect(Iterators.reverse(replace(x_reversed, reversed_pairs...))))
  else
    return String(collect(Iterators.reverse(replace(x_reversed, reversed_pairs...; count=count))))
  end
end

if Sys.isunix()
    const path_separator = '/'
elseif Sys.iswindows()
    const path_separator = '\\'
else
    error("path_separator for this OS need to be defined")
end

"""
    transform_page(filename)

Transform filename to title
seen in the navigation panel.
"""
function transform_page(filename)
  titlename = String(first(splitext(filename)))
  titlename = replace(titlename, "_" => " ")
  titlename = titlecase(titlename)
  return titlename
end

"""
    transform_subpage(local_path)

Transform local path to sub-page
title, seen in the navigation panel.
"""
function transform_subpage(local_path)
  return join(split_by_camel_case(local_path), " ")
end

"""
    split_by_camel_case(s::AbstractString)

Splits a string into an array of strings,
by its `CamelCase`. Examples:

```julia
julia> split_by_camel_case("CamelCase")
2-element Array{SubString{String},1}:
 "Camel"
 "Case"
julia> split_by_camel_case("PDESolver")
2-element Array{SubString{String},1}:
 "PDE"
 "Solver"
```
"""
split_by_camel_case(obj::AbstractString) =
    split(obj, r"((?<=\p{Ll})(?=\p{Lu})|(?<=\p{Lu})(?=\p{Lu}[^\p{Lu}]+))")

"""
    gather_pages!(array::Array,
                  filename,
                  folders_with_only_files,
                  fullpath = dirname(filename))

Construct a nested array of `Pair`s whose
keys mirror the given folder structure and
values which point to the files.
"""
function gather_pages!(
    array::Array,
    filename::S,
    folders_with_only_files::Vector{S},
    transform_page::Function,
    transform_subpage::Function,
    fullpath = dirname(filename),
) where {S<:AbstractString}
    paths = splitpath(filename)
    if length(paths) == 1
        file = paths[1]
        titlename = transform_page(file)
        push!(array, Pair(String(titlename), String(joinpath(fullpath, file))))
    else
        local_path = paths[1]
        key = transform_subpage(local_path)
        file_next = joinpath(paths[2:end]...)
        if !any(x.first == key for x in array)
            push!(array, Pair(key, Any[]))
        end
        arr_next = array[end].second
        gather_pages!(arr_next,
                      file_next,
                      folders_with_only_files,
                      transform_page,
                      transform_subpage,
                      fullpath)
    end
end

"""
    gather_pages(filenames::Array;
                 remove_first_level::Bool = false,
                 transform_page::Function=transform_page,
                 transform_subpage::Function=transform_subpage)

Construct a nested array of `Pair`s whose
keys mirror the given folder structure and
values which point to the files.

This was specifically designed for
Documenter.jl's `page` array.
See the tests for an example.

 - `directory` to get filenames from walk
 - `filenames` Array of `.md` files (using relative paths)
 - `remove_first_level` Bool indicating whether or not
                        to remove the root directory
                        (e.g., `generated/`)
 - `transform_page` function to transform the page (title) in navigation panel
 - `transform_subpage` function to transform the sub-page in navigation panel
 - `extension_filter` filter extensions
 - `transform_extension` transform extensions (e.g., ".jl" to ".md")
                         a useful one is `AutoPages.replace_reverse(x, ".jl" => ".md"; count=1)`
 - `transform_path` transform path
"""
function gather_pages(;
    directory::Union{AbstractString,Nothing}=nothing,
    filenames::Union{Array,Nothing}=nothing,
    remove_first_level::Bool = false,
    transform_page::Function=transform_page,
    transform_subpage::Function=transform_subpage,
    extension_filter=x->endswith(x, ".md"),
    transform_extension=x->x,
    transform_path=x->x,
    )

    if filenames === nothing && directory === nothing
      throw(ArgumentError("Need filenames or directory keyword"))
    end
    if filenames ≠ nothing && directory ≠ nothing
      throw(ArgumentError("Too many kwargs given"))
    end

    if directory ≠ nothing
      filenames = [String(joinpath(r, f)) for (r, _, files) in Base.Filesystem.walkdir(directory) for f in files]
      filenames = filter(extension_filter, filenames)
      filenames = map(x -> String(last(split(x, dirname(directory)))), filenames)
    end

    filenames = map(x -> transform_extension(x), filenames)
    filter!(x -> !(x == path_separator), filenames)
    filenames = map(x -> String(lstrip(x, path_separator)), filenames)
    filenames = map(x -> transform_path(x), filenames)
    dirnames = collect(Set(dirname.(filenames)))

    # TODO: Might be able to improve performance here
    dirnames = [x for x in dirnames if !any(occursin(x, y) && !(x == y) for y in dirnames)]
    folders_with_only_files = basename.(dirnames)
    array = Any[]
    for file in filenames
        gather_pages!(array,
                      file,
                      folders_with_only_files,
                      transform_page,
                      transform_subpage)
    end

    if remove_first_level
        array = array[1].second
    end

    return array, filenames
end

end


text:
module AtmosphericStates
using CUDA
using ..Device: array_type, array_device, CUDADevice, CPU
using DocStringExtensions
using Adapt
using GaussQuadrature
using CLIMAParameters
using CLIMAParameters.Planet: grav, R_d

using ..Vmrs

export AbstractAtmosphericState,
    AtmosphericState,
    ClearAtmosphericState,
    GrayAtmosphericState,
    setup_gray_as_pr_grid,
    setup_gray_as_alt_grid

abstract type AbstractAtmosphericState{FT,I,FTA1D} end

include("GrayAtmosphericStates.jl")

"""
    AtmosphericState{FT,FTA1D,FTA1DN,FTA2D,CLDP,CLDM,VMR,I} <:
        AbstractAtmosphericState{FT,I,FTA1D}

Atmospheric conditions, used to compute optical properties.

# Fields
$(DocStringExtensions.FIELDS)
"""
struct AtmosphericState{
    FT<:AbstractFloat,
    FTA1D<:AbstractArray{FT,1},
    FTA1DN<:Union{AbstractArray{FT,1},Nothing},
    FTA2D<:AbstractArray{FT,2},
    CLDP<:Union{AbstractArray{FT,2},Nothing},
    CLDM<:Union{AbstractArray{Bool,2},Nothing},
    VMR<:AbstractVmr{FT},
    I<:Int,
} <: AbstractAtmosphericState{FT,I,FTA1D}
    "longitude (`ncol`), optional"
    lon::FTA1DN
    "latitude (`ncol`), optional"
    lat::FTA1DN
    "Layer pressures `[Pa, mb]`; `(nlay,ncol)`"
    p_lay::FTA2D
    "Level pressures `[Pa, mb]`; `(nlay+1,ncol)`"
    p_lev::FTA2D
    "Layer temperatures `[K]`; `(nlay,ncol)`"
    t_lay::FTA2D
    "Level temperatures `[K]`; `(nlay+1,ncol)`"
    t_lev::FTA2D
    "Surface temperatures `[K]`; `(ncol)`"
    t_sfc::FTA1D
    "Number of molecules per cm^2 of dry air `(nlay, ncol)`"
    col_dry::FTA2D
    "volume mixing ratio of all relevant gases"
    vmr::VMR
    "effective radius of cloud liquid particles"
    cld_r_eff_liq::CLDP
    "effective radius of cloud ice particles"
    cld_r_eff_ice::CLDP
    "cloud water path"
    cld_path_liq::CLDP
    "cloud ice path"
    cld_path_ice::CLDP
    "cloud mask"
    cld_mask::CLDM
    "ice roughness, 1 = none, 2 = medium, 3 = rough"
    ice_rgh::I
    "Number of layers."
    nlay::I
    "Number of columns."
    ncol::I
    "Number of gases."
    ngas::I
end
Adapt.@adapt_structure AtmosphericState
end


text:
import IsApprox

@testset "FermiOp" begin
    @test FermiOp(0) === QuantumOps.FermiOps.I
    @test FermiOp(1) === QuantumOps.FermiOps.NumberOp
end

@testset "FermiTerm" begin
    t = FermiTerm("++-")
    @test reverse(t) != t
    @test reverse(reverse(t)) == t

    t1 = FermiTerm("+")
    t2 = FermiTerm("-")
    @test t1 * t2 == FermiTerm("N")

    t1 = FermiTerm("++")
    t2 = FermiTerm("--")
    @test t1 * t2 == FermiTerm("NN", -1)

    t1 = FermiTerm("+++")
    t2 = FermiTerm("---")
    @test t1 * t2 == FermiTerm("NNN", -1)

    t1 = FermiTerm("++++")
    t2 = FermiTerm("----")
    @test t1 * t2 == FermiTerm("NNNN", 1)

    t1 = FermiTerm("+E++")
    t2 = FermiTerm("-E--")
    @test t1 * t2 == FermiTerm("NENN", -1)

    @test (! IsApprox.ishermitian(FermiTerm("--++")))
    @test IsApprox.ishermitian(FermiTerm("IEINN"))
    @test (! IsApprox.ishermitian(FermiTerm("IEINN", 1 + 1.0im)))
end


text:
filename = replace(basename(@__FILE__), ".jl" => "")
filepath = joinpath(@__DIR__, "data/", filename, "input")
input = readlines(filepath)

diagnostic_report = transpose(parse.(Int, reduce(hcat, collect.(input)), base=2))
n, m = size(diagnostic_report)

# part 1

gamma_rate_bin = sum(diagnostic_report, dims=1) .> n/2
gamma_rate = sum(gamma_rate_bin[m-j]*2^j for j = 0:m-1)

epsilon_rate_bin = sum(diagnostic_report, dims=1) .< n/2
epsilon_rate = sum(epsilon_rate_bin[m-j]*2^j for j = 0:m-1)

part_one_answer = gamma_rate * epsilon_rate
println(part_one_answer)

# part 2

global oxygen_generator_rating, co2_scrubber_rating
let oxygen_generator_rating_bin = copy(diagnostic_report), co2_scrubber_rating_bin = copy(diagnostic_report)
    for j = 1:m
        count_one = sum(oxygen_generator_rating_bin[:, j])
        count_zero = sum(1 .- oxygen_generator_rating_bin[:, j])
        most_common = count_one .>= count_zero ? 1 : 0
        mask = oxygen_generator_rating_bin[:, j] .== most_common
        oxygen_generator_rating_bin = oxygen_generator_rating_bin[mask, :]
        if size(oxygen_generator_rating_bin, 1) == 1
            global oxygen_generator_rating = sum(oxygen_generator_rating_bin[m-j]*2^j for j = 0:m-1)
            break
        end
    end
    for j = 1:m
        count_one = sum(co2_scrubber_rating_bin[:, j])
        count_zero = sum(1 .- co2_scrubber_rating_bin[:, j])
        least_common = count_one .< count_zero ? 1 : 0
        mask = co2_scrubber_rating_bin[:, j] .== least_common
        co2_scrubber_rating_bin = co2_scrubber_rating_bin[mask, :]
        if size(co2_scrubber_rating_bin, 1) == 1
            global co2_scrubber_rating = sum(co2_scrubber_rating_bin[m-j]*2^j for j = 0:m-1)
            break
        end
    end
    part_two_answer = oxygen_generator_rating * co2_scrubber_rating
    println(part_two_answer)
end


text:
using DemoPackageNML
using Test

@testset "DemoPackageNML.jl" begin
    # Write your tests here.
end


text:
###################################################################################
# Antony Schutz 2015, ANR - POLCA - 2016
###################################################################################
# REF
###################################################################################
# [0] PAINTER
# Schutz, A., Ferrari, A., Mary, D., Soulez, F., Thiébaut, E., Vannier, M.
# Large scale 3D image reconstruction in optical interferometry
# EUSIPCO 2015
#
# [1] PAINTER
# Schutz, A., Ferrari, A., Mary, D., Soulez, F., Thiébaut, E., Vannier, M.
# PAINTER: a spatiospectral image reconstruction algorithm for optical interferometry
# J. Opt. Soc. Am. A, Vol. 31, N. 11, pp 2334--2345, 2014.
#
# [2] ADMM
# S. Boyd, N. Parikh, E. Chu, B. Peleato and J. Eckstein.
# “Distributed optimization and statistical learning
# via the alternating direction method of multipliers,”
# Found. Trends Mach. Learn., 3(1):1–122, January 2011.
###################################################################################
# # Operation involved in PAINTER
# ---------------------------------------------------------------------------------
# # # # # # # # # # #            W A V E L E T S
# ---------------------------------------------------------------------------------
# Object Estimation - Orthognal matrix - wavelet
# ---------------------------------------------------------------------------------
# function estimx_par{Tw<:WT.OrthoWaveletClass}(x::SharedArray{Float64,3},Fx::SharedArray{Complex{Float64},2},
function estimx_par{Tw<:WT.OrthoWaveletClass}( #x::Array{Float64,3},Fx::Array{Complex{Float64},2},
    rho_y::Float64,rho_spat::Float64,rho_spec::Float64,rho_ps::Float64,eta::Float64,
    yc::Array{Complex{Float64},2},z::Array{Float64,4},v::Array{Float64,3},w::Array{Float64,3},
    tau_xc::Array{Complex{Float64},2},tau_s::Array{Float64,4},tau_v::Array{Float64,3},tau_w::Array{Float64,3},
    nb::Int,nw::Int,nx::Int,NWvlt::Int,plan::Array{Any,1},Wvlt::Array{Tw,1},M::Array{Any,1},paral::Bool)
# Estimate the constrained, regularized 3D images from complexe visibilities
# step IV of PAINTER [0,1]
#
# x: images [(nx*nx)*nw], Fx is nb*nw is Non uniform Fourier transform of 3D images
# rho_x and tau_x: admm parameters and Lagrange mutlipliers
# yc: estimated complexe visibilities
# plan is the non uniform fft plan
# Wvlt is the list of used wavelets basis
# M is the used inverse matrix pre computed
    maty = (rho_y * yc) - tau_xc
    matz = (rho_spat * z) - tau_s
    matv = (rho_spec * v) - tau_v
    matw = (rho_ps * w) - tau_w
    Reg  = matv + matw
    matv = 0
    matw = 0
# parallel sum of wavelet basis
    wvd = SharedArray( Float64, (nx,nx,nw,NWvlt))
    @sync @parallel for ind in 1:nw*NWvlt
        n,b = ind2sub((nw,NWvlt),ind)
        wvd[:,:,n,b] = idwt(matz[:, :, n, b] , wavelet(Wvlt[b]) )
    end
    matz = 0
    wvd = sum(wvd,4)
# parallel image reconstruction
    x = SharedArray( Float64, (nx,nx,nw))
    Fx = SharedArray( Complex{Float64}, (nb, nw))
    @sync @parallel for n in 1:nw
        xtmp = (nfft_adjoint(plan[n], maty[:, n]) / nx) + Reg[:, :, n] + wvd[:, :, n] #+ wvd[n]
        xfst = nfft_adjoint(plan[n], M[n] * (nfft(plan[n], xtmp) / nx)) / nx
        xtmp = (xtmp - (rho_y / eta) * xfst) / eta
        x[:,:,n] = real(xtmp)
        Fx[:,n]  = nfft(plan[n], xtmp) / nx
    end
    return x, Fx
end
###################################################################################
# Regularization - constraint
# Section 4.B of PAINTER
###################################################################################
# Proximal operator for V2 and Cardano's formula
# ---------------------------------------------------------------------------------
# estimate the complexe visibilities from the squared modulus visibilities
# proximal operator
# P matrix of observed V2 [Nb*Nwvl]
# W matrix of observed V2 error variance [Nb*Nwvl]
# rho_y: Admm parameter (scalar>0)
# alpha: relative weight compared to phases difference: (scalar>0)
# (Global Cost) = alpha * (cost of V2) + Beta * (Cost of Phases diffrence)
# y_v2: auxiliary variable (last estimate of complexe visibility)
# nb is the number of base
# nw is the number of wavelength
#
# Equation 58-59 of PAINTER
#
function proxv2!(y_v2::Array,P::Array,W::Array,rho_y::Real,alpha::Real,nb::Int,nw::Int)
    mod_y = convert(SharedArray,abs(y_v2))
    ang_y = angle(y_v2)
    tmp1 = W.*rho_y /(4 * alpha)
    tmp2 = alpha ./ W
    @sync @parallel for z in 1:nb*nw
        m,n = ind2sub((nb,nw),z)
        sol = max(0., paintercubicroots( tmp1[m,n] - P[m,n], tmp1[m,n] .* mod_y[m,n]))
        cst = tmp2[m,n] .* (P[m,n] - sol.^2).^2 + .5 .* rho_y * (sol - mod_y[m,n] ).^2
        (a,b) = findmin(cst)
        mod_y[m,n] = sol[b]
    end

    y_v2[:] = copy(mod_y) .* exp(im .* ang_y)
end
# ---------------------------------------------------------------------------------
# ----- Cardano's formula
function paintercubicroots(c::Real,d::Real)
#-----   special  case  in  Painter -----
# a = 1.
# b = 0.
# d = -d
#-----  -----  -----  -----  -----  -----
# finds real valued roots of cubic
# Arguments:
#     a, b, c, d - coeffecients of cubic defined as
#                  ax^3 + bx^2 + cx + d = 0
# Returns:
# root   - an array of 1 or 3 real valued roots
# Reference:  mathworld.wolfram.com/CubicFormula.html
# Code follows Cardano's formula
# Copyright (c) 2008 Peter Kovesi
# School of Computer Science & Software Engineering
# The University of Western Australia
# pk at csse uwa edu au
# http://www.csse.uwa.edu.au/
# realcuberoot - computes real-valued cube root
function realcuberoot(x::Real)
    sign(x) .* abs(x).^(1 / 3)
end
    # Divide through by a to simplify things
    q = c / 3.
    r = d / 2.

    discriminant = q^3 + r^2
    if discriminant >= 0        # We have 1 real root and 2 imaginary
        s = realcuberoot(r + sqrt(discriminant))
        t = realcuberoot(r - sqrt(discriminant))
        root = s + t     # Just calculate the real root
    else                        # We have 3 real roots
        # In this case (r + sqrt(discriminate)) is complex so the following
        # code constructs the cube root of this complex quantity
        rho = sqrt(r^2 - discriminant)
        cubeRootrho = realcuberoot(rho)    # Cube root of complex magnitude
        thetaOn3 = acos(r / rho) / 3       # Complex angle/3
        crRhoCosThetaOn3 = cubeRootrho * cos(thetaOn3)
        crRhoSinThetaOn3 = cubeRootrho * sin(thetaOn3)
        root = zeros(3)
        root[1] = 2 * crRhoCosThetaOn3
        root[2] = -crRhoCosThetaOn3 - sqrt(3) * crRhoSinThetaOn3
        root[3] = -crRhoCosThetaOn3 + sqrt(3) * crRhoSinThetaOn3
    end
  return root
end
# ---------------------------------------------------------------------------------
# Proximal operator for phases difference (Section 5.2 PAINTER)
# for parallel
function proxphase(y_phi::Matrix,Xi::Vector,K::Vector,rho_y::Real,beta::Real
                  ,nb::Int,nw::Int,H::SparseMatrixCSC,pathoptpkpt::ASCIIString)
# estimate the complexe visibilities from the phases differences
# y_phi: vector of auxiliary variable \tilde{y} (last estimate of complexe visibility) is upadted
# Xi vector of observed Phases Difference [ {Nwvl*(Nb-1)*(Nb-2)/2 + Nb*(Nwvl-1)} *{1}]
# K the variance of phases must be corrected using Von Mises Ditribution (function EstKapVM)
# rho_y: Admm parameter (scalar>0)
# beta: relative weight compared to phases difference: (scalar>0)
# (Global Cost) = alpha * (cost of V2) + Beta * (Cost of Phases diffrence)
# nb is the number of base
# nw is the number of wavelength
# H: Phases difference to Phases matrix (function Ph2PhDiff)
# OPTOPT: structure of OptimPack vmlm option, see optimpack

    y_t = vec(y_phi)
    gam_t = abs(y_t)
    phi_t = angle(y_t)
    phi_0 = angle(y_t)

    function cost!{T<:Real}(x_phi::Array{T,1}, g_phi::Array{T,1})
        return costgradphi!(x_phi, g_phi, gam_t, phi_t, y_t, Xi, K, beta, rho_y, H)
    end

    if isempty(pathoptpkpt) # For travis
        phi = OptimPack.vmlm(cost!, phi_0, 100, verb = false
                            , grtol = 1e-3, gatol = 0, maxeval = 1000
                            , maxiter = 1000, stpmin = 1e-20, stpmax = e+20
                            , scaling = OptimPack.SCALING_OREN_SPEDICATO
                            , lnsrch = OptimPack.MoreThuenteLineSearch(ftol = 1e-8, gtol = 0.95))
    else
        include(pathoptpkpt)
        phi = OptimPack.vmlm(cost!, phi_0, memsize, verb = vt
                            , grtol = grt, gatol = gat, maxeval = mxvl
                            , maxiter = mxtr, stpmin = stpmn, stpmax = stpmx
                            , scaling = scl, lnsrch = ls)
    end

    if phi!=nothing
        Ek = phi_t - phi
        gam = max(0.0, gam_t .* cos(Ek))
        y_phi[:] = reshape(gam .* exp(im * phi), nb, nw)
    else
        y_phi[:] = reshape(gam_t .* exp(im * phi_t), nb, nw)
    end
end

function costgradphi!(x_phi::Vector,g_phi::Vector,gam_t::Vector,phi_t::Vector,y_t::Vector,Xi::Vector,K::Vector,beta::Real,rho_y::Real,H::SparseMatrixCSC)
# cost to minimize for phase difference minimization
# x_phi: Phases vector to estimate
# g_phi: gradiant of phase estimation cost function
# gam_t modulus of y_tilde
# phi_t phases of y_tilde
# Xi vector of observed Phases Difference [ {Nwvl*(Nb-1)*(Nb-2)/2 + Nb*(Nwvl-1)} *{1}]
# K the Von mises variance
# rho_y: Admm parameter (scalar>0)
# beta: relative weight compared to phases difference: (scalar>0)
# (Global Cost) = alpha * (cost of V2) + Beta * (Cost of Phases diffrence)
# H: Phases difference to Phases matrix (function Ph2PhDiff)
    Ek = phi_t - x_phi
    gam = max(0.0, gam_t .* cos(Ek))
    dphi = H * x_phi - Xi;
    yest = gam .* exp(im * x_phi)
    w1 = -sum(K .* cos(dphi))
    w2 = sum(abs(yest - y_t).^2)
    g_phi[:]= beta .* H' * (sin(dphi) ./ K) - rho_y .* gam .* gam_t .* sin(Ek)
    f = ((beta * w1) + (rho_y * w2)) / 2
    return f
end

###################################################################################
# MAIN ADMM LOOP
###################################################################################
# method 1, the algorithm is already initialized, structures are created and filed
# the admm can work with this informations
# function painteradmm(PDATA::PAINTER_Data,OIDATA::PAINTER_Input,OPTOPT::OptOptions,nbitermax::Int,aff::Bool)
#     painteradmm(PDATA,OIDATA,OPTOPT,nbitermax,aff)
# end
# for information about parameters read PAINTER [1] and admm [2]
function painteradmm(PDATA::PAINTER_Data,OIDATA::PAINTER_Input,nbitermax::Int,aff::Bool)
    const nx = OIDATA.nx
    const nb = OIDATA.nb
    const nw = OIDATA.nw
    const wvl = OIDATA.wvl
    const lambda_spat = OIDATA.lambda_spat
    const lambda_spec = OIDATA.lambda_spec
    const lambda_L1 = OIDATA.lambda_L1
    const epsilon = OIDATA.epsilon
    const rho_y = OIDATA.rho_y
    const rho_spat = OIDATA.rho_spat
    const rho_spec = OIDATA.rho_spec
    const rho_ps = OIDATA.rho_ps
    const alpha = OIDATA.alpha
    const beta = OIDATA.beta
    const eps1 = OIDATA.eps1
    const eps2 = OIDATA.eps2
    const Wvlt = OIDATA.Wvlt
    const mask3D = OIDATA.mask3D
    const P = OIDATA.P
    const W = OIDATA.W
    const Xi = OIDATA.Xi
    const K = OIDATA.K
    const paral = OIDATA.paral
    const eta = PDATA.eta
    const plan = PDATA.plan
    const F3D = PDATA.F3D
    const M = PDATA.M
    const NWvlt = length(Wvlt)
    const H = PDATA.H
    const baseNb = OIDATA.baseNb
    const Nt3indep = length(OIDATA.baseNb)
    const pathoptpkpt = OIDATA.pathoptpkpt
    yphidict = SharedArray( Complex128, (nb,nw) )
    Spcdct = convert( SharedArray, zeros( nx, nx, nw))
    vHt = convert( SharedArray, zeros(nx, nx, nw))
    Hx = convert( SharedArray, zeros( nx, nx, nw, NWvlt))

# ----------------------------------
    if !isempty(pathoptpkpt)
        include(pathoptpkpt)
        println(" ")
        println("------------------------ ")
        println("VMLM will run with file: ")
        println(pathoptpkpt)
        println("------------------------ ")
        println("memsize: ", memsize )
        println("verb: ", vt )
        println("grtol: ", grt )
        println("gatol: ", gat )
        println("maxeval: ", mxvl )
        println("maxiter: ", mxtr )
        println("stpmin: ", stpmn )
        println("stpmax: ", stpmx )
        println("scaling: ", scl )
        println("lnsrch: ", ls )
    end
# ----------------------------------
    println("")
    println("-----------------------------------------")
    println("| time |   primal   |    dual    |  It  |")
    println("-----------------------------------------")
    loop = true
    TiMe = zeros(nbitermax)
    while loop
        tic()
        PDATA.ind += 1

# update of yc from V2
        PDATA.y_v2 = PDATA.yc + PDATA.tau_pwc / rho_y
        proxv2!(PDATA.y_v2, P, W, rho_y, alpha, nb, nw)

# update of yc from phases difference
        y_phidat = PDATA.yc + PDATA.tau_xic / rho_y
        @sync @parallel for n in 1:Nt3indep
            yphidict[baseNb[n],:] = proxphase(y_phidat[baseNb[n],:], Xi[n], K[n]
                     , rho_y, beta, length( baseNb[n] ), nw, H[n], pathoptpkpt )
        end
        PDATA.y_phi = copy(yphidict)
        y_phidat = 0
# Consensus
        y_tmp = copy(PDATA.yc)
        PDATA.yc = ( PDATA.y_v2 + PDATA.y_phi + PDATA.Fx + (PDATA.tau_xc - (PDATA.tau_pwc + PDATA.tau_xic)) ./ rho_y ) ./ 3

# Object estimation
        x_tmp = copy(PDATA.x)
        PDATA.x,PDATA.Fx = estimx_par(rho_y, rho_spat, rho_spec, rho_ps,eta ,PDATA.yc, PDATA.z, PDATA.v, PDATA.w,
                                      PDATA.tau_xc, PDATA.tau_s, PDATA.tau_v, PDATA.tau_w, nb, nw, nx, NWvlt,
                                      plan, Wvlt, M, paral)
# update of auxiliary variables
        if rho_spat >0
            tmpx = copy(PDATA.x)
            @sync @parallel for ind in 1:nw*NWvlt
                n,b = ind2sub((nw,NWvlt),ind)
                Hx[:, :, n, b] = dwt(tmpx[:, :, n], wavelet(Wvlt[b]))
            end
            tmpx = 0
# update of z
            PDATA.z = Hx + (PDATA.tau_s / rho_spat)
            PDATA.z = max(1 - ((lambda_spat / rho_spat) ./ abs(PDATA.z)), 0.) .* PDATA.z
        end
# update of v
        if rho_spec >0
            tmp = permutedims(rho_spec * PDATA.r - PDATA.tau_r, [3, 1, 2])
            @sync @parallel for ind in 1:nx*nx
                m,n = ind2sub((nx,nx),ind)
                Spcdct[m, n, :]= idct(tmp[:, m, n] )
            end
            tmp = 0

            PDATA.v = ( Spcdct + (rho_spec * PDATA.x) + PDATA.tau_v) / (2 * rho_spec)
            vecv = permutedims(PDATA.v, [3, 1, 2])
            @sync @parallel for ind in 1:nx*nx
                m,n = ind2sub((nx,nx),ind)
                vHt[m, n, :] = dct(vecv[:, m, n] )
            end
            vecv = 0
# update of r
            PDATA.r = vHt + PDATA.tau_r / rho_spec
            PDATA.r = max(1 - (lambda_spec / rho_spec) ./ abs(PDATA.r), 0) .* PDATA.r
        end
# update of w
        if rho_ps>0
            u = PDATA.x + PDATA.tau_w ./ rho_ps
            PDATA.w = max(max(0.0, u) .* mask3D - lambda_L1, 0)
            PDATA.tau_w = PDATA.tau_w + rho_ps * (PDATA.x - PDATA.w)
        end
# update of Lagrange multipliers
        PDATA.tau_pwc = PDATA.tau_pwc + rho_y * (PDATA.yc - PDATA.y_v2)
        PDATA.tau_xic = PDATA.tau_xic + rho_y * (PDATA.yc - PDATA.y_phi)
        PDATA.tau_xc = PDATA.tau_xc + rho_y * (PDATA.Fx - PDATA.yc)
        if rho_spat >0
            PDATA.tau_s = PDATA.tau_s + rho_spat * (Hx - PDATA.z)
        end
        if rho_spec >0
            PDATA.tau_v = PDATA.tau_v + rho_spec * (PDATA.x - PDATA.v)
            PDATA.tau_r = PDATA.tau_r + rho_spec * (vHt- PDATA.r)
        end
# stopping criteria
        n1 = norm(vec(PDATA.x -x_tmp))
        n2 = norm(vec(PDATA.yc-y_tmp))
        push!(PDATA.crit1, n1)
        push!(PDATA.crit2, n2)
# Plot and verbose
        if aff&&(PDATA.ind - 1)==(PDATA.count * PDATA.CountPlot)
            OIDATA.PlotFct(PDATA,OIDATA)
            PDATA.count += 1
        end
        if (PDATA.ind >= nbitermax)||( (n1 < eps1)&&(n2 < eps2) )
            loop = false
        end
        @printf("| %02.02f | %02.04e | %02.04e | %04d |\n",toq(), PDATA.crit1[PDATA.ind], PDATA.crit2[PDATA.ind], PDATA.ind)
    end
    return PDATA
end
###################################################################################
# PAINTER MAIN FUNCTION
###################################################################################
function painter(;Folder = "", nbitermax = 1000, nx = 64, lambda_spat = 1/nx^2,
                 lambda_spec = 1/100, lambda_L1 = 0, epsilon = 1e-6,
                 rho_y = 1, rho_spat = 1, rho_spec = 1, rho_ps = 1, alpha = 1,
                 dptype = "all", dpprm = 0, pathoptpkpt = joinpath(dirname(@__FILE__),"optpckpt.jl"),
                 Wvlt  = [WT.db1, WT.db2, WT.db3, WT.db4, WT.db5, WT.db6, WT.db7, WT.db8, WT.haar],
                 beta = 1, eps1 = 1e-6, eps2 = 1e-6, FOV = 4e-2, mask3D = [], xinit3D = [], indfile = [], indwvl = [],
                 PlotFct = painterplotfct, aff = false, CountPlot = 10, admm = true, paral = true)
# Check if mandatory package are installed
    checkPack()
# PAINTER Data Type Creation
    OIDATA = painterinputinit()
    PDATA  = painterdatainit()
# PAINTER User parameter validation
    OIDATA = painterinit(OIDATA, Folder, nx, lambda_spat, lambda_spec, lambda_L1,
                         epsilon, rho_y, rho_spat, rho_spec, rho_ps, alpha, beta,
                         eps1, eps2, FOV, mask3D, xinit3D, Wvlt, paral,
                         dptype, dpprm, PlotFct,pathoptpkpt)
# OIFITS-FITS Data Read
    println("")
    OIDATA = readoifits(OIDATA, indfile, indwvl)
## PAINTER Matrices creation, Array Initialization
    println("")
    PDATA,OIDATA  = painterarrayinit(PDATA, OIDATA)
# Check, Create PAINTER object and mask initialisation from data or fits
    OIDATA.mask3D   = checkmask(OIDATA.mask3D, OIDATA.nx, OIDATA.nw)
    PDATA.CountPlot = CountPlot

# Main Loop ADMM
    if admm
        PDATA = painteradmm(PDATA, OIDATA, nbitermax, aff)
    end
    return OIDATA, PDATA
end
#########################################
function painter(OIDATA::PAINTER_Input,PDATA::PAINTER_Data,nbitermax::Int,aff::Bool;PlotFct = painterplotfct)
    OIDATA.PlotFct = PlotFct
    PDATA = painteradmm(PDATA, OIDATA, nbitermax, aff)
    return OIDATA, PDATA
end


text:
if !@isdefined ExoplanetsSysSim
    using ExoplanetsSysSim
end

import Compat: UTF8String, ASCIIString

## simulation_parameters
function setup_sim_param_model(args::Vector{String} = Array{String}(undef, 0) )   # allow this to take a list of parameter (e.g., from command line)
  sim_param = SimParam()
  # How many tatrges to generate
  #add_param_fixed(sim_param,"num_targets_sim_pass_one",150061)                      # Note this is used for the number of stars in the simulations, not necessarily related to number of Kepler targets
  add_param_fixed(sim_param,"num_targets_sim_pass_one",78005)                      # Note this is used for the number of stars in the simulations, not necessarily related to number of Kepler targets
  add_param_fixed(sim_param,"num_kepler_targets",78005)                            # Note this is used for the number of Kepler targets for the observational catalog

  # For generating target star properties
  add_param_fixed(sim_param,"generate_kepler_target",generate_kepler_target_from_table)
  add_param_fixed(sim_param,"star_table_setup",setup_star_table_christiansen)
  #add_param_fixed(sim_param,"stellar_catalog","q1_q17_dr25_stellar.jld") #"q1_q17_dr25_stellar.csv"
  add_param_fixed(sim_param,"stellar_catalog","q1q17_dr25_gaia_fgk.csv") #"q1q17_dr25_gaia_fgk.csv"
  #add_param_fixed(sim_param,"generate_kepler_target",ExoplanetsSysSim.generate_kepler_target_simple)  # An alternative that alternative can be used for testing if above breaks
  add_param_fixed(sim_param,"window_function","DR25topwinfuncs.jld2")

  # For generating planetary system properties
  add_param_fixed(sim_param,"generate_planetary_system", generate_planetary_system_clustered)

  add_param_fixed(sim_param,"generate_num_clusters",generate_num_clusters_poisson)
  add_param_fixed(sim_param,"generate_num_planets_in_cluster",generate_num_planets_in_cluster_poisson)
  add_param_active(sim_param,"log_rate_clusters",log(3.0))
  add_param_fixed(sim_param,"max_clusters_in_sys",10)
  add_param_active(sim_param,"log_rate_planets_per_cluster",log(3.0))
  add_param_fixed(sim_param,"max_planets_in_cluster",10)

  # generate_num_planets_in_cluster currently calls: generate_periods_power_law
  add_param_fixed(sim_param,"generate_sizes",ExoplanetsSysSim.generate_sizes_broken_power_law) # To choose the way we draw planetary radii; if "generate_sizes_power_law", then takes "power_law_r"; if "generate_sizes_broken_power_law", then takes "power_law_r1", "power_law_r2", and "break_radius"
  add_param_active(sim_param,"power_law_P",0.5)
  add_param_fixed(sim_param,"power_law_r",-2.5)
  add_param_active(sim_param,"power_law_r1",-1.5)
  add_param_active(sim_param,"power_law_r2",-6.0)
  add_param_fixed(sim_param,"min_period",3.0)
  add_param_fixed(sim_param,"max_period",300.0)
  add_param_fixed(sim_param,"min_radius",0.5*ExoplanetsSysSim.earth_radius)
  add_param_fixed(sim_param,"max_radius",10.0*ExoplanetsSysSim.earth_radius)
  add_param_fixed(sim_param,"break_radius",3.0*ExoplanetsSysSim.earth_radius)

  # generate_num_planets_in_cluster currently use these for the Inclination distribution
  add_param_active(sim_param,"sigma_incl",1.5) # degrees; 0 = coplanar w/ generate_kepler_target_simple; ignored by generate_planetary_system_uncorrelated_incl
  add_param_active(sim_param,"sigma_incl_near_mmr",1.5)

  add_param_fixed(sim_param,"max_incl_sys",0.0) #degrees; gives system inclinations from "max_incl_sys" (deg) to 90 (deg), so set to 0 for isotropic distribution of system inclinations; NOTE: make sure the difference between this and 90 (deg) is at least greater than "sigma_incl" and "sigma_incl_near_mmr"!

  # generate_num_planets_in_cluster currently use these for the Eccentricity distribution
  add_param_fixed(sim_param,"generate_e_omega",ExoplanetsSysSim.generate_e_omega_rayleigh)
  add_param_active(sim_param,"sigma_hk",0.05)
  #add_param_fixed(sim_param,"sigma_hk_one",0.1)
  #add_param_fixed(sim_param,"sigma_hk_multi",0.03)

  # generate_num_planets_in_cluster currently use these for the Stability tests
  add_param_fixed(sim_param,"num_mutual_hill_radii",8.0) #10.0
  add_param_fixed(sim_param,"generate_planet_mass_from_radius",generate_planet_mass_from_radius_Ning2018_table) # "ExoplanetsSysSim.generate_planet_mass_from_radius_powerlaw" or "generate_planet_mass_from_radius_Ning2018" or "generate_planet_mass_from_radius_Ning2018_table"
  add_param_fixed(sim_param,"mr_power_index",2.0)
  add_param_fixed(sim_param,"mr_const",1.0)
  add_param_fixed(sim_param,"mr_max_mass",1e3*ExoplanetsSysSim.earth_mass)
  add_param_active(sim_param,"sigma_log_radius_in_cluster",0.25)
  add_param_active(sim_param,"sigma_logperiod_per_pl_in_cluster",0.15)

  # Functions to calculate observables from physical system properties
  add_param_fixed(sim_param,"calc_target_obs_single_obs",ExoplanetsSysSim.calc_target_obs_single_obs)
  add_param_fixed(sim_param,"max_tranets_in_sys",8)     # SysSim ignores some planets in any systems with more than this many transiting planets to avoid wasting time on unphysical parameter values
  add_param_fixed(sim_param,"transit_noise_model",ExoplanetsSysSim.transit_noise_model_fixed_noise)
  #add_param_fixed(sim_param,"transit_noise_model",ExoplanetsSysSim.transit_noise_model_diagonal)
  # add_param_fixed(sim_param,"rng_seed",1234)   # If you want to be able to reproduce simulations
  add_param_fixed(sim_param,"read_target_obs",ExoplanetsSysSim.simulated_read_kepler_observations)  # Read Kepler observations to compare to from disk

  # Read any customizations, so its easy to keep them separate from the essentials
  #=if isfile("param_custom.jl")
     include("param_custom.jl")
  end
  =#
  if @isdefined add_param_custom
     sim_param = add_param_custom(sim_param)
  end
  return sim_param
end

function test_setup_sim_param()
  setup_sim_param_model()
end

function write_model_params(f, sim_param::SimParam)
    #This function writes all the model parameters to a file f as a header
    println(f, "# num_targets_sim_pass_one: ", get_int(sim_param,"num_targets_sim_pass_one"))
    println(f, "# max_incl_sys: ", get_real(sim_param,"max_incl_sys"))
    println(f, "# log_rate_clusters: ", get_real(sim_param,"log_rate_clusters"))
    println(f, "# max_clusters_in_sys: ", get_int(sim_param,"max_clusters_in_sys"))
    println(f, "# log_rate_planets_per_cluster: ", get_real(sim_param,"log_rate_planets_per_cluster"))
    println(f, "# max_planets_in_clusters: ", get_int(sim_param,"max_planets_in_cluster"))
    println(f, "# power_law_P: ", get_real(sim_param,"power_law_P"))
    println(f, "# min_period: ", get_real(sim_param,"min_period"))
    println(f, "# max_period: ", get_real(sim_param,"max_period"))

    if string(get_function(sim_param,"generate_sizes")) == "ExoplanetsSysSim.generate_sizes_power_law"
        println(f, "# power_law_r: ", get_real(sim_param,"power_law_r"))
    elseif string(get_function(sim_param,"generate_sizes")) == "ExoplanetsSysSim.generate_sizes_broken_power_law"
        println(f, "# power_law_r1: ", get_real(sim_param,"power_law_r1"))
        println(f, "# power_law_r2: ", get_real(sim_param,"power_law_r2"))
        println(f, "# break_radius (R_earth): ", get_real(sim_param,"break_radius")/ExoplanetsSysSim.earth_radius)
    end

    println(f, "# min_radius (R_earth): ", get_real(sim_param,"min_radius")/ExoplanetsSysSim.earth_radius)
    println(f, "# max_radius (R_earth): ", get_real(sim_param,"max_radius")/ExoplanetsSysSim.earth_radius)
    println(f, "# sigma_incl: ", get_real(sim_param,"sigma_incl"))
    println(f, "# sigma_incl_near_mmr: ", get_real(sim_param,"sigma_incl_near_mmr"))
    println(f, "# sigma_hk: ", get_real(sim_param,"sigma_hk"))
    println(f, "# num_mutual_hill_radii: ", get_real(sim_param,"num_mutual_hill_radii"))

    if string(get_function(sim_param,"generate_planet_mass_from_radius")) == "ExoplanetsSysSim.generate_planet_mass_from_radius_powerlaw"
        println(f, "# mr_power_index: ", get_real(sim_param,"mr_power_index"))
        println(f, "# mr_max_mass (M_earth): ", get_real(sim_param,"mr_max_mass")/ExoplanetsSysSim.earth_mass)
    elseif string(get_function(sim_param,"generate_planet_mass_from_radius")) == "generate_planet_mass_from_radius_Ning2018"
        println(f, "# mr_model: Ning2018")
    elseif string(get_function(sim_param,"generate_planet_mass_from_radius")) == "generate_planet_mass_from_radius_Ning2018_table"
        println(f, "# mr_model: Ning2018_table")
    end

    println(f, "# sigma_log_radius_in_cluster: ", get_real(sim_param,"sigma_log_radius_in_cluster"))
    println(f, "# sigma_logperiod_per_pl_in_cluster: ", get_real(sim_param,"sigma_logperiod_per_pl_in_cluster"))
    println(f, "#")
end


text:
function MagneticFieldDetector(arg0::jdouble, arg1::GeoMagneticFieldFactory_FieldModel, arg2::OneAxisEllipsoid)
    return MagneticFieldDetector((jdouble, GeoMagneticFieldFactory_FieldModel, OneAxisEllipsoid), arg0, arg1, arg2)
end

function MagneticFieldDetector(arg0::jdouble, arg1::GeoMagneticFieldFactory_FieldModel, arg2::OneAxisEllipsoid, arg3::jboolean)
    return MagneticFieldDetector((jdouble, GeoMagneticFieldFactory_FieldModel, OneAxisEllipsoid, jboolean), arg0, arg1, arg2, arg3)
end

function MagneticFieldDetector(arg0::jdouble, arg1::jdouble, arg2::jdouble, arg3::GeoMagneticFieldFactory_FieldModel, arg4::OneAxisEllipsoid, arg5::jboolean)
    return MagneticFieldDetector((jdouble, jdouble, jdouble, GeoMagneticFieldFactory_FieldModel, OneAxisEllipsoid, jboolean), arg0, arg1, arg2, arg3, arg4, arg5)
end

function MagneticFieldDetector(arg0::jdouble, arg1::jdouble, arg2::jdouble, arg3::GeoMagneticFieldFactory_FieldModel, arg4::OneAxisEllipsoid, arg5::jboolean, arg6::DataContext)
    return MagneticFieldDetector((jdouble, jdouble, jdouble, GeoMagneticFieldFactory_FieldModel, OneAxisEllipsoid, jboolean, DataContext), arg0, arg1, arg2, arg3, arg4, arg5, arg6)
end

function g(obj::MagneticFieldDetector, arg0::SpacecraftState)
    return jcall(obj, "g", jdouble, (SpacecraftState,), arg0)
end

function init(obj::MagneticFieldDetector, arg0::SpacecraftState, arg1::AbsoluteDate)
    return jcall(obj, "init", void, (SpacecraftState, AbsoluteDate), arg0, arg1)
end