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

using Markdown
using InteractiveUtils

# ╔═╡ fbce1e27-45db-465d-83f0-a0a7935abf14
using ActuaryUtilities

# ╔═╡ 2f700532-bc7b-4408-8268-f9f632a75451
using DataFrames

# ╔═╡ a3a5b911-1097-454a-b0ec-463874be92b4
using Setfield

# ╔═╡ ba85bf56-5401-4430-a0e2-eef7e06347a5
md" # Nested Projection Mechanics

A simple example of how one could define a nested projection system. Includes the following examples:
- Outer loop policy projections only
- Outer + Inner loop policy projections with padded cashflows determining reserves and capital
- Outer + Inner loop with a stochastic interest rate for the reserves

"

# ╔═╡ 1f4e28e0-c2f0-4c06-a9f8-d0d71999114a
md" ## Policy Data & Methods"

# ╔═╡ 001d049c-efd2-11ec-0083-cff70f84b072
abstract type Policy end

# ╔═╡ c7dbfd5b-a97e-4628-a2f8-8651f12c60f6
struct Term <: Policy
	inforce::Float64 # count of inforce
	term::Int        # length of benefit period (in months?)
	premium::Float64
	face::Float64
end

# ╔═╡ 50569a1e-5deb-48b6-be72-5be41dcb5b43
md"""
For consistency, use the function signature of:

```
(policy,assumptions,time) -> result
```
"""

# ╔═╡ 280c7ef0-61d5-4d28-87b2-d26855c76b6b
premiums(pol::Term,assumptions,time) = pol.inforce * pol.premium / 12

# ╔═╡ 741777ff-2d40-4393-a0d4-669f3460cc64
qx(pol::Term,assumptions,time)= assumptions.q / 12

# ╔═╡ b8fe0af1-29b7-419a-ad3e-7b120eb5cbe1
deaths(pol::Term,assumptions,time) = pol.inforce * qx(pol,assumptions,time)

# ╔═╡ 35bed858-6b76-4dbb-b49a-505f36b13434
claims(pol::Term,assumptions,time) = deaths(pol,assumptions,time) * pol.face

# ╔═╡ 1f59e516-4dc1-450e-87a4-9ce6ac64029d
md" ## Projection Assumptions"

# ╔═╡ dde19062-8a22-4caa-ae68-23895628dc3e
assumptions = (
	q = 0.012,
	int_reserve = 0.02,
	capital_factor = 0.1, # rate * reserves
)

# ╔═╡ 18e32b95-2d13-4502-8500-48b0c380291c
innerloop_assumption(outer_assump) = @set outer_assump.q *= 1.2

# ╔═╡ 10c9e253-2913-4f85-9989-8b3fffba423e
md"## Projection Logic"

# ╔═╡ ac656984-3bdd-4a33-90be-45042e0080a5
"""
	project(policy,assumptions;start_time=1,additional_processing=res)

The kwarg `additional_processing` defines an intermediate processing step where one can take the current model state and perform additional work, including nested projections. If left out of the arguments, the default for `additional_processing` is `res`, where `res` is (pol,assumpitons, result)->result (ie will just return the model point's results with no additional work being done).

"""
function project(
		pol::Term,
		assumptions;
		start_time=1,
		additional_processing=(pol,assumpitons, result)->result
	)
	A = assumptions
	map(start_time:pol.term) do t
		timestep = t
		premium = premiums(pol,A,t)
		q = qx(pol,A,t)
		death = deaths(pol,A,t)
		claim = claims(pol,A,t)
		net_cf = premium-claim
		inforce = pol.inforce - death
		pol = @set pol.inforce = inforce

		# return a vector of name tuples with the results
		result = (;
			timestep,
			premium,
			death,
			claim,
			net_cf,
			inforce,
			q,
		)

		# apply additional processing function
		additional_processing(pol,A,result)
	end
end

# ╔═╡ 1a903be4-57fb-466c-96e3-f2ad36acc77c
# Function signature: (policy, assumptions, result) -> updated result
function run_inner(policy,assumptions,result)
	A = innerloop_assumption(assumptions)
	p = project(policy,A;start_time=result.timestep+1)

	reserves = -pv(A.int_reserve,[modelpoint.net_cf for modelpoint in p])
	capital = reserves * A.capital_factor
	additional_results = (; reserves, capital)
	return merge(result,additional_results)

end

# ╔═╡ 324be56b-512b-4664-9f5b-943b6d9aeade
# Function signature: (policy, assumptions, result) -> updated result
function run_inner_stochastic(policy,assumptions,result)
	A = innerloop_assumption(assumptions)
	p = project(policy,A;start_time=result.timestep+1)

	# simple stochastic interest rate
	n = 100

	reserves = let
		i = A.int_reserve
		f = pv(i + 0.005 * randn(),[modelpoint.net_cf for modelpoint in p])

		-sum(f for _ in 1:n) / n

	end
	capital = reserves * A.capital_factor
	additional_results = (; reserves, capital)
	return merge(result,additional_results)

end

# ╔═╡ 5c13430a-96c4-447f-89a7-1fe99078feee
md" ## Projections

First, define a sample policy:"

# ╔═╡ be7bbb76-fc0f-40b5-982b-7e46126f9c28
p = Term(1.,120,1300.,100_000.)

# ╔═╡ b5013701-848a-4920-b4f5-360815a9f2e2
md" A projection without any additional processing:"

# ╔═╡ 1d59569f-c56b-45d0-ae57-ead57f68ed7b
project(p,outerloop_assumptions) |> DataFrame

# ╔═╡ 7dd4015f-6b58-4cf2-af3f-30bf3d43d468
md" And an example which uses a PADed inner loop to determine the resserves and capital:"

# ╔═╡ 2e7a438e-435f-4dc1-86c2-947ea8ceb8d4
project(p,assumptions;additional_processing = run_inner) |> DataFrame

# ╔═╡ 9e824c6d-a303-427f-973c-bc1b9658709d
project(p,assumptions;additional_processing = run_inner) |> DataFrame

# ╔═╡ 973a5a58-9b42-4532-be80-21f18fb6791e
project(p,assumptions;additional_processing = run_inner_stochastic) |> DataFrame

# ╔═╡ 98720096-e605-40e9-a32f-d1d5ade04be6
md"## Endnotes

### Further Work

This example is simple, but could be greatly optimized to reduce intermediate variable allocations, refine the timing of cashflows, add additional decrements, handle different types of `Policy`s, abstract some of the projection mechanics into an `Iterable` object, etc.

### Disclaimer
Created as a proof of concept and not indended to be interpreted as a meaningful projection."

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
ActuaryUtilities = "bdd23359-8b1c-4f88-b89b-d11982a786f4"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"

[compat]
ActuaryUtilities = "~3.4.0"
DataFrames = "~1.3.4"
Setfield = "~1.0.0"
"""

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

julia_version = "1.7.0"
manifest_format = "2.0"

[[deps.ActuaryUtilities]]
deps = ["Dates", "Distributions", "ForwardDiff", "QuadGK", "Roots", "StatsBase", "Yields"]
git-tree-sha1 = "4e3ff9120b210afe7b43ad8967c64d36e28f5a69"
uuid = "bdd23359-8b1c-4f88-b89b-d11982a786f4"
version = "3.4.0"

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

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

[[deps.ArrayInterfaceCore]]
deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
git-tree-sha1 = "d618d3cf75e8ed5064670e939289698ecf426c7f"
uuid = "30b0a656-2188-435a-8636-2ec0e6a096e2"
version = "0.1.12"

[[deps.ArrayLayouts]]
deps = ["FillArrays", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "ce5666a4081cfb66035a1a46e45a34ce889bac26"
uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
version = "0.8.7"

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

[[deps.AxisAlgorithms]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "WoodburyMatrices"]
git-tree-sha1 = "66771c8d21c8ff5e3a93379480a2307ac36863f7"
uuid = "13072b0f-2c55-5437-9ae7-d433b7a33950"
version = "1.0.1"

[[deps.BSplineKit]]
deps = ["ArrayLayouts", "BandedMatrices", "FastGaussQuadrature", "Interpolations", "LazyArrays", "LinearAlgebra", "Random", "Reexport", "SparseArrays", "StaticArrays"]
git-tree-sha1 = "651e79138b30580c5a1586b67c3298be6c7b8afe"
uuid = "093aae92-e908-43d7-9660-e50ee39d5a0a"
version = "0.8.4"

[[deps.BandedMatrices]]
deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra", "Random", "SparseArrays"]
git-tree-sha1 = "019aa88766e2493c59cbd0a9955e1bac683ffbcd"
uuid = "aae01518-5342-5314-be14-df237901396f"
version = "0.16.13"

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

[[deps.Bzip2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2"
uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
version = "1.0.8+0"

[[deps.CEnum]]
git-tree-sha1 = "eb4cb44a499229b3b8426dcfb5dd85333951ff90"
uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82"
version = "0.4.2"

[[deps.Calculus]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "f641eb0a4f00c343bbc32346e1217b86f3ce9dad"
uuid = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
version = "0.5.1"

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

[[deps.ChangesOfVariables]]
deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
git-tree-sha1 = "1e315e3f4b0b7ce40feded39c73049692126cf53"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.3"

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

[[deps.ColorVectorSpace]]
deps = ["ColorTypes", "FixedPointNumbers", "LinearAlgebra", "SpecialFunctions", "Statistics", "TensorCore"]
git-tree-sha1 = "d08c20eef1f2cbc6e60fd3612ac4340b89fea322"
uuid = "c3611d14-8923-5661-9e6a-0046d554d3a4"
version = "0.9.9"

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

[[deps.CommonSolve]]
git-tree-sha1 = "332a332c97c7071600984b3c31d9067e1a4e6e25"
uuid = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
version = "0.2.1"

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

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

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

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

[[deps.Contour]]
deps = ["StaticArrays"]
git-tree-sha1 = "9f02045d934dc030edad45944ea80dbd1f0ebea7"
uuid = "d38c429a-6771-53c6-b99e-75d170b6e991"
version = "0.5.7"

[[deps.Crayons]]
git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.1.1"

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

[[deps.DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "daa21eb85147f72e41f6352a57fccea377e310a9"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.3.4"

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

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

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

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

[[deps.DensityInterface]]
deps = ["InverseFunctions", "Test"]
git-tree-sha1 = "80c3e8639e3353e5d2912fb3a1916b8455e2494b"
uuid = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
version = "0.4.0"

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

[[deps.DiffRules]]
deps = ["IrrationalConstants", "LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"]
git-tree-sha1 = "28d605d9a0ac17118fe2c5e9ce0fbb76c3ceb120"
uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
version = "1.11.0"

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

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

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

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

[[deps.DualNumbers]]
deps = ["Calculus", "NaNMath", "SpecialFunctions"]
git-tree-sha1 = "5837a837389fccf076445fce071c8ddaea35a566"
uuid = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"
version = "0.6.8"

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

[[deps.FastGaussQuadrature]]
deps = ["LinearAlgebra", "SpecialFunctions", "StaticArrays"]
git-tree-sha1 = "58d83dd5a78a36205bdfddb82b1bb67682e64487"
uuid = "442a2c76-b920-505d-bb47-c5924d526838"
version = "0.4.9"

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

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

[[deps.FiniteDiff]]
deps = ["ArrayInterfaceCore", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays"]
git-tree-sha1 = "a0700c21266b55bf62c22e75af5668aa7841b500"
uuid = "6a86dc24-6348-571c-b903-95158fe2bd41"
version = "2.12.1"

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

[[deps.Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

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

[[deps.FreeType]]
deps = ["CEnum", "FreeType2_jll"]
git-tree-sha1 = "cabd77ab6a6fdff49bfd24af2ebe76e6e018a2b4"
uuid = "b38be410-82b0-50bf-ab77-7b57e271db43"
version = "4.0.0"

[[deps.FreeType2_jll]]
deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "87eb71354d8ec1a96d4a7636bd57a7347dde3ef9"
uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7"
version = "2.10.4+0"

[[deps.FreeTypeAbstraction]]
deps = ["ColorVectorSpace", "Colors", "FreeType", "GeometryBasics"]
git-tree-sha1 = "b5c7fe9cea653443736d264b85466bad8c574f4a"
uuid = "663a7486-cb36-511b-a19d-713bb74d65c9"
version = "0.9.9"

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

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

[[deps.HypergeometricFunctions]]
deps = ["DualNumbers", "LinearAlgebra", "SpecialFunctions", "Test"]
git-tree-sha1 = "cb7099a0109939f16a4d3b572ba8396b1f6c7c31"
uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
version = "0.3.10"

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

[[deps.Interpolations]]
deps = ["AxisAlgorithms", "ChainRulesCore", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "Requires", "SharedArrays", "SparseArrays", "StaticArrays", "WoodburyMatrices"]
git-tree-sha1 = "b7bc05649af456efc75d178846f47006c2c4c3c7"
uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
version = "0.13.6"

[[deps.InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "b3364212fb5d870f724876ffcd34dd8ec6d98918"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.7"

[[deps.InvertedIndices]]
git-tree-sha1 = "bee5f1ef5bf65df56bdd2e40447590b272a5471f"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.1.0"

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

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

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

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

[[deps.LazyArrays]]
deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra", "MacroTools", "MatrixFactorizations", "SparseArrays", "StaticArrays"]
git-tree-sha1 = "721bebe4d0f8581c18fccf272c62000e22a80a2d"
uuid = "5078a376-72f3-5289-bfd5-ec5146d43c02"
version = "0.22.10"

[[deps.LazyModules]]
git-tree-sha1 = "a560dd966b386ac9ae60bdd3a3d3a326062d3c3e"
uuid = "8cdb02fc-e678-4876-92c5-9defec4f444e"
version = "0.3.1"

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

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

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

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

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

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

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

[[deps.LogExpFunctions]]
deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "09e4b894ce6a976c354a69041a04748180d43637"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.15"

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

[[deps.LsqFit]]
deps = ["Distributions", "ForwardDiff", "LinearAlgebra", "NLSolversBase", "OptimBase", "Random", "StatsBase"]
git-tree-sha1 = "91aa1442e63a77f101aff01dec5a821a17f43922"
uuid = "2fda8390-95c7-5789-9bda-21331edee243"
version = "0.12.1"

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

[[deps.MarchingCubes]]
deps = ["StaticArrays"]
git-tree-sha1 = "3bf4baa9df7d1367168ebf60ed02b0379ea91099"
uuid = "299715c1-40a9-479a-aaf9-4a633d36f717"
version = "0.1.3"

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

[[deps.MatrixFactorizations]]
deps = ["ArrayLayouts", "LinearAlgebra", "Printf", "Random"]
git-tree-sha1 = "2212d36f97e01347adb1460a6914e20f2feee853"
uuid = "a3b82374-2e81-5b9e-98ce-41277c0e4c87"
version = "0.9.1"

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

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

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

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

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

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

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

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

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

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

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

[[deps.Optim]]
deps = ["Compat", "FillArrays", "ForwardDiff", "LineSearches", "LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "PositiveFactorizations", "Printf", "SparseArrays", "StatsBase"]
git-tree-sha1 = "7a28efc8e34d5df89fc87343318b0a8add2c4021"
uuid = "429524aa-4258-5aef-a3af-852621145aeb"
version = "1.7.0"

[[deps.OptimBase]]
deps = ["NLSolversBase", "Printf", "Reexport"]
git-tree-sha1 = "9cb1fee807b599b5f803809e85c81b582d2009d6"
uuid = "87e2bd06-a317-5318-96d9-3ecbac512eee"
version = "2.0.2"

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

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

[[deps.Parameters]]
deps = ["OrderedCollections", "UnPack"]
git-tree-sha1 = "34c0e9ad262e5f7fc75b10a9952ca7692cfc5fbe"
uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
version = "0.12.3"

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

[[deps.PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "a6062fe4063cdafe78f4a0a81cfffb89721b30e7"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.4.2"

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

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "47e5f437cc0e7ef2ce8406ce1e7e24d44915f88d"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.3.0"

[[deps.PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "dfb54c4e414caa595a1f2ed759b160f5a3ddcba5"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.3.1"

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

[[deps.QuadGK]]
deps = ["DataStructures", "LinearAlgebra"]
git-tree-sha1 = "78aadffb3efd2155af139781b8a8df1ef279ea39"
uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
version = "2.4.2"

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

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

[[deps.Ratios]]
deps = ["Requires"]
git-tree-sha1 = "dc84268fe0e3335a62e315a3a7cf2afa7178a734"
uuid = "c84ed2f1-dad5-54f0-aa8e-dbefe2724439"
version = "0.4.3"

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

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.3.0"

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

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

[[deps.Roots]]
deps = ["CommonSolve", "Printf"]
git-tree-sha1 = "06ba8114bf7fc4fd1688e2e4d2259d2000535985"
uuid = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
version = "1.2.0"

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

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

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

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

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

[[deps.SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

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

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

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

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

[[deps.StatsAPI]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "2c11d7290036fe7aac9038ff312d3b3a2a5bf89e"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.4.0"

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

[[deps.StatsFuns]]
deps = ["ChainRulesCore", "HypergeometricFunctions", "InverseFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
git-tree-sha1 = "5783b877201a82fc0014cbf381e7e6eb130473a4"
uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
version = "1.0.1"

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

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

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

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

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

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

[[deps.TensorCore]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "1feb45f88d133a655e001435632f019a9a1bcdb6"
uuid = "62fd8b95-f654-4bbd-a8a5-9c27f68ccd50"
version = "0.1.1"

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

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

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

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

[[deps.UnicodePlots]]
deps = ["ColorTypes", "Contour", "Crayons", "Dates", "FileIO", "FreeTypeAbstraction", "LazyModules", "LinearAlgebra", "MarchingCubes", "NaNMath", "Printf", "SparseArrays", "StaticArrays", "StatsBase", "Unitful"]
git-tree-sha1 = "05fa79b807e11f1a123ab636abcbc6cf4b6c51c4"
uuid = "b8865327-cd53-5732-bb35-84acbb429228"
version = "2.12.3"

[[deps.Unitful]]
deps = ["ConstructionBase", "Dates", "LinearAlgebra", "Random"]
git-tree-sha1 = "b649200e887a487468b71821e2644382699f1b0f"
uuid = "1986cc42-f94f-5a68-af5c-568840ba703d"
version = "1.11.0"

[[deps.WoodburyMatrices]]
deps = ["LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "de67fa59e33ad156a590055375a30b23c40299d3"
uuid = "efce3f68-66dc-5838-9240-27a6d6f5f9b6"
version = "0.5.5"

[[deps.Yields]]
deps = ["BSplineKit", "ForwardDiff", "LinearAlgebra", "LsqFit", "Optim", "Roots", "UnicodePlots"]
git-tree-sha1 = "e82a40304025d049e51defe7e6f9c29c650f5603"
uuid = "d7e99b2f-e7f3-4d9e-9f01-2338fc023ad3"
version = "3.0.1"

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

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

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

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

# ╔═╡ Cell order:
# ╟─ba85bf56-5401-4430-a0e2-eef7e06347a5
# ╠═fbce1e27-45db-465d-83f0-a0a7935abf14
# ╠═2f700532-bc7b-4408-8268-f9f632a75451
# ╠═1f4e28e0-c2f0-4c06-a9f8-d0d71999114a
# ╠═001d049c-efd2-11ec-0083-cff70f84b072
# ╠═c7dbfd5b-a97e-4628-a2f8-8651f12c60f6
# ╠═a3a5b911-1097-454a-b0ec-463874be92b4
# ╟─50569a1e-5deb-48b6-be72-5be41dcb5b43
# ╠═280c7ef0-61d5-4d28-87b2-d26855c76b6b
# ╠═b8fe0af1-29b7-419a-ad3e-7b120eb5cbe1
# ╠═35bed858-6b76-4dbb-b49a-505f36b13434
# ╠═741777ff-2d40-4393-a0d4-669f3460cc64
# ╠═1f59e516-4dc1-450e-87a4-9ce6ac64029d
# ╠═dde19062-8a22-4caa-ae68-23895628dc3e
# ╠═18e32b95-2d13-4502-8500-48b0c380291c
# ╠═10c9e253-2913-4f85-9989-8b3fffba423e
# ╠═ac656984-3bdd-4a33-90be-45042e0080a5
# ╠═1a903be4-57fb-466c-96e3-f2ad36acc77c
# ╠═324be56b-512b-4664-9f5b-943b6d9aeade
# ╟─5c13430a-96c4-447f-89a7-1fe99078feee
# ╠═be7bbb76-fc0f-40b5-982b-7e46126f9c28
# ╟─b5013701-848a-4920-b4f5-360815a9f2e2
# ╠═1d59569f-c56b-45d0-ae57-ead57f68ed7b
# ╠═7dd4015f-6b58-4cf2-af3f-30bf3d43d468
# ╠═2e7a438e-435f-4dc1-86c2-947ea8ceb8d4
# ╠═9e824c6d-a303-427f-973c-bc1b9658709d
# ╠═973a5a58-9b42-4532-be80-21f18fb6791e
# ╟─98720096-e605-40e9-a32f-d1d5ade04be6
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: JuliaActuary/JuliaActuaryLearn

max_issues_repo_name: JuliaActuary/JuliaActuaryLearn

max_forks_repo_name: JuliaActuary/JuliaActuaryLearn


text:
include("benchmark.jl")
include("p0.jl")
include("p1.jl")
include("p2.jl")

ext: jl

max_stars_repo_name: TuringLang/Libtask

max_issues_repo_name: TuringLang/Libtask

max_forks_repo_name: TuringLang/Libtask


text:
export ProbCircuit,
    multiply, summate, ismul, issum,
    num_parameters, num_parameters_node,
    mul_nodes, sum_nodes, params

using LogicCircuits

#####################
# Abstract probabilistic circuit nodes
#####################

"Root of the probabilistic circuit node hierarchy"
abstract type ProbCircuit <: LogicCircuit end

#####################
# node functions that need to be implemented for each type of probabilistic circuit
#####################

"Get the parameters associated with a sum node"
params(n::ProbCircuit) = n.log_probs

import LogicCircuits: children, compile # extend

"Multiply nodes into a single circuit"
function multiply end

"Sum nodes into a single circuit"
function summate end

#####################
# derived functions
#####################

"Is the node a multiplication?"
@inline ismul(n) = GateType(n) isa ⋀Gate
"Is the node a summation?"
@inline issum(n) = GateType(n) isa ⋁Gate

"Count the number of parameters in the circuit"
@inline num_parameters(c::ProbCircuit) =
    sum(n -> num_parameters_node(n), sum_nodes(c))

#####################
# methods to easily construct circuits
#####################

@inline multiply(xs::ProbCircuit...) = multiply(collect(xs))
@inline summate(xs::ProbCircuit...) = summate(collect(xs))

import LogicCircuits: conjoin, disjoin # make available for extension

# alias conjoin/disjoin using mul/sum terminology
@inline conjoin(args::Vector{<:ProbCircuit}; reuse=nothing) =
    multiply(args; reuse)
@inline disjoin(args::Vector{<:ProbCircuit}; reuse=nothing) =
    summate(args; reuse)

@inline Base.:*(x::ProbCircuit, y::ProbCircuit) = conjoin(x,y)
@inline Base.:*(xs::ProbCircuit...) = conjoin(xs...)
@inline Base.:+(x::ProbCircuit, y::ProbCircuit) = disjoin(x,y)
@inline Base.:+(xs::ProbCircuit...) = disjoin(xs...)

compile(::Type{<:ProbCircuit}, ::Bool) =
    error("Probabilistic circuits do not have constant leafs.")

struct WeightProbCircuit
    tmp_weight :: Float64
    circuit :: ProbCircuit
end

@inline Base.:*(w::Real, x::ProbCircuit) = WeightProbCircuit(w, x)
@inline Base.:*(x::ProbCircuit, w::Real) = w * x
@inline Base.:+(x::WeightProbCircuit...) = begin
    ch = collect(x)
    c = map(x -> x.circuit, ch)
    w = map(x -> x.tmp_weight, ch)
    pc = summate(c)
    pc.log_probs .= log.(w)
    pc
end

#####################
# circuit inspection
#####################

"Get the list of multiplication nodes in a given circuit"
mul_nodes(c::ProbCircuit) = ⋀_nodes(c)

"Get the list of summation nodes in a given circuit"
sum_nodes(c::ProbCircuit) = ⋁_nodes(c)

function check_parameter_integrity(circuit::ProbCircuit)
    for node in sum_nodes(circuit)
        @assert all(θ -> !isnan(θ), node.log_probs) "There is a NaN in one of the log_probs"
    end
    true
end


text:
 """
    UnionFind(n)

Classical Union-Find data structure with path compression as described
e.g. in CLRS in Chapter 21 "Disjooint Sets" (3rd edition).

## Example

```jldoctest
julia> import GraphIdx.Utils: UnionFind, find, unite!

julia> u = UnionFind(2);

julia> find(u, 1) == find(u, 2)
false

julia> unite!(u, find(u, 1), find(u, 2));

julia> find(u, 1) == find(u, 2)
true

```
"""
struct UnionFind
    p::Vector{Int}
    rank::Vector{Int}
    UnionFind(n::Int) = new(collect(1:n), zeros(Int, n))
end

function Base.show(io::IO, u::UnionFind)
    println(io, typeof(u), ": n=", length(u.p))
    println(io, "p:    ", u.p)
    println(io, "rank: ", u.rank)
end

function reset!(u::UnionFind)::UnionFind
    for i in 1:length(u)
        u.p[i] = i
    end
    u.rank .= 0
    u
end

"""
Representant of a partition in UnionFind.
"""
struct Rep
    i::Int
end

# make it possible to a[find(u, 4)], where a::Array
Base.to_index(r::Rep) = r.i

Base.getindex(u::UnionFind, x::Int)::Rep = find(u, x)

Base.length(u::UnionFind) = length(u.p)

function find(u::UnionFind, x::Int)::Rep
    if u.p[x] != x
        u.p[x] = find(u, u.p[x]).i
    end
    return Rep(u.p[x])
end

function unite!(u::UnionFind, fx::Rep, fy::Rep)::UnionFind
    let fx::Int = fx.i, fy::Int = fy.i
        if u.rank[fx] > u.rank[fy]
            u.p[fy] = fx
        else
            u.p[fx] = fy
            if u.rank[fx] == u.rank[fy]
                u.rank[fy] += 1
            end
        end
    end
    return u
end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
struct AllocationReplay{A}
    record::A
    step::Base.RefValue{Int}
end

AllocationReplay(record) = AllocationReplay(record, Ref(1))

Cassette.@context ReplayCtx
new_replay_ctx(record) = new_replay_ctx(AllocationReplay(record))
function new_replay_ctx(replay::AllocationReplay)
    reinitialize!(replay)
    return Cassette.disablehooks(ReplayCtx(metadata=replay))
end

"""
    reinitialize!(replay::AllocationReplay)

Return all recorded allocations to original sizes, and reset to first step.
"""
function reinitialize!(replay::AllocationReplay)
    reinitialize!(replay.record)
    replay.step[] = 1
    return replay
end

@inline function next_scheduled_alloc!(replay::AllocationReplay)
    step = replay.step[] :: Int
    alloc = replay.record.allocations[step] :: Array
    replay.step[] = step + 1
    return alloc
end
@inline next_scheduled_alloc!(ctx::ReplayCtx) = next_scheduled_alloc!(ctx.metadata)

@inline function Cassette.overdub(
    ctx::ReplayCtx, ::Type{Array{T,N}}, ::UndefInitializer, dims
)::Array{T,N} where {T,N}
    scheduled = next_scheduled_alloc!(ctx) :: Array{T,N}

    # Commented out until we can workout how to do this without allocations on the happy path
    # It seems like having any branch here makes it allocate
    # TODO: reenable this
    #==
    if  typeof(scheduled) !== Array{T,N} || size(scheduled) !== dims
        @warn "Allocation reuse failed. Indicates value dependent allocations." step=ctx.metadata.step[] expected_T=eltype(scheduled) actual_T=T expected_size=size actual_size=dims
        # Fallback to just doing the allocation
        return Array{T,N}(undef, dims)
    end
    ==#

    return scheduled
end

"""
    avoid_allocations(record, f, args...; kwargs...)

Avoid allocations stored in `record`. See also [`record_allocations`](@ref).
"""
function avoid_allocations(record, f, args...; kwargs...)
    ctx = new_replay_ctx(record)
    return Cassette.overdub(ctx, f, args...; kwargs...)
end


text:
"""
    MD_1C_TOKENS

Dictionary of single-char tokens for Markdown. Note that these characters are
exclusive, they cannot appear again in a larger token.
"""
const MD_1C_TOKENS = LittleDict{Char, Symbol}(
    '{'  => :LXB_OPEN,
    '}'  => :LXB_CLOSE,
    '\n' => :LINE_RETURN,
    )

"""
    MD_TOKENS_LX

Subset of `MD_1C_TOKENS` with only the latex tokens (for parsing what's in a math environment).
"""
const MD_1C_TOKENS_LX = LittleDict{Char, Symbol}(
    '{'  => :LXB_OPEN,
    '}'  => :LXB_CLOSE
    )

"""
    MD_TOKENS

Dictionary of tokens for Markdown. Note that for each, there may be several
possibilities to consider in which case the order is important: the first case
that works will be taken.
"""
const MD_TOKENS = LittleDict{Char, Vector{TokenFinder}}(
    '<'  => [ isexactly("<!--")  => :COMMENT_OPEN,     # <!-- ...
             ],
    '+'  => [ isexactly("+++", ('\n',)) => :MD_DEF_TOML,
             ],
    '-'  => [ isexactly("-->")   => :COMMENT_CLOSE,    #  ... -->
              incrlook(is_hr1)   => :HORIZONTAL_RULE,  # ---+
             ],
    '~'  => [ isexactly("~~~")   => :ESCAPE,           # ~~~  ... ~~~
             ],
    '['  => [ incrlook(is_footnote) => :FOOTNOTE_REF, # [^...](:)? defs will be separated after
             ],
    ']'  => [ isexactly("]: ") => :LINK_DEF,
             ],
    ':'  => [ incrlook(is_emoji) => :CAND_EMOJI,
             ],
    '\\' => [ # -- special characters, see `find_special_chars` in ocblocks
              isexactly("\\\\")       => :CHAR_LINEBREAK,   # --> <br/>
              isexactly("\\", (' ',)) => :CHAR_BACKSPACE,   # --> \
              isexactly("\\*")        => :CHAR_ASTERISK,    # --> *
              isexactly("\\_")        => :CHAR_UNDERSCORE,  # --> _
              isexactly("\\`")        => :CHAR_BACKTICK,    # --> `
              isexactly("\\@")        => :CHAR_ATSIGN,      # --> @
              # -- maths
              isexactly("\\{")        => :INACTIVE,         # See note [^1]
              isexactly("\\}")        => :INACTIVE,         # See note [^1]
              isexactly("\\\$")       => :INACTIVE,         # See note [^1]
              isexactly("\\[")        => :MATH_C_OPEN,      # \[ ...
              isexactly("\\]")        => :MATH_C_CLOSE,     #    ... \]
              # -- latex
              isexactly("\\newenvironment", ('{',)) => :LX_NEWENVIRONMENT,
              isexactly("\\newcommand", ('{',))     => :LX_NEWCOMMAND,
              isexactly("\\begin", ('{',)) => :CAND_LX_BEGIN,
              isexactly("\\end", ('{',))   => :CAND_LX_END,
              incrlook((_, c) -> α(c))     => :LX_COMMAND,  # \command⎵*
             ],
    '@'  => [ isexactly("@def", (' ',))   => :MD_DEF_OPEN,    # @def var = ...
              isexactly("@@", SPACE_CHAR) => :DIV_CLOSE,      # @@⎵*
              incrlook(is_div_open)       => :DIV_OPEN,       # @@dname
             ],
    '#'  => [ isexactly("#",      (' ',)) => :H1_OPEN, # see note [^2]
              isexactly("##",     (' ',)) => :H2_OPEN,
              isexactly("###",    (' ',)) => :H3_OPEN,
              isexactly("####",   (' ',)) => :H4_OPEN,
              isexactly("#####",  (' ',)) => :H5_OPEN,
              isexactly("######", (' ',)) => :H6_OPEN,
             ],
    '&'  => [ incrlook(is_html_entity) => :CHAR_HTML_ENTITY,
             ],
    '$'  => [ isexactly("\$", ('$',), false) => :MATH_A,  # $⎵*
              isexactly("\$\$")              => :MATH_B,  # $$⎵*
             ],
    '_'  => [ isexactly("_\$>_") => :MATH_I_OPEN,  # internal use when resolving a latex command
              isexactly("_\$<_") => :MATH_I_CLOSE, # within mathenv (e.g. \R <> \mathbb R)
              incrlook(is_hr2)   => :HORIZONTAL_RULE,
             ],
    '`'  => [ isexactly("`", ('`',), false)  => :CODE_SINGLE, # `⎵
              isexactly("``",('`',), false)  => :CODE_DOUBLE, # ``⎵*
              isexactly("```", SPACE_CHAR)   => :CODE_TRIPLE, # ```⎵*
              isexactly("`````", SPACE_CHAR) => :CODE_PENTA,  # `````⎵*
              is_language()                  => :CODE_LANG,   # ```lang*
              is_language2()                 => :CODE_LANG2,  # `````lang*
             ],
    '*'  => [ incrlook(is_hr3)   => :HORIZONTAL_RULE,
             ]
    ) # end dict
#= NOTE
[1] capturing \{ here will force the head to move after it thereby not
marking it as a potential open brace, same for the close brace.
[2] similar to @def except that it must be at the start of the line. =#

"""
    MD_TOKENS_LX

Subset of `MD_TOKENS` with only the latex tokens (for parsing what's in a math environment).
"""
const MD_TOKENS_LX = Dict{Char, Vector{TokenFinder}}(
    '\\' => [
        isexactly("\\{")         => :INACTIVE,
        isexactly("\\}")         => :INACTIVE,
        incrlook((_, c) -> α(c)) => :LX_COMMAND ]
    )

"""
    L_RETURNS

Convenience tuple containing the name for standard line returns and line
returns followed by an indentation (either a quadruple space or a tab).
"""
const L_RETURNS = (:LINE_RETURN, :LR_INDENT, :EOS)

"""
    MD_OCB

List of Open-Close Blocks whose content should be deactivated (any token within
their span should be marked as inactive) until further processing.
The keys are identifier for the type of block, the value is a pair with the
opening and closing tokens followed by a boolean indicating whether the content
of the block should be reprocessed.
The only `OCBlock` not in this dictionary is the brace block since it should
not deactivate its content which is needed to find latex definitions
(see parser/markdown/find_blocks/find_lxdefs).
"""
const MD_OCB = [
    # name                    opening token   closing token(s)
    # ---------------------------------------------------------------------
    OCProto(:COMMENT,         :COMMENT_OPEN, (:COMMENT_CLOSE,)),
    OCProto(:MD_DEF_BLOCK,    :MD_DEF_TOML,  (:MD_DEF_TOML,)  ),
    OCProto(:CODE_BLOCK_LANG, :CODE_LANG,    (:CODE_TRIPLE,)  ),
    OCProto(:CODE_BLOCK_LANG, :CODE_LANG2,   (:CODE_PENTA,)   ),
    OCProto(:CODE_BLOCK,      :CODE_TRIPLE,  (:CODE_TRIPLE,)  ),
    OCProto(:CODE_BLOCK,      :CODE_PENTA,   (:CODE_PENTA,)   ),
    OCProto(:CODE_INLINE,     :CODE_DOUBLE,  (:CODE_DOUBLE,)  ),
    OCProto(:CODE_INLINE,     :CODE_SINGLE,  (:CODE_SINGLE,)  ),
    OCProto(:MD_DEF,          :MD_DEF_OPEN,  L_RETURNS        ), # [^4]
    OCProto(:CODE_BLOCK_IND,  :LR_INDENT,    (:LINE_RETURN,)  ),
    OCProto(:ESCAPE,          :ESCAPE,       (:ESCAPE,)       ),
    OCProto(:FOOTNOTE_DEF,    :FOOTNOTE_DEF, L_RETURNS        ),
    OCProto(:LINK_DEF,        :LINK_DEF,     L_RETURNS        ),
    # ------------------------------------------------------------------
    OCProto(:H1,              :H1_OPEN,      L_RETURNS), # see [^3]
    OCProto(:H2,              :H2_OPEN,      L_RETURNS),
    OCProto(:H3,              :H3_OPEN,      L_RETURNS),
    OCProto(:H4,              :H4_OPEN,      L_RETURNS),
    OCProto(:H5,              :H5_OPEN,      L_RETURNS),
    OCProto(:H6,              :H6_OPEN,      L_RETURNS)
    ]
# the split is due to double brace blocks being allowed in markdown
const MD_OCB2 = [
    OCProto(:LXB,             :LXB_OPEN,     (:LXB_CLOSE,), nestable=true),
    OCProto(:DIV,             :DIV_OPEN,     (:DIV_CLOSE,), nestable=true),
    ]
#= NOTE:
* [3] a header can be closed by either a line return or an end of string (for
instance in the case where a user defines a latex command like so:
\newcommand{\section}{# blah} (no line return).)
* [4] MD_DEF take precedence over CODE_IND, note that if you have an indented
* block with @def in it, things may go bad.
* ordering matters!
=#

"""
    MD_HEADER

All header symbols.
"""
const MD_HEADER = (:H1, :H2, :H3, :H4, :H5, :H6)

"""
    MD_HEADER_OPEN

All header symbols (opening).
"""
const MD_HEADER_OPEN = (:H1_OPEN, :H2_OPEN, :H3_OPEN, :H4_OPEN, :H5_OPEN, :H6_OPEN)

"""
    MD_OCB_ESC

Blocks that will be escaped (tokens in their span will be ignored on the
current parsing round).
"""
const MD_OCB_ESC = [e.name for e ∈ MD_OCB if !e.nest]

"""
    MD_OCB_MATH

Same concept as `MD_OCB` but for math blocks, they can't be nested. Separating
them from the other dictionary makes their processing easier.
Dev note: order does not matter.
"""
const MD_OCB_MATH = [
    OCProto(:MATH_A,     :MATH_A,          (:MATH_A,)          ),
    OCProto(:MATH_B,     :MATH_B,          (:MATH_B,)          ),
    OCProto(:MATH_C,     :MATH_C_OPEN,     (:MATH_C_CLOSE,)    ),
    OCProto(:MATH_I,     :MATH_I_OPEN,     (:MATH_I_CLOSE,)    ),
    ]

"""
    MD_OCB_ALL

Combination of all `MD_OCB` in order.
DEV: only really used in tests.
"""
const MD_OCB_ALL = vcat(MD_OCB, MD_OCB2, MD_OCB_MATH)

"""
    MD_OCB_IGNORE

List of names of blocks that will need to be dropped at compile time.
"""
const MD_OCB_IGNORE = (:COMMENT, :MD_DEF)

"""
    MATH_DISPLAY_BLOCKS_NAMES

List of names of maths environments (display mode).
"""
const MATH_DISPLAY_BLOCKS_NAMES = collect(e.name for e ∈ MD_OCB_MATH if e.name != :MATH_A)

"""
    MATH_BLOCKS_NAMES

List of names of all maths environments.
"""
const MATH_BLOCKS_NAMES = tuple(:MATH_A, MATH_DISPLAY_BLOCKS_NAMES...)

"""
CODE_BLOCKS_NAMES

List of names of code blocks environments.
"""
const CODE_BLOCKS_NAMES = (:CODE_BLOCK_LANG, :CODE_BLOCK, :CODE_BLOCK_IND)

"""
    MD_CLOSEP

Blocks which, upon insertion, should close any open paragraph.
Order doesn't matter.
"""
const MD_CLOSEP = [MD_HEADER..., :DIV, CODE_BLOCKS_NAMES..., MATH_DISPLAY_BLOCKS_NAMES...]

"""
    MD_OCB_NO_INNER

List of names of blocks which will deactivate any block contained within them
as their content will be reprocessed later on.
See [`find_all_ocblocks`](@ref).
"""
const MD_OCB_NO_INNER = vcat(MD_OCB_ESC, MATH_BLOCKS_NAMES, :LXB, MD_HEADER)


text:
using Test, ParameterizedNotebooks

nb = ParameterizedNotebook("../examples/example.ipynb")
@test nb.parameters == [:seed]
@test_throws Exception nb()
@test [nb(;seed=seed) for seed in 1:4] == [1,4,9,16]


text:
@testset "1018.binary-prefix-divisible-by-5.jl" begin
    @test prefixes_div_by5([0, 1, 1]) == [true, false, false]
    @test prefixes_div_by5([1, 1, 1]) == [false, false, false]
    @test prefixes_div_by5([0, 1, 1, 1, 1, 1]) == [true, false, false, false, true, false]
    @test prefixes_div_by5([1, 1, 1, 0, 1]) == [false, false, false, false, false]
end


text:
# This file compares approaching multi-agent models in two ways:
# 1) using different Types to represent different agents. This leads to type
# stability within `step!`.
# 2) using a single type with some extra property `type` or `kind`, and then do
# an `if`-based branching on this type to dispatch to different functions.

# The result is that (2) is much faster.

################### DEFINITION 1 ###############
using Agents, Random, BenchmarkTools

mutable struct GridAgentOne <: AbstractAgent
    id::Int
    pos::Dims{2}
    one::Float64
    two::Bool
end

mutable struct GridAgentTwo <: AbstractAgent
    id::Int
    pos::Dims{2}
    one::Float64
    two::Bool
end

mutable struct GridAgentThree <: AbstractAgent
    id::Int
    pos::Dims{2}
    one::Float64
    two::Bool
end

mutable struct GridAgentFour <: AbstractAgent
    id::Int
    pos::Dims{2}
    one::Float64
    two::Bool
end

mutable struct GridAgentFive <: AbstractAgent
    id::Int
    pos::Dims{2}
    one::Float64
    two::Bool
end

model1 = ABM(
    Union{GridAgentOne,GridAgentTwo,GridAgentThree,GridAgentFour,GridAgentFive},
    GridSpace((15, 15));
    warn = false,
    rng = MersenneTwister(42),
    scheduler = schedule_randomly,
)

for _ in 1:50
    a = GridAgentOne(nextid(model1), (1,1), rand(model1.rng), rand(model1.rng, Bool))
    add_agent!(a, model1)
    a = GridAgentTwo(nextid(model1), (1,1), rand(model1.rng), rand(model1.rng, Bool))
    add_agent!(a, model1)
    a = GridAgentThree(nextid(model1), (1,1), rand(model1.rng), rand(model1.rng, Bool))
    add_agent!(a, model1)
    a = GridAgentFour(nextid(model1), (1,1), rand(model1.rng), rand(model1.rng, Bool))
    add_agent!(a, model1)
    a = GridAgentFive(nextid(model1), (1,1), rand(model1.rng), rand(model1.rng, Bool))
    add_agent!(a, model1)
end

agent_step!(agent::GridAgentOne, model1) = walk!(agent, rand, model1)
function agent_step!(agent::GridAgentTwo, model1)
    agent.one += rand(model1.rng)
    agent.two = rand(model1.rng, Bool)
end
function agent_step!(agent::GridAgentThree, model1)
    if any(a-> a isa GridAgentTwo, nearby_agents(agent, model1))
        agent.two = true
        walk!(agent, rand, model1)
    end
end
function agent_step!(agent::GridAgentFour, model1)
    agent.one += sum(a.one for a in nearby_agents(agent, model1))
end
function agent_step!(agent::GridAgentFive, model1)
    targets = filter!(a->a.one > 0.8, collect(nearby_agents(agent, model1, 3)))
    idx = argmax(map(t->edistance(agent, t, model1), targets))
    farthest = targets[idx]
    walk!(agent, sign.(farthest.pos .- agent.pos), model1)
end

################### DEFINITION 2 ###############

mutable struct GridAgentAll <: AbstractAgent
    id::Int
    pos::Dims{2}
    one::Float64
    two::Bool
    type::Symbol
end

model2 = ABM(
    GridAgentAll,
    GridSpace((15, 15));
    rng = MersenneTwister(42),
    scheduler = schedule_randomly,
)

for _ in 1:50
    a = GridAgentAll(nextid(model2), (1,1), rand(model2.rng), rand(model2.rng, Bool), :one)
    add_agent!(a, model2)
    a = GridAgentAll(nextid(model2), (1,1), rand(model2.rng), rand(model2.rng, Bool), :two)
    add_agent!(a, model2)
    a = GridAgentAll(nextid(model2), (1,1), rand(model2.rng), rand(model2.rng, Bool), :three)
    add_agent!(a, model2)
    a = GridAgentAll(nextid(model2), (1,1), rand(model2.rng), rand(model2.rng, Bool), :four)
    add_agent!(a, model2)
    a = GridAgentAll(nextid(model2), (1,1), rand(model2.rng), rand(model2.rng, Bool), :five)
    add_agent!(a, model2)
end

function agent_step!(agent::GridAgentAll, model2)
    if agent.type == :one
        agent_step_one!(agent, model2)
    elseif agent.type == :two
        agent_step_two!(agent, model2)
    elseif agent.type == :three
        agent_step_three!(agent, model2)
    elseif agent.type == :four
        agent_step_four!(agent, model2)
    else
        agent_step_five!(agent, model2)
    end
end
agent_step_one!(agent, model2) = walk!(agent, rand, model2)
function agent_step_two!(agent, model2)
    agent.one += rand(model2.rng)
    agent.two = rand(model2.rng, Bool)
end
function agent_step_three!(agent, model2)
    if any(a-> a.type == :two, nearby_agents(agent, model2))
        agent.two = true
        walk!(agent, rand, model2)
    end
end
function agent_step_four!(agent, model2)
    agent.one += sum(a.one for a in nearby_agents(agent, model2))
end
function agent_step_five!(agent, model2)
    targets = filter!(a->a.one > 1.0, collect(nearby_agents(agent, model2, 3)))
    if !isempty(targets)
        idx = argmax(map(t->edistance(agent, t, model2), targets))
        farthest = targets[idx]
        walk!(agent, sign.(farthest.pos .- agent.pos), model2)
    end
end

################### Benchmarks ###############

@btime step!($model1, agent_step!, dummystep, 500)
@btime step!($model2, agent_step!, dummystep, 500)

# Results:
# 718.589 ms (11581560 allocations: 704.26 MiB)
# 141.673 ms (2292318  allocations: 149.54 MiB)


text:
#!/usr/bin/env julia

# http://www.phpsimplex.com/en/simplex_method_example.htm

using Base.Test
using MathProgBase
using Logging
using Drs

A = Float64[2 1; 2 3; 3 1]
b = Float64[18, 42, 24]
c = Float64[-3, -2]

s = linprog(c, A, '<', b, -Inf, Inf, DrsMathProgSolver())

@test s.status == :Optimal
@test s.objval == -33
@test s.sol == [3, 12]

max_stars_repo_name: abogushov/lplib

max_issues_repo_name: abogushov/lplib

max_forks_repo_name: abogushov/lplib


text:
#!/usr/bin/env julia
using GraphIdx.Io.Dimacs10: open_dimacs10, parse_dimacs10_edges
using HDF5: h5open, attrs

fname = get(ARGS, 1, "belgium.bz2")
out = get(ARGS, 2, "out.h5")
output_bin = false

@time n, m, head, tail = open_dimacs10(parse_dimacs10_edges, fname)

isfile(out) && rm(out)
h5open(out, "w") do io
    attrs(io)["n"] = n
    attrs(io)["m"] = m
    write(io, "head", head .-1)
    write(io, "tail", tail .-1)
end

output_bin && open("bin", "w") do io
    write(io, head .-1);
    write(io, tail .-1)
end

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
"""
Tests for CR3BP dynamics.
"""
module CR3BPTests

    using ModelingToolkit, AstrodynamicalModels, Test

    @testset "CR3BP Model Constructors" begin
        model = CR3BP(; stm=false, structural_simplify=false)
        model = CR3BP(; stm=true, structural_simplify=true)
        @test true
    end

    @testset "CR3BP Model Calculations" begin
        vectorfield = CR3BPFunction()

        r = [0.9253021269565836, 0, 0]
        v = [0.0, 0.05852663414965813, 0.0]
        μ = 0.0009536838895767625

        @test isapprox(
                vectorfield(vcat(r,v), [μ], NaN),
                [0.0, 0.05852663414965813,  0.0, 0.053265045303684255, 0.0, -0.0];
                atol=1e-8
        )
    end

end


text:
using TextFileTools
using Test

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


text:
#--------------------------helper functions------------------------------------
using Statistics
function clean_data(data)

    t = data[:,1]
    Pa1 = data[:,2]
    Pa2 = data[:,3]

    Pa_se = std(hcat(Pa1,Pa2),dims=2)
    n1 = Pa2n(Pa1)
    n2 = Pa2n(Pa2)
    n_mean = mean(hcat(n1,n2),dims=2)
    n_se = std(hcat(n1,n2),dims=2)

    tstacked = vcat(t,t)
    Pa = vcat(Pa1,Pa2)

    p = sortperm(tstacked)
    tstacked = tstacked[p]
    Pa = Pa[p]
    return tstacked, Pa, hcat(t,n_mean,n_se)
end

function generate_data(sol)
  t = collect(range(0,stop=sol.t[end],length=101))

  solution = sum(sol(t),dims=1)
  # randomized = VectorOfArray([(sol(t[i]) + .02randn(size(sol,1))) for i in 1:length(t)])
  # dataset = convert(Array,randomized)
  randomized = solution .+ .01randn(size(solution))
  dataset = randomized
  display(scatter!(t,dataset'))
  return dataset
end

function generate_data(sol,t)

  solution = sum(sol(t),dims=1)
  # randomized = VectorOfArray([(sol(t[i]) + .02randn(size(sol,1))) for i in 1:length(t)])
  # dataset = convert(Array,randomized)
  randomized = solution .+ .01randn(size(solution))
  dataset = randomized
  display(scatter!(t,dataset'))
  return dataset
end

function n2Pa(n) # convert adhesion frequency to <n>
    Pa = 1 .- exp.(-n)
end

function Pa2n(Pa) # convert <n> to adhesion frequency
    n = log.(1 ./ (1 .- Pa))
end

using CSV
function clean_dissociation_data(data)
    n = size(data,2) ÷ 4 # 4 columns per force bin
    clean_data = []
    for i in 0:n-1
        ind_bin = i*4+1:i*4+4
        data_bin = data[:,ind_bin]
        mask_non_zeros = data_bin[:,1] .!= 0
        if sum(mask_non_zeros) > 1
             # clean_data[i+1] = data_bin[mask_non_zeros,:]
             push!(clean_data,data_bin[mask_non_zeros,:])
        end
    end
    return clean_data
end

function separate_diss_data(data)
    t = [bin[:,2] for bin in data]
    lnp = [bin[:,4] for bin in data]
    return t, lnp
end

function rescale_t(data)
    rescaled = data .- minimum(data)
end

function plot_dissociations(data)
    lb = 0
    for i in 1:length(data)
        F = mean(data[i][:,1])
        h = scatter!(Array(data[i][:,2]),Array(data[i][:,4]),
                        color=i,
                        markeralpha=0.5,
                        markerstrokealpha=1,
                        markersize=5,
                        label = string(F),
                        palette = :Dark2_5)
        if data[i][end,4] < lb
            lb = data[i][end,4]
        end
        # ylims! = (-10,0.)
        # display(h)
    end
end

function plot_F_vs_t(data)
    D = Ft_average_and_error(data)
    h = scatter!(D.f_mean,D.t_mean,
                    yerror=D.t_sem,
                    xerror=D.f_sem,
                    c=Color,
                    label="",
                    markerstrokecolor=:auto,
                    markersize=7)
    display(h)
end

function plot_F_vs_t(f,t)
    D = Ft_average_and_error(f,t)
    h = scatter!(D.f_mean,D.t_mean,
                    yerror=D.t_sem,
                    xerror=D.f_sem,
                    c=Color,
                    label="",
                    markerstrokecolor=:auto,
                    markersize=7)
end

function Ft_average_and_error(f,t)
    t_mean = []
    f_mean = []
    t_sem = []
    f_sem = []
    for (fi,ti) = zip(f,t)
        n = length(fi[:,1])
        push!(t_mean, mean(ti))
        push!(f_mean, mean(fi))
        push!(t_sem, std(ti)/sqrt(n) )
        push!(f_sem, std(fi)/sqrt(n) )
    end
    df = DataFrame(f_mean=f_mean,f_sem=f_sem,
                    t_mean=t_mean,t_sem=t_sem)
end
function Ft_average_and_error(data)
    t_mean = []
    f_mean = []
    t_sem = []
    f_sem = []
    for bin in data
        n = length(bin[:,1])
        push!(t_mean, mean(bin[:,2]))
        push!(f_mean, mean(bin[:,1]))
        push!(t_sem, std( bin[:,2])/sqrt(n) )
        push!(f_sem, std(bin[:,1])/sqrt(n) )
    end
    df = DataFrame(f_mean=f_mean,f_sem=f_sem,
                    t_mean=t_mean,t_sem=t_sem)
end

using Statistics
function extract_avg_forces(data)
    avg_fs = [mean(bin[:,1]) for bin in data]
end

function extract_max_lifetimes(data)
    max_t = [maximum(bin[:,2]) for bin in data]
end

# using Interpolations, Loess, ImageFiltering
#
# function max_lifetime_bounds(f,t,xs)
#     max_f = maximum([xs[end], f[end]])
#     fs = vcat(0., f, max_f)
#     ts = vcat(t[1], t, t[end])
#     # fs = vcat(f, max_f)
#     # ts = vcat(t, t[end])
#     interp_linear = LinearInterpolation(fs, ts)
#     ys = interp_linear(xs)
#     # h = plot!(xs,ys)
#
#     ker = ImageFiltering.Kernel.gaussian((3,))
#     smoothed_ys = imfilter(ys, ker)
#     # smoothed_ys = moving_average(ys,20)
#     # model = loess(collect(xs), ys)
#     # smoothed_ys = predict(model, xs)
#     # h = plot!(xs,smoothed_ys)
#     # display(h)
#     return smoothed_ys
# end

function lnp2p(lnp)
    p = exp.(lnp)
end

function plot_bell_equations(p)
    f = collect(0.0:1:35)
    plt = []
    if length(p) == 8
        y1 = log10.(bell_diss_kT.(p[1],p[2],f))
        y2 = log10.(bell_diss_kT.(p[3],p[4],f))
        y3 = log10.(bell_diss_kT.(p[5],p[6],f))
        y4 = log10.(bell_ass_kT.(p[7],p[8],f))
        # ytot = y1+y2
        rates = [y1,y2,y3,y4]
        plt = plot(f,rates,
                    ylims=(minimum(minimum.(rates)),maximum(maximum.(rates))),
                    labels=["k₋f" "k₋s" "ka" "k₋a"])
    elseif length(p) == 4
        y1 = bell_ass_kT.(p[1],p[2],f)
        y2 = bell_diss_kT.(p[3],p[4],f)
        # ytot = y1+y2
        plt = plot(f,[y1,y2],ylims=(0,30))
    end
    plt
end

function plot_diss_fit(model,Fs; generate_data = false)
    for (i,F) in enumerate(Fs)
        # Setup the ODE problem, then solve
        f = (du,u,p,t) -> model(du,u,p,t,dens,F)
        prob = ODEProblem(f,u₀,tspan,p)
        sol = solve(prob,AutoVern7(KenCarp4()),abstol=1e-8,reltol=1e-8)

        Σ = Complex.(sum(Array(sol),dims=1))
        # h = plot!(sol.t,real.(Σ)',ylims=(-5,0))
        h = plot!(sol.t,real.(log.(Σ)'),
            ylims=(-5.5,0),
            xlims=(0,6),
            color=i,
            label=string(F),
            palette = :Dark2_5)

        if generate_data
            push!(dataset,vec(generate_data(sol;logplot=true)))
        end

        # display(h)
    end
end

function show_diss_fit(model, data, forces)
  h = plot()
  plot_diss_fit(model,forces)
  plot_dissociations(data)
  return h
end
function show_diss_fit(model, t, lnp, forces)
    h = plot()
    plot_diss_fit(model,forces)
    scatter!(t, lnp)
    return h
end

using NumericalIntegration
function plot_force_lifetime_fit(model,Fs,data)

    mean_f  = extract_avg_forces(data)
    max_t  = extract_max_lifetimes(data)
    # mean_f = mean.(f)
    # max_t = maximum.(t)
    # scatter(mean_f,max_t)
    t_bounds = max_lifetime_bounds(mean_f,max_t,Fs)

    avg_lifetimes = zeros(length(Fs))
    for (i,F) in enumerate(Fs)
        # Setup the ODE problem, then solve
        f = (du,u,p,t) -> model(du,u,p,t,dens,F)
        prob = ODEProblem(f,u₀,tspan,p)
        sol = solve(prob,AutoVern7(KenCarp4()),abstol=1e-8,reltol=1e-8)

        t_query = collect(range(0.,stop=t_bounds[i],length=100))
        sol_at_bound_t = sol(t_query)

        Σ = real.(Complex.(sum(Array(sol_at_bound_t),dims=1)))
        # Σ_sol = real.(Complex.(sum(sol_at_bound_t)))
        # h = plot!(sol.t,real.(Σ)',ylims=(-5,0))
        # Σ = vec(Σ_sol)
        # t = sol.t[sol.t .<= t_bounds[i]]
        # Σ = Σ[sol.t .<= t_bounds[i]]
        avg_lifetime = integrate(t_query,vec(Σ))
        avg_lifetimes[i] = avg_lifetime
    end
    plot!(Fs,avg_lifetimes)
                # ylims=(0,Inf),
                # xlims=(0,Inf))
end

# given vector of lifetimes > compute lnp
function calculate_lnp(t)
    lnp = log.( (length(t):-1:1) ./ length(t) )
end

function measurements2Df(t,y)
    y_value = Measurements.value.(y)
    y_error = Measurements.uncertainty.(y)
    df = DataFrame(x = t, y = y_value, error = y_error)
end

function combine_bins(data)
    F = []
    t = []
    for bin in data
        F = vcat(F,bin[:,1])
        t = vcat(t,bin[:,2])
    end
    perm_id = sortperm(F)
    F = F[perm_id]
    t = t[perm_id]
    return Float64.(F), Float64.(t)
end


text:
struct NRTLParam <: EoSParam
    a::PairParam{Float64}
    b::PairParam{Float64}
    c::PairParam{Float64}
    Mw::SingleParam{Float64}
end

abstract type NRTLModel <: ActivityModel end

struct NRTL{c<:EoSModel} <: NRTLModel
    components::Array{String,1}
    icomponents::UnitRange{Int}
    params::NRTLParam
    puremodel::Vector{c}
    absolutetolerance::Float64
    references::Array{String,1}
end

@registermodel NRTL

export NRTL

function NRTL(components::Vector{String}; puremodel=PR,
    userlocations=String[],
     verbose=false)
    params = getparams(components, ["properties/critical.csv", "properties/molarmass.csv","Activity/NRTL/NRTL_unlike.csv"]; userlocations=userlocations, asymmetricparams=["a","b"], ignore_missing_singleparams=["a","b"], verbose=verbose)
    a  = params["a"]
    b  = params["b"]
    c  = params["c"]
    Mw  = params["Mw"]
    icomponents = 1:length(components)

    init_puremodel = [puremodel([components[i]]) for i in icomponents]
    packagedparams = NRTLParam(a,b,c,Mw)
    references = String[]
    model = NRTL(components,icomponents,packagedparams,init_puremodel,1e-12,references)
    return model
end

function activity_coefficient(model::NRTLModel,p,T,z)
    a = model.params.a.values
    b = model.params.b.values
    c = model.params.c.values

    x = z ./ sum(z)

    τ = @. a+b/T
    G = @. exp(-c*τ)
    lnγ = sum(x[j]*τ[j,:].*G[j,:] for j ∈ @comps)./sum(x[k]*G[k,:] for k ∈ @comps)+sum(x[j]*G[:,j]/sum(x[k]*G[k,j] for k ∈ @comps).*(τ[:,j] .-sum(x[m]*τ[m,j]*G[m,j] for m ∈ @comps)/sum(x[k]*G[k,j] for k ∈ @comps)) for j in @comps)
    return exp.(lnγ)
end


text:
#! /usr/bin/julia

# Rosetta Code, Command-line arguments

prog = basename(Base.source_path())

println(prog, "'s command-line arguments are:")
for s in ARGS
    println("    ", s)
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
push!(LOAD_PATH, "../src/")

using Documenter, GaussianProcesses

makedocs(
    modules = [GaussianProcesses],
    format = Documenter.HTML(),
    sitename = "GaussianProcesses.jl",
         pages = ["Introduction" => "index.md",
                     "Basic usage" => ["Regression.md", "plotting_gps.md"],
                     "Tutorials" => ["classification_example.md", "sparse_example.md", "mauna_loa.md", "poisson_regression.md"],
                     "Reference" => ["gp.md", "kernels.md", "mean.md","lik.md","optimization.md","sparse.md","crossvalidation.md"]
               ]
)

deploydocs(repo = "github.com/STOR-i/GaussianProcesses.jl.git")


text:
# Script to automatically insert Markdown footnotes for all [#xxxx] issue
# cross-references in the NEWS file.

NEWS = get(ARGS, 1, "NEWS.md")

s = read(NEWS, String)

m = match(r"\[#[0-9]+\]:", s)
if m !== nothing
    s = s[1:m.offset-1]
end

footnote(n) = "[#$n]: https://github.com/JuliaLang/julia/issues/$n"
N = map(m -> parse(Int,m.captures[1]), eachmatch(r"\[#([0-9]+)\]", s))
foots = join(map(footnote, sort!(unique(N))), "\n")

open(NEWS, "w") do f
    println(f, s, foots)
end

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
include("lm_util.jl")
# language model implementation
function initmodel(atype, hiddens, charhidden, chembedding, wordvocab, charvocab, init=xavier)
    model = Dict{Symbol, Any}()
    wordembedding = charhidden[1]

    model[:forw] = initweights(atype, hiddens, wordembedding, init)
    model[:back] = initweights(atype, hiddens, wordembedding, init)
    model[:char] = initweights(atype, charhidden, chembedding, init)
    model[:cembed] = atype(init(charvocab, chembedding))
    model[:soft] = [ atype(init(2hiddens[end], wordvocab)), atype(zeros(1, wordvocab)) ]
    return model
end

function chlstm(weight, bias, hidden, cell, input; mask=nothing)
    gates   = hcat(input,hidden) * weight .+ bias
    hsize   = size(hidden,2)
    forget  = sigm(gates[:,1:hsize])
    ingate  = sigm(gates[:,1+hsize:2hsize])
    outgate = sigm(gates[:,1+2hsize:3hsize])
    change  = tanh(gates[:,1+3hsize:end])
    cell    = cell .* forget + ingate .* change
    hidden  = outgate .* tanh(cell)

    # masking operation
    cell    = cell .* mask
    hidden  = hidden .* mask
    return (hidden,cell)
end

"""
Multilayer lstm forward function in single time step
"""
function chforw(weight, states, input; mask=nothing)
    x = input
    for i=1:2:length(states)
        (states[i], states[i+1]) = chlstm(weight[i], weight[i+1], states[i], states[i+1], x; mask=mask)
        x = states[i]
    end
    return x
end

"""
mchar, mcembed, states: model's character lstm, model's character embeddings, char-lstm states respectively.
wids, i2w_all, ch: batch of words, index to word array and character vocabulary respectively.
"""
function charembed(mchar, mcembed, states, wids, i2w_all, ch, atype)
    schar = copy(states)

    (data, masks) = charlup3(wids, i2w_all, ch)
    h = similar(schar[1])
    for (c, m) in zip(data, masks)
        embed = mcembed[c, :]
        mbon = convert(atype, m)
        h = chforw(mchar, schar, embed; mask=mbon)
    end
    return h
end

function charbilstm(model, chstates, states, sequence, i2w_all, chvocab, lval=[])
    total = 0.0
    count = 0
    atype = typeof(states[1])

    # extract the embeddings on character reading
    embeddings = Array(Any, length(sequence))
    for i=1:length(sequence)
        embeddings[i] = charembed(model[:char], model[:cembed], chstates, sequence[i], i2w_all, chvocab, atype)
    end

    # forward lstm
    fhiddens = Array(Any, length(sequence)-2)
    sf = copy(states)
    for i=1:length(sequence)-2
        x = embeddings[i]
        h = forward(model[:forw], sf, x)
        fhiddens[i] = copy(h)
    end

    # backward lstm
    bhiddens = Array(Any, length(sequence)-2)
    sb = copy(states)
    for i=length(sequence):-1:3
        x = embeddings[i]
        h = forward(model[:back], sb, x)
        bhiddens[i-2] = copy(h)
    end

    # concatenate layer
    for i=1:length(fhiddens)
        ypred = hcat(fhiddens[i], bhiddens[i]) * model[:soft][1] .+ model[:soft][2]
        ygold = map(x->x[1], sequence[i+1])
        total += logprob(ygold, ypred)
        count += length(ygold)
    end
    val = - total / count
    push!(lval, AutoGrad.getval(val))
    return val
end

gradcharbilstm = grad(charbilstm)

function train(model, chstates, states, sequence, i2w_all, chvocab, opts)
    lval = []
    gloss = gradcharbilstm(model, chstates, states, sequence, i2w_all, chvocab, lval)
    update!(model, gloss, opts)
    return lval[1]
end

function devperp(model, chstates, states, i2w_all, chvocab)
    devloss = []
    for d in dev
        charbilstm(model, chstates, states, d, i2w_all, chvocab, devloss)
    end
    return exp(mean(devloss))
end

max_stars_repo_name: kirnap/KU-dependencyParser

max_issues_repo_name: kirnap/KU-dependencyParser

max_forks_repo_name: kirnap/KU-dependencyParser


text:
module SolverTraces
using Crayons
using Formatting
using ProgressMeter
using Compat

"""
    TraceColumn

Base type for all columns that can appear in a
[`SolverTrace`](@ref). Each subtype is expected to have a property
named `fmt` containing a formatting string, and to be callable with
the current step number as argument, returning the necessary values to
render the formatting string, alternatively overload
`Formatting.format(c::C, i::Integer) where {C<:TraceColumn}`.
"""
abstract type TraceColumn end
export TraceColumn

Formatting.format(c::TraceColumn, i::Integer) =
    format(c.fmt, c(i)...)

include("colors.jl")
include("current_step.jl")
include("units.jl")
include("performance.jl")
include("eta.jl")
include("scalar_column.jl")
include("tolerance.jl")
include("solver_trace.jl")

end # module


text:
## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

module Gaston

export closefigure, closeall, clearfigure, figure, plot, histogram, imagesc,
surf, printfigure, gaston_demo, set_terminal, set_filename, set_default_legend,
set_default_plotstyle, set_default_color, set_default_marker,
set_default_linewidth, set_default_pointsize, set_default_title,
set_default_xlabel, set_default_ylabel, set_default_zlabel, set_default_fill,
set_default_grid, set_default_box, set_default_axis, set_default_xrange,
set_default_yrange, set_default_zrange, set_print_color, set_print_fontface,
set_print_fontsize, set_print_fontscale, set_print_linewidth, set_print_size,
gnuplot_send, addconf, addcoords, adderror, addfinancial, CurveConf, AxesConf,
llplot, meshgrid, gnuplot_exit

import Base.writemime

# before doing anything else, verify gnuplot is present on this system
if !success(`gnuplot --version`)
    error("Gaston cannot be loaded: gnuplot is not available on this system.")
end

# load files
include("gaston_types.jl")
include("gaston_aux.jl")
include("gaston_lowlvl.jl")
include("gaston_midlvl.jl")
include("gaston_hilvl.jl")
include("gaston_config.jl")
include("gaston_demo.jl")

# set up global variables
# global variable that stores gnuplot's state
gnuplot_state = GnuplotState(false,0,0,mktempdir(),"","","",false,Any[])

# when gnuplot_state goes out of scope, close the pipe
finalizer(gnuplot_state,gnuplot_exit)

# global variable that stores Gaston's configuration
gaston_config = GastonConfig()

end


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__()

module SynthesisFilters

using Compat
import Compat: view
using MelGeneralizedCepstrums

import SPTK
import Base: filt!
import MelGeneralizedCepstrums: allpass_alpha, glog_gamma

export
    # Synthesis filter types
    SynthesisFilter,
    MelGeneralizedCepstrumSynthesisFilter, # LMADF or MLSADF or MGLSADF
    LinearPredictionVariantsSynthesisFilter, # AllPoleDF

    # Speech waveform synthesis filters
    AllPoleDF,            # All-pole digital filter for synthesis from LPC
    AllPoleLatticeDF,     # All-pole lattice digital filter for synthesis from PARCOR
    LSPDF,                # LSP digital filter for synthesis from LSP
    LMADF,                # Log magnitude approximation digital filter for synthesis from cepstrum
    MLSADF,               # Mel-log spectrum approximation digital filter for synthesis from mel-cepstrum
    MGLSADF,              # Mel generalized-log spectrum approximation digital filter for synthesis from mel-generalized cepstrum

    # High-level interface for waveform synthesis
    synthesis,
    synthesis!,
    synthesis_one_frame!,

    filt!,                # filtering one sample
    to_filtcoef,          # spectral parameter to filter coefficients

    # Mel-generalized cepstrum synthesis filter properties
    allpass_alpha,        # all-pass constant (α)
    glog_gamma,           # parameter of generalized log function (γ)

    # MGLSADF property
    nstage                # -1/γ; number of stages for MGLSADF

for fname in [
              "const",
              "common",
              "poledf",
              "ltcdf",
              "lspdf",
              "lmadf",
              "mlsadf",
              "mglsadf",
              "helper"
    ]
    include(string(fname, ".jl"))
end

end # module


text:
module Util

import ..Clarus

function activity(;params...)
  return Clarus.api_request("Util","Activity";params...)
end

function tickers(;params...)
  return Clarus.api_request("Util","Tickers";params...)
end

function periodlength(;params...)
  return Clarus.api_request("Util","PeriodLength";params...)
end

function grid(;params...)
  return Clarus.api_request("Util","Grid";params...)
end

function domain(;params...)
  return Clarus.api_request("Util","Domain";params...)
end

end #MODULE-END


text:
module TestExamples2DEuler

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_2d_dgsem")

@testset "Compressible Euler" begin
  @trixi_testset "elixir_euler_source_terms.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
      l2   = [8.517808508019351e-7, 1.2350203856098537e-6, 1.2350203856728076e-6, 4.277886946638239e-6],
      linf = [8.357848139128876e-6, 1.0326302096741458e-5, 1.0326302101404394e-5, 4.496194024383726e-5])
  end

  @trixi_testset "elixir_euler_convergence_pure_fv.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_convergence_pure_fv.jl"),
      l2   = [0.026440292358506527, 0.013245905852168414, 0.013245905852168479, 0.03912520302609374],
      linf = [0.042130817806361964, 0.022685499230187034, 0.022685499230187922, 0.06999771202145322])
  end

  @trixi_testset "elixir_euler_density_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_density_wave.jl"),
      l2   = [0.0010600778457964775, 0.00010600778457634275, 0.00021201556915872665, 2.650194614399671e-5],
      linf = [0.006614198043413566, 0.0006614198043973507, 0.001322839608837334, 0.000165354951256802],
      tspan = (0.0, 0.5))
  end

  @trixi_testset "elixir_euler_source_terms_nonperiodic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonperiodic.jl"),
      l2   = [2.3653424742011065e-6, 2.1388875092652046e-6, 2.1388875093792834e-6, 6.010896863407163e-6],
      linf = [1.4080465934984687e-5, 1.7579850582816192e-5, 1.757985059525069e-5, 5.95689353266593e-5])
  end

  @trixi_testset "elixir_euler_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
      l2   = [0.061751715597716854, 0.05018223615408711, 0.05018989446443463, 0.225871559730513],
      linf = [0.29347582879608825, 0.31081249232844693, 0.3107380389947736, 1.0540358049885143])
  end

  @trixi_testset "elixir_euler_ec.jl with flux_kennedy_gruber" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
      l2   = [0.03481471610306124, 0.027694280613944234, 0.027697905866996532, 0.12932052501462554],
      linf = [0.31052098400669004, 0.3481295959664616, 0.34807152194137336, 1.1044947556170719],
      maxiters = 10,
      surface_flux = flux_kennedy_gruber,
      volume_flux = flux_kennedy_gruber)
  end

  @trixi_testset "elixir_euler_ec.jl with flux_chandrashekar" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
      l2   = [0.03481122603050542, 0.027662840593087695, 0.027665658732350273, 0.12927455860656786],
      linf = [0.3110089578739834, 0.34888111987218107, 0.3488278669826813, 1.1056349046774305],
      maxiters = 10,
      surface_flux = flux_chandrashekar,
      volume_flux = flux_chandrashekar)
  end

  @trixi_testset "elixir_euler_shockcapturing.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shockcapturing.jl"),
      l2   = [0.05380378209319121, 0.04696482040740498, 0.046967440141927635, 0.19686385916646665],
      linf = [0.1852693811520424, 0.24028641582658528, 0.2326441454434102, 0.6874069967414047])
  end

  @trixi_testset "elixir_euler_blast_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave.jl"),
      l2   = [0.1399704878290245, 0.11580885255188444, 0.115808885313645, 0.3391969791404581],
      linf = [1.4473441276982646, 1.3355435317505784, 1.335543531751034, 1.8203831482774346],
      maxiters = 30)
  end

  @trixi_testset "elixir_euler_blast_wave_pure_fv.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_pure_fv.jl"),
      l2   = [0.39957047631960346, 0.21006912294983154, 0.21006903549932, 0.6280328163981136],
      linf = [2.20417889887697, 1.5487238480003327, 1.5486788679247812, 2.4656795949035857],
      tspan = (0.0, 0.5))
  end

  @trixi_testset "elixir_euler_blast_wave_amr.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_amr.jl"),
      l2   = [0.6777351803934385, 0.28143110338521693, 0.28143035507881886, 0.7211949078634373],
      linf = [2.881942488461436, 1.8037294505391677, 1.803677581960154, 3.037915614619938],
      tspan = (0.0, 1.0))
  end

  @trixi_testset "elixir_euler_sedov_blast_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave.jl"),
      l2   = [0.4819231474652106, 0.1655453843075888, 0.16554538430758944, 0.6182451347074112],
      linf = [2.4849541392059202, 1.2816207838889486, 1.2816207838888944, 6.4743208586529875],
      tspan = (0.0, 1.0))
  end

  @trixi_testset "elixir_euler_positivity.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_positivity.jl"),
      l2   = [0.48401054859424775, 0.16613328471131006, 0.16613328471131006, 0.6181793382171535],
      linf = [2.4903327485761544, 1.2898210694161085, 1.2898210694161072, 6.4723873993158385],
      tspan = (0.0, 1.0))
  end

  @trixi_testset "elixir_euler_blob_mortar.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blob_mortar.jl"),
      l2   = [0.22206706482833133, 0.6296849068200803, 0.24351238064525166, 2.9118777552679345],
      linf = [9.584083314003589, 27.72635365994101, 9.343044146443756, 101.53747277352629],
      tspan = (0.0, 0.5))
  end

  @trixi_testset "elixir_euler_blob_amr.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blob_amr.jl"),
      l2   = [0.20150296718390992, 1.1816675092370863, 0.10129652123574924, 5.230357225596295],
      linf = [14.174958287930844, 71.16876942972455, 7.2799739019809415, 291.67608863630244],
      tspan = (0.0, 0.12))
  end

  @trixi_testset "elixir_euler_kelvin_helmholtz_instability.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_kelvin_helmholtz_instability.jl"),
      l2   = [0.055689867638859906, 0.032985369421826734, 0.052243738479375024, 0.08009259073134713],
      linf = [0.24051443859674326, 0.16611572070980668, 0.123559477474734, 0.2694484247345663],
      tspan = (0.0, 0.2))
  end

  @trixi_testset "elixir_euler_kelvin_helmholtz_instability_amr.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_kelvin_helmholtz_instability_amr.jl"),
      l2   = [0.05569140599245201, 0.033108757924402134, 0.05223526059785845, 0.08007239985228143],
      linf = [0.2540392117514094, 0.17432468822204936, 0.12323401271048683, 0.26897998287356195],
      tspan = (0.0, 0.2))
  end

  @trixi_testset "elixir_euler_colliding_flow.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_colliding_flow.jl"),
      l2   = [7.23705116e-03,   4.48861301e-02,   8.47273157e-07,   6.62667920e-01],
      linf = [1.93754293e-01,   5.52253475e-01,   4.98055004e-05,   1.50879674e+01],
      tspan = (0.0, 0.1))
  end

  @trixi_testset "elixir_euler_colliding_flow_amr.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_colliding_flow_amr.jl"),
      l2   = [6.76960244e-03,   3.21735698e-02,   2.63316434e-07,   6.78416591e-01],
      linf = [2.50119529e-01,   4.07493507e-01,   9.96993324e-05,   2.23204482e+01],
      tspan = (0.0, 0.1))
  end

  @trixi_testset "elixir_euler_astro_jet_amr.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_astro_jet_amr.jl"),
      l2   = [0.011338493575859507, 10.097552800226566, 0.0038670992882799088, 4031.828021246965],
      linf = [3.316670548834158, 2992.6286682433883, 7.919515633718346, 1.1914747265545786e6],
      tspan = (0.0, 1.0e-7))
  end

  @trixi_testset "elixir_euler_vortex.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex.jl"),
      l2   = [3.6343141788303172e-6, 0.0032111379945554378, 0.0032111482803927763, 0.004545715899533244],
      linf = [7.901851921288117e-5, 0.030561510550305093, 0.03050266851613559, 0.042876690674344076])
  end

  @trixi_testset "elixir_euler_vortex_mortar.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex_mortar.jl"),
      l2   = [2.1203067597905396e-6, 2.792922844344066e-5, 3.7593422075434476e-5, 8.813644239472934e-5],
      linf = [5.932046017209647e-5, 0.0007491265252463908, 0.0008165690091972433, 0.0022122638482677814])
  end

  @trixi_testset "elixir_euler_vortex_mortar_split.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex_mortar_split.jl"),
      l2   = [2.1201025902895562e-6, 2.8056524989293713e-5, 3.759500393896144e-5, 8.84137216637438e-5],
      linf = [5.934103282212444e-5, 0.0007552316670187409, 0.0008152449093751235, 0.0022069874183294758])
  end

  @trixi_testset "elixir_euler_vortex_shockcapturing.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex_shockcapturing.jl"),
      l2   = [3.7998676433961623e-6, 5.563687203373343e-5, 5.561628017766063e-5, 0.00015707075214260475],
      linf = [8.522682141076654e-5, 0.0009610342345369727, 0.0009656468947373265, 0.0030861481615822584])
  end

  @trixi_testset "elixir_euler_vortex_mortar_shockcapturing.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex_mortar_shockcapturing.jl"),
      l2   = [2.1203212494737245e-6, 2.8056721454951553e-5, 3.759510818832351e-5, 8.842116438932004e-5],
      linf = [5.9340649220973596e-5, 0.0007552376287922602, 0.0008152432005811283, 0.0022069702233196153])
  end

  @trixi_testset "elixir_euler_vortex_amr.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_vortex_amr.jl"),
      l2   = [2.1325938318497577e-6, 0.003281545066125044, 0.0032806321457482615, 0.00464605975938406],
      linf = [4.667437348149228e-5, 0.03175420871507906, 0.0318039789241531, 0.04561735256198318])
  end
end

end # module


text:
module SymbolicWedderburn

using LinearAlgebra
using SparseArrays
using Primes

using Cyclotomics
using GroupsCore
using PermutationGroups
using StarAlgebras

export symmetry_adapted_basis, WedderburnDecomposition
export basis,
    # degree, # too common name to export
    direct_summands,
    invariant_vectors,
    issimple,
    multiplicity

include("Characters/Characters.jl")
using .Characters
import .Characters: row_echelon_form!
import .Characters.FiniteFields

include("util.jl")
include("actions.jl")
include("action_characters.jl")
include("matrix_projections.jl")
include("minimal_projections.jl")
include("direct_summands.jl")
include("sa_basis.jl")
include("wedderburn_decomposition.jl")

end # module


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

module Reqs

import Base: ==
import ...Pkg.PkgError
using ..Types

# representing lines of REQUIRE files

abstract Line
immutable Comment <: Line
    content::AbstractString
end
immutable Requirement <: Line
    content::AbstractString
    package::AbstractString
    versions::VersionSet
    system::Vector{AbstractString}

    function Requirement(content::AbstractString)
        fields = split(replace(content, r"#.*$", ""))
        system = AbstractString[]
        while !isempty(fields) && fields[1][1] == '@'
            push!(system,shift!(fields)[2:end])
        end
        isempty(fields) && throw(PkgError("invalid requires entry: $content"))
        package = shift!(fields)
        all(field->ismatch(Base.VERSION_REGEX, field), fields) ||
            throw(PkgError("invalid requires entry for $package: $content"))
        versions = VersionNumber[fields...]
        issorted(versions) || throw(PkgError("invalid requires entry for $package: $content"))
        new(content, package, VersionSet(versions), system)
    end
    function Requirement(package::AbstractString, versions::VersionSet, system::Vector{AbstractString}=AbstractString[])
        content = ""
        for os in system
            content *= "@$os "
        end
        content *= package
        if versions != VersionSet()
            for ival in versions.intervals
                (content *= " $(ival.lower)")
                ival.upper < typemax(VersionNumber) &&
                (content *= " $(ival.upper)")
            end
        end
        new(content, package, versions, system)
    end
end

==(a::Line, b::Line) = a.content == b.content
hash(s::Line, h::UInt) = hash(s.content, h + (0x3f5a631add21cb1a % UInt))

# general machinery for parsing REQUIRE files

function read{T<:AbstractString}(readable::Vector{T})
    lines = Line[]
    for line in readable
        line = chomp(line)
        push!(lines, ismatch(r"^\s*(?:#|$)", line) ? Comment(line) : Requirement(line))
    end
    return lines
end

function read(readable::Union{IO,Base.AbstractCmd})
    lines = Line[]
    for line in eachline(readable)
        push!(lines, ismatch(r"^\s*(?:#|$)", line) ? Comment(line) : Requirement(line))
    end
    return lines
end
read(file::AbstractString) = isfile(file) ? open(read,file) : Line[]

function write(io::IO, lines::Vector{Line})
    for line in lines
        println(io, line.content)
    end
end
function write(io::IO, reqs::Requires)
    for pkg in sort!(collect(keys(reqs)), by=lowercase)
        println(io, Requirement(pkg, reqs[pkg]).content)
    end
end
write(file::AbstractString, r::Union{Vector{Line},Requires}) = open(io->write(io,r), file, "w")

function parse(lines::Vector{Line})
    reqs = Requires()
    for line in lines
        if isa(line,Requirement)
            if !isempty(line.system)
                applies = false
                if is_windows(); applies |=  ("windows"  in line.system); end
                if is_unix();    applies |=  ("unix"     in line.system); end
                if is_apple();   applies |=  ("osx"      in line.system); end
                if is_linux();   applies |=  ("linux"    in line.system); end
                if is_bsd();     applies |=  ("bsd"      in line.system); end
                if is_windows(); applies &= !("!windows" in line.system); end
                if is_unix();    applies &= !("!unix"    in line.system); end
                if is_apple();   applies &= !("!osx"     in line.system); end
                if is_linux();   applies &= !("!linux"   in line.system); end
                if is_bsd();     applies &= !("!bsd"     in line.system); end
                applies || continue
            end
            reqs[line.package] = haskey(reqs, line.package) ?
                intersect(reqs[line.package], line.versions) : line.versions
        end
    end
    return reqs
end
parse(x) = parse(read(x))

function dependents(packagename::AbstractString)
    pkgs = AbstractString[]
    cd(Pkg.dir()) do
        for (pkg,latest) in Pkg.Read.latest()
            if haskey(latest.requires, packagename)
                push!(pkgs, pkg)
            end
        end
    end
    pkgs
end

# add & rm – edit the content a requires file

function add(lines::Vector{Line}, pkg::AbstractString, versions::VersionSet=VersionSet())
    v = VersionSet[]
    filtered = filter(lines) do line
        if !isa(line,Comment) && line.package == pkg && isempty(line.system)
            push!(v, line.versions)
            return false
        end
        return true
    end
    length(v) == 1 && v[1] == intersect(v[1],versions) && return copy(lines)
    versions = reduce(intersect, versions, v)
    push!(filtered, Requirement(pkg, versions))
end

rm(lines::Vector{Line}, pkg::AbstractString) = filter(lines) do line
    isa(line,Comment) || line.package != pkg
end

end # module


text:
using RecipesBase
using ColorTypes
using LaTeXStrings
#import PlotUtils: cgrad

const mygreen = RGBA{Float64}(151/255,180/255,118/255,1)
const mygreen2 = RGBA{Float64}(113/255,161/255,103/255,1)
const myblue = RGBA{Float64}(74/255,144/255,226/255,1)

## Main scalar grid data plot recipe

@recipe function f(x::AbstractRange,y::AbstractRange,w::T;trim=0) where {T<:ScalarGridData}
    seriestype --> :contour
    aspect_ratio := 1
    legend --> :none
    grid --> :none
    linewidth --> 1
    framestyle --> :frame
    xlims --> (-Inf,Inf)
    ylims --> (-Inf,Inf)
    levels --> range(minimum(w.data),stop=maximum(w.data),length=16)
    x[1+trim:end-trim],y[1+trim:end-trim],transpose(w.data[1+trim:end-trim,1+trim:end-trim])
end

@recipe function f(w::T,g::PhysicalGrid) where {T<:ScalarGridData}
    xg, yg = coordinates(w,g)
    xg, yg, w
end

@recipe function f(w::T) where {T<:ScalarGridData}
    m,n = size(w.data)
    1:m, 1:n, w
end

## Vector grid data

@recipe function f(xu::AbstractRange,yu::AbstractRange,xv::AbstractRange,yv::AbstractRange,q::T) where {T<:VectorGridData}
    layout := (1,2)
    @series begin
      subplot := 1
      title --> L"u"
      xu, yu, q.u
    end

    @series begin
      subplot := 2
      title --> L"v"
      xv, yv, q.v
    end

end

@recipe function f(q::T,g::PhysicalGrid) where {T<:VectorGridData}
    xu, yu, xv, yv = coordinates(q,g)
    return xu, yu, xv, yv, q
end

@recipe function f(q::T) where {T<:VectorGridData}
    mu,nu = size(q.u)
    mv,nv = size(q.v)
    return 1:mu, 1:nu, 1:mv, 1:nv, q
end

## Tensor grid data

@recipe function f(xdudx::AbstractRange,ydudx::AbstractRange,
                   xdudy::AbstractRange,ydudy::AbstractRange,
                   xdvdx::AbstractRange,ydvdx::AbstractRange,
                   xdvdy::AbstractRange,ydvdy::AbstractRange,q::T) where {T<:TensorGridData}
    layout := (2,2)
    @series begin
      subplot := 1
      title --> L"du/dx"
      xdudx, ydudx, q.dudx
    end

    @series begin
      subplot := 2
      title --> L"dv/dx"
      xdvdx, ydvdx, q.dvdx
    end

    @series begin
      subplot := 3
      title --> L"du/dy"
      xdudy, ydudy, q.dudy
    end

    @series begin
      subplot := 4
      title --> L"dv/dy"
      xdvdy, ydvdy, q.dvdy
    end

end

@recipe function f(q::T,g::PhysicalGrid) where {T<:TensorGridData}
    xdudx, ydudx = coordinates(q.dudx,g)
    xdudy, ydudy = coordinates(q.dudy,g)
    xdvdx, ydvdx = coordinates(q.dvdx,g)
    xdvdy, ydvdy = coordinates(q.dvdy,g)
    return xdudx, ydudx, xdudy, ydudy, xdvdx, ydvdx, xdvdy, ydvdy, q
end

@recipe function f(q::T) where {T<:TensorGridData}
    mdudx,ndudx = size(q.dudx)
    mdudy,ndudy = size(q.dudy)
    mdvdx,ndvdx = size(q.dvdx)
    mdvdy,ndvdy = size(q.dvdy)
    return 1:mdudx, 1:ndudx, 1:mdudy, 1:ndudy, 1:mdvdx, 1:ndvdx, 1:mdvdy, 1:ndvdy, q
end

## Generated fields

@recipe function f(field::GeneratedField)
    field(), field.grid
end

max_stars_repo_name: JuliaIBPM/CartesianGrids

max_issues_repo_name: JuliaIBPM/CartesianGrids

max_forks_repo_name: JuliaIBPM/CartesianGrids


text:
# Copyright (c) 2014: SCS.jl contributors
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

struct ScaledPSDCone <: MOI.AbstractVectorSet
    side_dimension::Int
end

function MOI.Utilities.set_with_dimension(::Type{ScaledPSDCone}, dim)
    return ScaledPSDCone(div(-1 + isqrt(1 + 8 * dim), 2))
end

Base.copy(x::ScaledPSDCone) = ScaledPSDCone(x.side_dimension)

MOI.side_dimension(x::ScaledPSDCone) = x.side_dimension

function MOI.dimension(x::ScaledPSDCone)
    return div(x.side_dimension * (x.side_dimension + 1), 2)
end

struct ScaledPSDConeBridge{T,G} <: MOI.Bridges.Constraint.SetMapBridge{
    T,
    ScaledPSDCone,
    MOI.PositiveSemidefiniteConeTriangle,
    MOI.VectorAffineFunction{T},
    G,
}
    constraint::MOI.ConstraintIndex{MOI.VectorAffineFunction{T},ScaledPSDCone}
end

function MOI.Bridges.Constraint.concrete_bridge_type(
    ::Type{ScaledPSDConeBridge{T}},
    ::Type{G},
    ::Type{MOI.PositiveSemidefiniteConeTriangle},
) where {T,G<:Union{MOI.VectorOfVariables,MOI.VectorAffineFunction{T}}}
    return ScaledPSDConeBridge{T,G}
end

function MOI.Bridges.map_set(
    ::Type{<:ScaledPSDConeBridge},
    set::MOI.PositiveSemidefiniteConeTriangle,
)
    return ScaledPSDCone(set.side_dimension)
end

function MOI.Bridges.inverse_map_set(
    ::Type{<:ScaledPSDConeBridge},
    set::ScaledPSDCone,
)
    return MOI.PositiveSemidefiniteConeTriangle(set.side_dimension)
end

function _upper_to_lower_triangular_permutation(dim::Int)
    side_dimension = MOI.Utilities.side_dimension_for_vectorized_dimension(dim)
    permutation = zeros(Int, dim)
    i = 0
    for row in 1:side_dimension
        start = div(row * (row + 1), 2)
        for col in row:side_dimension
            i += 1
            permutation[i] = start
            start += col
        end
    end
    return sortperm(permutation), permutation
end

function _transform_function(
    func::MOI.VectorAffineFunction{T},
    scale,
    moi_to_scs::Bool,
) where {T}
    d = MOI.output_dimension(func)
    # upper_to_lower[i] maps the i'th element of the upper matrix to the linear
    #   index of the lower
    # lower_to_upper[i] maps the i'th element of the lower matrix to the linear
    #   index of the upper
    upper_to_lower, lower_to_upper = _upper_to_lower_triangular_permutation(d)
    # scale_factor[i] is the amount to scale the i'th linear index in the MOI
    # function by.
    scale_factor = fill(scale, d)
    for i in 1:d
        if MOI.Utilities.is_diagonal_vectorized_index(i)
            scale_factor[i] = 1.0
        end
    end
    scaled_constants = if moi_to_scs
        (func.constants.*scale_factor)[lower_to_upper]
    else
        func.constants[upper_to_lower] .* scale_factor
    end
    scaled_terms = MOI.VectorAffineTerm{T}[]
    for term in func.terms
        row = term.output_index
        i = moi_to_scs ? upper_to_lower[row] : lower_to_upper[row]
        scale_i = moi_to_scs ? scale_factor[row] : scale_factor[i]
        push!(
            scaled_terms,
            MOI.VectorAffineTerm(
                i,
                MOI.ScalarAffineTerm(
                    term.scalar_term.coefficient * scale_i,
                    term.scalar_term.variable,
                ),
            ),
        )
    end
    return MOI.VectorAffineFunction(scaled_terms, scaled_constants)
end

function _transform_function(func::MOI.VectorOfVariables, scale, moi_to_scs)
    new_f = MOI.Utilities.operate(*, Float64, 1.0, func)
    return _transform_function(new_f, scale, moi_to_scs)
end

function _transform_function(func::Vector{T}, scale, moi_to_scs::Bool) where {T}
    d = length(func)
    upper_to_lower, lower_to_upper = _upper_to_lower_triangular_permutation(d)
    scale_factor = fill(scale, d)
    for i in 1:d
        if MOI.Utilities.is_diagonal_vectorized_index(i)
            scale_factor[i] = 1.0
        end
    end
    if moi_to_scs
        return (func.*scale_factor)[lower_to_upper]
    else
        return func[upper_to_lower] .* scale_factor
    end
end

# Map ConstraintFunction from MOI -> SCS
function MOI.Bridges.map_function(::Type{<:ScaledPSDConeBridge}, f)
    return _transform_function(f, √2, true)
end

# Used to map the ConstraintPrimal from SCS -> MOI
function MOI.Bridges.inverse_map_function(::Type{<:ScaledPSDConeBridge}, f)
    return _transform_function(f, 1 / √2, false)
end

# Used to map the ConstraintDual from SCS -> MOI
function MOI.Bridges.adjoint_map_function(::Type{<:ScaledPSDConeBridge}, f)
    return _transform_function(f, 1 / √2, false)
end

# Used to set ConstraintDualStart
function MOI.Bridges.inverse_adjoint_map_function(
    ::Type{<:ScaledPSDConeBridge},
    f,
)
    return _transform_function(f, √2, true)
end