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

using Markdown
using InteractiveUtils

# ╔═╡ 6bf88b5c-3abc-4b03-90bb-a7eb8cb68577
using StochasticDiffEq, UnicodePlots, Statistics

# ╔═╡ 96c9b046-6e40-44b1-a19c-b7061f22470d
md"""## Eulerian Model

Starting from a step function we integrate a diffusion equation in one dimension, `z`, with closed boundary at the top and bottom.
"""

# ╔═╡ d78532b9-b17e-4a4a-bbab-e249bf852eba
begin
	## Eulerian

	function EulerianModel(nt=1)
	    N=20
	    dt=1e-4
	    dx=1.0/2/N
	    T=[zeros(N);ones(N)]
	    T0=deepcopy(T)
	    for tt in 1:nt
	        dTr=(circshift(T,-1)-T); dTr[end]=0;
	        dTl=(T-circshift(T,+1)); dTl[1]=0;
	        T.+=(dTr-dTl)*dt/dx/dx
	    end
	    return T,T0
	end

	function plot_EulerianModel(T,T0)
	    plt=lineplot(T0); lineplot!(plt,T)
	    return plt
	end

	T,T0=EulerianModel(1000);
	f=plot_EulerianModel(T,T0)

	"Eulerian Model Has Completed"
end

# ╔═╡ 41487611-deba-49fc-a5f7-f4f63f7935c6
md"""## Lagrangian Model

Starting from a step-like distribution of two particle groups over the `z=(0,1)` interval. Particle positions evolve according to a standard `Wiener` process. They bounce at the top and bottom boundaries. The charge of particles in one group is 0 at the start; 1 in the other group. Afterwards particles can interact with one another within each model layer.

_Credits: the presented model is a simplified version of the `aquacosm` model from Paparella & Vichi M (2020) Stirring, Mixing, Growing: Microscale Processes Change Larger Scale Phytoplankton Dynamics. doi: 10.3389/fmars.2020.00654_
"""

# ╔═╡ 32a692ab-4835-4943-8b46-8345ce881eb5
begin
	## Lagrangian

	function solve_paths(u₀)
	    dt=1e-3
	    f(u,p,t) = 0.0
	    g(u,p,t) = 0.1
	    #dt = 1//2^(4)
	    tspan = (0.0,1.0)
	    prob = SDEProblem(f,g,u₀,(0.0,1.0))
	    solve(prob,EM(),dt=dt)
	end

	function fold_tails(z)
	    while !isempty(findall( xor.(z.>1.0,z.<0.0) ))
	        z[findall(z.<0.0)].=-z[findall(z.<0.0)]
	        z[findall(z.>1.0)].=(2.0 .- z[findall(z.>1.0)])
	    end
	end

	## mix between cells that are at the same level
	function mix_neighbors(za,ca,zb,cb,p)
	    dz=0.1
	    for i0 in 1:10
	        z0=dz*(i0-1)
	        z1=dz*i0
	        ia=findall((za.>z0).*(za.<=z1))
	        ib=findall((zb.>z0).*(zb.<=z1))
	        tmp=mean([ca[ia];cb[ib]])
	        ca[ia].=(1-p)*ca[ia] .+ p*tmp
	        cb[ib].=(1-p)*cb[ib] .+ p*tmp
	    end
	end

	"""
		main_loop(;p=0.5,nt=5)
	```
	p=0.5 # fraction of mass exchanged with neighbors every time step
	nt=5 # number of time steps
	```
	"""
	function main_loop(;p=0.5,nt=5)
		for tt in 1:nt
			u₀=u₀a
			sol=solve_paths(u₀)
			za=sol(0:0.01:1)[:,:]
			fold_tails(za)
			u₀a[:]=za[:,end]

			u₀=u₀b
			sol=solve_paths(u₀)
			zb=sol(0:0.01:1)[:,:]
			fold_tails(zb)
			u₀b[:]=zb[:,end]

			mix_neighbors(u₀a,ca,u₀b,cb,p)
		end
		return "done with model run"
	end

	"Lagrangian model formulated"
end

# ╔═╡ c86cf1e9-d0b6-48e8-8a01-0dcc591afc85
begin
	## Plotting

	function plot_paths(z)
		np=size(z,1)
		plt=lineplot(z[1,:],ylim=(-0.1,1.1))
		np>1 ? lineplot!(plt,z[2,:]) : nothing
		plt
	end

	function plot_paths(sol::RODESolution)
		np=size(sol,1)
		plt=lineplot(sol(0:0.01:1)[1,:],ylim=(-0.1,1.1))
		np>1 ? lineplot!(plt,sol(0:0.01:1)[2,:]) : nothing
		plt
	end

	function gridded_stats(za,ca,zb,cb)
		out=zeros(10,2)
		dz=0.1
		t=size(za,2)
		for i0=1:10
			z0=0+dz*(i0-1)
			ia=findall( (za[:,t].>z0).*(za[:,t].<=z0+dz) );
			ib=findall( (zb[:,t].>z0).*(zb[:,t].<=z0+dz) );
			tmp=[ca[ia,t];cb[ib,t]]
			out[i0,1]=mean(tmp)
			out[i0,2]=std(tmp)
		end
		out
	end

	function plot_stats(st)
		plt=lineplot(st[:,1],ylim=(-0.1,1.1))
		lineplot!(plt,st[:,1].+st[:,2])
		lineplot!(plt,st[:,1].-st[:,2])
		plt
	end

	"Plotting functions"
end

# ╔═╡ 2195a7a4-30d6-47fe-92db-d7ebef67ef03
begin
	# initial conditions
	np=10000
	u₀a=0.5*rand(np)
	ca=zeros(np)
	u₀b=0.5 .+ 0.5*rand(np)
	cb=ones(np)
	"done with initialization"
end

# ╔═╡ fd2db517-2282-43da-a543-a1ef59a439c2
# Main model run
main_loop(p=0.1,nt=10)

# ╔═╡ 5ab5321a-15b9-11ec-3321-b7d7cab1ba35
begin
	# Another run of the dispersion model just for plotting

	u₀=u₀a
	sol=solve_paths(u₀)
	za=sol(0:0.01:1)[:,:]
	fold_tails(za)
	fa=plot_paths(za)

	u₀=u₀b
	sol=solve_paths(u₀)
	zb=sol(0:0.01:1)[:,:]
	fold_tails(zb)
	fb=plot_paths(zb)

	##

	st=gridded_stats(u₀a,ca,u₀b,cb)
	fs=plot_stats(st);

	## various

	#εf(z)C(1−C)

	"figures generated : f, fs, fa, fb"
end

# ╔═╡ aa7aa58b-e2e7-4b40-aeab-38137f8d4e2d
md"""## Summary / Plot Collection

From left to right:

- solution of the Eulerian model
- solution of the Lagrangian model (with irreversible mixing)
- pair or Lagrangian trajectories
- pair or Lagrangian trajectories
"""

# ╔═╡ 1df22a2f-c344-4d71-b44d-9a821c14e548
(f,fs,fa,fb)

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"

[compat]
StochasticDiffEq = "~6.37.1"
UnicodePlots = "~2.4.0"
"""

# ╔═╡ 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-beta3.0"
manifest_format = "2.0"

[[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 = ["IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"]
git-tree-sha1 = "d84c956c4c0548b4caf0e4e96cf5b6494b5b1529"
uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
version = "3.1.32"

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

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

[[deps.BitTwiddlingConvenienceFunctions]]
deps = ["Static"]
git-tree-sha1 = "652aab0fc0d6d4db4cc726425cadf700e9f473f1"
uuid = "62783981-4cbd-42fc-bca8-16325de8dc4b"
version = "0.1.0"

[[deps.CPUSummary]]
deps = ["Hwloc", "IfElse", "Static"]
git-tree-sha1 = "ed720e2622820bf584d4ad90e6fcb93d95170b44"
uuid = "2a0fbf3d-bb9c-48f3-b0a9-814d99fd7ab9"
version = "0.1.3"

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

[[deps.CloseOpenIntervals]]
deps = ["ArrayInterface", "Static"]
git-tree-sha1 = "ce9c0d07ed6e1a4fecd2df6ace144cbd29ba6f37"
uuid = "fb6a15b2-703c-40df-9091-08a04967cfa9"
version = "0.1.2"

[[deps.CommonSolve]]
git-tree-sha1 = "68a0743f578349ada8bc911a5cbd5a2ef6ed6d1f"
uuid = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
version = "0.2.0"

[[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 = "4866e381721b30fac8dda4c8cb1d9db45c8d2994"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.37.0"

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

[[deps.ConstructionBase]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "f74e9d5388b8620b4cee35d4c5a618dd4dc547f4"
uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
version = "1.3.0"

[[deps.Crayons]]
git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"

[[deps.DEDataArrays]]
deps = ["ArrayInterface", "DocStringExtensions", "LinearAlgebra", "RecursiveArrayTools", "SciMLBase", "StaticArrays"]
git-tree-sha1 = "31186e61936fbbccb41d809ad4338c9f7addf7ae"
uuid = "754358af-613d-5f8d-9788-280bf1605d4c"
version = "0.2.0"

[[deps.DataAPI]]
git-tree-sha1 = "bec2532f8adb82005476c141ec23e921fc20971b"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.8.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.DiffEqBase]]
deps = ["ArrayInterface", "ChainRulesCore", "DEDataArrays", "DataStructures", "Distributions", "DocStringExtensions", "FastBroadcast", "ForwardDiff", "FunctionWrappers", "IterativeSolvers", "LabelledArrays", "LinearAlgebra", "Logging", "MuladdMacro", "NonlinearSolve", "Parameters", "PreallocationTools", "Printf", "RecursiveArrayTools", "RecursiveFactorization", "Reexport", "Requires", "SciMLBase", "Setfield", "SparseArrays", "StaticArrays", "Statistics", "SuiteSparse", "ZygoteRules"]
git-tree-sha1 = "420ad175d5e420e2c55a0ed8a9c18556e6735f80"
uuid = "2b5f629d-d688-5b77-993f-72d75c75574e"
version = "6.73.2"

[[deps.DiffEqJump]]
deps = ["ArrayInterface", "Compat", "DataStructures", "DiffEqBase", "FunctionWrappers", "LightGraphs", "LinearAlgebra", "PoissonRandom", "Random", "RandomNumbers", "RecursiveArrayTools", "Reexport", "StaticArrays", "TreeViews", "UnPack"]
git-tree-sha1 = "99a65172d95e5ccc016c9be0542fa3858cf97a18"
uuid = "c894b116-72e5-5b58-be3c-e6d8d4ac2b12"
version = "7.3.0"

[[deps.DiffEqNoiseProcess]]
deps = ["DiffEqBase", "Distributions", "LinearAlgebra", "Optim", "PoissonRandom", "QuadGK", "Random", "Random123", "RandomNumbers", "RecipesBase", "RecursiveArrayTools", "Requires", "ResettableStacks", "SciMLBase", "StaticArrays", "Statistics"]
git-tree-sha1 = "d6839a44a268c69ef0ed927b22a6f43c8a4c2e73"
uuid = "77a26b50-5914-5dd7-bc55-306e6241c503"
version = "5.9.0"

[[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 = "3ed8fa7178a10d1cd0f1ca524f249ba6937490c0"
uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
version = "1.3.0"

[[deps.Distances]]
deps = ["LinearAlgebra", "Statistics", "StatsAPI"]
git-tree-sha1 = "9f46deb4d4ee4494ffb5a40a27a2aced67bdd838"
uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
version = "0.10.4"

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

[[deps.Distributions]]
deps = ["ChainRulesCore", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns"]
git-tree-sha1 = "f4efaa4b5157e0cdb8283ae0b5428bc9208436ed"
uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
version = "0.25.16"

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

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

[[deps.ExponentialUtilities]]
deps = ["ArrayInterface", "LinearAlgebra", "Printf", "Requires", "SparseArrays"]
git-tree-sha1 = "7a541ee92e2f8b16356ed6066d0c44b85984b780"
uuid = "d4d017d3-3776-5f7e-afef-a10c40355c18"
version = "1.9.0"

[[deps.FastBroadcast]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "26be48918640ce002f5833e8fc537b2ba7ed0234"
uuid = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
version = "0.1.8"

[[deps.FastClosures]]
git-tree-sha1 = "acebe244d53ee1b461970f8910c235b259e772ef"
uuid = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a"
version = "0.3.2"

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

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

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

[[deps.FunctionWrappers]]
git-tree-sha1 = "241552bc2209f0fa068b6415b1942cc0aa486bcc"
uuid = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e"
version = "1.1.2"

[[deps.Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[deps.HostCPUFeatures]]
deps = ["BitTwiddlingConvenienceFunctions", "IfElse", "Libdl", "Static"]
git-tree-sha1 = "3169c8b31863f9a409be1d17693751314241e3eb"
uuid = "3e5b6fbb-0976-4d2c-9146-d79de83f2fb0"
version = "0.1.4"

[[deps.Hwloc]]
deps = ["Hwloc_jll"]
git-tree-sha1 = "92d99146066c5c6888d5a3abc871e6a214388b91"
uuid = "0e44f5e4-bd66-52a0-8798-143a42290a1d"
version = "2.0.0"

[[deps.Hwloc_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "3395d4d4aeb3c9d31f5929d32760d8baeee88aaf"
uuid = "e33a78d0-f292-5ffc-b300-72abe9b543c8"
version = "2.5.0+0"

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

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

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

[[deps.IrrationalConstants]]
git-tree-sha1 = "f76424439413893a832026ca355fe273e93bce94"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.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.JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[deps.LabelledArrays]]
deps = ["ArrayInterface", "LinearAlgebra", "MacroTools", "StaticArrays"]
git-tree-sha1 = "bdde43e002847c34c206735b1cf860bc3abd35e7"
uuid = "2ee39098-c373-598a-b85f-a56591580800"
version = "1.6.4"

[[deps.LayoutPointers]]
deps = ["ArrayInterface", "LinearAlgebra", "ManualMemory", "SIMDTypes", "Static"]
git-tree-sha1 = "d2bda6aa0b03ce6f141a2dc73d0bcb7070131adc"
uuid = "10f19ff3-798f-405d-979b-55457f8fc047"
version = "0.1.3"

[[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.LineSearches]]
deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "Printf"]
git-tree-sha1 = "f27132e551e959b3667d8c93eae90973225032dd"
uuid = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
version = "7.1.1"

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

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

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

[[deps.LoopVectorization]]
deps = ["ArrayInterface", "CPUSummary", "CloseOpenIntervals", "DocStringExtensions", "HostCPUFeatures", "IfElse", "LayoutPointers", "LinearAlgebra", "OffsetArrays", "PolyesterWeave", "Requires", "SLEEFPirates", "Static", "ThreadingUtilities", "UnPack", "VectorizationBase"]
git-tree-sha1 = "d469fcf148475a74c221f14d42ee75da7ccb3b4e"
uuid = "bdcacae8-1622-11e9-2a5c-532679323890"
version = "0.12.73"

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

[[deps.ManualMemory]]
git-tree-sha1 = "9cb207b18148b2199db259adfa923b45593fe08e"
uuid = "d125e4d3-2237-4719-b19c-fa641b8a4667"
version = "0.1.6"

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

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

[[deps.Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

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

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

[[deps.MuladdMacro]]
git-tree-sha1 = "c6190f9a7fc5d9d5915ab29f2134421b12d24a68"
uuid = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
version = "0.2.2"

[[deps.NLSolversBase]]
deps = ["DiffResults", "Distributed", "FiniteDiff", "ForwardDiff"]
git-tree-sha1 = "144bab5b1443545bc4e791536c9f1eacb4eed06a"
uuid = "d41bc354-129a-5804-8e4c-c37616107c6c"
version = "7.8.1"

[[deps.NLsolve]]
deps = ["Distances", "LineSearches", "LinearAlgebra", "NLSolversBase", "Printf", "Reexport"]
git-tree-sha1 = "019f12e9a1a7880459d0173c182e6a99365d7ac1"
uuid = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
version = "4.5.1"

[[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.NonlinearSolve]]
deps = ["ArrayInterface", "FiniteDiff", "ForwardDiff", "IterativeSolvers", "LinearAlgebra", "RecursiveArrayTools", "RecursiveFactorization", "Reexport", "SciMLBase", "Setfield", "StaticArrays", "UnPack"]
git-tree-sha1 = "35585534c0c79c161241f2e65e759a11a79d25d0"
uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
version = "0.3.10"

[[deps.OffsetArrays]]
deps = ["Adapt"]
git-tree-sha1 = "c870a0d713b51e4b49be6432eff0e26a4325afee"
uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
version = "1.10.6"

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

[[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.Optim]]
deps = ["Compat", "FillArrays", "LineSearches", "LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "PositiveFactorizations", "Printf", "SparseArrays", "StatsBase"]
git-tree-sha1 = "7863df65dbb2a0fa8f85fcaf0a41167640d2ebed"
uuid = "429524aa-4258-5aef-a3af-852621145aeb"
version = "1.4.1"

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

[[deps.OrdinaryDiffEq]]
deps = ["Adapt", "ArrayInterface", "DataStructures", "DiffEqBase", "DocStringExtensions", "ExponentialUtilities", "FastClosures", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "Logging", "LoopVectorization", "MacroTools", "MuladdMacro", "NLsolve", "Polyester", "RecursiveArrayTools", "Reexport", "SparseArrays", "SparseDiffTools", "StaticArrays", "UnPack"]
git-tree-sha1 = "66816b9f09c3925ae3b072eb9ef1835c68282ed1"
uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
version = "5.64.0"

[[deps.PDMats]]
deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
git-tree-sha1 = "4dd403333bcf0909341cfe57ec115152f937d7d8"
uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
version = "0.11.1"

[[deps.Parameters]]
deps = ["OrderedCollections", "UnPack"]
git-tree-sha1 = "2276ac65f1e236e0a6ea70baff3f62ad4c625345"
uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
version = "0.12.2"

[[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.PoissonRandom]]
deps = ["Random", "Statistics", "Test"]
git-tree-sha1 = "44d018211a56626288b5d3f8c6497d28c26dc850"
uuid = "e409e4f3-bfea-5376-8464-e040bb5c01ab"
version = "0.4.0"

[[deps.Polyester]]
deps = ["ArrayInterface", "BitTwiddlingConvenienceFunctions", "CPUSummary", "IfElse", "ManualMemory", "PolyesterWeave", "Requires", "Static", "StrideArraysCore", "ThreadingUtilities"]
git-tree-sha1 = "4a67fa65b250256e214a6a2482c87c835ec0fc02"
uuid = "f517fe37-dbe3-4b94-8317-1923a5111588"
version = "0.5.0"

[[deps.PolyesterWeave]]
deps = ["BitTwiddlingConvenienceFunctions", "CPUSummary", "IfElse", "Static", "ThreadingUtilities"]
git-tree-sha1 = "371a19bb801c1b420b29141750f3a34d6c6634b9"
uuid = "1d0040c9-8b98-4ee7-8388-3f51789ca0ad"
version = "0.1.0"

[[deps.PositiveFactorizations]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "17275485f373e6673f7e7f97051f703ed5b15b20"
uuid = "85a6dd25-e78a-55b7-8502-1745935b8125"
version = "0.2.4"

[[deps.PreallocationTools]]
deps = ["ArrayInterface", "ForwardDiff", "LabelledArrays"]
git-tree-sha1 = "9e917b108c4aaf47e8606542325bd2ccbcac7ca4"
uuid = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
version = "0.1.0"

[[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.QuadGK]]
deps = ["DataStructures", "LinearAlgebra"]
git-tree-sha1 = "12fbe86da16df6679be7521dfb39fbc861e1dc7b"
uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
version = "2.4.1"

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

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

[[deps.Random123]]
deps = ["Libdl", "Random", "RandomNumbers"]
git-tree-sha1 = "0e8b146557ad1c6deb1367655e052276690e71a3"
uuid = "74087812-796a-5b5d-8853-05524746bad3"
version = "1.4.2"

[[deps.RandomNumbers]]
deps = ["Random", "Requires"]
git-tree-sha1 = "441e6fc35597524ada7f85e13df1f4e10137d16f"
uuid = "e6cf234a-135c-5ec9-84dd-332b85af5143"
version = "1.4.0"

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

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

[[deps.RecursiveFactorization]]
deps = ["LinearAlgebra", "LoopVectorization", "Polyester", "StrideArraysCore", "TriangularSolve"]
git-tree-sha1 = "575c18c6b00ce409f75d96fefe33ebe01575457a"
uuid = "f2c3362d-daeb-58d1-803e-2bc74f2840b4"
version = "0.2.4"

[[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.ResettableStacks]]
deps = ["StaticArrays"]
git-tree-sha1 = "256eeeec186fa7f26f2801732774ccf277f05db9"
uuid = "ae5879a3-cd67-5da8-be7f-38c6eb64a37b"
version = "1.1.1"

[[deps.Rmath]]
deps = ["Random", "Rmath_jll"]
git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f"
uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
version = "0.7.0"

[[deps.Rmath_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7"
uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
version = "0.3.0+0"

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

[[deps.SIMDTypes]]
git-tree-sha1 = "330289636fb8107c5f32088d2741e9fd7a061a5c"
uuid = "94e857df-77ce-4151-89e5-788b33177be4"
version = "0.1.0"

[[deps.SLEEFPirates]]
deps = ["IfElse", "Static", "VectorizationBase"]
git-tree-sha1 = "947491c30d4293bebb00781bcaf787ba09e7c20d"
uuid = "476501e8-09a2-5ece-8869-fb82de89a1fa"
version = "0.6.26"

[[deps.SciMLBase]]
deps = ["ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "RecipesBase", "RecursiveArrayTools", "StaticArrays", "Statistics", "Tables", "TreeViews"]
git-tree-sha1 = "ff686e0c79dbe91767f4c1e44257621a5455b1c6"
uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
version = "1.18.7"

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

[[deps.Setfield]]
deps = ["ConstructionBase", "Future", "MacroTools", "Requires"]
git-tree-sha1 = "fca29e68c5062722b5b4435594c3d1ba557072a3"
uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46"
version = "0.7.1"

[[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.SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[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 = "aebcead0644d3b3396c205a09544590b5115e282"
uuid = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
version = "1.16.4"

[[deps.SpecialFunctions]]
deps = ["ChainRulesCore", "LogExpFunctions", "OpenSpecFun_jll"]
git-tree-sha1 = "a322a9493e49c5f3a10b50df3aedaf1cdb3244b7"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "1.6.1"

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

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

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

[[deps.StatsAPI]]
git-tree-sha1 = "1958272568dc176a1d881acb797beb909c785510"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.0.0"

[[deps.StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "8cbbc098554648c84f79a463c9ff0fd277144b6c"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.10"

[[deps.StatsFuns]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
git-tree-sha1 = "46d7ccc7104860c38b11966dd1f72ff042f382e4"
uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
version = "0.9.10"

[[deps.StochasticDiffEq]]
deps = ["ArrayInterface", "DataStructures", "DiffEqBase", "DiffEqJump", "DiffEqNoiseProcess", "DocStringExtensions", "FillArrays", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "Logging", "MuladdMacro", "NLsolve", "OrdinaryDiffEq", "Random", "RandomNumbers", "RecursiveArrayTools", "Reexport", "SparseArrays", "SparseDiffTools", "StaticArrays", "UnPack"]
git-tree-sha1 = "1d5dff61a15532dd6c67fafebc52f2feaa2b3b82"
uuid = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
version = "6.37.1"

[[deps.StrideArraysCore]]
deps = ["ArrayInterface", "CloseOpenIntervals", "IfElse", "LayoutPointers", "ManualMemory", "Requires", "SIMDTypes", "Static", "ThreadingUtilities"]
git-tree-sha1 = "1258e25e171aec339866f283a11e7d75867e77d7"
uuid = "7792a7ef-975c-4747-a70f-980b88e8d1da"
version = "0.2.4"

[[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 = "1162ce4a6c4b7e31e0e6b14486a6986951c73be9"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.5.2"

[[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.ThreadingUtilities]]
deps = ["ManualMemory"]
git-tree-sha1 = "03013c6ae7f1824131b2ae2fc1d49793b51e8394"
uuid = "8290d209-cae3-49c0-8002-c8c24d57dab5"
version = "0.4.6"

[[deps.TreeViews]]
deps = ["Test"]
git-tree-sha1 = "8d0d7a3fe2f30d6a7f833a5f19f7c7a5b396eae6"
uuid = "a2a6695c-b41b-5b7d-aed9-dbfdeacea5d7"
version = "0.3.0"

[[deps.TriangularSolve]]
deps = ["CloseOpenIntervals", "IfElse", "LayoutPointers", "LinearAlgebra", "LoopVectorization", "Polyester", "Static", "VectorizationBase"]
git-tree-sha1 = "ed55426a514db35f58d36c3812aae89cfc057401"
uuid = "d5829a12-d9aa-46ab-831f-fb7c9ab06edf"
version = "0.1.6"

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

[[deps.UnPack]]
git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b"
uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
version = "1.0.2"

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

[[deps.UnicodePlots]]
deps = ["Crayons", "Dates", "SparseArrays", "StatsBase"]
git-tree-sha1 = "7ee3014ae511a5ef11523c38c216e7ad59650ce7"
uuid = "b8865327-cd53-5732-bb35-84acbb429228"
version = "2.4.0"

[[deps.VectorizationBase]]
deps = ["ArrayInterface", "CPUSummary", "HostCPUFeatures", "Hwloc", "IfElse", "LayoutPointers", "Libdl", "LinearAlgebra", "SIMDTypes", "Static"]
git-tree-sha1 = "43c605e008ac67adb672ef08721d4720dfe2ad41"
uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f"
version = "0.21.7"

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

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

[[deps.ZygoteRules]]
deps = ["MacroTools"]
git-tree-sha1 = "9e7a1e8ca60b742e508a315c17eef5211e7fbfd7"
uuid = "700de1a5-db45-46bc-99cf-38207098b444"
version = "0.2.1"

[[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:
# ╠═6bf88b5c-3abc-4b03-90bb-a7eb8cb68577
# ╟─96c9b046-6e40-44b1-a19c-b7061f22470d
# ╟─d78532b9-b17e-4a4a-bbab-e249bf852eba
# ╟─41487611-deba-49fc-a5f7-f4f63f7935c6
# ╟─32a692ab-4835-4943-8b46-8345ce881eb5
# ╟─c86cf1e9-d0b6-48e8-8a01-0dcc591afc85
# ╟─2195a7a4-30d6-47fe-92db-d7ebef67ef03
# ╟─fd2db517-2282-43da-a543-a1ef59a439c2
# ╟─5ab5321a-15b9-11ec-3321-b7d7cab1ba35
# ╟─aa7aa58b-e2e7-4b40-aeab-38137f8d4e2d
# ╠═1df22a2f-c344-4d71-b44d-9a821c14e548
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
import Base: ==
function solve_hfg(domain_name, problem_name; max_levels=20)
    graph = create_funnel_graph(domain_name, problem_name; max_levels=max_levels)

    model = Model(Gurobi.Optimizer)
    set_silent(model)

    T = graph.num_levels
    max_nf = max([length(get_fluents(graph, i)) for i=1:T]...)
    # max_na = max([length(graph.acts[i]) for i=1:T-1]...)
    max_na = max([length([act for act in graph.acts[i] if act.name != :noop && act.name != :maintain]) for i=1:T-1]...)

    @variable(model, p[i=1:max_nf, j=1:T], Bin)
    @variable(model, a[i=1:max_na, j=1:T-1], Bin)

    @variable(model, xᵣ[i=1:2, j=1:T])
    @variable(model, xₒ[i=1:2, j=1:T])

    #inits
    num_init_props = length(get_fluents(graph, 1))
    [@constraint(model, p[i, 1] == 1) for i = 1 : num_init_props]
    [@constraint(model, p[i, 1] == 0) for i = num_init_props+1 : max_nf]
    @constraint(model, xᵣ[1, 1] == 0.25)
    @constraint(model, xᵣ[2, 1] == 0.25)
    @constraint(model, xₒ[1, 1] == 10.5)
    @constraint(model, xₒ[2, 1] == 10.5)

    #goals
    goalinds = get_fluent_terminal_indices(graph; loc=:goals)
    println("GOAL: ",goalinds)
    [@constraint(model, p[f, T] == 1) for f in goalinds]
    # [@constraint(model, p[f, T] == 0) for f=1:max_nf if !(f in goalinds)]
    @constraint(model, xᵣ[1, T] == 40.5)
    @constraint(model, xᵣ[2, T] == 40.5)
    @constraint(model, xₒ[1, T] == 20.5)
    @constraint(model, xₒ[2, T] == 20.5)

    #constraints
    #precondition constraint
    for t=1:T-1
        fluents = get_fluents(graph, t)
        for f=1:length(fluents)
            fluent = fluents[f]
            act_inds = get_actions_precondition_sat(fluent, t, graph)

            [@constraint(model, a[i,t] ≤ p[f,t]) for i in act_inds]
        end
    end

    #effect constraint
    for t=1:T-1
        fluents = get_fluents(graph, t+1)
        for f=1:length(fluents)
            fluent = fluents[f]
            act_inds = get_actions_effect_sat(fluent, t, graph)
            fnext = get_fluent_index(fluent, t+1, graph)
            # if t==8 println("sats ",fluent," ",act_inds) end
            if !isnothing(fnext) && !isempty(act_inds)
                @constraint(model, p[fnext, t+1] ≤ sum([a[i,t] for i in act_inds]))
            end
        end
        na = length(graph.acts[t])
        @constraint(model, a[na+1:max_na, t] .== 0)
    end

    #mutex constraint

    # for t=1:T-1
    #     fluents = get_fluents(graph, t)
    #     for f=1:length(fluents)
    #         fluent = fluents[f]
    #         act_pairs = get_mutex_actions(fluent, t, graph)
    #         if t==8 println(act_pairs) end
    #         [@constraint(model, a[i[1], t] + a[i[2], t] ≤ 1) for i in act_pairs]
    #     end
    # end
    for t=1:T-2
        act_mutexes = get_mutex_action_indices(t, graph)
        println(act_mutexes)
        [@constraint(model, a[i[1], t] + a[i[2], t] ≤ 1) for i in act_mutexes]
    end

    for t=1:T-1
        prop_pairs = get_mutex_indices(t, graph)
        # println("pairs ",prop_pairs)
        [@constraint(model, p[i[1], t] + p[i[2], t] ≤ 1) for i in prop_pairs]
    end
    #=

    #continuous constraints

    for t=1:T
        fluents = get_fluents(graph, t)
        nd = length(graph.props[t][:discrete])
        for f = nd+1:length(fluents)
            fluent = fluents[f]
            if fluent.name == :robot
                @constraint(model, p[f, t] => {fluent.r[1].þ ≤xᵣ[1, t]≤ fluent.r[2].þ})
                @constraint(model, p[f, t] => {fluent.r[3].þ ≤xᵣ[2, t]≤ fluent.r[4].þ})
            elseif fluent.name == :b1
                @constraint(model, p[f, t] => {fluent.r[1].þ ≤xₒ[1, t]≤ fluent.r[2].þ})
                @constraint(model, p[f, t] => {fluent.r[3].þ ≤xₒ[2, t]≤ fluent.r[4].þ})
            end
        end
        # for f = 1:length(fluents)
        #     if fluents[f].name == :robot
        #         @constraint(model, p[f,t] == 1)
        #     end
        # end
    end
    for t=1:T-1
        act_inds = get_pick_place_inds(graph, t)
        for i in act_inds
            act = graph.acts[t][i]
            for reg in act.continuous_prec
                if reg.name == :robot
                    @constraint(model, a[i, t] => {reg.r[1].þ ≤xᵣ[1, t]≤ reg.r[2].þ})
                    @constraint(model, a[i, t] => {reg.r[3].þ ≤xᵣ[2, t]≤ reg.r[4].þ})
                elseif reg.name == :b1
                    @constraint(model, a[i, t] => {reg.r[1].þ ≤xₒ[1, t]≤ reg.r[2].þ})
                    @constraint(model, a[i, t] => {reg.r[3].þ ≤xₒ[2, t]≤ reg.r[4].þ})
                end
            end
        end
        # nempty = length([act for act in graph.acts[t] if act.name != :noop && act.name != :maintain])
        # println("level ",t," num acts ",nempty)
        # [@constraint(model, a[i,t] == 0) for i = nempty+1:max_na]
        # @constraint(model, sum([a[i, t] for i=1:max_na])==1 )

    end

    =#
    t=@variable(model)
    @constraint(model, [t; sum([xᵣ[:,t]-xᵣ[:,t+1] for t=1:T-2])] in SecondOrderCone())
    @objective(model, Min, t)
    @time optimize!(model)
    sol = value.(a)
    # sol = value.(p)
    # sol = value.(xᵣ)
    render_action(sol, graph)
    # graph
    # value.(xₒ)
end

==(reg1::Region, reg2::Region) = reg1.r[1].ϕ₁ == reg2.r[1].ϕ₁ && reg1.r[1].ϕ₂ == reg2.r[1].ϕ₂ && reg1.r[1].ϕ₃ == reg2.r[1].ϕ₃ && reg1.r[1].þ == reg2.r[1].þ &&
reg1.r[2].ϕ₁ == reg2.r[2].ϕ₁ && reg1.r[2].ϕ₂ == reg2.r[2].ϕ₂ && reg1.r[2].ϕ₃ == reg2.r[2].ϕ₃ && reg1.r[2].þ == reg2.r[2].þ &&
reg1.r[3].ϕ₁ == reg2.r[3].ϕ₁ && reg1.r[3].ϕ₂ == reg2.r[3].ϕ₂ && reg1.r[3].ϕ₃ == reg2.r[3].ϕ₃ && reg1.r[3].þ == reg2.r[3].þ &&
reg1.r[4].ϕ₁ == reg2.r[4].ϕ₁ && reg1.r[4].ϕ₂ == reg2.r[4].ϕ₂ && reg1.r[4].ϕ₃ == reg2.r[4].ϕ₃ && reg1.r[4].þ == reg2.r[4].þ

function get_fluent_terminal_indices(graph; loc=:init)
    indices = Set()
    T = graph.num_levels
    if loc == :init
        for i=1:length(graph.initprops[:discrete])
            ind = findall(x->x == graph.initprops[:discrete], graph.props[0][:discrete])[1]
            push!(indices, ind)
        end
    else
        for i=1:length(graph.goalprops[:discrete])
            ind = findall(x->x == graph.goalprops[:discrete][i], graph.props[T][:discrete])[1]
            push!(indices, ind)
        end
        nf = length(graph.props[T][:discrete])

        for i=1:length(graph.goalprops[:continuous])
            ind = findall(x-> intersects(x, graph.goalprops[:continuous][i]), graph.props[T][:continuous])[1]
            push!(indices, ind+nf)
        end
    end
    return collect(indices)
end

# function get_fluents(graph, level)
#     fluents=[graph.props[level][:discrete]...,graph.props[level][:continuous]...]
#     return fluents
# end

function get_fluents(graph, level)
    fluents=Any[graph.props[level][:discrete]...]
    for f in graph.props[level][:continuous]
        if f.name != :b1
            push!(fluents, f)
        else
            if !any([f.name == :b1 for f in fluents])
                push!(fluents, f)
            end
        end
    end
    return fluents
end

function get_fluent_index(fluent, t, graph)
    fluents = get_fluents(graph, t)
    for i=1:length(fluents)
        if fluents[i] == fluent
            return i
        end
    end
    return nothing
end

#gets discrete actions as well as funnels too
function get_actions_precondition_sat(fluent, level, graph)
    # actions = graph.acts[level]
    actions = [act for act in graph.acts[level] if act.name != :noop && act.name != :maintain]
    act_inds = []
    for i=1:length(actions)
        act = graph.acts[level][i]
        if fluent in act.pos_prec
            push!(act_inds, i)
        end
        if typeof(fluent) == Region
            if !isempty(act.continuous_prec)
                for j = 1:length(act.continuous_prec)
                    if act.continuous_prec[j].name == fluent.name
                        if intersects(act.continuous_prec[j], fluent)
                            # if level == 8 && act.name == :place println("sat") end
                            push!(act_inds, i)
                        end
                    end
                end
            end
        end
    end
    return act_inds
end

#gets discrete actions as well as funnels too
function get_actions_effect_sat(fluent, level, graph)
    # actions = graph.acts[level]
    actions = [act for act in graph.acts[level] if act.name != :noop && act.name != :maintain]

    act_inds = []
    for i=1:length(actions)
        act = graph.acts[level][i]
        if fluent in act.pos_eff
            push!(act_inds, i)
        end
        if typeof(fluent) == Region
            for er in act.end_region
                if er.name == fluent.name
                    if intersects(er, fluent)
                        push!(act_inds, i)
                        # if level == 9 println(act.name," intersects") end
                        break
                    end
                end
            end
        end
    end
    return act_inds
end

function get_mutex_actions(fluent, level, graph)
    pref = get_actions_precondition_sat(fluent, level, graph)
    preneff = union(pref, get_actions_effect_sat(fluent, level, graph))
    act_mutexes = []
    [push!(act_mutexes, [i,j]) for i in pref for j in preneff if i!=j && !([j,i] in act_mutexes)]

    return act_mutexes
end

function get_mutex_action_indices(level, graph)
    act_mutexes = []
    for i=1:length(graph.acts[level])
        a=graph.acts[level][i]
        for j=1:length(graph.acts[level])
            b=graph.acts[level][j]
            if i!=j && a.name != :maintain && a.name != :noop  && b.name != :maintain && b.name != :noop
                if !isempty(intersect(a.neg_eff, union(b.pos_prec, b.pos_eff))) ||
                    !isempty(intersect(b.neg_eff, union(a.pos_prec, a.pos_eff)))
                    if !([j,i] in act_mutexes) push!(act_mutexes, [i,j]) end
                end
            end
        end
    end
    return act_mutexes
end

function get_mutex_indices(level, graph)
    prop_mutexes = []
    pairs = [pair for pair in graph.μprops[level] if typeof(pair[1])!= Region]
    props = get_fluents(graph, level)
    for pair in pairs
        i = findall(x->x == pair[1], props)[1]
        j = findall(x->x == pair[2], props)[1]

        if !([j,i] in prop_mutexes) && [i,j]!=[1,2] push!(prop_mutexes, [i,j]) end
    end
    return prop_mutexes
end

function get_pick_place_inds(graph, level)
    act_inds = []
    for i=1:length(graph.acts[level])
        if graph.acts[level][i].name == :pick || graph.acts[level][i].name == :place
            push!(act_inds, i)
        end
    end
    return act_inds
end

function get_nonactinds(graph, level)
    nonactinds = []
    for i=1:length(graph.acts[level])
        if graph.acts[level][i].name == :noop || graph.acts[level][i].name == :maintain
            push!(nonactinds, i)
        end
    end
    return nonactinds
end

function render_action(y, graph)
    T = graph.num_levels-1
    actions = []
    for t=1:T
        # try
        ind = findall(x->x == 1.0, y[1:end,t])[1]
        act  = graph.acts[t][ind]
        # act = all_actions[ind]
        push!(actions, (act.name, act.params))
        # catch
        # end
    end
    # println(actions)
    return actions

end

function get_same_region_types(graph, level, type)
    prop_inds = []
    fluents = get_fluents(graph, level)
    for i = 1:length(fluents)
        if typeof(fluents[i]) == Region
            if fluents[i].name == type
                push!(prop_inds, i)
            end
        end
    end
    return prop_inds
end

ext: jl


text:
using Documenter, BinaryTwoStageDesigns

makedocs(
  format = :html,
  sitename = "BinaryTwoStageDesigns",
  pages = [
    "index.md",
    "Sample Spaces" => "sample_space.md",
    "Design Parameters" => "parameters.md",
    "Two-Stage Designs" => "designs.md",
    "Finding Optimal Designs" => "optimal_designs.md",
    "Inference" => "inference.md"
  ],
  repo = "https://github.com/imbi-heidelberg/BinaryTwoStageDesigns/{commit}"
)

max_stars_repo_name: kkmann/binaryTwoStageDesigns

max_issues_repo_name: kkmann/binaryTwoStageDesigns

max_forks_repo_name: kkmann/binaryTwoStageDesigns


text:
"""
    execute_sts_daily_api_call(user_model::PSUserModel, ticker_symbol::String;
        data_type::Symbol = :json, outputsize::Symbol = :compact, logger::Union{Nothing,AbstractLogger} = nothing) -> PSResult

Stuff will go here. Awesome stuff, the most beautiful stuff ever.
"""
function execute_sts_daily_api_call(user_model::PSUserModel, ticker_symbol::String;
    data_type::Symbol = :json, outputsize::Symbol = :compact, logger::Union{Nothing,AbstractLogger} = nothing)::PSResult

    try

        # get key's and id -
        api_key = user_model.alphavantage_api_key

        # use the alpha_vantage_api to download the data -
        url = "$(alphavantage_api_url_string)?function=TIME_SERIES_DAILY&symbol=$(ticker_symbol)&apikey=$(api_key)&datatype=$(string(data_type))&outputsize=$(string(outputsize))"
        api_call_raw_data = http_get_call_with_url(url) |> checkresult

        # log if we have one -
        if (isnothing(logger) == false)
            log_api_call(logger, user_model, url)
        end

        # make the handler calls, depending upon the data type -
        if (data_type == :json)

            # process -
            data_series_key = "Time Series (Daily)"
            return process_raw_json_api_data_sts(api_call_raw_data, data_series_key)

        elseif (data_type == :csv)

            # return the data back to the caller -
            return process_raw_csv_api_data(api_call_raw_data)
        else
            # formulate the error message -
            error_message = "$(data_type) is not supported. Supported values are {:json,:csv}"

            # throw -
            return PSResult{PSError}(PSError(error_message))
        end

    catch error
        return PSResult(error)
    end
end

"""
    execute_sts_adjusted_daily_api_call(user_model::PSUserModel, ticker_symbol::String;
        data_type::Symbol = :json, outputsize::Symbol = :compact, logger::Union{Nothing,AbstractLogger} = nothing) -> PSResult

Stuff will go here. Awesome stuff, the most beautiful stuff ever.
"""
function execute_sts_adjusted_daily_api_call(user_model::PSUserModel, ticker_symbol::String;
    data_type::Symbol = :json, outputsize::Symbol = :compact, logger::Union{Nothing,AbstractLogger} = nothing)::PSResult

    try

        # get key's and id -
        api_key = user_model.alphavantage_api_key

        # call to alpha_vantage_api to get data
        url = "$(alphavantage_api_url_string)?function=TIME_SERIES_DAILY_ADJUSTED&symbol=$(ticker_symbol)&apikey=$(api_key)&datatype=$(string(data_type))&outputsize=$(string(outputsize))"
        api_call_raw_data = http_get_call_with_url(url) |> checkresult

        # call to logger
        if (isnothing(logger) == false)
            log_api_call(logger, user_model, url)
        end

        # parse json if data called as a .json
        if (data_type == :json)

            # process json
            data_series_key = "Time Series (Daily)"
            return process_raw_json_api_data_sts_daily_adjusted(api_call_raw_data, data_series_key)

        elseif (data_type == :csv)

            # return process .csv
            return process_raw_csv_api_data(api_call_raw_data)

        else

            # tell user they requested an unsupported type of data
            error_message = "$(data_type) isn't supported by AlphaVantage. Supported values are {:json, :csv}"
            return PSResult{PSError}(PSError(error_message))
        end

    catch error
        return PSResult(error)
    end
end


text:
module Cthulhu

using TerminalMenus
using InteractiveUtils

export descend, @descend, descend_code_typed, @descend_code_typed

struct Callsite
    id::Int # ssa-id
    f
    tt
    rt
end

function Base.string(c::Callsite)
    io = IOBuffer()
    print(io, "%", c.id, " = invoke ", c.f, "(")
    TT = c.tt.parameters
    for T in TT
        print(io, "::", T, ",")
    end
    !isempty(TT) && seek(io, position(io)-1)
    print(io, ")", "::", c.rt)
    String(take!(io))
end

function Callsite(id, mi, rt)
    f = getfield(mi.def.module, mi.def.name)
    tt = Tuple{mi.specTypes.parameters[2:end]...}
    return Callsite(id, f, tt, rt)
end

"""
  @descend_code_typed

  Evaluates the arguments to the function or macro call, determines their
types, and calls `code_typed on the resulting expression.
"""
macro descend_code_typed(ex0...)
    InteractiveUtils.gen_call_with_extracted_types_and_kwargs(__module__, :descend_code_typed, ex0)
end

"""
  @descend

  Shortcut for [`@descend_code_typed`](@ref).
"""
macro descend(ex0...)
    esc(:(@descend_code_typed($(ex0...))))
end

"""
    descend_code_typed(f, tt; kwargs...)

Given a function and a tuple-type, interactively explore the output of
`code_typed` by descending into `invoke` statements. Type enter to select an
`invoke` to descend into, select ↩  to ascend, and press q or control-c to
quit.

# Usage:
```julia
function foo()
    T = rand() > 0.5 ? Int64 : Float64
    sum(rand(T, 100))
end

descend_code_typed(foo, Tuple{})
```
"""
function descend_code_typed(f, @nospecialize(tt); kwargs...)
    try
        _descend(f, tt; kwargs...)
    catch x
        if x isa InterruptException
            return nothing
        else
            rethrow(x)
        end
    end
    return nothing
end

"""
  descend

  Shortcut for [`descend_code_typed`](@ref).
"""
const descend = descend_code_typed

function _descend(@nospecialize(F), @nospecialize(TT); kwargs...)
    methods = code_typed(F, TT; kwargs...)
    if isempty(methods)
        println("$(string(Callsite(-1 ,F, TT, Any))) has no methods")
        return
    end
    CI, rt = first(methods)
    callsites = Callsite[]
    for (id, c) in enumerate(CI.code)
        if c isa Expr
            if c.head === :invoke
                rt = CI.ssavaluetypes[id]
                push!(callsites, Callsite(id, c.args[1], rt))
            elseif c.head === :call
                rt = CI.ssavaluetypes[id]
                if c.args[1] isa Function
                    f = c.args[1]
                elseif c.args[1] isa GlobalRef
                    mod = c.args[1].mod
                    name = c.args[1].name
                    f = getfield(mod, name)
                elseif c.args[1] isa Core.SSAValue
                    # probably somthing of form
                    # %1 = Base.Broadcast.materialize::Const(materialize, false)
                    # ...
                    # %9 = (%1)(%8)::Any
                    _T = CI.ssavaluetypes[c.args[1].id]
                    if _T isa Core.Compiler.Const && _T.val isa Function
                        f = _T.val
                    elseif _T isa Type
                        continue
                    else
                        @warn "Don't know how to handle call: " c
                        dump(c.args[1])
                        continue
                    end
                else
                    @warn "Don't know how to handle call: " c
                    dump(c.args[1])
                    continue
                end
                args = c.args[2:end]
                types = Any[]
                for arg in args
                    if arg isa Core.SSAValue
                        T = CI.ssavaluetypes[arg.id]
                    elseif arg isa Core.SlotNumber
                        T = TT.parameters[arg.id - 1]
                    elseif arg isa Expr # arrrgh
                        @assert arg.head === :static_parameter
                        T = typeof(args.args[1])
                    else
                        T = typeof(arg)
                    end
                    if T isa Core.Compiler.Const
                        T = typeof(T.val)
                    end
                    push!(types, T)
                end

                # Filter out builtin functions and intrinsic function
                if f isa Core.Builtin || f isa Core.IntrinsicFunction
                    continue
                end

                # Filter out abstract signatures
                if any(isabstracttype, types)
                    continue
                end

                push!(callsites, Callsite(id, f, Tuple{types...}, rt))
            end
        end
    end
    while true
        println()
        println("│ ─ $(string(Callsite(-1, F, TT, rt)))")
        display(CI=>rt)
        println()
        TerminalMenus.config(cursor = '•', scroll = :wrap)
        menu = RadioMenu(vcat(map(string, callsites), ["↩ "]))
        println("In `$F` select a call to descend into or ↩ to ascend. [q] to quit.")
        cid = request(menu)
        if cid == length(callsites) + 1
            break
        end
        if cid == -1
            throw(InterruptException())
        end
        callsite = callsites[cid]
        _descend(callsite.f, callsite.tt)
    end
end

end


text:
abstract type VariableType end
abstract type Quantitative <: VariableType end
struct Categorical <: VariableType end
struct Continuous <: Quantitative end
struct Discrete <: Quantitative end

# Variable type of a single value
function vartype(v)
    if is_quantitative(v)
        if is_continuous(v)
            return Continuous
        elseif is_discrete(v)
            return Discrete
        else
            return Quantitative
        end
    elseif is_categorical(v)
        # Not currently making a distinction between nominal/dichotomous or ordinal
        return Categorical
    else
        error("Could not classify variable")
    end
end

# Variable type for a vector
# will check the lesser of the vector length or `searchlen`
function vartype(v::T where T <: AbstractVector, searchlen = 100)
    subarr = view(v, 1:min(length(v), searchlen))

    # Assess
    if any(is_categorical.(subarr))
        return Categorical
    else
        if any(is_continuous.(subarr))
            return Continuous
        elseif any(is_discrete.(subarr))
            return Discrete
        else
            return Quantitative
        end
    end
end

# Variable types for a matrix
function vartype(v::T where T <: AbstractMatrix, dim=1, searchlen = 100)
    if dim !== 1 && dim !== 2
        error("dim must equal 1 or 2")
    end

    len = size(v, dim == 1 ? 2 : 1)
    vartypes = Array{Any, 1}(len)

    @inbounds for i = 1:len
        if dim == 1
            vartypes[i] = vartype(v[:, i], searchlen)
        else
            vartypes[i] = vartype(v[i, :], searchlen)
        end
    end

    return vartypes
end

# Type => Variable Type logic
is_quantitative(v) = typeof(v) <: Number && typeof(v) !== Bool
is_categorical(v) = !is_quantitative(v)
is_continuous(v) = typeof(v) <: AbstractFloat || typeof(v) <: Complex || typeof(v) <: Irrational
is_discrete(v) = typeof(v) <: Int || typeof(v) <: Rational


text:
# ExchangeMarket constructor

"""
    AssetsPair

Type for asset pairs exchange/trade.

Asset pairs are made of two String fields, a `base_asset` and a `quote_asset`.

# Examples

```jldoctest
julia> AssetsPair("EUR", "BRL")
AssetsPair("EUR", "BRL")

julia> AssetsPair("euro", "BRL")
ERROR: ArgumentError: The given code symbol pair ("euro", "BRL") is not allowed, both should be all in ascii uppercase letters and at least three-character long.
Stacktrace:
  ...
```
"""
struct AssetsPair
    base_asset::String
    quote_asset::String
end

"""
    Rate

Type for exchange/stock rates.

A rate is simply a positive Number.

The structure's constructor checks whether this requirement is met,
otherwise an ArgumentError is thrown.

# Examples

```jldoctest
julia> Rate(1.2)
Rate(1.2)

julia> Rate(-2)
ERROR: ArgumentError: The exchange rate must be a positive number
Stacktrace:
  ...
```
"""
struct Rate
    value::Number
    Rate(r) = r > Unitful.zero(r) ? new(r) :
        throw(ArgumentError("The rate must be positive"))
end

"""
    ExchangeMarket

Type used for a dictionary of `AssetsPair() => Rate()` pair quotes.

It is given as a Dict{AssetsPair, Rate}, where the keys are
the pairs with the base and quote assets and the value
is the exchange/trade rate for this pair (i..e. how much in quote units
is needed to buy/trande one unit of the base asset).

For instance, the market

    mkt = ExchangeMarket(AssetsPair("EUR", "USD") => Rate(1.164151))

contains the pair `AssetsPair("EUR", "USD")` and the exchange rate
`Rate(1.164151)`, which means that one can buy 1 EUR with 1.164151 USD.
"""
ExchangeMarket = Dict{AssetsPair, Rate}

"""
    get_rate(u::String, v::String, rate_value::Number)

Return the exchange rate as a Unitful.Quantity in proper units.
"""
function get_rate(u::String, v::String, rate_value::Number)
    return Main.eval(Meta.parse("(" * string(rate_value) * ")u\"" * u * "/" * v * "\""))
end

"""
    is_asset(code_abbr::String)

Check whether `code_abbr` refers to a registered asset.
"""
function is_asset(code_abbr)
    has_unit = m->(isdefined(m,Symbol(code_abbr)))
    return length(findall(has_unit, Unitful.unitmodules)) > 0
end

"""
    generate_exchmkt(d::Dict{Tuple{String,String},T}) where {T<:Number}

Generate an instance of `ExchangeMarket` from a dictionary of base-quote-value rates.

# Examples

```jldoctest
julia> generate_exchmkt(Dict(("EUR", "USD") => 1.164151))
Dict{AssetsPair,Float64} with 1 entry:
  AssetsPair("EUR", "USD") => Rate(1.16415)
```
"""
function generate_exchmkt(d::Dict{Tuple{String,String},T}) where {T<:Number}
    valid_d = Dict(key => value for (key, value) in d if is_asset(key[1]) && is_asset(key[2]))
    return Dict([AssetsPair(key[1], key[2]) => Rate(get_rate(key[2], key[1], value)) for (key,value) in valid_d])
end

"""
    generate_exchmkt(a::Array{Pair{Tuple{String,String},Float64},1})

Generate an instance of `ExchangeMarket` from an array of base-quote-value rates.

# Examples

```jldoctest
julia> generate_exchmkt([("EUR","USD") => 1.19536, ("USD","EUR") => 0.836570])
Dict{AssetsPair,Float64} with 2 entries:
  AssetsPair("EUR", "USD") => Rate(1.19536)
  AssetsPair("USD", "EUR") => Rate(0.83657)
```
"""
function generate_exchmkt(a::Array{Pair{Tuple{String,String},T},1}) where {T<:Number}
    return generate_exchmkt(Dict(a))
end

"""
    generate_exchmkt(a::Array{Pair{Tuple{String,String},Float64},1})

Generates an instance of `ExchangeMarket` from a single of base-quote-value rate.

# Examples

```jldoctest
julia> generate_exchmkt(("EUR", "USD") => 1.164151)
Dict{AssetsPair,Float64} with 1 entry:
  AssetsPair("EUR", "USD") => Rate(1.16415)
```
"""
function generate_exchmkt(p::Pair{Tuple{String,String},T}) where {T<:Number}
    return generate_exchmkt([p])
end

"""
    uconvert(u::Units, x::Quantity, e::ExchangeMarket; mode::Int=1)

Convert between assets, allowing for inverse and secondary rates.

If mode=1, which is the default, a direct conversion is attempted, i.e.
if the given market includes the conversion rate from `unit(x)` to `u`,
then the conversion takes place with this rate.

If mode=-1 and the given market includes the rate from `u` to `unit(x)`,
then the conversion of `x` to `u` is achieved with the rate which is the
multiplicative inverse of the exchange rate from `u` to `unit(x)`.

If mode=2, and the given market includes the rate from
`unit(x)` to an intermediate asset `v` and from  `v` to `u`, then
the exchange/trade takes place with the product of these two rates.
If there is more than one intermediate rates available, then the first
one encountered in a nested loop in which the second pair is in the
inner loop is the one chosen.

If mode=-2, a combination of `-1` and `2` is used, i.e. an intermediate
asset is used for the inverse rate from `u` to `unit(x)`.

An `ArgumentError` is thrown if mode is none of the above or if `u` or `x`
are not assets, or if the necessary rates cannot be found.

# Examples

Assuming `forex_exchmkt["2020-11-01"]` ExchangeMarket contains the key-value
pair `("EUR","BRL") => 6.685598`, then the following exchange takes place:

```jldoctest
julia> uconvert(u"BRL", 1u"EUR", forex_exchmkt["2020-11-01"])
6.685598 BRL
julia> uconvert(u"BRL", 1u"BRL", forex_exchmkt["2020-11-01"], mode=-1)
0.149575251159283 EUR
```
"""
function uconvert(u::Unitful.Units, x::Unitful.Quantity, e::ExchangeMarket; mode::Int=1)
    u_match = match(r".+?\{(.+)\}", string(Unitful.dimension(u)))
    x_match = match(r".+?\{(.+)\}", string(Unitful.dimension(x)))

    if (u_match !== nothing) && (x_match !== nothing)
        u_asset = deboldify(u_match.captures[1])
        x_asset = deboldify(x_match.captures[1])
        pair = AssetsPair(x_asset, u_asset)
        pairinv = AssetsPair(u_asset, x_asset)
        if mode == 1 && pair in keys(e)
            return Unitful.uconvert(u, e[pair].value * x)
        elseif mode == -1 && pairinv in keys(e)
            return Unitful.uconvert(u, x / e[pairinv].value)
        elseif mode == 2
            for (pair1, rate1) in e
                for (pair2, rate2) in e
                    if pair1.base_asset == x_asset && pair2.quote_asset == u_asset && pair1.quote_asset == pair2.base_asset
                        return Unitful.uconvert(u, rate1.value * rate2.value * x)
                    end
                end
            end
        elseif mode == -2
            for (pair1, rate1) in e
                for (pair2, rate2) in e
                    if pair1.base_asset == u_asset && pair2.quote_asset == x_asset && pair1.quote_asset == pair2.base_asset
                        return Unitful.uconvert(u, x / (rate1.value * rate2.value))
                    end
                end
            end
        end
        throw(ArgumentError(
            "No such rate available in the given market" *
            "for the conversion from $(Unitful.unit(x)) to $u."
            )
        )
    else
        throw(ArgumentError("$u and $x must be proper assets"))
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct DialogApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_dialogOpen(_api::DialogApi, dialog::String, trigger_id::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/dialog.open", ["slackAuth"])
    Swagger.set_param(_ctx.query, "dialog", dialog)  # type String
    Swagger.set_param(_ctx.query, "trigger_id", trigger_id)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Open a dialog with a user
Param: dialog::String (required)
Param: trigger_id::String (required)
Return: Dict{String, Any}
"""
function dialogOpen(_api::DialogApi, dialog::String, trigger_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_dialogOpen(_api, dialog, trigger_id; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function dialogOpen(_api::DialogApi, response_stream::Channel, dialog::String, trigger_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_dialogOpen(_api, dialog, trigger_id; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export dialogOpen


text:
function ground(args...)
    background("white")
    sethue("black")
end

@testset "Javis Viewer" begin
    astar(args...) = star(O, 50)
    acirc(args...) = circle(Point(100, 100), 50)

    vid = Video(500, 500)
    action_list = [
        BackgroundAction(1:100, ground),
        Action(1:100, morph(astar, acirc; action = :fill)),
    ]

    javis(vid, action_list, pathname = "")

    viewer_win, frame_dims, r_slide, tbox, canvas, actions, total_frames, video =
        Javis._javis_viewer(vid, 100, action_list, false)
    visible(viewer_win, false)

    @test get_gtk_property(viewer_win, :title, String) == "Javis Viewer"

    Javis._increment(video, [r_slide, tbox], actions, frame_dims, canvas, total_frames)
    sleep(0.1)
    curr_frame = Reactive.value(r_slide)
    second_frame = Javis.get_javis_frame(video, actions, curr_frame)
    @test Reactive.value(r_slide) == 2

    Javis._decrement(video, [r_slide, tbox], actions, frame_dims, canvas, total_frames)
    sleep(0.1)
    curr_frame = Reactive.value(r_slide)
    first_frame = Javis.get_javis_frame(video, actions, curr_frame)
    @test Reactive.value(r_slide) == 1

    @test first_frame != second_frame

    Javis._decrement(video, [r_slide, tbox], actions, frame_dims, canvas, total_frames)
    sleep(0.1)
    curr_frame = Reactive.value(r_slide)
    last_frame = Javis.get_javis_frame(video, actions, curr_frame)
    @test curr_frame == total_frames

    Javis._increment(video, [r_slide, tbox], actions, frame_dims, canvas, total_frames)
    sleep(0.1)
    curr_frame = Reactive.value(r_slide)
    first_frame = Javis.get_javis_frame(video, actions, curr_frame)
    @test curr_frame == 1

    @test last_frame != first_frame
end


text:
function polynomial_t(p::AbstractPolynomial{<:Number}, timestamp=0)
    @assert nvariables(p) == 1
    d = maxdegree(p)
    coeffs = zeros(d + 1)
    for term in terms(p)
        @assert nvariables(term) == 1  # should be redundant with the nvariables call above, can probably be removed
        coeffs[degree(term) + 1] = coefficient(term)
    end
    polynomial_t(timestamp, length(coeffs), coeffs)
end

function polynomial_matrix_t(p::AbstractMatrix{<:AbstractPolynomial}, timestamp=0)
    msg = polynomial_matrix_t()
    msg.timestamp = timestamp
    msg.rows = size(p, 1)
    msg.cols = size(p, 2)
    msg.polynomials = [[polynomial_t(p[i, j], timestamp) for j in 1:size(p, 2)] for i in 1:size(p, 1)]
    msg
end

polynomial_matrix_t(p::AbstractVector{<:AbstractPolynomial}, timestamp=0) = polynomial_matrix_t(reshape(p, :, 1), timestamp)

function piecewise_polynomial_t(pf::PiecewiseFunction{<:AbstractPolynomial}, args...)
    piecewise_polynomial_t(fill(pf, 1, 1), args...)
end

function piecewise_polynomial_t(pfs::AbstractVecOrMat{<:PiecewiseFunction{<:AbstractPolynomial}}, timestamp=0)
    pf1 = first(pfs)
    b = breaks(pf1)
    num_pieces = length(b) - 1
    for pf in pfs
        if pf !== pf1
            @assert breaks(pf) == b
            @assert length(pieces(pf)) == num_pieces
        end
    end
    msg = piecewise_polynomial_t(
        timestamp,
        length(b),
        b,
        num_pieces,
        Vector{polynomial_matrix_t}(num_pieces)
    )
    for k in 1:num_pieces
        mat = polynomial_matrix_t()
        mat.rows = size(pfs, 1)
        mat.cols = size(pfs, 2)  # returns 1 when pfs is a Vector
        resize!(mat)
        msg.polynomial_matrices[k] = mat
    end
    for I in eachindex(pfs)
        pf = pfs[I]
        for k in 1:num_pieces
            piece = pieces(pf)[k]
            msg.polynomial_matrices[k].polynomials[I] = polynomial_t(piece)
        end
    end
    msg
end


text:
using Statistics

struct Interval{T}
    previous::Union{Event{T}, Nothing}
    main::Event{T}
    next::Event{T}
end
duration(I::Interval) = I.next.time - I.main.time

intervals(ts::EventTS) =  (Interval(i==1 ? nothing : ts[i-1], ts[i], ts[i+1]) for i in 1:(length(ts)-1))
intervals(ts::EventTS, value) = (Interval(i==1 ? nothing : ts[i-1], ts[i], ts[i+1]) for i in 1:(length(ts)-1) if ts[i].val == value)

previous_current_next_value(I::Interval) = (I.previous == nothing ? nothing : I.previous.val, I.main.val, I.next.val)

function interval_duration_stats(ts, val; groupby=previous_current_next_value)
    groups = Dict()
    for I in intervals(ts, val)
        k,v = groupby(I), duration(I)
        if haskey(groups,k)
            push!(groups[k],v)
        else
            groups[k]=[v]
        end
    end
    Dict(k=>(count=length(v), duration_mean=mean(v), duration_std=std(v, corrected=false)) for (k,v) in groups)
end

export interval_duration_stats


text:
# Model Lion & Zebra (predetor vs. prey) over the course of 100 years.

# Using Lotka-Volterra equation, we have:
# dx/dt=x(α−βy)  dy/dt=−y(γ−δx)

start_time = 0;
end_time = 100;

# How much time passes between each successive calculation
time_step = 1/4 # In years
end_step = int(((end_time - start_time)/ time_step));

# The number of zebras when the simulation is started
initial_zebras = 30000;

# The number of lions when the simulation is started
initial_lions = 15;

# The number of zebras that must be killed for a lion to be born. (Zebras / Lion)

zebras_per_lion = 1000;

# The chance a zebra will die when a zebra lion cross paths
fighting_chance = 0.50;

meeting_chance = 0.02;

zebra_growth_rate = 0.20;

lion_death_rate = 0.10;

# Initialization

zebras_over_time = fill(0.0, end_step+1);
lions_over_time = fill(0.0, end_step+1);
model_time = fill(0.0, end_step+1);

zebras = initial_zebras;
lions = initial_lions;

zebras_over_time[1]	= zebras;
lions_over_time[1] = lions;

# Run the model

for sim_step = 1:end_step

	sim_time = start_time + sim_step * time_step;
	model_time[sim_step] = sim_time;

	# First we must calculate our flows (our rates)
	zebra_births = zebras * zebra_growth_rate;
	zebra_deaths  = min(zebras, meeting_chance*fighting_chance*zebras*lions);

	lion_births = 1/zebras_per_lion * zebra_deaths;
	lion_deaths = lions * lion_death_rate;

	# Update the stock levels
	zebras = zebras + zebra_births - zebra_deaths;
	lions = lions + lion_births - lion_deaths

	# Stock values always update in the next time step
	zebras_over_time[sim_step+1] = zebras;
	lions_over_time[sim_step+1] = lions;
end

println("Time, Zebras (Thousands), Lions");

for i = 1:end_step
	print(model_time[i]);
	print(",");
	print(zebras_over_time[i]/1000);
	print(",");
	println(lions_over_time[i]);
end

max_stars_repo_name: neilpanchal/Julia-Algorithms

max_issues_repo_name: neilpanchal/Julia-Algorithms

max_forks_repo_name: neilpanchal/Julia-Algorithms


text:
export flux_density, specific_power_loss
export turns, copper_weight_to_meters, winding_resistance
export volt_seconds_per_turn, volt_seconds
export chan_inductor, leakage_inductance, capacitance

function interpolate_third_point(x1,y1,x2,y2, x3)
  m = (y2-y1)/(x2-x1)
  b = y1-x1*m
  return m*x3+b
end
function find_nearest_spl_frequency_indices(spl::SpecificPowerLossData,f)
  l = length(spl.frequency)
  for i in 2:l # there will never be many of these
    if spl.frequency[i]>f
      return i-1:i
    end
  end
  return l-1:l
end

"""
    specific_power_loss(spl::SpecificPowerLossData, flux_density, frequency)
    specific_power_loss(fp::FerriteProperties, flux_density, frequency)
    specific_power_loss(t::Transformer, flux_density, frequency)
    specific_power_loss(tpa::TransformerPowerAnalysis)

Returns specific power loss.
"""
specific_power_loss

function specific_power_loss(spl::SpecificPowerLossData, flux_density, f)
  # spl = tabulated specific power loss data from graph on datasheet
  # flux_density = magnetic field strength in Tesla
  # f = frequency in Hz
  # retruns specfic power loss (W/m^3) at flux_density (Tesla), f (Hz)
  indexrange = find_nearest_spl_frequency_indices(spl,f)
  f_array = [spl.frequency[i] for i in indexrange]
  pv_array = [spl.mb[i][1]*log10(flux_density)+spl.mb[i][2] for i in indexrange]
  pv = 10^interpolate_third_point(f_array[1],pv_array[1],
                                    f_array[2],pv_array[2], f)
  return pv # specicic power loss (W/m^3) at flux_density, frequency f
end
specific_power_loss(fp::FerriteProperties, flux_density, f) =
  specific_power_loss(fp.spl_hot,flux_density,f)
specific_power_loss(t::Transformer, flux_density::Float64, f::Float64)=
  specific_power_loss(ferrite(t),flux_density,f)
specific_power_loss(tpa::TransformerPowerAnalysis) = tpa.core_specific_power

"""
    flux_density(spl::SpecificPowerLossData, coreloss, frequency)
    flux_density(fp::FerriteProperties, coreloss, frequency)
    flux_density(t::Transformer, coreloss, frequency)
    flux_density(tpa::TransformerPowerAnalysis)

Returns magnetic field strength in Tesla.
"""
flux_density

function flux_density(spl::SpecificPowerLossData, coreloss::Float64, f::Float64)
  # spl = tabulated specific power loss data from graph on datasheet
  # coreloss = specfic power loss at b,f in W/m^3
  # f = frequency in Hz
  # returns magametic field strength in Tesla
  indexrange = find_nearest_spl_frequency_indices(spl,f)
  f_array = [spl.frequency[i] for i in indexrange]
  b_array = [(log10(coreloss)-spl.mb[i][2])/spl.mb[i][1] for i in indexrange]
  b = 10^interpolate_third_point(f_array[1],b_array[1],
                              f_array[2],b_array[2], f)
  return b # magnetic field strength in Tesla
end
flux_density(fp::FerriteProperties, coreloss::Float64, f::Float64) =
  flux_density(fp.spl_hot,coreloss,f)
flux_density(t::Transformer, coreloss::Float64, f::Float64) =
  flux_density(ferriteproperties(t),coreloss,f)

"""
    copper_weight_to_meters(oz)

PCB copper thickness is typicaly given in ounces.  This function multiplies
by 35e-6 to give thickness in meters.
"""
copper_weight_to_meters(oz) = 35e-6*oz

"""
    turns(x :: WindingGeometry)
    turns(x :: Windings)
    turns(x :: Transformer)
    turns(x :: TransformerPowerAnalysis)

Number of turns, Tuple for objects with two windings.
"""
function turns(w::Windings)
  primary_turns = turns(primarywindinggeometry(w))
  primary_turns *= isprimaryseries(w) ? count(isprimary(w)) : 1
  secondary_turns = turns(secondarywindinggeometry(w))
  secondary_turns *= issecondaryseries(w) ? length(isprimary(w)) - count(isprimary(w)) : 1
  (primary_turns, secondary_turns)
end
turns(t::Transformer) = turns(windings(t))
turns(tpa::TransformerPowerAnalysis) = turns(transformer(tpa))

"""
    volt_seconds_per_turn(effective_area, flux_density_pp)
    volt_seconds_per_turn(x ::CoreGeometry, flux_density_pp)
    volt_seconds_per_turn(x ::Windings, flux_density_pp)
    volt_seconds_per_turn(x ::Transformer, flux_density_pp)
    volt_seconds_per_turn(x ::TransformerPowerAnalysis, flux_density_pp)
    volt_seconds_per_turn(x ::TransformerPowerAnalysis)

Volt seconds per turn at `flux_density_pp`.  The first form is just an alias for multiply.
"""
volt_seconds_per_turn(effective_area, flux_density_pp) =
  effective_area * flux_density_pp
volt_seconds_per_turn(cg::CoreGeometry, flux_density_pp) =
  volt_seconds_per_turn(effective_area(cg), flux_density_pp)
volt_seconds_per_turn(w::Windings, flux_density_pp) =
  volt_seconds_per_turn(core(w), flux_density_pp)
volt_seconds_per_turn(t::Transformer, flux_density_pp) =
  volt_seconds_per_turn(windings(t), flux_density_pp)
volt_seconds_per_turn(tpa::TransformerPowerAnalysis, flux_density_pp)=
  volt_seconds_per_turn(transformer(tpa), flux_density_pp)
volt_seconds_per_turn(tpa::TransformerPowerAnalysis) =
  volt_seconds_per_turn(transformer(tpa), flux_density(tpa))
"""
    volt_seconds(x :: Windings, flux_density_pp)
    volt_seconds(x :: Transformer, flux_density_pp)
    volt_seconds(x :: TransformerPowerAnalysis, flux_density_pp)
    volt_seconds(x :: TransformerPowerAnalysis)

Volt seconds at `flux_density_pp`.
"""
volt_seconds(w::Windings, flux_density_pp) =
  volt_seconds_per_turn(core(w), flux_density_pp).*turns(w)
volt_seconds(t::Transformer, flux_density_pp) =
  volt_seconds(windings(t), flux_density_pp)
volt_seconds(tpa::TransformerPowerAnalysis, flux_density_pp) =
  volt_seconds(transformer(tpa), flux_density_pp)
volt_seconds(tpa::TransformerPowerAnalysis) =
  volt_seconds(transformer(tpa), flux_density(tpa))

const μ0 =4π*1e-7
skin_depth(ρ,f)=√(ρ/(π*f*μ0)) # skin depth
function conductivity(ρ_20, temperature_coefficient, temperature)
  ρ_20*(1 + temperature_coefficient*(temperature-20.0))
end

layer_resistance(wg::WindingGeometry, thickness, ρ) =
  ρ*length(wg)/(width(wg)*thickness)
function layer_resistance(w::Windings, i, frequency=0.0, temperature=100.0)
  winding_geometry = isprimary(w)[i] ? primarywindinggeometry(w) : secondarywindinggeometry(w)
  ρ = conductivity(ρ_20(material(stackup(pcb(w)))[i*2-1]), tc(material(stackup(pcb(w)))[i*2-1]), temperature)
  effective_thickness = thickness(stackup(pcb(w)))[i*2-1]
  if frequency>0.0
    δ = skin_depth(ρ, frequency)
    effective_thickness = minimum((sides(w)[i]*δ, effective_thickness))
  end
  layer_resistance(winding_geometry, effective_thickness, ρ)
end
function layer_resistance_tuple(w::Windings, frequency=0.0, temperature=100.0)
  ntuple(i->layer_resistance(w, i, frequency, temperature), length(isprimary(w)))
end
"""
    winding_resistance(x :: Windings, frequency=0.0, temperature=100.0)
    winding_resistance(x :: Transformer, frequency=0.0, temperature=100.0)
    winding_resistance(x :: TransformerPowerAnalysis, frequency=0.0)

Returns tuple `(primary_resistance, secondary_resistance)`
"""
function winding_resistance(w::Windings, frequency=0.0, temperature=100.0)
  primary = 0.0
  secondary = 0.0
  for i in eachindex(isprimary(w))
    x = layer_resistance(w,i,frequency,temperature)
    if isprimary(w)[i]
      primary += isprimaryseries(w) ? x : 1/x
    else
      secondary += issecondaryseries(w) ? x : 1/x
    end
  end
  (isprimaryseries(w) ? primary : 1/primary, issecondaryseries(w) ? secondary : 1/secondary)
end
winding_resistance(t::Transformer, frequency=0.0, temperature=100.0) =
  winding_resistance(windings(t),frequency,temperature)
winding_resistance(tpa::TransformerPowerAnalysis, frequency=0.0) =
  winding_resistance(transformer(tpa), frequency)

"""
    center_frequency(x :: FerriteProperties)
    center_frequency(x :: Transformer)

Middle of the minimum and maximum recomended operating frequencys from datasheet.
"""
center_frequency(fp::FerriteProperties) = (fmin(fp)+fmax(fp))/2.0
center_frequency(t::Transformer) = center_frequency(ferrite(t))

"""
    chan_inductor(ferriteproperties, effective_area, effective_length,ishot=true)
    chan_inductor(transformer, ishot=true)
    chan_inductor(transformer_power_analysis, ishot=true)

Parameters for LTspice [Chan inductor](http://ltwiki.org/?title=The_Chan_model)
to be used for magnetizing inductance.
"""
chan_inductor
function chan_inductor(fp::FerriteProperties,
                       effective_area,effective_length,ishot::Bool=true)
  bh = ishot ? bh_hot(fp) : bh_room(fp)
  ChanInductor(hc(bh_hot(fp)), bs(bh_hot(fp)), br(bh_hot(fp)), effective_area, effective_length, 0.0, 1.0)
end
chan_inductor(t::Transformer, ishot::Bool=true) =
  chan_inductor(ferrite(t),
                effective_area(t),
                effective_length(t),
                ishot)
chan_inductor(tpa::TransformerPowerAnalysis, ishot::Bool=true) =
  chan_inductor(transformer(tpa), ishot)

"""
    winding_area(x)

Area the windings occupy.  For internal use.
"""
function winding_area(cg::CoreGeometry)
  π*(half_center_width(cg)+winding_aperture(cg))^2-π*half_center_width(cg)^2+
  2.0*winding_aperture(cg)*center_length(cg)
end
winding_area(w::Windings) = winding_area(core(w))
winding_area(t::Transformer) = winding_area(windings(t))
winding_area(tpa::TransformerPowerAnalysis) = winding_area(transformer(tpa))
"""
    winding_breadth_volume(windings)

Tuple (breadth, volume).  For internal use calculating leakage inductance.
"""
function winding_breadth_volume(w::Windings)
  breadth = 0.0
  volume = 0.0
  wa = winding_area(core(w))
  for i in 2:length(isprimary(w))
    if isprimary(w)[i] != isprimary(w)[i-1]
      breadth += winding_aperture(w)
      volume += thickness(stackup(pcb(w)))[2*(i-1)]*wa
    end
  end
  (breadth, volume)
end
"""
    leakage_inductance(x)

Leakage inductance from volume between primary and secondary, referenced to
primary.  Fringing fields are not taken into account.
"""
function leakage_inductance(w::Windings)
  (breadth, volume) = winding_breadth_volume(w)
  μ0*volume*(turns(w)[1]/breadth)^2
end
leakage_inductance(t::Transformer) = leakage_inductance(windings(t))
leakage_inductance(tpa::TransformerPowerAnalysis) = leakage_inductance(transformer(tpa))

"""
    capacitance(x)

Capacitance between primary and secondary.
"""
function capacitance(w::Windings)
  a = winding_area(w)
  c = 0.0
  thick = thickness(stackup(pcb(w)))
  mat = material(stackup(pcb(w)))
  for i in 2:length(isprimary(w))
    if isprimary(w)[i] != isprimary(w)[i-1]
      c += ϵ(mat[2*(i-1)])*a/thick[2*(i-1)]
    end
  end
  return c
end
capacitance(t::Transformer) = capacitance(windings(t))
capacitance(tpa::TransformerPowerAnalysis) = capacitance(transformer(tpa))

"""
    pcb_thickness(x)

Overall thickness of PCB.
"""
pcb_thickness(x::Stackup) = sum(thickness(x))
pcb_thickness(x::PCB_Specification) = pcb_thickness(stackup(x))
pcb_thickness(x::Windings) = pcb_thickness(pcb(x))
pcb_thickness(x::Transformer) = pcb_thickness(windings(x))
pcb_thickness(x::TransformerPowerAnalysis) = pcb_thickness(transformer(x))


text:
### Solver options ###
opts = TrajectoryOptimization.SolverOptions()
opts.square_root = false
opts.verbose = false
opts.cost_tolerance_intermediate = 1e-6
opts.constraint_tolerance = 1e-6
opts.cost_tolerance = 1e-6
opts.constraint_decrease_ratio = 0.25
opts.penalty_scaling = 10.0
opts.bp_reg_initial = 0.0
opts.outer_loop_update_type = :default
opts.use_static = false
opts.resolve_feasible = false
######################

### Set up model, objective, solver ###
# Model, objective (unconstrained)
model_pendulum, obj_uncon_pendulum = TrajectoryOptimization.Dynamics.pendulum!
model_dubins, obj_uncon_dubins = TrajectoryOptimization.Dynamics.dubinscar!
model_cartpole, obj_uncon_cartpole = TrajectoryOptimization.Dynamics.cartpole_udp
model_quadrotor, obj_uncon_quadrotor = TrajectoryOptimization.Dynamics.quadrotor

## Constraints

# pendulum
u_min_pendulum = -2
u_max_pendulum = 2
x_min_pendulum = [-20;-20]
x_max_pendulum = [20; 20]

# dubins car
u_min_dubins = [-1; -1]
u_max_dubins = [1; 1]
x_min_dubins = [0; -100; -100]
x_max_dubins = [1.0; 100; 100]

# cartpole
u_min_cartpole = -20
u_max_cartpole = 40
x_min_cartpole = [-10; -1000; -1000; -1000]
x_max_cartpole = [10; 1000; 1000; 1000]

# quadrotor
u_min = -100.0
u_max = 100.0
# -constraint that quaternion should be unit
function cE(c,x,u)
    c[1] = x[4]^2 + x[5]^2 + x[6]^2 + x[7]^2 - 1.0
    return nothing
end

# -Constrained objective
obj_con_pendulum = ConstrainedObjective(obj_uncon_pendulum, u_min=u_min_pendulum, u_max=u_max_pendulum, x_min=x_min_pendulum, x_max=x_max_pendulum)
obj_con_dubins = ConstrainedObjective(obj_uncon_dubins, u_min=u_min_dubins, u_max=u_max_dubins, x_min=x_min_dubins, x_max=x_max_dubins)
obj_con_cartpole = ConstrainedObjective(obj_uncon_cartpole, u_min=u_min_cartpole, u_max=u_max_cartpole, x_min=x_min_cartpole, x_max=x_max_cartpole)
obj_con_quadrotor = TrajectoryOptimization.ConstrainedObjective(obj_uncon_quadrotor, u_min=u_min, u_max=u_max)#, cE=cE)

# System selection
model = model_pendulum
obj = obj_uncon_pendulum
u_max = u_max_pendulum
u_min = u_min_pendulum

# Solver
intergrator_zoh = :rk3
intergrator_foh = :rk3_foh

# dt = 0.1
N = 100
solver_zoh = Solver(model,obj,integration=intergrator_zoh,N=N,opts=opts)
solver_foh = Solver(model,obj,integration=intergrator_foh,N=N,opts=opts)

# -Initial state and control trajectories
X_interp = line_trajectory(solver_zoh.obj.x0,solver_zoh.obj.xf,solver_zoh.N)
U = zeros(solver_foh.model.m,solver_foh.N)
#######################################

### Solve ###
@time results_zoh, stats_zoh = solve(solver_zoh,U)
@time results_foh, stats_foh = solve(solver_foh,U)

# plot(log.(stats_foh["cost"]),title="Unconstrained Pendulum",ylabel="log(cost)",xlabel="iteration",label=["foh","zoh"])
# plot!(log.(stats_zoh["cost"]))

### Results ###
println("Final state (zoh)-> res: $(results_zoh.X[end]), goal: $(solver_zoh.obj.xf)\n Iterations: $(stats_zoh["iterations"])\n Outer loop iterations: $(stats_zoh["major iterations"])\n Max violation: $(stats_zoh["c_max"][end])\n Max μ: $(maximum([to_array(results_zoh.μ)[:]; results_zoh.μN[:]]))\n Max abs(λ): $(maximum(abs.([to_array(results_zoh.λ)[:]; results_zoh.λN[:]])))\n")
println("Final state (foh)-> res: $(results_foh.X[end]), goal: $(solver_foh.obj.xf)\n Iterations: $(stats_foh["iterations"])\n Outer loop iterations: $(stats_foh["major iterations"])\n Max violation: $(stats_foh["c_max"][end])\n Max μ: $(maximum([to_array(results_foh.μ)[:]; results_foh.μN[:]]))\n Max abs(λ): $(maximum(abs.([to_array(results_foh.λ)[:]; results_foh.λN[:]])))\n")

# # Controllers
controller_zoh = generate_controller(to_array(results_zoh.X),to_array(results_zoh.U),to_array(results_zoh.K),solver_zoh.N,solver_zoh.dt,:zoh,u_min,u_max)
controller_foh = generate_controller(to_array(results_foh.X),to_array(results_foh.U),to_array(results_foh.K),solver_foh.N,solver_foh.dt,:foh,u_min,u_max)

# Simulate
dt_sim = 0.001
tf = obj.tf
t_solve = collect(get_time(solver_foh))
t_sim = collect(0:dt_sim:tf)
x0 = obj.x0
f = model.f
n = model.n
m = model.m

using ODE
integrator_sim = :ode45 # note that ode45 can be used with 'simulate_controller' but the time indicies must be altered when plotting

X_zoh_sim, U_zoh_sim = simulate_controller(f,integrator_sim,controller_zoh,n,m,dt_sim,x0,tf,u_min,u_max)
X_foh_sim, U_foh_sim = simulate_controller(f,integrator_sim,controller_foh,n,m,dt_sim,x0,tf,u_min,u_max)

# RMS error
xf_rms_zoh = sqrt(mean((results_zoh.X[end] - X_zoh_sim[:,end]).^2))
xf_rms_foh = sqrt(mean((results_foh.X[end] - X_foh_sim[:,end]).^2))

plot(range(0,stop=solver_zoh.obj.tf,length=size(X_zoh_sim,2)),X_zoh_sim')
plot!(range(0,stop=solver_foh.obj.tf,length=size(X_foh_sim,2)),X_foh_sim')
length(t_sim)

# ZOH plotting
X_zoh_interp, U_zoh_interp = interpolate_trajectory(solver_zoh, to_array(results_zoh.X), to_array(results_zoh.U), t_sim)
p_zoh = plot(t_sim,X_zoh_interp',color="black",title="zoh (N = $(solver_zoh.N),tf=$(solver_zoh.obj.tf)) sim @ $(convert(Int64,1/dt_sim))Hz",labels="")
p_zoh = plot!(range(0,stop=solver_zoh.obj.tf,length=size(X_zoh_sim,2)),X_zoh_sim',color="green",labels="",xlabel="xf RMS error: $(round(xf_rms_zoh,digits=5))")
display(p_zoh)
# savefig("knotpointtest_zoh.png")

p_zoh = plot(t_sim,U_zoh_interp',color="black",title="zoh (N = $(solver_zoh.N),tf=$(solver_zoh.obj.tf)) sim @ $(convert(Int64,1/dt_sim))Hz",labels="")
p_zoh = plot!(range(0,stop=solver_zoh.obj.tf,length=size(U_zoh_sim,2)),U_zoh_sim',color="green",labels="",xlabel="xf RMS error: $(round(xf_rms_zoh,digits=5))")
display(p_zoh)
# savefig("knotpointtest_zoh_control.png")

# FOH plotting
X_foh_interp, U_foh_interp = interpolate_trajectory(solver_foh, to_array(results_foh.X), to_array(results_foh.U), t_sim)
p_foh = plot(t_sim,X_foh_interp',color="black",title="foh (N = $(solver_foh.N),tf=$(solver_foh.obj.tf)) sim @ $(convert(Int64,1/dt_sim))Hz",labels="")
p_foh = plot!(range(0,stop=solver_foh.obj.tf,length=size(X_foh_sim,2)),X_foh_sim',color="purple",labels="",xlabel="xf RMS error: $(round(xf_rms_foh,digits=5))")
display(p_foh)
# savefig("knotpointtest_foh.png")

p_foh = plot(t_sim, U_foh_interp',color="black",title="foh (N = $(solver_foh.N),tf=$(solver_foh.obj.tf)) sim @ $(convert(Int64,1/dt_sim))Hz",labels="")
p_foh = plot!(range(0,stop=solver_foh.obj.tf,length=size(U_foh_sim,2)),U_foh_sim',color="purple",labels="",xlabel="xf RMS error: $(round(xf_rms_foh,digits=5))")
display(p_foh)
# savefig("knotpointtest_foh_control.png")


text:
using LinearAlgebra
using Random
using Polynomials
using SemiLagrangian:
    get_kl_ku, getbspline, LuSpline, B_SplineLU, sol, decal, get_order, interpolate!
import SemiLagrangian: sol, sol!

Random.seed!(5431221)
# create a band or circular matrix from a vector of non-zero data
function topl(n, t, iscirc = true)
    res = zeros(Rational{BigInt}, n, n)
    kl, ku = get_kl_ku(size(t, 1))
    for i = 1:n
        for (j, v) in enumerate(t)
            ind = i + j - kl - 1
            if 1 <= ind <= n
                res[i, ind] = v
            elseif iscirc
                if ind < 1
                    ind += n
                else
                    ind -= n
                end
                res[i, ind] = v
            end
        end
    end
    return res
end
# inplace lu decomposition
function decLU(A::Matrix{T}) where {T}
    n = size(A, 1)
    for k = 1:n
        pivot = A[k, k]
        for i = k+1:n
            A[i, k] /= pivot
        end
        for i = k+1:n, j = k+1:n
            # if A[i,k] != 0 && A[k,j] != 0
            #     println("trace avant k=$k i=$i j=$j A[i,j]=$(A[i,j]) A[i,k]=$(A[i,k]) A[k,j]=$(A[k,j])")
            # end
            A[i, j] -= A[i, k] * A[k, j]
            # if A[i,k] != 0 && A[k,j] != 0
            #     println("trace après k=$k i=$i j=$j A[i,j]=$(A[i,j])")
            # end
        end
    end
    return A
end
function getLU(A::Matrix{T}) where {T}
    n = size(A, 1)
    L = zeros(T, n, n)
    U = zeros(T, n, n)
    for i = 1:n
        L[i, i] = 1
        for j = 1:i-1
            L[i, j] = A[i, j]
        end
        for j = i:n
            U[i, j] = A[i, j]
        end
    end
    return L, U
end

# function testTools()

function sol(A::Matrix{T}, Y::Vector{T}) where {T}
    L, U = getLU(A)
    n = size(A, 1)
    Y1 = zeros(T, n)
    for i = 1:n
        Y1[i] = Y[i] - sum(Y1[1:i-1] .* A[i, 1:i-1])
    end
    #   @assert Y1 == (L^(-1))*Y "Erreur 1"
    #   @assert isapprox(Y1, (L^(-1))*Y, atol=1e-60) "Erreur 1"
    X = zeros(T, n)
    for i = n:-1:1
        X[i] = (Y1[i] - sum(X[i+1:n] .* A[i, i+1:n])) / A[i, i]
    end
    #    @assert X == (U^(-1))*Y1 "Erreur 2"
    #    @assert isapprox(X,(U^(-1))*Y1,atol=1e-60) "Erreur 2"
    #    @assert X == ((L*U)^(-1))*Y "Erreur3"
    #    @assert isapprox(X, ((L*U)^(-1))*Y, atol=1e-60) "Erreur3"
    return X, Y1
end

function test_decLU(n)
    dividende = 10000
    A = big.((rand(Int, n, n) .% dividende) .// dividende)
    B = copy(A)
    decLU(A)
    L, U = getLU(A)
    @test L * U == B
    b = big.((rand(Int, n) .% dividende) .// dividende)
    x, _ = sol(A, b)
    @test B * x == b
end

function test_sol(T, n,order)
    A = topl(n, getbspline(order, 0).(1:order))
    bsp = B_SplineLU(order,n,Float64)
    b = rand(n)
#    y = zeros(n)
    y = sol(bsp, b)
    diff = norm(A*y-b)
    @show diff
    return diff
end

function test_splu(n, order, iscirc, isLU; type = Rational{BigInt}, tol = NaN, perf = false)
    t = convert.(type, getbspline(big(order), 0).(1:order))
    A = convert.(type, topl(n, t, iscirc))
    Acp = copy(A)
    ku = div(order, 2)
    kl = order - 1 - ku
    if isLU
        perf && @time decLU(A)
        !perf && decLU(A)
    end
    spA = LuSpline(A, ku, kl; iscirc = iscirc, isLU = isLU)
    perf && @time spB = LuSpline(n, t; iscirc = iscirc, isLU = isLU)
    !perf && (spB = LuSpline(n, t; iscirc = iscirc, isLU = isLU))
    if isnan(tol)
        @test spA == spB
    end

    dividende = 10000
    b = big.((rand(Int, n) .% dividende) .// dividende)
    b = convert.(type, b)
    if isLU
        x, y = sol(A, b)
        x2, y2 = sol(spB, b)

        if isnan(tol)
            #            @show norm(x-x2,Inf)
            @test x == x2
        else
            println("norm=$(norm(x-x2))")
            res = Acp * x
            println("normres=$(norm(res-b))")
            @test isapprox(res, b, atol = tol)
            @test isapprox(x, x2, atol = tol)
        end
    end
end

function test_perf(n, order, iscirc)
    t = getbspline(big(order), 0).(1:order)
    A = topl(n, t, iscirc)
    ku = div(order, 2)
    kl = order - 1 - ku
    @time decLU(A)

    spA = LuSpline(A, ku, kl; iscirc = iscirc, isLU = true)
    @time spB = LuSpline(n, t; iscirc = iscirc, isLU = true)

    dividende = 10000
    b = big.((rand(Int, n) .% dividende) .// dividende)
    x, y = sol(A, b)
    x2, y2 = sol(spB, b)
    @test y == y2
    @test x == x2

end

function test_interface()
    getpar(bsp::B_SplineLU{T,tedge,order}) where {T,tedge,order} =
        (T, tedge, order)
    for i = 3:2:29
        bsp = B_SplineLU(i, 104, big"0.")
        (_1, _2, order) = getpar(bsp)
        @test i == get_order(bsp)
        @test i == order
        @test "B_SplineLU{BigFloat,CircEdge,$order}" == replace("$bsp", " " => "")
    end
end

@testset "test decLU that is a tool for test" begin
    test_decLU(30)
end

@testset "test B_SplineLU sol" begin
    @test test_sol(Float64, 128, 9) < 1e-10
end

@testset "test LuSpline" begin
    Random.seed!(1899221)
    test_splu(30, 9, true, false)
    test_splu(30, 10, true, false)
    test_splu(30, 9, false, false)
    test_splu(30, 10, false, false)
    test_splu(30, 9, false, true)
    test_splu(30, 10, false, true)
    test_splu(30, 9, true, true)
    test_splu(31, 10, true, true)
    test_splu(30, 9, true, true, type = BigFloat, tol = 1e-70)
    test_splu(31, 10, true, true, type = BigFloat, tol = 1e-70)
    test_splu(30, 11, false, true, type = BigFloat, tol = 1e-70)
    test_splu(31, 10, false, true, type = BigFloat, tol = 1e-70)
    test_splu(30, 3, true, false)
    test_splu(30, 4, true, false)
    test_splu(30, 3, false, false)
    test_splu(30, 4, false, false)
    test_splu(30, 3, false, true)
    test_splu(30, 4, false, true)
    test_splu(30, 3, true, true)
    test_splu(31, 4, true, true)
    test_splu(30, 3, true, true, type = BigFloat, tol = 1e-70)
    test_splu(31, 4, true, true, type = BigFloat, tol = 1e-70)
    test_splu(30, 5, false, true, type = BigFloat, tol = 1e-70)
    test_splu(31, 4, false, true, type = BigFloat, tol = 1e-70)
    test_interface()
end

@testset "test perf" begin
    test_splu(1000, 9, true, true, type = Float64, tol = 1e-12, perf = true)
end