text:
using JuLIP.FIO: load_dict
using JuLIP: chemical_symbol, atomic_number
using JuLIP.Potentials: WrappedAnalyticFunction, F64fun, @D

import JuLIP.Potentials: evaluate!, evaluate_d!

export EMT

"""
`EMT`: a re-implementation of the `EMT` calculator (a variant of EAM) of ASE
in Julia, largely for fun and comparison with Python, but also to demonstrate
how to implement a multi-component calculator in JuLIP
"""
struct EMT <: SitePotential
   pair::Vector{WrappedPairPotential}
   Cpair::Vector{Float64}
   rho::Vector{WrappedPairPotential}
   embed::Vector{WrappedAnalyticFunction}
   z2ind::Dict{AtomicNumber, Int}
   Z::Vector{AtomicNumber}
   rc::Float64
end

cutoff(calc::EMT) = calc.rc + 0.5

# ========================== Initialisation  ============================

function _load_emt()
   D = load_dict( joinpath(@__DIR__(), "..", "data", "emt.json") )
   return D["params"], D["acut"], D["rc"], D["beta"]
end

EMT(Z::Symbol) = EMT([Z])
EMT(Z::AbstractVector{AtomicNumber}) = EMT(chemical_symbol.(Z))

function EMT(symbols = nothing)
   # load all the parameters
   params, acut, rc, β = _load_emt()
   if symbols == nothing
      symbols = Symbol.(collect(keys(params)))
   end
   Z = atomic_number.(symbols)
   rcplus = rc +  0.5   # the actual cutoff
   # get unique Z numbers and use them to allocate the EMT calculator
   Z = unique(Z)
   nZ = length(Z)
   emt = EMT(Vector{WrappedPairPotential}(undef, nZ),  # pair
             Vector{Float64}(undef, nZ),   # Cpair
             Vector{WrappedPairPotential}(undef, nZ),  # rho
             Vector{WrappedAnalyticFunction}(undef, nZ), # embed
             Dict{AtomicNumber,Int}(),
             Z,
             rc)
   # loop through unique symbols/Z and for each compute the relevant potentials
   for (i, z) in enumerate(Z)
      # store Z -> i maps
      emt.z2ind[z] = i
      sym = chemical_symbol(z)

      p = params[String(sym)]
      # cut-off: this is a FD function; this is AWFUL! replace with spline???
      # Θ = 1.0 / (1.0 + exp($acut * (r - $rc) ))
      # pair potential
      n0, κ, s0 = p["n0"], p["κ"], p["s0"]
      emt.pair[i] = WrappedPairPotential(
            @analytic(r -> n0 * exp( -κ * (r / β - s0) ) * θ,
                      θ = 1.0 / (1.0 + exp(acut * (r - rc) )) ), rcplus)
      emt.Cpair[i] = p["Cpair"]
      # radial electron density function
      η2, γ1 = p["η2"], p["γ1"]
      emt.rho[i] = WrappedPairPotential(
            @analytic(
            r -> n0 * exp( -η2 * (r - (β*s0)) ) * θ,
            θ = 1.0 / (1.0 + exp(acut * (r - rc) )) ), rcplus)
      # embedding function
      E0, V0, λ = p["E0"], p["V0"], p["λ"]
      emt.embed[i] = F64fun( @analytic(
         ρ̄-> E0 * ((1.0 + λ * DS) * exp(-λ * DS) - 1.0) + 6.0 * V0 * exp(-κ * DS),
         DS = - log( (1.0/(12.0*γ1*n0)) * ρ̄ ) / (β * η2) ) )
   end
   return emt
end

# ========================== Main Functionality ============================

function evaluate!(tmp, emt::EMT, Rs, Zs, z0)
   ρ̄ = 0.0
   Es = 0.0
   i0 = emt.z2ind[z0]
   for (R, Z) in zip(Rs, Zs)
      i = emt.z2ind[Z]
      r = norm(R)
      Es += emt.Cpair[i0] * emt.pair[i](r)
      ρ̄ +=  emt.rho[i](r)
   end
   Es += emt.embed[i0](ρ̄)
   return Es
end

function evaluate_d!(dEs, tmp, emt::EMT, Rs, Zs, z0)
   ρ̄ = 0.0
   i0 = emt.z2ind[z0]
   for (R, Z) in zip(Rs, Zs)
      i = emt.z2ind[Z]
      r = norm(R)
      ρ̄ += emt.rho[i](r)
   end
   dF = @D emt.embed[i0](ρ̄)
   for (j, (R, Z)) in enumerate(zip(Rs, Zs))
      i = emt.z2ind[Z]
      r = norm(R)
      R̂ = R/r
      dpair = emt.Cpair[i0] * (@D emt.pair[i](r))
      drho = @D emt.rho[i](r)
      dEs[j] = (dpair + dF * drho) * R̂
   end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 26aaebd9-aaa1-4b07-9075-1e20c0b7090a
# hideall

using Genie, Genie.Router, Genie.Requests, Genie.Renderer.Json

# ╔═╡ 1b4522b2-39b7-11ec-3946-e5aaf1cc1f01
md"""
# Force compiling Genie routes

Because Julia JIT-compiles the code, upon starting up the application, at the time of the first to a certain route, Julia will need to compile the function(s) which are responsible for handling the request. This means that the first request will potentially take many seconds to resolve, which might not be acceptable in a production environment. For such cases, once can use Genie itself to visit each defined route and trigger the compilation.

Here is a sample script which defines two `GET` routes and then automatically triggers them upon starting up the application:
"""

# ╔═╡ 309f98f2-8056-41a5-ad4d-a6e8ec093760
md"""
If you are using Julia REPL instead of Pluto notebooks. Include the below code:

```julia
import Pkg
Pkg.add("Genie")
```
"""

# ╔═╡ 32e73ac8-26e1-4908-b37a-e2a75e2102ec
# hideall

md"""
```julia
using Genie, Genie.Router, Genie.Requests, Genie.Renderer.Json

route("/foo") do
  json(:foo => "Foo")
end

route("/bar") do
  json(:bar => "Bar")
end

function force_compile()
  sleep(5)

  for (name, r) in Router.named_routes()
    Genie.Requests.HTTP.request(r.method, "http://localhost:8000" * tolink(name))
  end
end

@async force_compile()

up(async = false)
```
"""

# ╔═╡ cd5198ce-5792-42e3-92b7-d16f206f936b
# hideall

route("/foo") do
  json(:foo => "Foo")
end;

# ╔═╡ 227dec9e-c836-40be-9811-50da04bacc8d
# hideall

route("/bar") do
  json(:bar => "Bar")
end;

# ╔═╡ a9e41693-0178-4cac-88b9-70abc44b8fca
# hideall

function force_compile()
  sleep(5)

  for (name, r) in Router.named_routes()
    Genie.Requests.HTTP.request(r.method, "http://localhost:8000" * tolink(name))
  end
end;

# ╔═╡ 6ceb043b-c6c8-45cb-a257-289dfb37d276
# hideall

@async force_compile();

# ╔═╡ 936d966c-b61c-4e48-a4f5-6ae6f9c7e880
md"""
In the snippet we define two routes. Then we add a function, `force_compile` which iterates over the routes and hits them through the web server. We then invoke the function, which will be executed with a 5 seconds delay, enough to allow the web server to start up.
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Genie = "c43c736e-a2d1-11e8-161f-af95117fbd1e"

[compat]
Genie = "~3.3.0"
"""

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

[[ArgParse]]
deps = ["Logging", "TextWrap"]
git-tree-sha1 = "3102bce13da501c9104df33549f511cd25264d7d"
uuid = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
version = "1.1.4"

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

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

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

[[CSTParser]]
deps = ["Tokenize"]
git-tree-sha1 = "b2667530e42347b10c10ba6623cfebc09ac5c7b6"
uuid = "00ebfdb7-1f24-5e51-bd34-a7502290713f"
version = "3.2.4"

[[CodeTracking]]
deps = ["InteractiveUtils", "UUIDs"]
git-tree-sha1 = "9aa8a5ebb6b5bf469a7e0e2b5202cf6f8c291104"
uuid = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
version = "1.0.6"

[[CommonMark]]
deps = ["Crayons", "JSON", "URIs"]
git-tree-sha1 = "393ac9df4eb085c2ab12005fc496dae2e1da344e"
uuid = "a80b9123-70ca-4bc0-993e-6e3bcb318db6"
version = "0.8.3"

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

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

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

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

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

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

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

[[EzXML]]
deps = ["Printf", "XML2_jll"]
git-tree-sha1 = "0fa3b52a04a4e210aeb1626def9c90df3ae65268"
uuid = "8f5d6c58-4d21-5cfd-889c-e3ad7ee6a615"
version = "1.1.0"

[[FilePathsBase]]
deps = ["Dates", "Mmap", "Printf", "Test", "UUIDs"]
git-tree-sha1 = "d962b5a47b6d191dbcd8ae0db841bc70a05a3f5b"
uuid = "48062228-2e41-5def-b9a4-89aafe57970f"
version = "0.9.13"

[[FileWatching]]
uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"

[[GMP_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "781609d7-10c4-51f6-84f2-b8444358ff6d"

[[Genie]]
deps = ["ArgParse", "Dates", "Distributed", "EzXML", "FilePathsBase", "HTTP", "HttpCommon", "Inflector", "JSON3", "JuliaFormatter", "Logging", "Markdown", "MbedTLS", "Millboard", "Nettle", "OrderedCollections", "Pkg", "REPL", "Reexport", "Revise", "SHA", "Serialization", "Sockets", "UUIDs", "Unicode", "YAML"]
git-tree-sha1 = "4f5526913b239ece648a73d4445c7978f1251935"
uuid = "c43c736e-a2d1-11e8-161f-af95117fbd1e"
version = "3.3.0"

[[HTTP]]
deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
git-tree-sha1 = "14eece7a3308b4d8be910e265c724a6ba51a9798"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.9.16"

[[HttpCommon]]
deps = ["Dates", "Nullables", "Test", "URIParser"]
git-tree-sha1 = "46313284237aa6ca67a6bce6d6fbd323d19cff59"
uuid = "77172c1b-203f-54ac-aa54-3f1198fe9f90"
version = "0.5.0"

[[Inflector]]
deps = ["Unicode"]
git-tree-sha1 = "8555b54ddf27806b070ce1d1cf623e1feb13750c"
uuid = "6d011eab-0732-4556-8808-e463c76bf3b6"
version = "1.0.1"

[[IniFile]]
deps = ["Test"]
git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.0"

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

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

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

[[JSON3]]
deps = ["Dates", "Mmap", "Parsers", "StructTypes", "UUIDs"]
git-tree-sha1 = "7d58534ffb62cd947950b3aa9b993e63307a6125"
uuid = "0f8b85d8-7281-11e9-16c2-39a750bddbf1"
version = "1.9.2"

[[JuliaFormatter]]
deps = ["CSTParser", "CommonMark", "DataStructures", "Pkg", "Tokenize"]
git-tree-sha1 = "10c95cebcfa37c1f510a726c90886db4745e1238"
uuid = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
version = "0.15.11"

[[JuliaInterpreter]]
deps = ["CodeTracking", "InteractiveUtils", "Random", "UUIDs"]
git-tree-sha1 = "e273807f38074f033d94207a201e6e827d8417db"
uuid = "aa1ae85d-cabe-5617-a682-6adf51b2e16a"
version = "0.8.21"

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

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

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

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

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

[[Libiconv_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "42b62845d70a619f063a7da093d995ec8e15e778"
uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
version = "1.16.1+1"

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

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

[[LoweredCodeUtils]]
deps = ["JuliaInterpreter"]
git-tree-sha1 = "491a883c4fef1103077a7f648961adbf9c8dd933"
uuid = "6f1432cf-f94c-5a45-995e-cdbf5db27b0b"
version = "2.1.2"

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

[[MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

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

[[Millboard]]
git-tree-sha1 = "ea6a5b7e56e76d8051023faaa11d91d1d881dac3"
uuid = "39ec1447-df44-5f4c-beaa-866f30b4d3b2"
version = "0.2.5"

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

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

[[Nettle]]
deps = ["Libdl", "Nettle_jll"]
git-tree-sha1 = "a68340b9edfd98d0ed96aee8137cb716ea3b6dea"
uuid = "49dea1ee-f6fa-5aa6-9a11-8816cee7d4b9"
version = "0.5.1"

[[Nettle_jll]]
deps = ["Artifacts", "GMP_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "eca63e3847dad608cfa6a3329b95ef674c7160b4"
uuid = "4c82536e-c426-54e4-b420-14f461c4ed8b"
version = "3.7.2+0"

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

[[Nullables]]
git-tree-sha1 = "8f87854cc8f3685a60689d8edecaa29d2251979b"
uuid = "4d1e1d77-625e-5b40-9113-a560ec7a8ecd"
version = "1.0.0"

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

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "d911b6a12ba974dabe2291c6d450094a7226b372"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.1.1"

[[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"

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

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

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

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

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

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

[[Revise]]
deps = ["CodeTracking", "Distributed", "FileWatching", "JuliaInterpreter", "LibGit2", "LoweredCodeUtils", "OrderedCollections", "Pkg", "REPL", "Requires", "UUIDs", "Unicode"]
git-tree-sha1 = "41deb3df28ecf75307b6e492a738821b031f8425"
uuid = "295af30f-e4ad-537b-8983-00126c2a3abe"
version = "3.1.20"

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

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

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

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

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

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

[[StringEncodings]]
deps = ["Libiconv_jll"]
git-tree-sha1 = "50ccd5ddb00d19392577902f0079267a72c5ab04"
uuid = "69024149-9ee7-55f6-a4c4-859efe599b68"
version = "0.3.5"

[[StructTypes]]
deps = ["Dates", "UUIDs"]
git-tree-sha1 = "d24a825a95a6d98c385001212dc9020d609f2d4f"
uuid = "856f2bd8-1eba-4b0a-8007-ebc267875bd4"
version = "1.8.1"

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

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

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

[[TextWrap]]
git-tree-sha1 = "9250ef9b01b66667380cf3275b3f7488d0e25faf"
uuid = "b718987f-49a8-5099-9789-dcd902bef87d"
version = "1.0.1"

[[Tokenize]]
git-tree-sha1 = "0952c9cee34988092d73a5708780b3917166a0dd"
uuid = "0796e94c-ce3b-5d07-9a54-7f471281c624"
version = "0.5.21"

[[URIParser]]
deps = ["Unicode"]
git-tree-sha1 = "53a9f49546b8d2dd2e688d216421d050c9a31d0d"
uuid = "30578b45-9adc-5946-b283-645ec420af67"
version = "0.4.1"

[[URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

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

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

[[XML2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "1acf5bdf07aa0907e0a37d3718bb88d4b687b74a"
uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
version = "2.9.12+0"

[[YAML]]
deps = ["Base64", "Dates", "Printf", "StringEncodings"]
git-tree-sha1 = "3c6e8b9f5cdaaa21340f841653942e1a6b6561e5"
uuid = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
version = "0.4.7"

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

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

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

# ╔═╡ Cell order:
# ╟─1b4522b2-39b7-11ec-3946-e5aaf1cc1f01
# ╟─309f98f2-8056-41a5-ad4d-a6e8ec093760
# ╟─32e73ac8-26e1-4908-b37a-e2a75e2102ec
# ╠═26aaebd9-aaa1-4b07-9075-1e20c0b7090a
# ╠═cd5198ce-5792-42e3-92b7-d16f206f936b
# ╠═227dec9e-c836-40be-9811-50da04bacc8d
# ╠═a9e41693-0178-4cac-88b9-70abc44b8fca
# ╠═6ceb043b-c6c8-45cb-a257-289dfb37d276
# ╟─936d966c-b61c-4e48-a4f5-6ae6f9c7e880
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

ext: jl


text:
"""
FinEtoolsDeforLinear (C) 2017-2019, Petr Krysl

Finite Element tools.  Julia implementation  of the finite element method
for continuum mechanics. Package for linear stress analysis problems.
"""
module FinEtoolsDeforLinear

__precompile__(true)

include("allmodules.jl")

# Exports follow:

###########################################################################
# Linear deformation functionality
###########################################################################
using .DeforModelRedModule: AbstractDeforModelRed, DeforModelRed1D, DeforModelRed1DStress, DeforModelRed1DStrain, DeforModelRed2DStrain, DeforModelRed2DStress, DeforModelRed2DAxisymm, DeforModelRed3D, nstressstrain, nthermstrain, stresscomponentmap, Blmat!, divmat, vgradmat
# Exported: types  for model reduction in stress analysis
export AbstractDeforModelRed, DeforModelRed1D, DeforModelRed1DStress, DeforModelRed1DStrain, DeforModelRed2DStrain, DeforModelRed2DStress, DeforModelRed2DAxisymm, DeforModelRed3D
# Exported: num stresses/strains,  number of thermal strains, and map of  the numbering of stress components
export nstressstrain, nthermstrain, stresscomponentmap
# Exported: strain-displacement matrix and divergence matrix for all model-reduction types
export Blmat!, divmat, vgradmat

using .MatDeforModule: AbstractMatDefor
# using .MatDeforModule: AbstractMatDefor, strain2x2tto3v!, strain3vto2x2t!, strain3x3tto6v!, strain6vto3x3t!, strain9vto3x3t!, strain3x3tto9v!, strain9vto6v!, strain6vto9v!, stress2x2to3v!,  stress3vto2x2t!, stress3vto3x3t!, stress4vto3x3t!, stress6vto3x3t!, stress3x3tto6v!, stress9vto6v!,  stress6vto9v!
using .MatDeforModule: rotstressvec!
using .MatDeforModule: dett, strainvdet, strainvtr, strainttov!, strainvtot!, stressttov!, stressvtot!
using .MatDeforModule: tens4checksymmetry, tens4symmtto6x6t!, tens4symmt6x6tot!, tens4dot2!, tens4ijkl!, tens4ikjl!, tens4iljk!, tens4identity!, tens4transposor!, tens4tracor!, tens4symmetrizor!, tens4skewor!, tens4deviator!
# Exported: abstract type for  models of deformation,  conversion methods  for strain and stress, transformations  of strain and stress
export AbstractMatDefor
export rotstressvec!
export dett, strainvdet, strainvtr, strainttov!, strainvtot!, stressttov!, stressvtot!
export tens4checksymmetry, tens4symmtto6x6t!, tens4symmt6x6tot!, tens4dot2!, tens4ijkl!, tens4ikjl!, tens4iljk!, tens4identity!, tens4transposor!, tens4tracor!, tens4symmetrizor!, tens4skewor!, tens4deviator!

using .MatDeforLinearElasticModule: AbstractMatDeforLinearElastic
# Exported: type of  isotropic elastic material
export AbstractMatDeforLinearElastic

using .MatDeforElastIsoModule: MatDeforElastIso
# Exported: type of  isotropic elastic material
export MatDeforElastIso

using .MatDeforElastOrthoModule: MatDeforElastOrtho
# Exported: type of orthotropic elastic material
export MatDeforElastOrtho

using .FEMMDeforLinearBaseModule: AbstractFEMMDeforLinear, stiffness, nzebcloadsstiffness, thermalstrainloads, mass, inspectintegpoints
# Exported: abstract type for linear information, discretization methods for the abstract type
export AbstractFEMMDeforLinear, stiffness, nzebcloadsstiffness, thermalstrainloads, mass, inspectintegpoints

using .FEMMDeforLinearModule: FEMMDeforLinear
# Exported: type for linear deformation
export FEMMDeforLinear

using .FEMMDeforWinklerModule: FEMMDeforWinkler, surfacenormalspringstiffness
# Exported: type for distributed-spring support, discretization method
export FEMMDeforWinkler, surfacenormalspringstiffness

using .FEMMDeforLinearMSModule: FEMMDeforLinearMSH8, FEMMDeforLinearMST10, stiffness, nzebcloadsstiffness, thermalstrainloads, inspectintegpoints
# Exported: type for mean-strain solid elements, discretization methods
export FEMMDeforLinearMSH8, FEMMDeforLinearMST10, stiffness, nzebcloadsstiffness, thermalstrainloads, inspectintegpoints

using .FEMMDeforLinearIMModule: FEMMDeforLinearIMH8, stiffness, nzebcloadsstiffness, thermalstrainloads, inspectintegpoints
# Exported: type for mean-strain solid elements, discretization methods
export FEMMDeforLinearIMH8, stiffness, nzebcloadsstiffness, thermalstrainloads, inspectintegpoints

using .FEMMDeforSurfaceDampingModule: FEMMDeforSurfaceDamping, dampingABC
#Exported: type for surface damping (absorbing boundary conditions)
export FEMMDeforSurfaceDamping, dampingABC

using .FEMMDeforLinearNICEModule: FEMMDeforLinearNICEH8, FEMMDeforLinearNICET4, stiffness, nzebcloadsstiffness, thermalstrainloads, inspectintegpoints
# Exported: type for NICE (Nodally-integrated continuum elements) solid elements, discretization methods
export FEMMDeforLinearNICEH8, FEMMDeforLinearNICET4, stiffness, nzebcloadsstiffness, thermalstrainloads, inspectintegpoints

using .FEMMDeforLinearESNICEModule: FEMMDeforLinearESNICET4, FEMMDeforLinearESNICEH8, stiffness, nzebcloadsstiffness, thermalstrainloads, inspectintegpoints
# Exported: type for ESICE (Energy-sampling stabilized nodally-integrated continuum elements) solid elements, discretization methods
export FEMMDeforLinearESNICET4, FEMMDeforLinearESNICEH8, stiffness, nzebcloadsstiffness, thermalstrainloads, inspectintegpoints

end # module


text:
module MSTM

include("v3.jl")
include("v4.jl")

end


text:
using Documenter, MAFE

makedocs(
    modules = [MAFE],
    format = :html,
    checkdocs = :exports,
    sitename = "MAFE.jl",
    pages = Any["index.md"]
)

deploydocs(
    repo = "github.com/hendri54/MAFE.jl.git",
)


text:
"Connection points at which energy balance is upheld, has `axis` and `addto!`"
struct Bus <: Component
    model::AbstractEnergyModel
    class::Symbol
end

"Represents a synchronuous zone. Induced by determine_subnetworks!(model)"
struct SubNetwork <: Component
    model::AbstractEnergyModel
    class::Symbol
    buses::Axis
end

mutable struct EnergyModel{MT <: ModelType, TF <: PM.AbstractPowerFormulation} <: AbstractEnergyModel
    devices::Dict{Symbol,Device}
    subnetworks::Dict{Symbol,SubNetwork}
    buses::Dict{Symbol,Bus}
    axes::Dict{Symbol,Axis}
    data::AbstractData
    parent::Union{AbstractEnergyModel,Nothing}
    jumpmodel::Union{JuMP.AbstractModel,Nothing}
    jumpnames::Bool
    jumpobjects::Dict{Symbol,Dict{Symbol}{Any}}
end

ensure_optimizerfactory(opt::JuMP.OptimizerFactory) = opt
ensure_optimizerfactory(::Type{T}) where T <: MOI.AbstractOptimizer = with_optimizer(T)

function EnergyModel(::Type{MT}, ::Type{TF}, data::AbstractData;
                     parent=nothing, jumpmodel=nothing, jumpnames=true, optimizer=nothing) where
        {MT <: ModelType, TF <: PM.AbstractPowerFormulation}

    if isnothing(jumpmodel)
        jumpmodel = isnothing(optimizer) ? JuMP.Model() : JuMP.Model(ensure_optimizerfactory(optimizer))
    end

    EnergyModel{MT,TF}(Dict{Symbol,Device}(),
                       Dict{Symbol,SubNetwork}(),
                       Dict{Symbol,Bus}(),
                       Dict{Symbol,Axis}(),
                       data,
                       parent,
                       jumpmodel,
                       jumpnames,
                       Dict{Symbol, Dict{Symbol}{Any}}())
end

function EnergyModel(other::EnergyModel{MT,TF};
                     modeltype=nothing, transmissionform=nothing,
                     jumpnames=nothing, kwargs...) where
        {MT <: ModelType, TF <: PM.AbstractPowerFormulation}

    if isnothing(modeltype) modeltype = MT end
    if isnothing(transmissionform) transmissionform = TF end

    kwargs = Dict{Symbol,Any}(kwargs)
    get!(kwargs, :jumpnames, other.jumpnames)

    model = EnergyModel(modeltype, transmissionform, other.data; kwargs...)

    merge!(model.buses, Dict(k=>withmodel(v, model) for (k,v) in other.buses))
    merge!(model.devices, Dict(k=>withmodel(v, model) for (k,v) in other.devices))
    merge!(model.axes, other.axes)

    determine_subnetworks!(model)

    model
end

EnergyModel(; kwargs...) = EnergyModel(Data(); kwargs...)
EnergyModel(filename; kwargs...) = EnergyModel(load(filename); kwargs...)
EnergyModel(data::AbstractData; modeltype=ExpansionModel, transmissionform=PM.DCPlosslessForm, kwargs...) = load(EnergyModel(modeltype, transmissionform, data; kwargs...))

function load(m::EnergyModel)
    for (class, T) in components(m.data)
        push!(m, T(m, class))
    end
    determine_subnetworks!(m)
    m
end

Base.push!(m::EnergyModel, d::Device) = (m.devices[d.class] = d)
Base.push!(m::EnergyModel, c::SubNetwork) = (m.subnetworks[c.class] = c)
Base.push!(m::EnergyModel, c::Bus) = (m.buses[c.class] = c)
Base.push!(m::EnergyModel, ax::Axis{name}) where name = (m.axes[name] = ax)

set_snapshots!(m::EnergyModel, snapshots::AbstractArray) =
    push!(m, Axis{:snapshots}(collect(snapshots)))

with_model(d::Device, DF::DeviceFormulation) = with_formulation(d, DF)
with_model(d::Device, ::Type{T}) where T <: Device = model()
function set_formulation!(m::EnergyModel, changes::Pair...)
    for (selector, newmodel) in changes
        for d in devices(m, selector)
            push!(m, with_formulation(d, newmodel))
        end
    end
end

store_results!(m::EnergyModel) = store_results!(m.data, m)

devices(m::EnergyModel) = values(m.devices)
devices(sn::SubNetwork) = devices(model(sn))

devices(m::EnergyModel, T::Type{<:Device}) = (d for d = devices(m) if isa(d, T))
devices(sn::SubNetwork, T::Type{<:Device}) = (d for d = devices(sn) if isa(d, T))

buses(m::EnergyModel) = values(m.buses)
subnetworks(m::EnergyModel) = values(m.subnetworks)

JuMP.optimize!(m::EnergyModel; kwargs...) = optimize!(m.jumpmodel; kwargs...)

Base.show(io::IO, m::EnergyModel) = print(io, typeof(m), " with ", length(axis(m, Bus)), " buses and ", length(m.devices), " devices")
Base.show(io::IO, c::Component) = print(io, typeof(c), "($(c.class))")
function Base.show(io::IO, c::Device{DF}) where DF <: DeviceFormulation
    print(io, typeof(c), "($(c.class))")
    demoted_DF = demote_formulation(modeltype(model(c)), DF)
    demoted_DF != DF && print(io, " (demoted to $demoted_DF)")
end

# function Base.show(io::IO, ::MIME"text/plain", d::Device)
#     println(io, d, " with ")
#     println(io, "* ", length(d.vars), " variables")
#     print(io, "* ", length(d.constrs), " constraints")
# end

model(m::EnergyModel) = m
model(c::Component) = c.model
modeltype(m::EnergyModel{MT,TF}) where {MT <: ModelType, TF <: PM.AbstractPowerFormulation} = MT
transmissionform(m::EnergyModel{MT,TF}) where {MT <: ModelType, TF <: PM.AbstractPowerFormulation} = TF

naming(c::Component) = c.class
naming(c::Component, args...) = Symbol(naming(c), flatten((:(::), a) for a=args)...)

Base.findall(pred::Base.Fix2{typeof(in), <:Axis}, d::Device) = intersect((findall(pred, get(d, attr, axis(d))) for attr = busattributes(d))...)
Base.findall(pred::Base.Fix2{typeof(in), <:Axis}, c::Bus) = findall(pred, axis(c).val)

function Base.getindex(m::EnergyModel, class::Symbol)
    for t = (:devices, :buses, :subnetworks)
        haskey(getfield(m, t), class) && return getfield(m, t)[class]
    end
    throw(KeyError(class))
end
Base.getindex(m::EnergyModel, ::Type{T}) where T <: Device = ContainerView(m, Dict{Symbol,T}(d.class=>d for d=devices(m, T)))
Base.getindex(m::EnergyModel, ::Type{Bus}) = ContainerView(m, Dict{Symbol,Bus}(d.class=>d for d=buses(m)))

Base.getindex(sn::SubNetwork, ::Type{T}) where T <: Device = SubContainerView(model(sn), Dict{Symbol,T}(d.class=>d for d=devices(sn, T)), sn.buses)
Base.getindex(sn::SubNetwork, ::Type{Bus}) = SubContainerView(model(sn), Dict{Symbol,Bus}(d.class=>d for d=buses(model(sn))), sn.buses)

"""
    Base.get(c::Component, attr::Symbol)
    Base.get(c::Component, attr::Symbol, axes...)

For the component `c` get the attribute `attr`, which is either the

- JuMP variable or constraint, or the
- parameter data from the `Data` object for the `class` of the component

If `axes` are provided they must be `Symbol`s or AxisArrays.Axis objects, to
specify which dimensions need to be added flexibly, ie by wrapping in a
WrappedArray.

If you want a JuMP object or a parameter specifically, use `getjump` or
`getparam` instead.

# Examples

```julia
c = m[:onwind]

get(c, :p)   # gets the JuMP variable

get(c, :p_max_pu) # gets the plant availability
                  # (which might be static or a timeseries)

get(c, :p_max_pu, axis(c), :snapshots) # plant availability as timeseries
```
"""
function Base.get(c::Component, attr::Symbol)
    ret = getjump(c, attr)
    !isnothing(ret) ? ret : getparam(c, attr)
end
Base.get(c::Component, attr::Symbol, axes...) = WrappedArray(get(c, attr), axes...)

Base.getindex(c::Component, attr::Symbol) = get(c, attr)

getjump(m::EnergyModel, c::Component, attr::Symbol) = getjump(m, c.class, attr)
function getjump(m::EnergyModel, class::Symbol, attr::Symbol)
    ret = get(get!(m.jumpobjects, class, Dict{Symbol}{Any}()), attr, nothing)
    if !isnothing(ret)
        AxisArray(ret)
    elseif !isnothing(m.parent)
        getjump(m.parent, class, attr)
    else
        nothing
    end
end
getjump(c::Component, attr::Symbol) = getjump(model(c), c, attr)
JuMP.getvalue(c::Component, attr::Symbol) = getvalue.(getjump(c, attr))
JuMP.getdual(c::Component, attr::Symbol) = getdual.(getjump(c, attr))
getparam(c::Component, attr::Symbol) = get(model(c).data, c, attr)

function axis(m::EnergyModel, name::Symbol)
    ax = get(m.axes, name, nothing)
    !isnothing(ax) ? ax : axis(m.data, name)
end
axis(m::EnergyModel, args...) = axis(m.data, args...)
axis(c::Component) = axis(model(c), c)
axis(c::Component, attr) = axis(model(c), attr)

axis(m::EnergyModel, T::Type{<:Component}) = axis(m[T])
axis(m::SubNetwork, T::Type{<:Component}) = axis(m[T])

# Could be specialized to not have to retrieve the whole axis (on the other
# hand, the axis should be cached, anyway)
Base.length(c::Component) = length(axis(c))

issimple(c::Component) = issimple(model(c).data, naming(c))
issimple(::SubNetwork) = false

withmodel(c::T, m::EnergyModel) where T <: Component = T(m, c.class)

add!(m::EnergyModel, ::Type{T}, class::Symbol; kwargs...) where T <: Component = add!(m, T, Axis{class}([class]); kwargs...)
add!(m::EnergyModel, ::Type{T}, class::Symbol, names; kwargs...) where T <: Component = add!(m, T, Axis{class}(names); kwargs...)
function add!(m::EnergyModel, ::Type{T}, ax::Axis{class}; suffix=nothing, parameters...) where {class, T <: Component}
    if !isnothing(suffix)
        ax = Axis{class}(Symbol.(ax.val, suffix))
    end

    push!(m.data, T, ax; parameters...)
    push!(m, T(m, class))
end

max_stars_repo_name: coroa/EnergyModels

max_issues_repo_name: coroa/EnergyModels

max_forks_repo_name: coroa/EnergyModels


text:
mutable struct houseRent <: Rent

end

max_stars_repo_name: adamwillisXanax/accountingJulia

max_issues_repo_name: adamwillisXanax/accountingJulia

max_forks_repo_name: adamwillisXanax/accountingJulia


text:
using Documenter

include("../src/Metatheory.jl")

makedocs(modules = [Metatheory], sitename = "Metatheory.jl",
    pages = [
        "index.md"
        "theories.md"
        "egraphs.md"
        "analysis.md"
        "classic.md"
        "api.md"
    ])

deploydocs(repo = "github.com/0x0f0f0f/Metatheory.jl.git")


text:
using AssemblyCalculus: NeuralArea, neuron_attrib_graph, hebb_update!
using AssemblyCalculus: fire!
using Graphs: adjacency_matrix

@testset "Hebbian Updates" begin
    # Make an area
    num_neurons = 10
    T = Float32
    area1 = NeuralArea(num_neurons, 3, 1.f0)
    area2 = NeuralArea(num_neurons, 3, 1.f0)
    current = ones(T, 10)
    ag = neuron_attrib_graph(T, NeuralArea{T})
    # Make synapses
    for i in 1:num_neurons
        source = area1[i]
        target = area2[i]
        source[area2, target] = 1.f0
        if i > (num_neurons / 2)
            contrib = get!(ag, target)
            push!(contrib, source)
        end
    end

    # Ion current array update
    hebb_update!(area1, current)
    @test all(current .== 1.f0)
    random_firing!(area1)
    hebb_update!(area1, current)
    @test length(findall(current .> 1f0)) == 3
    area2.firing =[1, 9, 10]
    # Synapse Update
    hebb_update!(area2, ag)
    syn = map(x -> collect(values(x[area2]))[1], area1)
    @test all(syn[1:8] .== 1f0)
    @test all(syn[9:10] .== 2f0)

    init_currents = zeros(T, 10)
    init_currents[1:5] .= rand(T, 5)
    currents = zeros(T, 10)
    pows = [1., 2, 4, 1, 8]
    currents[1:5] .= init_currents[1:5] .* (2f0 .^ pows)
    ag = neuron_attrib_graph(T, NeuralArea{T})
    ic_currents = zeros(T, 10)
    ic_currents[1:5] .= init_currents[1:5]
    ic = Stimulus(area2, ic_currents)
    hebb_update!(ag, ic, currents, init_currents)
    @test all(ic.currents .== currents)

    ic_currents = zeros(T, 10)
    ic_currents[1:3] .= init_currents[1:3]
    ic = Stimulus(area2, ic_currents)
    for i in 1:3
        contrib = get!(ag, area2[i])
        push!(contrib, area1[i])
    end
    hebb_update!(ag, ic, currents, init_currents)

    @test all(ic.currents[1:3] .== currents[1:3])
    @test all(ic.currents[4:5] .!= currents[4:5])
    syn = map(x -> collect(values(x[area2]))[1], area1)
    @test all(syn[1:3] .== (2f0 .^ pows[1:3]))
    @test all(syn[4:5] .== 1f0)
end

@testset "Firing" begin
    T = Float32
    num_neurons = 10
    areas = [NeuralArea(num_neurons, 3, 1.f0) for i in 1:3]
    a1, a2, a3 = areas
    # Make synapses
    for i in 1:num_neurons
        source = a1[i]
        target = a2[i]
        target2 = a3[i]
        source[a2, target] = 1.1f0
        source[a3, target2] = 1.2f0
        target[a3, target2] = 1.3f0
    end

    # Assembly firing
    currents = zeros(Float32, num_neurons)
    random_firing!.(areas)
    assem = Assembly(a2, Stimulus{T}[], Assembly{T}[])
    ag = neuron_attrib_graph(Float32, typeof(a2))
    fire!(assem, a1, currents, ag)
    @test all(currents .== 0f0)
    @test length(ag) == 0
    fire!(assem, a3, currents, ag)
    @test all(currents[assem.neurons] .== 1.3f0)
    currents[assem.neurons] .= 0f0
    @test all(currents .== 0f0)
    @test length(ag) .== assem.area.assembly_size
    @test all([length(v) for v in values(ag)] .== 1)

    # Active Area firing
    active_areas = [a1, a3]
    area_currents  = [zeros(T, num_neurons) for i in 1:2]
    ag = fire!(active_areas, area_currents)
    @test all(area_currents[1] .== 0f0)
    @test all(area_currents[2][a1.firing] .== 1.2f0)
    @test length(ag) == a1.assembly_size
    @test all([length(v) for v in values(ag)] .== 1)
    @test length(intersect(collect(keys(ag)), a3.neurons)) == a1.assembly_size
end

@testset "Full Simulation" begin
# From scratch simulation
T = Float32
adj = T.([
    0 1 0 0 1 0;
    1 0 0 0 0 1;
    0 0 0 1 0 0;
    0 0 1 0 0 0;
    0 0 0 0 0 1;
    0 0 0 0 1 0;
])

plasticities = [1f0, 1f0, 2f0]
assem_current = [1f0, 1f0, 0f0, 0f0]
ion_current = [1f0, 1f0, 0f0, 0f0]
max_iters = 5
firing = [1, 0, 1, 0]
true_spikes = [ones(1, 1), ones(1, 1)]
init_currents = ion_current + assem_current
currents = deepcopy(init_currents)
true_adj = deepcopy(adj)
sub = true_adj[1:4, 1:4]
for i in 1:max_iters
    total_current = sub * firing + currents
    firing_next = zeros(4)
    spike1 = argmax(total_current[1:2])
    spike2 = argmax(total_current[3:4])
    firing_next[spike1] = 1
    firing_next[spike2 + 2] = 1
    sub[firing_next * transpose(firing) .== 1] .*= (1 + plasticities[1])
    currents[firing_next .== 1] .*= (1 + plasticities[1])
    firing = firing_next
    true_spikes = map(hcat, true_spikes, [[spike1], [spike2]])
end
# Adj update
true_adj[1:4, 1:4] .= sub
# Ion current update
nz = ion_current .!= 0
times_grown = log.(currents[nz] ./ init_currents[nz]) ./ log(2f0)
true_ion_current = zeros(4)
true_ion_current[nz] = (1 + plasticities[1]) .^ times_grown
# Assembly edge update
true_adj[1, 5] *= (1f0 + plasticities[3]) ^ times_grown[1]
true_adj[2, 6] *= (1f0 + plasticities[3]) ^ times_grown[2]

true_distances = vcat(zeros(T, 1, 6), ones(T, 1, 6))
true_distances[1, 1:2] .= 1

# Run sim
area_sizes = [2, 2, 2]
assembly_sizes =[1, 1, 2]
plasticities = [1f0, 1f0, 2f0]
ba = Brain(adj, area_sizes, assembly_sizes, plasticities)
# Initial firing
a1, a2, a3 = ba.areas
a1.firing = [1]
a2.firing = [1]
a3.firing = [1, 2]
# IonCurrents and Assemblies
inputs = [Stimulus(a1, ones(T, 2)), zero_stim(a2)]
assemblies = [Assembly(a3)]
# Simulate
new_assems, spikes, distances = simulate!(
    inputs,
    assemblies,
    tol=0,
    max_iters=5,
    random_initial=false
    )
# Test
@test all(true_adj .== adjacency_matrix(ba))
@test all([all(true_spikes[i] .== spikes[i]) for i in 1:2])
@test all(true_distances .== distances)
@test all([new_assems[i].neurons[1] == true_spikes[i][end] for i in 1:2])

# Rerun with MaxIters stop criteria
timesteps = size(distances)[2] - 1
ba = Brain(adj, area_sizes, assembly_sizes, plasticities)
# Initial firing
a1, a2, a3 = ba.areas
a1.firing = [1]
a2.firing = [1]
a3.firing = [1, 2]
# IonCurrents and Assemblies
inputs = [Stimulus(a1, ones(T, 2)), zero_stim(a2)]
assemblies = [Assembly(a3)]
new_assems, spikes = simulate!(inputs, assemblies, timesteps, random_initial=false)
@test all(true_adj .== adjacency_matrix(ba))
@test all([all(true_spikes[i] .== spikes[i]) for i in 1:2])
@test all([new_assems[i].neurons[1] == true_spikes[i][end] for i in 1:2])

# Simulate w/o recording
ba = Brain(adj, area_sizes, assembly_sizes, plasticities)
# Initial firing
a1, a2, a3 = ba.areas
a1.firing = [1]
a2.firing = [1]
a3.firing = [1, 2]
# IonCurrents and Assemblies
inputs = [Stimulus(a1, ones(T, 2)), zero_stim(a2)]
assemblies = [Assembly(a3)]
new_assems = simulate!(
    inputs,
    assemblies,
    tol=0,
    max_iters=5,
    random_initial=false,
    record=false
    )
@test all(true_adj .== adjacency_matrix(ba))
@test all([new_assems[i].neurons[1] == true_spikes[i][end] for i in 1:2])

end


text:
using SpectralKit: __block_length, __cumulative_block_length, SmolyakIndices,
    __smolyak_length, SmolyakGridShuffle, SmolyakProduct

####
#### blocks
####

@testset "block length" begin
    @test __cumulative_block_length(0) == 1
    @test __cumulative_block_length(1) == 1 + 2
    @test __cumulative_block_length(2) == 1 + 2 + 2
    @test __cumulative_block_length(3) == 1 + 2 + 2 + 4
    @test __cumulative_block_length(4) == 1 + 2 + 2 + 4 + 8
    for i in 1:10
        c = __cumulative_block_length(i)
        cprev = i == 0 ? 0 : __cumulative_block_length(i - 1)
        @test __block_length(i) == c - cprev
    end
end

@testset "block shuffle" begin
    results = Dict([0 => [1],
                    1 => [2, 1, 3],
                    2 => [3, 1, 5, 2, 4],
                    3 => vcat([5, 1, 9, 3, 7], 2:2:8),
                    4 => vcat([9, 1, 17, 5, 13], 3:4:17, 2:2:17)])
    for (b, i) in pairs(results)
        len = length(i)
        @test len == __cumulative_block_length(b)
        ι = SmolyakGridShuffle(len)
        @test length(ι) == len
        @test @inferred eltype(ι) == Int
        @test collect(ι) == i
    end
end

####
#### traversal
####

"""
Naive implementation of Smolyan index iteration, traversing a `CartesianIndices` and keeping
valid indexes. For testing/comparison. Returns a vector of `indexes => blocks` pairs.
"""
function smolyak_indices_check(N, B, M)
    m = __cumulative_block_length(M)
    b_table = fill(M, m)
    for b in (M-1):(-1):0
        b_table[1:__cumulative_block_length(b)] .= b
    end
    T = NTuple{N,Int}
    result = Vector{Pair{T,T}}()
    for ι in CartesianIndices(ntuple(_ -> 1:m, N))
        ix = Tuple(ι)
        blocks = map(i -> b_table[i], ix)
        if sum(blocks) ≤ B
            push!(result, ix => blocks)
        end
    end
    result
end

@testset "Smolyak indices" begin
    for B in 0:3
        for M in 0:B
            for N in 1:4
                ι = SmolyakIndices{N}(SmolyakParameters(B, M))
                x1 = @inferred collect(ι)
                x2 = first.(smolyak_indices_check(N, B, M))
                len = @inferred __smolyak_length(Val(N), Val(B), M)
                @test x1 == x2
                @test len == length(x1) == length(ι) == length(x2)
            end
        end
    end
end

@testset "Smolyak product primitives" begin
    for B in 0:3
        for M in 0:B
            for N in 1:4
                ι = SmolyakIndices{N}(SmolyakParameters(B, M))
                ℓ = __cumulative_block_length(min(B,M))
                sources = SVector{N}([rand(SVector{ℓ, Float64}) for _ in 1:N])
                P = SmolyakProduct(ι, sources)
                @test length(ι) == length(P)
                @test eltype(P) == Float64
                for (i, p) in zip(ι, P)
                    @test prod(getindex.(sources, i)) ≈ p
                end
            end
        end
    end
end

####
#### api
####

@testset "Smolyak API sanity checks" begin
    f(x) = (x[1] - 3) * (x[2] + 5) # linear function, just a sanity check
    basis = smolyak_basis(Chebyshev, InteriorGrid(), SmolyakParameters(3),
                          (BoundedLinear(0, 4), BoundedLinear(0, 3)))
    @test @inferred(domain(basis)) == ((0, 4), (0, 3))
    x = @inferred grid(Float64, basis)
    M = @inferred collocation_matrix(basis, x)
    θ = M \ f.(x)
    @test sum(abs.(θ) .> 1e-8) == 4
    y1, y2 = range(domain(basis)[1]...; length = 100), range(domain(basis)[2]...; length = 100)
    for y1 in y1
        for y2 in y2
            y = SVector(y1, y2)
            @test linear_combination(basis, θ, y) ≈ f(y)
        end
    end
end

@testset "Smolyak API allocations" begin
    basis = smolyak_basis(Chebyshev, InteriorGrid(), SmolyakParameters(3),
                          (BoundedLinear(0, 4), BoundedLinear(0, 3)))
    y = SVector(1.0, 2.0)
    θ = randn(dimension(basis))
    @inferred linear_combination(basis, θ, y)
    @test @ballocated(linear_combination($basis, $θ, $y)) == 0
end

@testset "Smolyak approximation" begin
    function f(x)               # cross-exp is particularly nasty for Smolyak
        x1, x2, x3, x4 = x
        exp(0.3 * (x1 - 2.0) * (x2 - 1.0)) + exp(-x3 * (x2 + 5)) - 1 / log(x4^2 + 10)
    end
    basis = smolyak_basis(Chebyshev, InteriorGrid(), SmolyakParameters(5),
                          (BoundedLinear(0, 4), BoundedLinear(-2.0, 2.0),
                           SemiInfRational(1.0, 2.0), InfRational(0, 4)))
    x = grid(basis)
    θ = collocation_matrix(basis, x) \ f.(x)
    for x in x
        @test linear_combination(basis, θ, x) ≈ f(x) atol = 1e-10
    end
    s = SobolSeq([0, -2, 1, -10], [4, 2, 5, 10])
    skip(s, 1 << 5)
    Δabs, Δrel = mapreduce((x, y) -> max.(x, y), Iterators.take(s, 10^5)) do x
        fx = f(x)
        Δabs = abs(fx - linear_combination(basis, θ, x))
        Δrel = Δabs / (1 + abs(fx))
        Δabs, Δrel
    end
    @test Δabs ≤ 1e-3
    @test Δrel ≤ 5e-4
end

@testset "Smolyak derivatives" begin
    function f(x)
        x1, x2 = x
        7.0 + x1 + 2*x2 + 3*x2*x1
    end
    basis = smolyak_basis(Chebyshev, InteriorGrid(), SmolyakParameters(2),
                          (BoundedLinear(0, 4), BoundedLinear(-5.0, 5.0)))
    x = grid(basis)
    θ = collocation_matrix(basis, x) \ f.(x)
    F = linear_combination(basis, θ)
    for x in x
        @test F(x) ≈ f(x) atol = 1e-10
        @test gradient(F, x) ≈ gradient(f, x)
    end
end


text:
using StructuralCausalModels, Test

ProjDir = @__DIR__

amat_data = transpose(reshape([
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,
  0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,
  0,0,0,0,1,0,1,0,1,1,0,0,0,0,0,0,
  1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  1,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,
  0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0
], (16,16)));

vars = [Symbol("n$i") for i in 1:size(amat_data, 1)]
a = NamedArray(Int.(amat_data), (vars, vars), ("Rows", "Cols"));

dag = DAG("ex_11", a)

@testset "set_dag_df" begin

  @test dag.df == nothing

  try
    set_dag_df!(dag, DataFrame())
  catch ex
    @test ex == AssertionError("DataFrame has different number of columns")
  end

  @test dag.df == nothing
  set_dag_df!(dag, nothing)
  @test dag.df == nothing
  set_dag_df!(dag, DataFrame(); force=true)
  @test dag.df == DataFrame()
  set_dag_df!(dag, nothing)
  @test dag.df == nothing

end


text:
module Translation

using Random

const filename = joinpath(tempdir(), "translation-$(randstring()).jld")

module Writing

using JLD
import ..Translation: filename

mutable struct MyType
    a::Int
end

jldopen(filename, "w") do file
    truncate_module_path(file, Writing)
    write(file, "x", MyType(3))
end

end # Writing

module Reading

using JLD, Compat.Test
import ..Translation: filename

mutable struct MyType
    a::Int
    b::Float32
end

mutable struct MyOldType
    a::Int
end

translate("MyType", "Translation.Reading.MyOldType")

t = jldopen(filename) do file
    read(file, "x")
end

@test isa(t, MyOldType)

JLD.readas(x::MyOldType) = MyType(x.a, 0f0)

t = jldopen(filename) do file
    read(file, "x")
end

@test isa(t, MyType)

end # Reading

end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export spoiler

"""
    spoiler(;kwargs...)
    spoiler(children::Any;kwargs...)
    spoiler(children_maker::Function;kwargs...)

A Spoiler component.
Hide long sections of content under spoiler. For more information, see: https://mantine.dev/core/spoiler/
Keyword arguments:
- `children` (String; optional): Primary content
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `hideLabel` (String; optional): Label for close spoiler action
- `initialState` (Bool; optional): Initial spoiler state, true to wrap content in spoiler, false to show content without spoiler, opened state will be updated on mount
- `loading_state` (optional): Object that holds the loading state object coming from dash-renderer. loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `maxHeight` (Real; optional): Max height of visible content, when this point is reached spoiler appears
- `showLabel` (String; optional): Label for open spoiler action
- `style` (Dict; optional): Inline style override
"""
function spoiler(; kwargs...)
        available_props = Symbol[:children, :id, :hideLabel, :initialState, :loading_state, :maxHeight, :showLabel, :style]
        wild_props = Symbol[]
        return Component("spoiler", "Spoiler", "dash_mantine_components", available_props, wild_props; kwargs...)
end

spoiler(children::Any; kwargs...) = spoiler(;kwargs..., children = children)
spoiler(children_maker::Function; kwargs...) = spoiler(children_maker(); kwargs...)


text:
import Term.renderables: Renderable, RenderableText, AbstractRenderable
import Term.segment: Segment
import Term: fillin

@testset "\e[34mSegment" begin
    seg = Segment("test", "default")
    @test seg.text == "\e[22mtest\e[22m"
    @test seg.measure.w == 4

    @test_nothrow println(seg)

    seg = Segment("test", "red")
    @test seg.text == "\e[31mtest\e[39m"
    @test seg.measure.w == 4

    seg = Segment("aa\n[blue]123[/blue]")
    @test seg.measure.w == 3
    @test seg.measure.h == 2

    seg = Segment("test")
    @test seg * "test" isa Segment
    @test (seg*"t2").text == "testt2"
    @test ("t2"*seg).text == "t2test"
    @test (seg*seg).text == "testtest"
end

@testset "\e[34mRenderables - Renderable" begin
    lorem = "Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore"
    r = RenderableText(lorem)

    @test r.measure.w == 99
    @test r.measure.h == 1

    r2 = RenderableText(r)
    @test r.measure.w == 99
    @test r.measure.h == 1

    r = RenderableText(lorem; width=22)
    @test string(r) == "Lorem ipsum dolor    \nsit amet, consectetur\nadipiscing elit,     \nsed do eiusmod tempor\nincididunt ut labore "
    @test r.measure.w == 21

    r = RenderableText(lorem; width=22, style="red")
    @test string(r) == "\e[31mLorem ipsum dolor    \e[39m\n\e[31msit amet, consectetur\e[39m\n\e[31madipiscing elit,     \e[39m\n\e[31msed do eiusmod tempor\e[39m\n\e[31mincididunt ut labore \e[39m"
    @test r.measure.w == 21
end


text:
#=##############################################################################
# DESCRIPTION
    Testing of an isolated, planar, 45-deg swept-back wing (Bertin's planar wing
    in Example 7.2, pp. 343 of Bertin's Aerodynamics for Engineers). Validated
    with experimental data from Weber and Brebner (1958), Low-speed tests on
    45-deg swept-back wings, part I, Tables 3 and 4.

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Dec 2019
  * License   : MIT
=###############################################################################

"""
    Test FLOWVLM solver on kinematics of an isolated, planar, swept wing.
"""
function isolatedwing(; xfoil=true,
                        # OUTPUT OPTIONS
                        save_path=nothing,
                        run_name="isolatedwing",
                        prompt=true,
                        verbose=true, v_lvl=0
                        )

    # ------------ PARAMETERS --------------------------------------------------
    if verbose; println("\t"^(v_lvl+1)*"Defining parameters..."); end;

    # Rotor parameters
    rotor_file = "apc10x7.csv"          # Rotor geometry
    data_path = uns.def_data_path       # Path to rotor database
    pitch = 0.0                         # (deg) collective pitch of blades
    n_r = 10                            # Number of blade elements
    # xfoil = false                     # Whether to run XFOIL
    pos_bs = [-0.5, 0.5]                # Position on the semispan of every rotor
    CWs = [false, false]                # Clock-wise rotation of every rotor

    # Read radius of this rotor and number of blades
    R, B = uns.read_rotor(rotor_file; data_path=data_path)[[1,3]]

    # Simulation parameters
    J = 0.6                             # Advance ratio Vinf/(nD)
    ReD07 = 1.5e6                       # Diameter-based Reynolds at 70% span
    ReD = ReD07/0.7                     # Diameter-based Reynolds
    rho = 1.225                         # (kg/m^3) air density
    mu = 1.81e-5                        # (kg/ms) air dynamic viscosity
    nu = mu/rho
    sound_spd = 343                     # (m/s) speed of sound
    RPM = uns.calc_RPM(ReD, J, R, 2*R, nu) # RPM
    magVinf = J*RPM/60*2*R              # (m/s) freestream velocity
    Minf = magVinf / sound_spd
    Mtip = 2*pi*RPM/60*R / sound_spd
    Rob = 0.075                         # Ratio of rotor radius over span
    b = R/Rob                           # (m) wing span
    Reb = magVinf*b/nu

    # magVinf = 163*0.3048              # (m/s) freestream
    rhoinf = rho                        # (kg/m^3) air density
    alpha = 4.2                         # (deg) angle of attack
    qinf = 0.5*rhoinf*magVinf^2         # (Pa) static pressure

    # Wing geometry
    twist = 0.0                         # (deg) root twist
    lambda = 45.0                       # (deg) sweep
    gamma = 0.0                         # (deg) Dihedral
    # b = 98*0.0254                     # (m) span
    ar = 5.0                            # Aspect ratio
    tr = 1.0                            # Taper ratio
    n_w = 4*2^3                         # Number of horseshoes
    r_w = 3.0                           # Expansion ratio of horeshoes
    n_rot = 1.75                        # Factor for HS density behind rotors

    chord = b/ar
    Rec = magVinf*chord/nu

    # Freestream function
    # Vinf(X, t) = magVinf*[cos(alpha*pi/180), 0.0, sin(alpha*pi/180)]
    # Vinf(X, t) = 1e-12*ones(3)      # (Don't make this zero or things will break)
    # Here I had to give it an initial freestream or the unsteady shedding would
    # in the first step
    Vinf(X, t) = t==0 ? magVinf*[1,0,0] : 1e-12*ones(3)

    # Solver parameters
    nrevs = 10                          # Number of revolutions in simulation
    nsteps_per_rev = 72                 # Time steps per revolution
    # p_per_step = 2                    # Sheds per time step
    p_per_step = 1
    ttot = nrevs/(RPM/60)               # (s) total simulation time
    nsteps = nrevs*nsteps_per_rev       # Number of time steps
    lambda_vpm = 2.125                  # Core overlap
    overwrite_sigma = lambda_vpm * 2*pi*R/(nsteps_per_rev*p_per_step) # Smoothing core size
    surf_sigma = R/10                   # Smoothing radius of lifting surface
    # vlm_sigma = surf_sigma              # Smoothing radius of VLM
    vlm_sigma = -1                      # NOTE: keep this disabled or the VLM will go awry
    shed_unsteady = true                # Shed particles from unsteady loading
                                        # Max particles for memory pre-allocation
    max_particles = ((2*n_r+1)*B + 2*n_w+1)*nrevs*nsteps_per_rev*p_per_step
    plot_disc = true                    # Plot blade discretization for debugging
    wake_coupled = true                 # Coupled VLM solver with VPM
    vlm_rlx = 0.75                      # VLM relaxation (deactivated with -1)
    # vlm_rlx = 0.0
    vlm_init = true                     # Initialize with the VLM semi-infinite wake solution

    if verbose
        println("\t"^(v_lvl+2)*"J:\t\t$(J)")
        println("\t"^(v_lvl+2)*"ReD07:\t\t$(ReD07)")
        println("\t"^(v_lvl+2)*"RPM:\t\t$(ceil(Int, RPM))")
        println("\t"^(v_lvl+2)*"Mtip:\t\t$(round(Mtip, digits=3))")
        println("\t"^(v_lvl+2)*"Minf:\t\t$(round(Minf, digits=3))")
        println("\t"^(v_lvl+2)*"R/b:\t\t$(round(Rob, digits=3))")
        println("\t"^(v_lvl+2)*"b:\t\t$(round(b, digits=2)) m")
        println("\t"^(v_lvl+2)*"c:\t\t$(round(chord, digits=2)) m")
        println("\t"^(v_lvl+2)*"Reb:\t\t$(ceil(Int, Reb))")
        println("\t"^(v_lvl+2)*"Rec:\t\t$(ceil(Int, Rec))")
        println("\t"^(v_lvl+2)*"Solution time: "*
            "$(round(magVinf*nrevs/(RPM/60)/b, digits=3)) spans ($nrevs revolutions)")
        println("\t"^(v_lvl+2)*"vlm_fsgm:\t\t$(round(vlm_sigma/b, digits=6))")
    end

    # ------------ GENERATING GEOMETRY -----------------------------------------
    if verbose; println("\t"^(v_lvl+1)*"Generating geometry..."); end;

    # # Generate rotor
    # rotor = uns.generate_rotor(rotor_file; pitch=pitch,
    #                                         n=n_r, CW=CW, ReD=ReD,
    #                                         verbose=verbose, xfoil=xfoil,
    #                                         data_path=data_path,
    #                                         plot_disc=plot_disc,
    #                                         verbose=verbose, v_lvl=v_lvl+1)

    # Wing discretization
    _pos_bs = unique(abs.(pos_bs))
    if length(pos_bs)!=0 && length(pos_bs)!=1 && length(pos_bs)/length(_pos_bs)!=2
        error("Logic error: Position of rotors is expected to be symmetric.")
    end
    refinement = []
    for pos_b in _pos_bs
        # Add wing section before rotor
        pos = pos_b - R/(b/2)                   # Position along semi span
        if pos > 0                              # Length of segment
            len = pos - (length(refinement)!=0 ? refinement[end][1] : 0)
            nfrac = len                         # fraction of horseshoes
            if length(refinement) != 0
                push!(refinement, [len/2, nfrac/2, r_w])
                push!(refinement, [len/2, nfrac/2, 1/r_w])
            else
                push!(refinement, [len, nfrac, 1/r_w])
            end
        end

        # Add rotor section
        pos = min(pos_b + R/(b/2), 1.0)
        len = pos - (length(refinement)!=0 ? refinement[end][1] : 0)
        nfrac = len
        push!(refinement, [len, n_rot*nfrac, 1.01])
    end
    if length(refinement) == 0
        push!(refinement, [1.0, 1.0, 1/r_w])
    end
    # Add tip
    sumc = sum(ref[1] for ref in refinement)
    if abs(sumc-1.0) > 0.0001
        len = abs(sumc-1)
        nfrac = 1.5*len
        push!(refinement, [len/2, nfrac/2, r_w])
        push!(refinement, [len/2, nfrac/2, 1/r_w])
    end

    # Generate wing
    wing = vlm.simpleWing(b, ar, tr, twist, lambda, gamma;
                                                n=n_w, refinement=refinement)

    # Pitch wing to corresponding angle of attack
    O = zeros(3)                                    # Coordinate system origin
    Oaxis = gt.rotation_matrix2(0.0, -alpha, 0.0)   # Coordinate system axes
    vlm.setcoordsystem(wing, O, Oaxis)

    # ------------- SIMULATION SETUP -------------------------------------------
    if verbose; println("\t"^(v_lvl+1)*"Simulation setup..."); end;

    # ----- VEHICLE DEFINITION
    # System of all FLOWVLM objects
    system = vlm.WingSystem()
    vlm.addwing(system, run_name*"_wing", wing)
    # vlm.addwing(system, run_name*"_rotor", rotor)

    # Systems of rotors
    # rotors = vlm.Rotor[rotor]   # Defining this rotor as its own system
    rotor_systems = ()

    # System solved through VLM solver
    vlm_system = system

    # Wake-shedding system
    wake_system = system

    # FVS's Vehicle object
    vehicle = uns.VLMVehicle(   system;
                                vlm_system=vlm_system,
                                rotor_systems=rotor_systems,
                                wake_system=wake_system
                             )

    # ----- MANEUVER DEFINITION
    RPM_fun(t) = 1.0                # RPM (normalized by reference RPM) as a
                                    # function of normalized time

    angle = ()                      # Angle of each tilting system (none in this case)
    # sysRPM = (RPM_fun, )          # RPM of each rotor system
    sysRPM = ()
    Vvehicle(t) = [-1,0,0]          # Translational velocity of vehicle over Vref
    anglevehicle(t) = zeros(3)      # (deg) angle of the vehicle

    Vref = magVinf                  # (m/s) aircraft velocity during cruise
    # RPMref = RPM                  # Reference RPM
    RPMref = 0.0

    # FVS's Maneuver object
    maneuver = uns.KinematicManeuver(angle, sysRPM, Vvehicle, anglevehicle)

    # Plot maneuver path and controls
    uns.plot_maneuver(maneuver; vis_nsteps=nsteps)

    # ----- SIMULATION DEFINITION
    Vinit = Vref*Vvehicle(0)       # Initial vehicle velocity
    simulation = uns.Simulation(vehicle, maneuver, Vref, RPMref, ttot; Vinit=Vinit)

    # monitor = generate_monitor_rotor(J, rho, RPM, nsteps; save_path=save_path,
    #                                                         run_name=run_name)
    monitor = generate_monitor_wing(wing, b, ar, nsteps, Vinf, rhoinf, qinf,
                                                        magVinf, wake_coupled)

    # ------------- RUN SIMULATION ---------------------------------------------
    if verbose; println("\t"^(v_lvl+1)*"Running simulation..."); end;
    pfield = uns.run_simulation(simulation, nsteps;
                                      # SIMULATION OPTIONS
                                      Vinf=Vinf,
                                      # SOLVERS OPTIONS
                                      p_per_step=p_per_step,
                                      overwrite_sigma=overwrite_sigma,
                                      vlm_sigma=vlm_sigma,
                                      surf_sigma=surf_sigma,
                                      vlm_init=vlm_init,
                                      vlm_rlx=vlm_rlx,
                                      max_particles=max_particles,
                                      wake_coupled=wake_coupled,
                                      shed_unsteady=shed_unsteady,
                                      extra_runtime_function=monitor,
                                      # OUTPUT OPTIONS
                                      save_path=save_path,
                                      run_name=run_name,
                                      prompt=prompt,
                                      verbose=verbose, v_lvl=v_lvl+1,
                                      save_code=splitdir(@__FILE__)[1],
                                      save_horseshoes=false
                                      )
    return pfield, wing#, rotor
end

function generate_monitor_wing(wing, b, ar, nsteps, Vinf, rhoinf, qinf, magVinf,
                                wake_coupled;
                                figname="monitor_wing", nsteps_plot=1,
                                disp_plot=true, figsize_factor=5/6,
                                save_path=nothing, run_name="wing",
                                nsteps_savefig=10,
                                extra_plots=true)

    # ------------- SIMULATION MONITOR -----------------------------------------
    y2b = 2*wing._ym/b

    # Weber's lift distribution data (Table 3)
    web_2yb = [0.0, 0.041, 0.082, 0.163, 0.245, 0.367, 0.510, 0.653, 0.898, 0.949]
    web_Cl = [0.235, 0.241, 0.248, 0.253, 0.251, 0.251, 0.251, 0.246, 0.192, 0.171]
    web_CL = 0.238
    web_ClCL = web_Cl/web_CL

    # Weber's drag distribution data (Table 3)
    web_Cd = [0.059, 0.025, 0.016, 0.009, 0.007, 0.006, 0.006, 0.004, -0.002, -0.007]
    web_CD = 0.005
    web_CdCD = web_Cd/web_CD

    prev_wing = nothing

    # Critical length to ignore horseshoe forces
    meanchord = b/ar
    lencrit = 0.5*meanchord/vlm.get_m(wing)

    # Name of convergence file
    if save_path!=nothing
        fname = joinpath(save_path, run_name*"_convergence.csv")
    end

    function extra_runtime_function(sim, PFIELD, T, DT)

        aux = PFIELD.nt/nsteps
        clr = (1-aux, 0, aux)

        if PFIELD.nt==0 && (disp_plot || save_path!=nothing)
            figure(figname, figsize=[7*2, 5*2]*figsize_factor)
            subplot(221)
            xlim([-1,1])
            xlabel(L"$\frac{2y}{b}$")
            ylabel(L"$\frac{Cl}{CL}$")
            title("Spanwise lift distribution")

            subplot(222)
            xlim([-1,1])
            xlabel(L"$\frac{2y}{b}$")
            ylabel(L"$\frac{Cd}{CD}$")
            title("Spanwise drag distribution")

            subplot(223)
            xlabel("Simulation time (s)")
            ylabel(L"Lift Coefficient $C_L$")

            subplot(224)
            xlabel("Simulation time (s)")
            ylabel(L"Drag Coefficient $C_D$")

            figure(figname*"_2", figsize=[7*2, 5*1]*figsize_factor)
            subplot(121)
            xlabel(L"$\frac{2y}{b}$")
            ylabel(L"Circulation $\Gamma$")
            subplot(122)
            xlabel(L"$\frac{2y}{b}$")
            ylabel(L"Effective velocity $V_\infty$")

            if extra_plots
                for num in 1:3 # (1==ApA, 2==AB, 3==BBp)
                    figure(figname*"_3_$num", figsize=[7*3, 5*3]*figsize_factor)
                    suptitle("Velocitiy at $(num==1 ? "ApA" : num==2 ? "AB" : "BBp")")
                    for i in 7:9; subplot(330+i); xlabel(L"$\frac{2y}{b}$"); end
                    subplot(331)
                    ylabel(L"$V_\mathrm{VPM}$ Velocity")
                    subplot(334)
                    ylabel(L"$V_\mathrm{VLM}$ Velocity")
                    subplot(337)
                    ylabel(L"$V_\mathrm{kin}$ and $V_\infty$ Velocity")
                    subplot(331)
                    title(L"$x$-component")
                    subplot(332)
                    title(L"$y$-component")
                    subplot(333)
                    title(L"$z$-component")
                end
                figure(figname*"_4", figsize=[7*3, 5*1]*figsize_factor)
                for num in 1:3
                    subplot(130+num)
                    title("Length at $(num==1 ? "ApA" : num==2 ? "AB" : "BBp")")
                    xlabel(L"$\frac{2y}{b}$")
                    if num==1; ylabel("Bound-vortex length"); end;
                end
            end

            # Convergence file header
            if save_path!=nothing
                f = open(fname, "w")
                print(f, "T,CL,CD\n")
                close(f)
            end
        end

        if PFIELD.nt!=0 && PFIELD.nt%nsteps_plot==0 && (disp_plot || save_path!=nothing)
            figure(figname)

            # Force at each VLM element
            Ftot = uns.calc_aerodynamicforce(wing, prev_wing, PFIELD, Vinf, DT,
                                                            rhoinf; t=PFIELD.t,
                                                            lencrit=lencrit)
            L, D, S = uns.decompose(Ftot, [0,0,1], [-1,0,0])
            vlm._addsolution(wing, "L", L)
            vlm._addsolution(wing, "D", D)
            vlm._addsolution(wing, "S", S)

            # Force per unit span at each VLM element
            Vout, lenout = extra_plots ? ([], []) : (nothing, nothing)
            ftot = uns.calc_aerodynamicforce(wing, prev_wing, PFIELD, Vinf, DT,
                                        rhoinf; t=PFIELD.t, per_unit_span=true,
                                        Vout=Vout, lenout=lenout,
                                        lencrit=lencrit)
            l, d, s = uns.decompose(ftot, [0,0,1], [-1,0,0])

            # Lift of the wing
            Lwing = norm(sum(L))
            CLwing = Lwing/(qinf*b^2/ar)
            ClCL = norm.(l) / (Lwing/b)

            # Drag of the wing
            Dwing = norm(sum(D))
            CDwing = Dwing/(qinf*b^2/ar)
            CdCD = [sign(dot(this_d, [1,0,0])) for this_d in d].*norm.(d) / (Dwing/b) # Preserves the sign of drag

            vlm._addsolution(wing, "Cl/CL", ClCL)
            vlm._addsolution(wing, "Cd/CD", CdCD)

            subplot(221)
            plot(web_2yb, web_ClCL, "ok", label="Weber's experimental data")
            plot(-web_2yb, web_ClCL, "ok")
            plot(y2b, ClCL, "-", label="FLOWVLM", alpha=0.5, color=clr)

            subplot(222)
            plot(web_2yb, web_CdCD, "ok", label="Weber's experimental data")
            plot(-web_2yb, web_CdCD, "ok")
            plot(y2b, CdCD, "-", label="FLOWVLM", alpha=0.5, color=clr)

            subplot(223)
            plot([0, T], web_CL*ones(2), ":k", label="Weber's experimental data")
            plot([T], [CLwing], "o", label="FLOWVLM", alpha=0.5, color=clr)

            subplot(224)
            plot([0, T], web_CD*ones(2), ":k", label="Weber's experimental data")
            plot([T], [CDwing], "o", label="FLOWVLM", alpha=0.5, color=clr)

            figure(figname*"_2")
            subplot(121)
            plot(y2b, wing.sol["Gamma"], "-", label="FLOWVLM", alpha=0.5, color=clr)
            if wake_coupled && PFIELD.nt!=0
                subplot(122)
                plot(y2b, norm.(wing.sol["Vkin"])/magVinf, "-", label="FLOWVLM", alpha=0.5, color=[clr[1], 1, clr[3]])
                if "Vvpm" in keys(wing.sol)
                    plot(y2b, norm.(wing.sol["Vvpm"]), "-", label="FLOWVLM", alpha=0.5, color=clr)
                end
                plot(y2b, [norm(Vinf(vlm.getControlPoint(wing, i), T)) for i in 1:vlm.get_m(wing)],
                                                            "-k", label="FLOWVLM", alpha=0.5)
            end

            if extra_plots
                m = vlm.get_m(wing)
                for num in 1:3               # (1==ApA, 2==AB, 3==BBp)
                    figure(figname*"_3_$num")
                    for Vi in 1:3           # (1==Vvpm, 2==Vvlm, 3==Vinf && Vkin)
                        for xi in 1:3       # (1==Vx, 2==Vy, 3==Vz)
                            subplot(330 + (Vi-1)*3 + xi)
                            if Vi!=3
                                plot(y2b, [Vout[(i-1)*3 + num][Vi][xi] for i in 1:m], color=clr, alpha=0.5)
                            else
                                plot(y2b, [Vout[(i-1)*3 + num][Vi][xi] for i in 1:m], "k", alpha=0.5)
                                plot(y2b, [Vout[(i-1)*3 + num][Vi+1][xi] for i in 1:m], color=clr, alpha=0.5)
                            end
                        end
                    end
                end
                figure(figname*"_4")
                for num in 1:3
                    subplot(130+num)
                    plot(y2b, [lenout[(i-1)*3 + num] for i in 1:m], color=clr, alpha=0.5)
                end
            end
            if save_path!=nothing

                # Write rotor position and time on convergence file
                f = open(fname, "a")
                print(f, T, ",", CLwing, ",", CDwing, "\n")
                close(f)

                # Save figures
                if PFIELD.nt%nsteps_savefig==0
                    figure(figname)
                    savefig(joinpath(save_path, run_name*"_convergence.png"),
                                                            transparent=false)
                    figure(figname*"_2")
                    savefig(joinpath(save_path, run_name*"_convergence2.png"),
                                                            transparent=false)
                end
            end
        end

        prev_wing = deepcopy(wing)

        return false
    end
end


text:
# Official spec : http://www.libpng.org/pub/png/spec/1.2/png-1.2.pdf

"""
    PNG

A wrapper around the binary encoding of a PNG image, holding it's attributes
"""
struct PNG
    attr::Dict
    data::Vector{UInt8}
end

function PNG(buffer::Base.IOBuffer)
    !check_valid_png(buffer.data) && return false
    attr  = read_info(buffer.data)
    return PNG(attr, buffer.data)
end

"""
    check_valid_png(data)

Checks if the header of the provided buffer is a valid png header.
Also checks if the first chunk is valid according to the PNG standard

Refs:
[1] http://www.libpng.org/pub/png/spec/1.2/png-1.2.pdf
[2] https://en.wikipedia.org/wiki/Portable_Network_Graphics#File_header
"""
function check_valid_png(data)
    #Header https://en.wikipedia.org/wiki/Portable_Network_Graphics#File_header
    header_value = UInt8[0x89, 0x50, 0x4e, 0x47, 0x0d, 0x0a, 0x1a, 0x0a]

    # Avoid accessing memory that is not there
    length(data) < 12 && return false
    # Check the header
    data[1:8] != header_value && return false

    # First chunk: check length
    cnk_start = 8;
    cnk_len = reinterpret(Int32, data[cnk_start.+(1:4)]) |> first |> ntoh
    String(data[cnk_start.+(5:8)]) != "IHDR" && return false
    cnk_len <13 && return false

    return true
end

"""
    read_info(data)

Reads the content of the first chunk of a PNG, containing information such as
height, width and color type.
"""
function read_info(data)
    cnk_start = 8;
    cnk_data = cnk_start + 8;

    # Read chunk data
    width, height = ntoh.(reinterpret(Int32, data[cnk_data.+(1:8)]))
    bit_depth     = ntoh.(reinterpret(Int8,  data[cnk_data+9]))
    col_type      = ntoh.(reinterpret(Int8,  data[cnk_data+10]))
    compression   = ntoh.(reinterpret(Int8,  data[cnk_data+11]))
    filter        = ntoh.(reinterpret(Int8,  data[cnk_data+12]))
    interlace     = ntoh.(reinterpret(Int8,  data[cnk_data+12]))

    attr = Dict(
        :width=>width,
        :height=>height,
        :bitdepth=>bit_depth,
        :colorspace=>col_type,
        :compression=>compression,
        :filter=>filter,
        :interlace=>interlace
    )
    return attr
end

function Base.convert(::Type{PNG}, img::AbstractArray{<:Colorant})
    pb = PipeBuffer()
    save(Stream(format"PNG", pb), img)
    return PNG(pb)
end


text:
# Use
#
#     DOCUMENTER_DEBUG=true julia --color=yes make.jl local [nonstrict] [fixdoctests]
#
# for local builds.

using Documenter
using H3DDetectorSystems

makedocs(
    sitename = "H3DDetectorSystems",
    modules = [H3DDetectorSystems],
    format = Documenter.HTML(
        prettyurls = !("local" in ARGS),
        canonical = "https://oschulz.github.io/H3DDetectorSystems.jl/stable/"
    ),
    pages=[
        "Home" => "index.md",
        "API" => "api.md",
        "LICENSE" => "LICENSE.md",
    ],
    doctest = ("fixdoctests" in ARGS) ? :fix : true,
    linkcheck = ("linkcheck" in ARGS),
    strict = !("nonstrict" in ARGS),
)

deploydocs(
    repo = "github.com/oschulz/H3DDetectorSystems.jl.git",
    forcepush = true
)


text:
module FrostHelperEntityBatcherController

using ..Ahorn, Maple

@mapdef Entity "FrostHelper/EntityBatcher" EntityBatcher(x::Integer, y::Integer, types::String="", effect::String="", depth::Integer=0, parameters::String="", dynamicDepthBatchSplitField::String="", flag::String="", flagInverted::Bool=false)

const placements = Ahorn.PlacementDict(
	"Entity Batcher (Frost Helper + Shader Helper)" => Ahorn.EntityPlacement(
		EntityBatcher
	)
)

const sprite = "editor/FrostHelper/RainbowTilesetController"

function Ahorn.selection(entity::EntityBatcher)
	x, y = Ahorn.position(entity)

	return Ahorn.getSpriteRectangle(sprite, x, y)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::EntityBatcher, room::Maple.Room)
	x, y = Ahorn.position(entity)

	Ahorn.drawSprite(ctx, sprite, x, y)
end

end

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper


text:
##############################
# BUBBLE SORT
##############################

function bubble_sort(arr)
    for i = 1:length(arr), j = 1:(length(arr)-i)
        if arr[j] > arr[j+1]
            arr[j+1], arr[j] = arr[j], arr[j+1]
        end
    end
end

# Get the array from user input
input = map(x->parse(Int, x), split(readline()))

# Sort the array
bubble_sort(input)

# Print the array
for x in input
    print(x," ")
end

max_stars_repo_name: Nikola-Mircic/julia-cheatsheet

max_issues_repo_name: Nikola-Mircic/julia-cheatsheet

max_forks_repo_name: Nikola-Mircic/julia-cheatsheet


text:
using Test

@testset "test object" begin
    @test true
    @test true
    @test true
end


text:
function ImpulseManeuver(arg0::EventDetector, arg1::AttitudeProvider, arg2::Vector3D, arg3::jdouble)
    return ImpulseManeuver((EventDetector, AttitudeProvider, Vector3D, jdouble), arg0, arg1, arg2, arg3)
end

function ImpulseManeuver(arg0::EventDetector, arg1::Vector3D, arg2::jdouble)
    return ImpulseManeuver((EventDetector, Vector3D, jdouble), arg0, arg1, arg2)
end

function event_occurred(obj::AbstractDetector, arg0::SpacecraftState, arg1::jboolean)
    return jcall(obj, "eventOccurred", Action, (SpacecraftState, jboolean), arg0, arg1)
end

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

function get_attitude_override(obj::ImpulseManeuver)
    return jcall(obj, "getAttitudeOverride", AttitudeProvider, ())
end

function get_delta_v_sat(obj::ImpulseManeuver)
    return jcall(obj, "getDeltaVSat", Vector3D, ())
end

function get_handler(obj::AbstractDetector)
    return jcall(obj, "getHandler", EventHandler, ())
end

function get_isp(obj::ImpulseManeuver)
    return jcall(obj, "getIsp", jdouble, ())
end

function get_max_check_interval(obj::AbstractDetector)
    return jcall(obj, "getMaxCheckInterval", jdouble, ())
end

function get_max_iteration_count(obj::AbstractDetector)
    return jcall(obj, "getMaxIterationCount", jint, ())
end

function get_threshold(obj::AbstractDetector)
    return jcall(obj, "getThreshold", jdouble, ())
end

function get_trigger(obj::ImpulseManeuver)
    return jcall(obj, "getTrigger", EventDetector, ())
end

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

function is_forward(obj::AbstractDetector)
    return jcall(obj, "isForward", jboolean, ())
end

function reset_state(obj::AbstractDetector, arg0::SpacecraftState)
    return jcall(obj, "resetState", SpacecraftState, (SpacecraftState,), arg0)
end

function with_handler(obj::AbstractDetector, arg0::EventHandler)
    return jcall(obj, "withHandler", AbstractDetector, (EventHandler,), arg0)
end

function with_max_check(obj::AbstractDetector, arg0::jdouble)
    return jcall(obj, "withMaxCheck", AbstractDetector, (jdouble,), arg0)
end

function with_max_iter(obj::AbstractDetector, arg0::jint)
    return jcall(obj, "withMaxIter", AbstractDetector, (jint,), arg0)
end

function with_threshold(obj::AbstractDetector, arg0::jdouble)
    return jcall(obj, "withThreshold", AbstractDetector, (jdouble,), arg0)
end


text:
module SSLTestGBIF
using SimpleSDMLayers
using GBIF
using Test

temperature = worldclim(1)

kingfisher = GBIF.taxon("Megaceryle alcyon", strict=true)

o = GBIF.occurrences(kingfisher, "hasCoordinate" => "true")

# Extract from a single record
for oc in o
    @test typeof(temperature[oc]) <: Number
end

# Modify a value
numboc = SimpleSDMResponse(zeros(Int64, (20, 20)))
for oc in o
    numboc[oc] += 1
end

@test typeof(convert(Matrix, numboc)) == Matrix{Union{Int64,Nothing}}

@test typeof(clip(temperature, o)) == typeof(temperature)

end


text:
###############################################################################
#
# Equilibrium pressure from volume
#
###############################################################################
"""
    p_from_volume(
                pv::AbstractCapacity{FT},
                rvol::FT,
                T::FT
    ) where {FT<:AbstractFloat}

Calculate equilibrium pressure from relative volume, given
- `pv` [`AbstractCapacity`](@ref) type of struct
"""
function p_from_volume(
            pv::PVCurveLinear{FT},
            rvol::FT,
            T::FT
) where {FT<:AbstractFloat}
    return (rvol - 1) / pv.slope
end

function p_from_volume(
            pv::PVCurveSegmented{FT},
            rvol::FT,
            T::FT
) where {FT<:AbstractFloat}
    @unpack c_all, RWC_apo, RWC_TLP, ϵ_bulk = pv;

    if rvol > RWC_TLP
        return -c_all * R_GAS(FT) * T / (rvol - RWC_apo) * FT(1e-6) +
                ϵ_bulk * (rvol - RWC_TLP)
    elseif rvol > RWC_apo
        return -c_all * R_GAS(FT) * T / (rvol - RWC_apo) * FT(1e-6)
    else
        return FT(-100)
    end
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
#' ---
#' title: Simulation-based calibration of FirstPassageTools.jl
#' author: Garrett Smith
#' date: 07.04.2022
#' ---
#'
#' # Simulation-based Calibration
#'
#' Simulation-based calibration [@talts2018validating] is used to ensure that a Bayesian
#' sampling algorithm is functioning properly. This script does a reasonably thorough job of
#' testing the NUTS sampler available in `Turing.jl` for a simple, one-parameter,
#' non-hierarchical first-passage time model.
#'
#' Note that this script should be run with `julia -t X SimulationBasedCalibration.jl` in
#' order to run with X threads

using Distributions
using Turing
using Plots
using Pkg
Pkg.activate("../../FirstPassageTools.jl/")
using FirstPassageTools
using LinearAlgebra
LinearAlgebra.BLAS.set_num_threads(1)  # Necessary to prevent crashing with the mulithreaded code below.

#' # Defining things
#'
#' Here, we define a simple model with two transient states and a single absorbing state.

T = [-3.0 0; 3.0 -3.0]
A = [0.0 3]
p0 = [1.0, 0]

#' Here, I define the prior on τ, the scaling parameter to be fit. If the sampling algorithm
#' works as it's supposed to, then the posterior should eventually approach the prior when
#' repeatedly sampling τs from the prior, simulating data, and then fitting the model to
#' that data.

#prior = truncated(Normal(1.0, 0.5), lower=0.0)
prior = Normal(0.0, 0.5)

#' Now, we define the model, a combination of the prior and a likelihood.

@model function mod(y)
    # Define prior
    τ ~ prior

    # Define likelihood
    y ~ filldist(fpdistribution(exp(τ)*T, exp(τ)*A, p0), length(y))
end

#' Now, we define the SBC function, which takes a single argument that is the number of
#' times to run the prior-generate-fit sequence, the number of data points to generate, and
#' the number of posterior samples to sample.

function sbc(nranks=1000, ndata=100, nposterior=500)
    ranks = zeros(nranks)
    taus = rand(prior, nranks)
    Threads.@threads for i = 1:nranks
        # Sample a τ from the prior
        #curr = rand(prior)
        #curr = exp(rand(prior))
        curr = exp(taus[i])
        # Generate data from the fpdistribution
        dat = rand(fpdistribution(curr*T, curr*A, p0), ndata)
        # Fit the fpdistribution to the simulated data
        #posterior = sample(mod(dat), NUTS(0, 0.65), nposterior, progress=false)
        #posterior = sample(mod(dat), SMC(25), nposterior, progress=false)
        posterior = sample(mod(dat), PG(25), nposterior, progress=false)
        # Get the rank of curr in the posterior
        #ranks[i] = count(x -> x < curr, posterior[:τ])
        ranks[i] = count(x -> exp(x) < curr, posterior[:τ])
    end
    return ranks
end

#' # Running the calibration
#'
#' Now we run the calibration:

nr = 1000
nd = 50
np = 100
@time rks = sbc(nr, nd, np)

#' # Visualizing the results
#'
#' If the sampler is well-calibrated, the whole histogram should fall within the confidence
#' interval in the plot.

function plot_sbc(ranks, nposterior)
    bn = Binomial(length(ranks), (nposterior+1)^-1)
    fig = histogram(ranks, bins=0:nposterior)
    hspan!(fig, quantile(bn, [0.005, 0.995]), color=:grey, alpha=0.25)
    savefig(fig, "SBCResults.pdf")
end

plot_sbc(rks, np)


text:
import Base: @kwdef
Pkg.add("BitemporalPostgres")
Pkg.add("HTTP")
Pkg.add("JSON")
Pkg.add("TimeZones")
Pkg.add("ToStruct")
Pkg.add("Intervals")
Pkg.add("SearchLight")
Pkg.add("SearchLightPostgreSQL")
Pkg.add("Test")
using Test
import InsuranceContractsController
using InsuranceContractsController.InsuranceContracts
using InsuranceContractsController.InsurancePartners
using InsuranceContractsController.InsuranceTariffs
using BitemporalPostgres
using SearchLight
using TimeZones
using ToStruct
using JSON
using HTTP
using Genie, Genie.Router, Genie.Requests

@testset "CreateContract" begin

    SearchLight.Configuration.load() |> SearchLight.connect
    SearchLight.Migrations.create_migrations_table()
    BitemporalPostgres.up()
    SearchLight.Migrations.up()
# create Partner
    p = Partner()
    pr = PartnerRevision(description = "blue")
    w = Workflow(
        tsw_validfrom = ZonedDateTime(2014, 5, 30, 21, 0, 1, 1, tz"Africa/Porto-Novo"),
    )
    create_entity!(w)
    create_component!(p, pr, w)
    commit_workflow!(w)

# create Tariff
t = Tariff()
tr = TariffRevision(description = "blue")
w0 = Workflow(
    tsw_validfrom = ZonedDateTime(2014, 5, 30, 21, 0, 1, 1, tz"Africa/Porto-Novo"),
)
create_entity!(w0)
create_component!(t, tr, w0)
commit_workflow!(w0)

# create Contract
    c = Contract()
    cr = ContractRevision(description = "blue")
    cpr = ContractPartnerRef(ref_super = c.id)
    cprr = ContractPartnerRefRevision(ref_partner = p.id, description = "blue")

    cpi = ProductItem(ref_super = c.id)
    cpir = ProductItemRevision(position=1,description="blue")

    pitr = ProductItemTariffRef(ref_super = cpi.id)
    pitrr = ProductItemTariffRefRevision(ref_tariff = t.id, description = "blue")

    pipr = ProductItemPartnerRef(ref_super = cpi.id)
    piprr = ProductItemPartnerRefRevision(ref_partner = p.id, description = "blue")

    w1 = Workflow(
        tsw_validfrom = ZonedDateTime(2014, 5, 30, 21, 0, 1, 1, tz"Africa/Porto-Novo"),
    )
    create_entity!(w1)
    create_component!(c, cr, w1)
    create_subcomponent!(c, cpr, cprr, w1)
    create_subcomponent!(c, cpi, cpir, w1)
    create_subcomponent!(cpi, pitr, pitrr, w1)
    create_subcomponent!(cpi, pipr, piprr, w1)
    commit_workflow!(w1)

# end

# update Contract yellow
# @testset "UpdateContractYellow" begin

    cr1 = ContractRevision(ref_component = c.id, description = "yellow")
    w2 = Workflow(
        ref_history = w1.ref_history,
        tsw_validfrom = ZonedDateTime(2016, 5, 30, 21, 0, 1, 1, tz"Africa/Porto-Novo"),
    )
    update_entity!(w2)
    update_component!(cr, cr1, w2)
    commit_workflow!(w2)
    @test w2.ref_history == w1.ref_history

# nd

# update Contract red
# @testset "UpdateContractRed" begin
    cr2 = ContractRevision(ref_component = c.id, description = "red")
    w3 = Workflow(
        ref_history = w2.ref_history,
        tsw_validfrom = ZonedDateTime(2015, 5, 30, 21, 0, 1, 1, tz"Africa/Porto-Novo"),
    )
    update_entity!(w3)
    update_component!(cr1, cr2, w3)
    commit_workflow!(w3)
    @test w3.ref_history == w2.ref_history

    # end of mutations
end

hforest = mkforest(DbId(3), MaxDate, ZonedDateTime(1900, 1, 1, 0, 0, 0, 0, tz"UTC"), MaxDate)
print_tree(hforest)

max_stars_repo_name: michaelfliegner/BitempWebApp

max_issues_repo_name: michaelfliegner/BitempWebApp


text:
cd(@__DIR__)
using Pkg
Pkg.activate(".")

using F16Model, LinearAlgebra, Plots
using BSON: @load
using Flux
using Statistics

#loading neural networks
@load "Models/Pressure_overall_best.bson" Pressure_overall_best
pressuremodel = Pressure_overall_best
@load "Models/Fault_Model.bson" Fault_Model
faultmodel = Fault_Model

#TODO: are these values still correct?
#just copied and pasted these... might want to automate in future
fault_std = 0.310076;fault_mean = 0.741667;
alt_std = 3345.34; alt_mean = 6762.5;
mach_std = 0.123024;mach_mean = 0.488333;
pressure_std = 24622.9; pressure_mean = 53718.4;

#=
Aircraft states are:
npos = North position in ft
epos = East position in ft
alt = height/altitude in ft
phi = Euler angle in rad TODO is this equivalent to Yaw,Pitch,Roll?
theta = Euler angle in rad
psi = Euler angle in rad
Vt = Flight velocity in ft/s
alp = Angle of attack in rad
bet = Side slip angle in rad
p = roll rate
q = pitch rate
r = yaw rate

Aircraft controls are:
T = Thrust in lbs
dele = Elevator angle in deg
dail = Aileron angle in deg
drud = Rudder angle in deg
dlef = Leading edge flap angle in deg
=#

# ===================================================
# Trim the aircraft for steady-level flight at h0,V0
# xbar = trim state
# ubar = trim control
# status = status of the optimization, status = 0 means optimization found solution and (xbar, ubar) defines a valid trim  point/
# prob = data structure from IpOpt.

h0 = 15000; # ft
Vt0 = 500;   # ft/s
xbar, ubar, status, prob = F16Model.Trim(h0,Vt0); # Default is steady-level

# ===================================================
# Linearize about some trim point (x0,u0)
Ac, Bc = F16Model.Linearize(xbar,ubar); # Get continuous-time model.
nx,nu = size(Bc);

# to indclude gps we want Npos and Epos (x[1,2])
ix = [1,2,3,5,7,8,11]; # States I want...  now including Npos (1) and Epos (2)
iu = [2]; # Control variables I want, this how the disturbance comes in.

# Pull out longitudinal model
Al = Ac[ix,ix];
Bl = Bc[ix,iu];

# ============= TO DO #1 =========
# 1. Define the sensor model and write continuous time dynamics w.r.t z (the augmented state var)
# 2. Discretize the model using code below.
Vbar = xbar[7]; # trim velocity
αbar = xbar[8]; # trim angle of attack

 # Mach, dynamic pressure, static pressure from atmospheric values
Mach,qbar,ps = F16Model.atmos(h0,Vt0);
ρ = 2*qbar/Vt0^2; # calculating density in backwards way

# Matrix of what? TODO figure out what this matrix is doing
M = [cos(αbar) -Vbar*sin(αbar);
     sin(αbar)  Vbar*cos(αbar)];

# Sensor 1 -- Dynamic Pressure
C1 = [0 0 0 0 ρ*Vbar 0 0 ];

# Sensor 2,3 - Accelerometer
C23 = M*Al[5:6,:]; # need to change to corresponding values from Al and Bl
D23 = M*Bl[5:6,:];

# Sensor 4 -- Gyro
C4 = [0 0 0 0 0 0 1];

#sensor 5 - angle of attack sensor
#=
TODO figure out proper value... Assuming 1 for now...
most likely b/c sensor is directly measuring the state used in system
dynamic pressure had to convert pressure perturbation to velocity perturbation
=#
C5 = [0 0 0 0 0 1 0];

#=
 sensor 6 - gps sensor; Provides Npos and Epos
=#
C6 = [1 0 0 0 0 0 0]; # Npos
C7 = [0 1 0 0 0 0 0]; # Epos

#=
sensor 7 - barometer; Assume this outputs altitude directly
=#
C8 = [0 0 1 0 0 0 0]; # altitude

# Combined sensor y = C*x + D*d + n
C = [C6;C7;C8;C1;C23;C4;C5;];
#TODO understand why some are getting the D matrix while others are zero
D = [0;0;0;0;D23;0;0;];

# Add Filter to the disturbance
# F(s) = 1/(s/ωc + 1)
# (s/ωc + 1)do = di
# xfdot = -ωc*xf + ωc*di, do = xf
ωc = 1; #1 rads/
af = -ωc; bf = ωc; cf = 1;

# Augmented System Dynamics
Aa = [Al Bl; zeros(1,7) af]; #size changed b/c added Epos and Npos states
Ba = [zeros(7,1);bf];
Ca = [C D]; # Augmented measurement model to be used in Kalman filtering.

nz = size(Aa,1)
nu = size(Ba,2);

# Convert continuous time dynamics to discrete-time dynamics
S = [Aa Ba;
    zeros(nu,nz) zeros(nu,nu)];

dt = 1/50;
Md = exp(S*dt); # Matrix exponential.

# --- Discrete-time Dynamics ----
Ad = Md[1:nz,1:nz];
Bd =  Md[1:nz,(nz+1):end];

# Check stability
println("Absolute eigen values:\n", abs.(eigen(Ad).values))

#= Propagate uncertainty with AAd and BBd.
Define how disturbance is coming in.
Assume it is coming in from dele.
=#

d2r = pi/180;
#this increases to amount of states we have + 1 (added 2 from Npos and Epos)
xTrue = [0.,0.,0.,0.,0.,0.,0.,0.]; #TODO will Npos and Epos xtrue change over timesteps?

#=
ymeasured = g(x) + n
ybar + ytil = g(x + xtil)
ytil = C*xtil + n -- Model
ytil = ymeasured - ybar; the measurement we feed to Kalman filter.
ybar = g(xbar)
=#

##
μ_post = zeros(nz,1);
Σ_post = zeros(nz,nz);

μ_prior = zeros(nz,1);
Σ_prior = zeros(nz,nz);

 # aoa wil be on magnitude of +-0.08726656 rads (5deg) : |1%| is about .000873
# gps sensor is accurate to w/in 2.3ft TODO is this correct way of implementing?
#SSM need to be some small % of trim values (simulates small errors)
SensorScalingMatrix = Diagonal([2.3,         # Npos, Unsure about value used
                                2.3,         # Epos, Unsure about value used
                                5,           # altitude
                                2,           # qbar scaling
                                1,           # u acceleration
                                1,           # w acceleration
                                0.1,         # pitch rate, q?
                                0.000873,    # Angle of Attack

]);
DisturbanceScalingMatrix = sqrt(5);

Q = DisturbanceScalingMatrix*DisturbanceScalingMatrix'; #  degree^2 covariance
R = SensorScalingMatrix*SensorScalingMatrix'; # Noises of the sensors - diagonal matrix.

ny = size(C,1);

nSteps = 200; # Propagate nSteps
Sig = [];faultlist = []; machlist = []; altlist = [];
basefault = collect(range(0,1,length=nSteps)) .^ 2
for i in 1:nSteps
    μ_prior .= Ad*μ_post; # Propagate μ
    Σ_prior .= Ad*Σ_post*Ad' + Bd*Q*Bd'; # Propagate Σ
    # println(" Diagonal(Σprior): ", Diagonal(Σ_prior).diag);

    # Get sensor data
    n = SensorScalingMatrix*randn(ny); # E(n*n') := R
    d = DisturbanceScalingMatrix*rand(1); # E(d*d') := Q;

    # Simulating the measurement.
    y = Ca*xTrue + n; # Measurement data
    xTrue .= Ad*xTrue + Bd*d; # Dynamics of the true system

    # Compute Kalman gain
    K = Σ_prior*Ca'*inv(Ca*Σ_prior*Ca' + R);

    # Compute posterior with Kalman filteringqbar
    μ_post .= μ_prior + K*(y-Ca*μ_prior);
    Σ_post .= (I(nz) - K*Ca)*Σ_prior;

    #exact state
    altexact = (h0+xTrue[3])*0.3048  #converted ft to meters
    altexact_norm = (altexact - alt_mean)/alt_std #normalizing
    machexact,qbar,ps = F16Model.atmos(h0+xTrue[3],Vt0+xTrue[5]);
    machexact_norm = (machexact - mach_mean)/mach_std #normalizing
    #saving mach and alt steps for graphing later
    push!(machlist, machexact)
    push!(altlist, altexact)

    #kalman state
    altKF = (h0+μ_post[3])*0.3048  #converted feet to meters
    altKF_norm = (altKF - alt_mean)/alt_std  #normalizing
    machKF,qbar,ps = F16Model.atmos(h0+μ_post[3],Vt0+μ_post[5]);
    machKF_norm = (machKF - mach_mean)/mach_std  #normalizing

    #pressure prediction
    global fault = basefault[i]
    global fault_norm = (fault - fault_mean)/fault_std #set to desired fault
    Pstag_norm = pressuremodel([altexact_norm,machexact_norm,fault_norm])[1]
    #fault prediction
    faultpredict_norm = faultmodel([altKF_norm,machKF_norm,Pstag_norm])[1]
    faultpredict = (faultpredict_norm * fault_std) + fault_mean
    push!(faultlist, faultpredict)

    #print(faultpredict,"\n") #should desired fault parameter
    #print(μ_post,"\n")
    # No filtering.
    # μ_post = μ_prior;
    # Σ_post = Σ_prior;

    # Saving the data the open-loop uncertainty propagation.
    # print("μ: ",μ);
    # println(" Diagonal(Σ): ", Diagonal(Σ_post).diag);
    #push!(Sig,Diagonal(Σ_post).diag)
    #
    # # Save the  data with Kalman filtering switch ed on.
end

#TODO: Error Analysis do std and mean of error
#not sure if this is allowed... maybe try different approach to representing error
errorlist = abs.(faultlist .- basefault)
errormean = Statistics.mean(errorlist)
errorstd = Statistics.std(errorlist)
confidence = errormean + 2 * errorstd
print("Error Mean = $errormean\nError std = $errorstd\n")
print("this path predicted the fault parameter to +/-$confidence with a 95% confidence\n")

#creating 3d line plot to display results
plot(machlist,altlist,faultlist,
        title="Fault Prediction over Mach/Altitude Flight Path",
        xlabel="mach",
        ylabel="alt (m)",
        zlabel="Fault Parameter",
        label="Neural Network",
        zlims=[-0.2,1],
        lw=2,
        camera=(15,30),
        legend=:bottomleft)   #actual results
plot!(machlist,altlist,basefault,
        label="Base Fault",
        l2=2)  #base values to compare

plot!(machlist,altlist,zeros(nSteps),
        label="Path Trace")  #trace on bottom

for i in 1:nSteps÷5
    plot!([machlist[5*i],machlist[5*i]],
            [altlist[5*i],altlist[5*i]],
            [basefault[5*i],faultlist[5*i]],
            linecolor=:grey,label=nothing)
    if i == 1
        plot!([machlist[5*i],machlist[5*i]],
                [altlist[5*i],altlist[5*i]],
                [basefault[5*i],faultlist[5*i]],
                linecolor=:grey,label="Difference")
    end
end
current()

max_stars_repo_name: NoahLopezGit/PitotTubeNeuralNetwork

max_issues_repo_name: NoahLopezGit/PitotTubeNeuralNetwork

max_forks_repo_name: NoahLopezGit/PitotTubeNeuralNetwork


text:
# Autogenerated wrapper script for WCS_jll for aarch64-apple-darwin
export libwcs

JLLWrappers.@generate_wrapper_header("WCS")
JLLWrappers.@declare_library_product(libwcs, "@rpath/libwcs.7.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libwcs,
        "lib/libwcs.7.7.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()