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

using Markdown
using InteractiveUtils

# ╔═╡ ca80f5a1-e389-4522-91a5-b2bf35da8f56
begin
	using PlutoUI
	using Optim
end

# ╔═╡ 16cd6686-6b4b-40c4-92a0-dd38b4df4b2b
md"""
##### Example: Equilibrium Extent of Reaction for Ideal Aqueous Reaction Mixtures

Consider the reversible aqueous phase reaction A$_{aq}$ $\rightleftharpoons$ B$_{aq}$ (reaction 1) in which we assume _ideal_ liquid phase behavior. Compute the equlibrium extent $\epsilon$ of reaction as a function of different $A$ and $B$ compounds.

__Background__:
We derived the expression relating the equilibrium constant to the Gibbs reaction energy change:

$$\prod_{i=1}^{\mathcal{M}}\left(x_{i}\gamma_{i}\right)^{\sigma_{i1}} = \exp\left(-\frac{\Delta{G}^{\circ}}{RT}\right)$$

where the product (the $\prod\star$ term) is the liquid phase equilibrium constant $K$, $x_{i}$ denotes the mol fraction of component $i$, and $\sigma_{i1}$ denotes the stoichiometric coefficient governing component $i$ in reaction $1$ (the only reaction that we have). We have assumed the reaction is _ideal_, meaning the activity coefficient $\gamma_{i}\rightarrow{1},\forall{i}$. This allows us to expand the product (where A is species 1, B is species 2, and $\sigma_{11}$ = -1 and $\sigma_{21}$ = 1):

$$\frac{x_{2}}{x_{1}} = \exp\left(-\frac{\Delta{G}^{\circ}}{RT}\right)$$

or (after substitution of the definition of mole fraction):

$$\frac{n^{\circ}_{2}+\epsilon}{n^{\circ}_{1}-\epsilon} = \exp\left(-\frac{\Delta{G}^{\circ}}{RT}\right)$$

In general, we need to estimate $\epsilon$ numerically, but un this case, we can determine an analytical solution for the equilibrium extent $\epsilon$:

$$\epsilon = \frac{bn_{1}^{\circ}-n_{2}^{\circ}}{b+1}$$

where $b\equiv\exp\left(-\Delta{G}^{\circ}/RT\right)$.

"""

# ╔═╡ 1987041a-3cb4-4d8f-b262-e1eeed63cfe5
md"""
##### Example problem setup
"""

# ╔═╡ 1ca273c0-014b-4da8-abb9-df3747e07332
begin

	# array for Gibbs formation data -
	gibbs_formation_data_array = [

		-1304.7 -1302.1 	; # G6P and F6P 	units: kJ/mol
		-1097.2 -1091.5 	; # DHAP and GAP 	units: kJ/mol
		-1356.7 -1352.1 	; # G3P and G2P 	units: kJ/mol

	]; # dim: N_examples rows x 2 cols (A and B)

	# convert to kJ/mmol -
	gibbs_formation_data_array = (1/1000)*gibbs_formation_data_array

	# setup the state array -
	initial_mol_array = [

		10.0 	0.0 	; # G6P and F6P units: mmol
		10.0 	0.0 	; # DHAP and GAP units: mmol
		10.0 	0.0 	; # G3P and G2P units: mmol
	];

	# show -
	nothing
end

# ╔═╡ 68f9d52f-96ca-42e0-ba7d-e3268da9e6ea
# R constant in kJ/mmol-K
R = 8.314*(1/1000)*(1/1000)

# ╔═╡ 26bc5887-04af-457c-be76-e1d9a836129c
# Temperature in K -
T = 25+273.15

# ╔═╡ 11ed237f-9739-4294-ba37-546238ad3c5e
md"""
##### Approach 1: Use the analytical solution to estimate the equilibrium ϵ
"""

# ╔═╡ af3d64a0-7553-4bad-a122-179febd11773
begin

	# how many examples do we have?
	number_of_examples = size(gibbs_formation_data_array,1)

	# initialize some space -
	ϵ_array = Array{Float64,1}(undef,number_of_examples)

	# alias -
	RT = R*T

	# compute the eq extent for each example -
	for example_index = 1:number_of_examples

		# compute the ΔG term -
		ΔG = gibbs_formation_data_array[example_index,2] - gibbs_formation_data_array[example_index,1]

		# get initial mols -
		A_initial = initial_mol_array[example_index,1]
		B_initial = initial_mol_array[example_index,2]

		# compute b -
		b = exp(-ΔG/RT)

		# compute ϵ -
		ϵ = (b*A_initial - B_initial)/(b+1)

		# grab ϵ and store -
		ϵ_array[example_index] = ϵ*(1/A_initial) 	# return a scaled ϵ -
	end
end

# ╔═╡ 36a2c2a1-2598-4fc0-adf4-23fb108de1f0
ϵ_array

# ╔═╡ 8755b296-9d1a-45c7-a66a-d3315a921b75
md"""
##### Approach 2: Numerically estimate the equilibrium ϵ
In the general case, we may not be able to solve for ϵ at equilibrium. In these cases, we need to estimate the solution numerically using one of many possible software tools e.g., MATLAB, Excel, Python or our favorite programming language [Julia](https://julialang.org).

Imagine that we have the general problem (single reaction, not necessarily ideal):

$$f\left(ϵ\right) = \exp\left(-\frac{\Delta{G}^{\circ}}{RT}\right)$$

where $f\left(ϵ\right)$ is the equilibrium constant expression written in terms of the extent of reaction $\epsilon$. We need to find an $\epsilon$, call it $\epsilon^{\star}$, such that:

$$f\left(ϵ^{\star}\right) - \exp\left(-\frac{\Delta{G}^{\circ}}{RT}\right) = 0$$

To do that we solve the following optimization problem:

$$\min_\epsilon \left(f\left(ϵ^{\star}\right) - \exp\left(-\frac{\Delta{G}^{\circ}}{RT}\right)\right)^{2}$$

subject to the constraint $\epsilon\in(0,\min\left(n^{\circ}_{1},n^{\circ}_{2},\dots\right))$ or for an extent scaled by the initial number of mols of the limiting reactant $\epsilon\in(0,1)$.

"""

# ╔═╡ 9dd10d12-32f1-4f9c-acb1-e5b951d0cbb0
begin

	# put problem stuff into a dictionary -
	parameters_dictionary = Dict{String,Any}()
	parameters_dictionary["gibbs_formation_data_array"] = gibbs_formation_data_array
	parameters_dictionary["initial_mol_array"] = initial_mol_array

	# show -
	nothing
end

# ╔═╡ 84460784-4c05-4217-baed-4c420c39822d
function objective_function(epsilon, parameters, example_index)

	# grab stuff from the parameters dict -
	gibbs_formation_data_array = parameters["gibbs_formation_data_array"]
	initial_mol_array = parameters["initial_mol_array"]

	# compute the ΔG term -
	ΔG = gibbs_formation_data_array[example_index,2] - gibbs_formation_data_array[example_index,1]

	# get initial mols -
	A_initial = initial_mol_array[example_index,1]
	B_initial = initial_mol_array[example_index,2]

	# compute b -
	b = exp(-ΔG/RT)

	# unscale ϵ -
	ϵ = (A_initial)*epsilon[1]

	# compute the residual -
	δ = (B_initial + ϵ)/(A_initial - ϵ) - b

	# return -
	return δ^2
end

# ╔═╡ 0576ca4a-f918-4332-a178-e74c68314f1f
begin

	# which example are we doing?
	example_index = 1

	# Setup an initial guess for ϵ -
	ϵ_initial_guess = [0.9999]

	# setup bounds -
	L = 0.00001
	U = 0.99999

	# setup the objective function -
	OF(p) = objective_function(p, parameters_dictionary, example_index)

    # call the optimizer -
    opt_result = optimize(OF,L, U, ϵ_initial_guess, Fminbox(BFGS()))
end

# ╔═╡ 667e2d50-3f29-41cf-9bea-3488a4c90d16
with_terminal() do

	bgfs_soln = Optim.minimizer(opt_result)[1]
	println("Optim found ϵ = $(bgfs_soln) for example index = $(example_index)")
end

# ╔═╡ 64d7437c-35b8-11ec-2ba3-13ef42e58e20
html"""
<style>
main {
    max-width: 1200px;
    width: 90%;
    margin: auto;
    font-family: "Roboto, monospace";
}

a {
    color: blue;
    text-decoration: none;
}

.H1 {
    padding: 0px 30px;
}
</style>"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Optim = "429524aa-4258-5aef-a3af-852621145aeb"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
Optim = "~1.4.1"
PlutoUI = "~0.7.16"
"""

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

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

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

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

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

[[ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "0541d306de71e267c1a724f84d44bbc981f287b4"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.10.2"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[LogExpFunctions]]
deps = ["ChainRulesCore", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "6193c3815f13ba1b78a51ce391db8be016ae9214"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.4"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

# ╔═╡ Cell order:
# ╟─16cd6686-6b4b-40c4-92a0-dd38b4df4b2b
# ╟─1987041a-3cb4-4d8f-b262-e1eeed63cfe5
# ╠═1ca273c0-014b-4da8-abb9-df3747e07332
# ╠═68f9d52f-96ca-42e0-ba7d-e3268da9e6ea
# ╠═26bc5887-04af-457c-be76-e1d9a836129c
# ╟─11ed237f-9739-4294-ba37-546238ad3c5e
# ╠═af3d64a0-7553-4bad-a122-179febd11773
# ╠═36a2c2a1-2598-4fc0-adf4-23fb108de1f0
# ╟─8755b296-9d1a-45c7-a66a-d3315a921b75
# ╠═9dd10d12-32f1-4f9c-acb1-e5b951d0cbb0
# ╠═0576ca4a-f918-4332-a178-e74c68314f1f
# ╠═667e2d50-3f29-41cf-9bea-3488a4c90d16
# ╠═84460784-4c05-4217-baed-4c420c39822d
# ╠═ca80f5a1-e389-4522-91a5-b2bf35da8f56
# ╟─64d7437c-35b8-11ec-2ba3-13ef42e58e20
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using Clang.cindex
using Clang.wrap_c

JULIAHOME=EnvHash()["JULIAHOME"]
SODIUM_PATH = "/usr/local/include/sodium/"
SODIUM_TL = "/usr/local/include/sodium.h"
OUT_DIR = "../src/"

clang_includes = map(x->joinpath(JULIAHOME, x), [
  "deps/llvm-3.2/build/Release/lib/clang/3.2/include",
  "deps/llvm-3.2/include",
  "deps/llvm-3.2/include",
  "deps/llvm-3.2/build/include/",
  "deps/llvm-3.2/include/"
  ])
clang_extraargs = ["-D", "__STDC_LIMIT_MACROS", "-D", "__STDC_CONSTANT_MACROS"]

wc = wrap_c.init(
    ".",
    "$(OUT_DIR)ls_common_h.jl",
    clang_includes,
    clang_extraargs,
    (th, h) ->
    begin
        if search(h, "sodium") == 0:-1
            return false
        else
#            println("th : $th, h : $h")
            return true
        end
        true
    end,
    h -> "libsodium",
    h ->
    begin
#        println("filename for header : $h")
        "$(OUT_DIR)ls_" * replace(last(split(h, "/")), ".", "_")  * ".jl"
    end)

wc.options.wrap_structs = true
wrap_c.wrap_c_headers(wc, [SODIUM_TL])

f = open("$(OUT_DIR)ls_defines_h.jl", "w+")
println(f, "#   Generating #define constants")

fe = open("$(OUT_DIR)ls_exports_h.jl", "w+")
println(fe, "#   Generating exports")
for e in split(open(readall, "$(OUT_DIR)ls_sodium_h.jl"), "\n")
  m = match(r"^\s*\@c\s+[\w\:\{\}\_]+\s+(\w+)", e)
  if (m != nothing)
#    println(m)
    @printf fe "export %s\n"  m.captures[1]
  end
end

for e in split(open(readall, "$(OUT_DIR)ls_common_h.jl"), "\n")
  m = match(r"^\s*\@ctypedef\s+(\w+)", e)
  if (m != nothing)
#   println(m.captures[1])
    @printf fe "export %s\n"  m.captures[1]
  else
    m = match(r"^\s*const\s+(\w+)", e)
    if (m != nothing)
#       println(m.captures[1])
        @printf fe "export %s\n"  m.captures[1]
    end
  end
end

g_consts = Set()

for fn in split(readall(`ls $(SODIUM_PATH)`))
    hashdefs = split(readall(`gcc -E -dD -P $(SODIUM_PATH)$fn`), "\n")
    for e in hashdefs
        m = match(r"^\s*#define\s+crypto_(\w+)\s+(.+)", e)
        if (m != nothing)
    #        println(m.captures[2])
            if !(contains(g_consts, m.captures[1]))
                m2 = match(r"^\s*(\d+)([a-zA-Z]+)", m.captures[2])
                if m2 != nothing
                    println ("replacing $(m.captures[2]) with $(m2.captures[1])")
                    c2 = m2.captures[1]
                else
                    c2 = m.captures[2]
                end

                @printf f "const crypto_%-30s = %s\n"  m.captures[1]  c2
                @printf fe "export crypto_%s\n"  m.captures[1]

                add!(g_consts, m.captures[1])
            end
        end
    end
end

close(f)
close(fe)

ext: jl


text:
propositions = [T, F]
p = Proposition(rand(propositions))
q = Proposition(rand(propositions))
r = Proposition(rand(propositions))

## Axiom 1.1 [Commutativity] ##

@test (p * q ↔ q * p) == T
@test (p + q ↔ q + p) == T
@test ((p ↔ q) ↔ (q ↔ p)) == T

## Axiom 1.2 [Associativity] ##

@test (p * (q * r) ↔ (p * q) * r) == T
@test (p + (q + r) ↔ (p + q) + r) == T

## Axiom 1.3 [Distributivity] ##

@test (p + (q * r) ↔ (p + q) * (p + r)) == T
@test (p * (q + r) ↔ (p * q) + (p * r)) == T

## Axiom 1.4 [De Morgan]

@test (!(p * q) ↔ !p + !q) == T
@test (!(p + q) ↔ !p * !q) == T

## Axiom 1.5 [Negation] ##

@test (!!p ↔ p) == T

## Axiom 1.6 [Excluded Middle] ##

@test (p + !p ↔ T) == T

## Axiom 1.7 [Contradiction] ##

@test (p * !p ↔ F) == T

## Axiom 1.8 [Implication] ##

@test (p ⟹ q ↔ !p + q) == T

## Axiom 1.9 [Equality] ##

@test ((p ↔ q) ↔ (p ⟹ q) * (q ⟹ p)) == T

## Axiom 1.10 [or-simplification] ##

@test (p + p ↔ p) == T
@test (p + T ↔ T) == T
@test (p + F ↔ p) == T
@test (p + (p * q) ↔ p) == T

## Axiom 1.11 [and-simplification] ##

@test (p * p ↔ p) == T
@test (p * T ↔ p) == T
@test (p * F ↔ F) == T
@test (p * (p + q) ↔ p) == T

## Axiom 1.12 [Identity] ##

@test (p ↔ p) == T

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
using ExcelReaders
using Statistics
import Dates
using StatsBase
using LinearAlgebra

abstract type DataFrequency end
struct MonthlyData <: DataFrequency
    nobs::Int # number of observations
    ns::Int # number of series
end

struct QuarterlyData <: DataFrequency
    nobs::Int # number of observations
    ns::Int # number of series
end

abstract type DetrendMethod end
struct BiWeight{TR <: Real} <: DetrendMethod
    weight::TR
end
struct Mean <: DetrendMethod end
struct NoDetrend <: DetrendMethod end

get_sample_periods(initvec::Array, lastvec::Array, ::MonthlyData) =
    12*(lastvec[1]-initvec[1]-1)+lastvec[2]+(12-initvec[2]+1)
get_sample_periods(initvec::Array, lastvec::Array, ::QuarterlyData) =
    4*(lastvec[1]-initvec[1]-1)+lastvec[2]+(4-initvec[2]+1)
function MonthlyData(initvec, lastvec, ns)
    nobs = 12*(lastvec[1]-initvec[1]-1)+lastvec[2]+(12-initvec[2]+1)
    return MonthlyData(nobs, ns)
end
function QuarterlyData(initvec, lastvec, ns)
    nobs = 4*(lastvec[1]-initvec[1]-1)+lastvec[2]+(4-initvec[2]+1)
    return QuarterlyData(nobs, ns)
end
"""
deflate series by PCE
"""
function deflate_series!(labvec_long::AbstractVector, labvec_short::AbstractVector,
                        x::AbstractVector, price_def::AbstractVector,
                        price_def_lfe::AbstractVector, price_def_pgdp,
                        is::Integer, defcode::Val{1})
    x .= x./price_def
    labvec_long[is] = labvec_long[is] * " Defl by PCE Def"
    labvec_short[is] = "Real_" * labvec_short[is]
    return x
end
"""
deflate series by PCE(LFE)
"""
function deflate_series!(labvec_long::AbstractVector, labvec_short::AbstractVector,
                         x::AbstractVector, price_def::AbstractVector,
                         price_def_lfe::AbstractVector, price_def_pgdp,
                         is::Integer, defcode::Val{2})
    x .= x./price_def_lfe
    labvec_long[is] = labvec_long[is] * " Defl by PCE(LFE) Def"
    labvec_short[is] = "Real_" * labvec_short[is]
    return x
end
"""
deflate series by GDP
"""
function deflate_series!(labvec_long::AbstractVector, labvec_short::AbstractVector,
                         x::AbstractVector, price_def::AbstractVector,
                         price_def_lfe::AbstractVector, price_def_pgdp,
                         is::Integer, defcode::Val{3})
    x .= x./price_def_pgdp
    labvec_long[is] = labvec_long[is] * " Defl by GDP Def"
    labvec_short[is] = "Real_" * labvec_short[is]
    return x
end
deflate_series!(labvec_long::AbstractVector, labvec_short::AbstractVector,
                x::AbstractVector, price_def::AbstractVector,
                price_def_lfe::AbstractVector, price_def_pgdp,
                is::Integer, defcode::Val{0}) = x
"""
transform monthly data to quarterly data by averaging
##### Arguments
- `data_m::AbstractArray`: monthly data. `size(y, 1) == length(date)`
- `date_m::AbstractVector{Date}`: Vector of Date.
"""
function monthly_to_quarterly(data_m::AbstractArray, date_m::AbstractVector{Dates.Date}, ::MonthlyData)
    if size(data_m, 1) != length(date_m)
        @error "number of rows of data_m must be same as the length of date"
    end
    years = Dates.year.(date_m)
    date_m_quarter = transform_date.(date_m)
    T_q = length(unique(date_m_quarter))
    data_q = Array{Union{Missing, Float64}}(undef, T_q, size(data_m, 2))
    for t = 1:T_q
        id = (date_m_quarter .== [unique(date_m_quarter)[t]])
        data_q[t, :] = mean(data_m[id, :], dims=1)
    end
    return data_q, unique(date_m_quarter)
end
function monthly_to_quarterly(data_q::AbstractArray, date_q::AbstractVector{Dates.Date}, ::QuarterlyData)
    date_q_quarter = transform_date.(date_q)
    return data_q, date_q_quarter
end
transform_date(date::Dates.Date) = (Dates.year(date), get_quarter(date))
get_quarter(date::Dates.Date) = cld(Dates.month(date), 3)

# Level
transform(x::AbstractVector, ::Val{1}) = x
# First Difference
transform(x::AbstractVector, ::Val{2}) = vcat(missing, x[2:end].-x[1:end-1])
# Second Difference
transform(x::AbstractVector, ::Val{3}) = vcat(missing, missing, x[3:end]-2*x[2:end-1].+x[1:end-2])
# Log-level
transform(x::AbstractVector, ::Val{4}) = log.(x)
# Log-First-Difference
transform(x::AbstractVector, ::Val{5}) = transform(log.(x), Val(2))
# Log-Second-Difference
transform(x::AbstractVector, ::Val{6}) = transform(log.(x), Val(3))

function transform!(data::AbstractMatrix, tcode::AbstractVector)
    if size(data, 2) != length(tcode)
        @error "size(data, 2) must equal to length(tcode)"
    end
    for i = 1:size(data, 2)
        data[:, i] = transform(data[:, i], Val(tcode[i]))
    end
    return nothing
end
# Threshold multiple for IQR is 4.5
adjust_outlier!(x::AbstractVector, ::Val{1}, io_method::Integer) =
    adjust_outlier!(x, 4.5, io_method)
# Threshold multiple for IQR is 3
adjust_outlier!(x::AbstractVector, ::Val{2}, io_method::Integer) =
    adjust_outlier!(x, 3, io_method)
adjust_outlier!(x::AbstractVector, ::Val{0}, io_method::Integer) = nothing

function adjust_outlier!(x::AbstractVector, thr::Real, tflag::Integer)
    # compute median and IQR
    zm = median(skipmissing(x))
    iqr = quantile(skipmissing(x), 0.75) - quantile(skipmissing(x), 0.25)
    z = x[.!ismissing.(x)]
    # sort!(z)
    # tmp = z[ceil.(Int, length(z)*[1/4, 1/2, 3/4])]
    # zm = tmp[2]
    # iqr = tmp[3]-tmp[1]
    (iqr >= 1e-6) || @error "error in adjusting outlier"

    ya = abs.(x .- zm)
    i_outlier = ya .> (thr*iqr)
    adjust_x!(x, i_outlier, zm, thr, iqr, Val(tflag))
    return nothing
end
"""
replace with missing value
"""
function adjust_x!(x::AbstractVector, i_outlier, zm, thr, iqr, ::Val{0})
     x[i_outlier] = missing
     return nothing
end
"""
replace with maximum or minimum value
"""
function adjust_x!(x, i_outlier, zm, thr, iqr, ::Val{1})
    isign = Int.(x[i_outlier] .> 0) - Int.(x[i_outlier] .< 0)
    yt = zm .+ isign.*(thr*iqr)
    x[i_outlier] = yt
    return nothing
end
"""
replace with median value
"""
function adjust_x!(x, i_outlier, zm, thr, iqr, ::Val{2})
    x[i_outlier] = zm
    return nothing
end
"""
replace with local median (obs + or - 3 on each side)
"""
function adjust_x!(x, i_outlier, zm, thr, iqr, ::Val{3})
    # Compute rolling median
    iwin = 3  # Window on either side
    for i in findall(i_outlier)
        j1 = max(1, i-iwin)
        j2 = min(length(x), i+iwin)
        x[i] = median(skipmissing(x[j1:j2]))
    end
    return nothing
end
"""
replace with one-sided median (5 preceding obs)
"""
function adjust_x!(x, i_outlier, zm, thr, iqr, ::Val{4})
    # Compute rolling median
    iwin=5  # Window on one side
    for i in findall(i_outlier .& .!ismissing.(i_outlier))
        j1 = max(1, i-iwin)
        j2 = i
        x[i] = median(skipmissing(x[j1:j2]))
    end
    return nothing
end

readin_data(frequency::MonthlyData, datatype) =
    readin_monthly_data(true, 4, 2, 6, [1, 2, 3, 5], frequency, datatype)
readin_data(frequency::QuarterlyData, datatype) =
    readin_monthly_data(true, 4, 2, 5, [1, 2, 3, 5], frequency, datatype)
readdata(::MonthlyData) = readxlsheet("data/hom_fac_1.xlsx", "Monthly")
readdata(::QuarterlyData) = readxlsheet("data/hom_fac_1.xlsx", "Quarterly")
function readin_monthly_data(correct_outlier::Bool,
                             io_method::Integer, # Replacement of outliers,
                             ndesc::Integer,     # number of "description" rows in Excel file
                             ncodes::Integer,    # number of rows of "codes" in Excel file
                             cat_include::AbstractVector,
                             frequency::DataFrequency,
                             datatype::Symbol
                             )
    dnobs = frequency.nobs
    ns_m = frequency.ns
    ## read monthly data
    wholedata = readdata(frequency)

    maindata = wholedata[1:1+ndesc+ncodes+dnobs, 2:ns_m+1]
    date = Dates.Date.(wholedata[1+ndesc+ncodes+1:1+ndesc+ncodes+dnobs, 1])

    namevec, labvec_long, labvec_short, aggcode, tcode, defcode,
           outliercode, includecode, catcode = get_headers(maindata, frequency)
    datamat_m = maindata[1+ndesc+ncodes+1:end, :]
    datamat_m[.!isa.(datamat_m, Float64)] .= missing # set missing values
    datamat_m = convert(Array{Union{Missing, Float64}, 2}, datamat_m)

    # Price Deflators
    price_def, price_def_lfe, price_def_pgdp =
        get_deflators(namevec, datamat_m, frequency)

    # Standardize Killian Activity Index
    standardize_killian!(datamat_m, namevec, frequency)

    usedcols_id = (includecode .!=0) .& [in(cat, cat_include) for cat in floor.(catcode)]
    usedcols_id = get_usedcols_id(includecode, catcode, cat_include, Val(datatype))
    data_m = datamat_m[:, usedcols_id]
    bpdefcode = defcode[usedcols_id]
    bpoutlier = outliercode[usedcols_id]
    bplab_long = labvec_long[usedcols_id]
    bplab_short = labvec_short[usedcols_id]
    ns = size(data_m, 2)

    [deflate_series!(bplab_long, bplab_short, view(data_m, :, i),
                     price_def, price_def_lfe, price_def_pgdp, i, Val(bpdefcode[i])) for i =1:ns]
    data_q, date_q = monthly_to_quarterly(data_m, date, frequency)  # Temporally aggregated to quarterly
    bpdata_raw = copy(data_q)
    transform!(data_q, tcode[usedcols_id]) # Transform .. log, first difference, etc.
    bpdata_noa = copy(data_q)
    !correct_outlier || [adjust_outlier!(view(data_q, :, i), Val(bpoutlier[i]), io_method) for i = 1:ns]
    return data_q, bpdata_raw, bpdata_noa, date_q, catcode[usedcols_id],
           includecode[usedcols_id], namevec[usedcols_id]
end
get_usedcols_id(includecode, catcode, cat_include, ::Val{:Real}) =
     (includecode .!=0) .& [in(cat, cat_include) for cat in floor.(catcode)]
get_usedcols_id(includecode, catcode, cat_include, ::Val{:All}) = includecode .!=0

function get_headers(maindata, ::MonthlyData)
    namevec = uppercase.(convert(Vector{String}, maindata[1, :]))
    labvec_long = convert(Vector{String}, maindata[2, :])   # Vector of "long" labels
    labvec_short = convert(Vector{String}, maindata[3, :])  # Vector of "short" labels
    aggcode = Int.(maindata[4, :]) # Temporal aggregation code
    tcode = Int.(maindata[5, :])   # transformation code
    defcode = Int.(maindata[6, :]) # code for price deflation (nominal to real)
    outliercode = Int.(maindata[7, :]) # code for outlier adjustment
    includecode = Int.(maindata[8, :]) # code for use in factor estimation
    catcode = maindata[9, :]     # category code for ordering variables
    return namevec, labvec_long, labvec_short, aggcode, tcode, defcode,
           outliercode, includecode, catcode
end

function get_headers(maindata, ::QuarterlyData)
    namevec = uppercase.(convert(Vector{String}, maindata[1, :]))
    labvec_long = convert(Vector{String}, maindata[2, :])   # Vector of "long" labels
    labvec_short = convert(Vector{String}, maindata[3, :])  # Vector of "short" labels
    tcode = Int.(maindata[4, :])   # transformation code
    defcode = Int.(maindata[5, :]) # code for price deflation (nominal to real)
    outliercode = Int.(maindata[6, :]) # code for outlier adjustment
    includecode = Int.(maindata[7, :]) # code for use in factor estimation
    catcode = maindata[8, :]     # category code for ordering variables
    return namevec, labvec_long, labvec_short, nothing, tcode, defcode,
           outliercode, includecode, catcode
end

function get_deflators(namevec, datamat, ::MonthlyData)
    j = findall(namevec.=="PCEPI")[1] # PCE Price Deflator
    price_def = datamat[:, j]
    j = findall(namevec.=="PCEPILFE")[1] # PCE-xFE Price Deflator
    price_def_lfe = datamat[:, j]
    return price_def, price_def_lfe, NaN
end

function get_deflators(namevec, datamat, ::QuarterlyData)
    j = findall(namevec.=="PCECTPI")[1] # PCE Price Deflator
    price_def = datamat[:, j]
    j = findall(namevec.=="JCXFE")[1] # PCE Excl. food and energy
    price_def_lfe = datamat[:, j]
    j = findall(namevec.=="GDPCTPI")[1] # GDP Deflator
    price_def_pgdp = datamat[:, j]
    return price_def, price_def_lfe, price_def_pgdp
end

"""
Standardize Killian Activity Index
"""
function standardize_killian!(datamat_m, namevec, ::MonthlyData)
    j = findall(namevec.=="GLOBAL_ACT")[1] # Killian Index
    tmp = datamat_m[:, j]
    tmp1 = tmp[.!ismissing.(tmp)]
    tmp2 = (tmp1.-mean(tmp1))./std(tmp1)
    datamat_m[.!ismissing.(tmp), j] = tmp2
    return nothing
end

standardize_killian!(datamat_m, namevec, ::QuarterlyData) = nothing

function detrend_var!(data::AbstractArray, bw::BiWeight)
    trend = Array{Union{Missing, Float64}}(undef, size(data))
    for i = 1:size(data, 2)
        trend[:, i] = bi_weight_filter(data[:, i], bw.weight)
        data[:, i] = data[:, i] - trend[:, i]
    end
    return trend
end

function detrend_var!(data::AbstractArray, ::Mean)
    trend = repmat(mean(skipmissing(data), dims=1), size(data, 1), 1)
    data .= data .- trend
    return trend
end

detrend_var!(data::AbstractArray, ::NoDetrend) = fill(missing, size(data))

function bi_weight_filter(y::AbstractVector, weight::Real)
    T = length(y)
    ytrend = Array{Union{Missing, Float64}}(undef, T)
    trend = 1:T
    for t = findall(.!ismissing.(y))
        dt = (trend .- t)/weight
        bw_weight = 15/16*((1 .-dt.^2).^2) # bi-weight
        bw_weight[abs.(dt).>=1] .= 0
        bw_weight_ignore_missing = bw_weight[.!ismissing.(y)]
        bw_weight_ignore_missing = bw_weight_ignore_missing ./ sum(bw_weight_ignore_missing)
        ytrend[t] = dot(bw_weight_ignore_missing, y[.!ismissing.(y)])
    end
    return ytrend
end

"""
##### Arguments
- `detrend_method::DetrendMethod`: method for detrending.
- `datatype::Symbol`: type of read data, `:Real` or `:All`
"""
function readin_data(md::MonthlyData, qd::QuarterlyData,
                     detrend_method::DetrendMethod, datatype::Symbol)
    data_m, bpdata_raw_m, bpdata_noa_m, date_m, bpcatcode_m, inclcode_m, namevec_m =
        readin_data(md, datatype)

    data_q, bpdata_raw_q, bpdata_noa_q, date_q, bpcatcode_q, inclcode_q, namevec_q =
        readin_data(qd, datatype)

    all(date_m .== date_q) ||
        @error "inconsistent sample size for monthly and quarterly data"

    id = sortperm(vcat(bpcatcode_m, bpcatcode_q))
    bpdata = hcat(data_m, data_q)[:, id]
    bpdata_unfiltered = copy(bpdata)
    bpdata_trend = detrend_var!(bpdata, detrend_method)

    dataset = (bpdata_raw = hcat(bpdata_raw_m, bpdata_raw_q)[:, id],
               bpcatcode =vcat(bpcatcode_m, bpcatcode_q)[id],
               bpdata = bpdata,
               bpdata_unfiltered =bpdata_unfiltered,
               bpdata_noa = hcat(bpdata_noa_m, bpdata_noa_q)[:, id],
               bpdata_trend = bpdata_trend,
               inclcode = vcat(inclcode_m, inclcode_q)[id],
               bpnamevec = vcat(namevec_m, namevec_q)[id],
               calvec = calender_value.(date_q),
               calds = date_q)
    return dataset
end

calender_value(year::Integer, quarter::Integer) = year+(quarter-1)/4
calender_value(calvec::Tuple{Int, Int}) = calender_value(calvec[1], calvec[2])


text:
module TestDTM
    using Base.Test
    using Languages
    using TextAnalysis

    sample_file = joinpath(dirname(@__FILE__), "data", "poem.txt")

    fd = FileDocument(sample_file)
    sd = StringDocument(text(fd))

    crps = Corpus(Any[fd, sd])

    m = DocumentTermMatrix(crps)
    dtm(m)
    dtm(m, :dense)

    update_lexicon!(crps)

    m = DocumentTermMatrix(crps)
    dtm(m)
    dtm(m, :dense)

    tf_idf(dtm(m, :dense))

    dtv(crps[1], lexicon(crps))

    hash_dtv(crps[1], TextHashFunction())
    hash_dtv(crps[1])

    dtm(crps)
    sparse(dtm(crps))
    hash_dtm(crps)

    tdm(crps)
    hash_tdm(crps)
end


text:
using BinaryBuilder

name = "FLAC"
version = v"1.3.2"

# Collection of sources required to build Ipopt
sources = [
    "https://downloads.xiph.org/releases/flac/flac-1.3.2.tar.xz" =>
    "91cfc3ed61dc40f47f050a109b08610667d73477af6ef36dcad31c31a4a8d53f",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/flac-*

# For programs with older configure steps, we often need to update and regenerate them
update_configure_scripts
autoreconf -i -f

# Do the building dance
./configure --prefix=${prefix} --host=${target} --disable-avx
make -j${nproc} V=1
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line.
platforms = supported_platforms()

dependencies = [
    # We want libogg to power up our FLAC
    "https://github.com/staticfloat/OggBuilder/releases/download/v1.3.3-7/build_Ogg.v1.3.3.jl",
]

# The products that we will ensure are always built
products = prefix -> [
    LibraryProduct(prefix, "libFLAC", :libflac),
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
using Compose, Colors, CSV
using DataFrames

"""
    add_acquisition_points(imgpath::AbstractString, mime::AbstractString, x::Real, y::Real, fovx::Real, fovy::Real, pts::DataFrame)

A generic function for adding marker points to an image.  The image at `imgpath` of MIME-type `mime` centered at stage position `x`,`y`
with X and Y field-of-views `fovx`, `fovy` is annotated.  The point markers are placed at the coordinates pts[:,:X], pts[:,:Y].
"""
function add_acquisition_points(
    img::AbstractString,
    mime::AbstractString,
    x::Real,
    y::Real,
    fovx::Real,
    fovy::Real,
    pts::DataFrame,
)
    col = RGBA(colorant"yellow", 0.4) # Transparent yellow
    x0, y0 = x + 0.5 * fovx, y - 0.5 * fovy
    xx, yy = ((x0 .- pts[:, :X]) ./ fovx), (pts[:, :Y] .- y0) ./ fovy
    onimg =
        filter(i -> xx[i] > 0.0 && xx[i] < 1.0 && yy[i] > 0.0 && yy[i] < 1.0, eachindex(xx))
    compose(
        context(units = UnitBox(0.0, 0.0, 1.0, 1.0)),
        (
            context(),
            circle(xx[onimg], yy[onimg], [0.005]),
            Compose.stroke(col),
            linewidth(0.003w),
            fill("transparent"),
        ),
        (context(), bitmap(mime, read(img), 0.0, 0.0, 1.0, 1.0)),
    )
end

"""
    add_annotations(imgpath::AbstractString, modes::Symbol...; color = RGBA(colorant"yellow", 0.4), mime::AbstractString = "image/png")

Adds annotations to a SEMAntics image from the `images.txt` and `spectra.txt` meta-data files.

  Modes
  * `:ScaleBar`           Add a scale bar
  * `:AcquisitionPoints`  Mark all spectrum acquisition points.
"""
function add_annotations(
    imgpath::AbstractString,
    modes::Symbol...;
    color = RGBA(colorant"yellow", 0.4),
    mime::AbstractString = "image/png",
)
    imgdir, imgfn = splitdir(imgpath)
    df = CSV.read(joinpath(imgdir, "images.txt")) |> DataFrame
    img = splitext(imgfn)[1][1:end-3]
    rws = filter(r -> r[:Name] == img, df)
    if nrow(rws) > 0
        r = last(rws)
        aspect, fov = r[:YDim] / r[:XDim], r[:FOV]
        contexts = Any[context(units = UnitBox(0.0, 0.0, 1.0, aspect))]
        if :ScaleBar in modes
            sc = 10.0^floor(log10(fov))
            scale =
                fov / sc > 8.0 ? 5.0 * sc :
                (fov / sc > 4.0 ? 2.0 * sc : (fov / sc > 2.0 ? 1.0 * sc : 0.5 * sc))
            backcol = RGBA(0.8, 0.8, 0.8, 0.1)
            unit = (sc < 0.001 ? "nm" : (sc < 1.0 ? "μm" : "mm"))
            val = sc < 0.001 ? scale * 1.0e6 : (sc < 1.0 ? scale * 1.0e3 : scale)
            push!( # Scale bar
                contexts,
                (
                    context(),
                    rectangle(0.96 - scale / fov, 0.036, scale / fov, 0.006),
                    fill(color),
                ),
            )
            push!( # Text for scale bar
                contexts,
                (
                    context(),
                    text(0.96 - 0.5 * scale / fov, 0.06, "$val $unit", hcenter, vtop),
                    fontsize(0.04h),
                    font("Helvetica"),
                    fill(color),
                ),
            )
            push!( # Transparent gray background
                contexts,
                (
                    context(),
                    rectangle(0.954 - scale / fov, 0.03, scale / fov + 0.012, 0.07),
                    fill(backcol),
                ),
            )
        end
        if :AcquisitionPoints in modes
            x0, y0 = r[:X] + 0.5 * fov, r[:Y] - 0.5 * fov * aspect
            fovx, fovy = r[:FOV], aspect * r[:FOV]
            pts = CSV.read(joinpath(imgdir, "spectra.txt")) |> DataFrame
            xx, yy = ((x0 .- pts[:, :X]) ./ fovx), (pts[:, :Y] .- y0) ./ fovy
            onimg = filter(
                i -> xx[i] > 0.0 && xx[i] < 1.0 && yy[i] > 0.0 && yy[i] < 1.0,
                eachindex(xx),
            )
            push!( # X markers
                contexts,
                (
                    context(),
                    xgon(xx[onimg], yy[onimg], [0.02 * w], [4]),
                    stroke(color),
                    linewidth(0.003 * w),
                    fill("transparent"),
                ),
            )
        end
        # Image
        push!(contexts, (context(), bitmap(mime, read(imgpath), 0.0, 0.0, 1.0, aspect)))
        return compose(contexts...)
    end
end

max_stars_repo_name: NicholasWMRitchie/NeXLSpectrum

max_issues_repo_name: NicholasWMRitchie/NeXLSpectrum

max_forks_repo_name: NicholasWMRitchie/NeXLSpectrum


text:
using ..Ocean: surface_flux

"""
    ocean_temperature_boundary_state!(::Union{NumericalFluxFirstOrder, NumericalFluxGradient}, ::Insulating, ::HBModel)

apply insulating boundary condition for temperature
sets transmissive ghost point
"""
function ocean_temperature_boundary_state!(
    nf::Union{NumericalFluxFirstOrder, NumericalFluxGradient},
    bc_temperature::Insulating,
    ocean,
    Q⁺,
    A⁺,
    n,
    Q⁻,
    A⁻,
    t,
)
    Q⁺.θ = Q⁻.θ

    return nothing
end

"""
    ocean_temperature_boundary_state!(::NumericalFluxSecondOrder, ::Insulating, ::HBModel)

apply insulating boundary condition for velocity
sets ghost point to have no numerical flux on the boundary for κ∇θ
"""
@inline function ocean_temperature_boundary_state!(
    nf::NumericalFluxSecondOrder,
    bc_temperature::Insulating,
    ocean,
    Q⁺,
    D⁺,
    A⁺,
    n⁻,
    Q⁻,
    D⁻,
    A⁻,
    t,
)
    Q⁺.θ = Q⁻.θ
    D⁺.κ∇θ = n⁻ * -0

    return nothing
end

"""
    ocean_temperature_boundary_state!(::Union{NumericalFluxFirstOrder, NumericalFluxGradient}, ::TemperatureFlux, ::HBModel)

apply temperature flux boundary condition for velocity
applies insulating conditions for first-order and gradient fluxes
"""
function ocean_temperature_boundary_state!(
    nf::Union{NumericalFluxFirstOrder, NumericalFluxGradient},
    bc_velocity::TemperatureFlux,
    ocean,
    args...,
)
    return ocean_temperature_boundary_state!(nf, Insulating(), ocean, args...)
end

"""
    ocean_temperature_boundary_state!(::NumericalFluxSecondOrder, ::TemperatureFlux, ::HBModel)

apply insulating boundary condition for velocity
sets ghost point to have specified flux on the boundary for κ∇θ
"""
@inline function ocean_temperature_boundary_state!(
    nf::NumericalFluxSecondOrder,
    bc_temperature::TemperatureFlux,
    ocean,
    Q⁺,
    D⁺,
    A⁺,
    n⁻,
    Q⁻,
    D⁻,
    A⁻,
    t,
)
    Q⁺.θ = Q⁻.θ
    D⁺.κ∇θ = n⁻ * surface_flux(ocean.problem, A⁻.y, Q⁻.θ)

    return nothing
end


text:
using Test, LinearAlgebra, Statistics, DiscriminantAnalysis

const DA = DiscriminantAnalysis

function random_centroids(T::Type{<:AbstractFloat}, m::Int, p::Int)
    M = zeros(T, p, m)
    for k = 1:m
        M[:, k] = T[rand() < 0.5 ? -2k : 2k for i = 1:p]
    end
    return M
end

function random_data(T::Type{<:AbstractFloat}, nₘ::Vector{Int}, p::Int)
    n = sum(nₘ)
    m = length(nₘ)

    X = zeros(T, p, n)
    M = zeros(T, p, m)

    ix = sortperm(rand(n))

    y = [k for (k, nₖ) in enumerate(nₘ) for i = 1:nₖ][ix]
    M = random_centroids(T, m, p)

    for (k, nₖ) in enumerate(nₘ)
        Xₖ = randn(T, p, nₖ)

        Xₖ .-= mean(Xₖ, dims=2)
        Xₖ .+= view(M, :, k)

        X[:, y .== k] = Xₖ
    end

    return (X, y, M)
end

function random_cov(T::Type{<:AbstractFloat}, p::Integer)
    Q = qr(randn(T, p, p)).Q
    D = Diagonal(2p*rand(T, p))

    Σ = Symmetric(transpose(Q)*D*Q)

    W_svd = inv(√(D))*Q
    W_chol = copy(transpose(inv(cholesky(Σ).U)))

    return (Σ, W_svd, W_chol)
end

function perturb(M::Matrix{T}, tol::Real=convert(T,0.05)) where {T <: Real}
    ϵ = convert(T, tol)
    jitter = range(one(T)-ϵ, stop=one(T)+ϵ, length=10000)
    return broadcast(x -> x*rand(jitter), M)
end

@testset "DiscriminantAnalysis.jl" begin
    include("test_common.jl")
    include("test_whiten.jl")
    include("test_parameters.jl")
    include("test_classifiers_linear.jl")
end


text:
using Random
include("boardrobots.jl")
include("solver.jl")

"""
    write_board(board::Matrix{Int}, filename::String)

Write the challenge `board` in the file `filename`.
"""
function write_board(board::Matrix{Int}, filename::String)
    open(filename, "w") do io
        for i in 1:5
            for j in 1:5
                if board[i,j] != 0
                    write(io, "$(board[i,j])")
                else
                    write(io, " ")
                end
            end
            write(io, "\n")
        end
    end
    return nothing
end

"""
    generate_random_instance(nb_move::Int)

Return an instance that is doable in `nb_move` moves.
"""
function generate_random_instance(nb_move::Int)
    instance_found = false
    board = initialize_board()
    while !instance_found
        board = initialize_board()
        i = rand(1:5)
        j = rand(1:5)
        while i == 3 && j == 3
            i = rand(1:5)
            j = rand(1:5)
        end
        board[i,j] = 1
        for robot in 2:6
            i = rand(1:5)
            j = rand(1:5)
            while board[i,j] != 0
                i = rand(1:5)
                j = rand(1:5)
            end
            board[i,j] = robot
            if !isempty(find_solution(board)) && length(find_solution(board)[1]) == nb_move
                instance_found = true
                break
            end
        end
    end
    robot = 2
    for i in 1:5
        for j in 1:5
            if !(board[i,j] in [0,1])
                board[i,j] = robot
                robot += 1
            end
        end
    end
    return board
end

max_stars_repo_name: remi-garcia/LunarLockout

max_issues_repo_name: remi-garcia/LunarLockout

max_forks_repo_name: remi-garcia/LunarLockout


text:
"""
```
cleared_img = clearborder(img)
cleared_img = clearborder(img, width)
cleared_img = clearborder(img, width, background)
```
Returns a copy of the original image after clearing objects connected to the border of the image.
Parameters:
 -  img          = Binary/Grayscale input image
 -  width        = Width of the border examined (Default value is 1)
 -  background   = Value to be given to pixels/elements that are cleared (Default value is 0)
"""
function clearborder(img::AbstractArray, width::Integer=1, background::Integer=0)

    for i in size(img)
        if(width > i)
            throw(ArgumentError("Border width must not be greater than size of the image."))
        end
    end

    connectivity = ntuple(i -> 3, ndims(img))
    labels = label_components(img,trues(connectivity))
    number = maximum(labels) + 1

    dimensions = size(img)
    outerrange = CartesianIndices(map(i -> 1:i, dimensions))
    innerrange = CartesianIndices(map(i -> 1+width:i-width, dimensions))

    border_labels = Set{Int64}()
    for i in EdgeIterator(outerrange,innerrange)
        push!(border_labels, labels[i])
    end

    new_img = similar(img)
    for itr in eachindex(labels)
        if labels[itr] in border_labels
            new_img[itr] = background
        else
        	new_img[itr] = img[itr]
        end
    end

    return new_img

end


text:
## ---
using Transformers: Vocabulary, encode
using DataStructures: SortedSet

function openDataset(f)
    s = open(f) do file
        read(file, String)
    end
    asDataset(s)
end

function asDataset(s)
    @show datasize=length(s)
    chars=[SortedSet(s)...]
    unksym = '□'
    vocab = Vocabulary(chars, unksym)

    getBlockFrom(start, size)  =vocab(collect(s[start:start+size-1]))
    function sampledata(blocksize=100)
        #TODO if blocksize is greater that the data, then throw exception
        lastbyte = datasize-blocksize   # last byte that can be the start for this sample
        getBlockFrom(rand(1:lastbyte), blocksize)  #and encoded block of data
    end
    datasize, vocab, sampledata
end

max_stars_repo_name: daryoush/minGPT

max_issues_repo_name: daryoush/minGPT

max_forks_repo_name: daryoush/minGPT


text:
# This file includes an example file for 1D integrate.

using FractalTools
using Makie

# Generate data
xi, dx, xf = 0., 0.01, 10.
x = collect(xi : dx : xf)
y = sin.(x)
pts = collect.(zip(x, y))
res = integrate(pts, Interp1D(0.0001))


text:
flatten{T}(x::Vector{T}) = x
flatten(X) = reshape(X, prod(size(X)))

# convert vector of data samples into a single data matrix, flatten samples if needed
datamatrix{T,N}(samples::Vector{Array{T,N}}) =
    hcat([flatten(sample) for sample in samples]...)

# convert single data matrix into a vector of samples
# samplevec{T}(datamat::Matrix{T}, sampledims::(Int...)) =
#     Matrix{T}[reshape(datamat[:, i], sampledims) for i=1:size(datamat, 2)]

function gs_orthonorm(M)
    O = zeros(eltype(M), size(M))
    ncol = size(M, 2)
    tol = 100*eps(Float64)
    k = 1
    for i=1:ncol
        v = M[:, i]   # column to orthonormalize
        # subtract projections over previous vectors
        for j=1:k-1
            v = v - dot(O[:,j], v) * O[:,j]
        end
        # only keep nonzero vectors
        n = norm(v)
        if n > tol
            O[:, k] = v / n
            k = k + 1
        end
    end
    return O[:, 1:k-1]
end

function StatsBase.describe{T,N}(A::Array{T,N})
    println("size   : $(size(A))")
    println("nzs    : $(countnz(A) / length(A))")
    println("min    : $(minimum(A))")
    println("mean   : $(mean(A))")
    println("median : $(median(A))")
    println("max    : $(maximum(A))")
    println("sum    : $(sum(A))")
end


text:
export CyclotomicExt, cyclotomic_extension

################################################################################
#
#  Type definition
#
################################################################################
mutable struct CyclotomicExt
  k::AnticNumberField
  n::Int
  Kr::Hecke.NfRel{nf_elem}
  Ka::AnticNumberField
  mp::Tuple{NfToNfRel, NfToNfMor}

  kummer_exts::Dict{Set{fmpz}, Tuple{Vector{NfOrdIdl}, KummerExt}}
                      #I save the kummer extensions used in the class field construction
                      #The keys are the factors of the minimum of the conductor
  function CyclotomicExt()
    return new()
  end
end

function Base.show(io::IO, c::CyclotomicExt)
  print(io, "Cyclotomic extension by zeta_$(c.n) of degree $(degree(c.Ka))")
end

################################################################################
#
#  Interface and creation
#
################################################################################

@doc Markdown.doc"""
    cyclotomic_extension(k::AnticNumberField, n::Int) -> CyclotomicExt
Computes $k(\zeta_n)$, in particular, a structure containing $k(\zeta_n)$
both as an absolute extension, as a relative extension (of $k$) and the maps
between them.
"""
function cyclotomic_extension(k::AnticNumberField, n::Int; cached::Bool = true, compute_maximal_order::Bool = true, compute_LLL_basis::Bool = true)
  Ac = CyclotomicExt[]
  if cached
    try
      Ac = Hecke._get_cyclotomic_ext_nf(k)::Vector{CyclotomicExt}
      for i = Ac
        if i.n == n
          return i
        end
      end
    catch e
      if !(e isa AbstractAlgebra.AccessorNotSetError)
        rethrow(e)
      end
    end
  end

  kt, t = PolynomialRing(k, "t", cached = false)
  c = CyclotomicExt()
  c.kummer_exts = Dict{Set{fmpz}, Tuple{Vector{NfOrdIdl}, KummerExt}}()
  c.k = k
  c.n = n

  if n <= 2
    #Easy, just return the field
    Kr = number_field(t+1, cached = false, check = false)[1]
    if compute_maximal_order
      Ok = maximal_order(k)
      if compute_LLL_basis
        lll(Ok)
      end
    end
    abs2rel = NfToNfRel(k, Kr, gen(k), k(-1), Kr(gen(k)))
    small2abs = id_hom(k)
    c.Kr = Kr
    c.Ka = k
    c.mp = (abs2rel, small2abs)
    if cached
      push!(Ac, c)
      Hecke._set_cyclotomic_ext_nf(k, Ac)
    end
    return c
  end

  ZX, X = PolynomialRing(FlintZZ, cached = false)
  f = cyclotomic(n, X)
  fk = change_ring(f, kt)
  if n < 5
    #For n = 3, 4 the cyclotomic polynomial has degree 2,
    #so we can just ask for the roots.
    if !isone(gcd(degree(fk), degree(k))) && !isone(gcd(discriminant(maximal_order(k)), n)) && !istotally_real(k)
      rt = _roots_hensel(fk, max_roots = 1, isnormal = true)
    else
      rt = nf_elem[]
    end
    if length(rt) == 1
      #The polynomial splits completely!
      Kr, gKr = number_field(t - rt[1], cached = false, check = false)
      abs2rel = NfToNfRel(k, Kr, gen(k), rt[1], Kr(gen(k)))
      small2abs = id_hom(k)
      if compute_maximal_order
        Ok = maximal_order(k)
        if compute_LLL_basis
          lll(Ok)
        end
      end
      c.Kr = Kr
      c.Ka = k
      c.mp = (abs2rel, small2abs)
    else
      Kr, Kr_gen = number_field(fk, "z_$n", cached = false, check = false)
      Ka, abs2rel, small2abs = Hecke.absolute_field(Kr, false)
      if compute_maximal_order
        Zk = maximal_order(k)
        b_k = basis(Zk, k)
        B_k = Vector{nf_elem}(undef, degree(Ka))
        for i = 1:length(b_k)
          B_k[i] = small2abs(b_k[i])
        end
        g = abs2rel\(Kr_gen)
        for j = 1:length(b_k)
          B_k[j+length(b_k)] = B_k[j]*g
        end
        ZKa = Hecke.NfOrd(B_k)
        ZKa.disc = discriminant(f)^degree(k)*discriminant(Zk)^degree(fk)
        ZKa.index = root(divexact(numerator(discriminant(Ka)), ZKa.disc), 2)
        ZKa.gen_index = fmpq(ZKa.index)
        for (p,v) = factor(gcd(discriminant(Zk), fmpz(n))).fac
          ZKa = pmaximal_overorder(ZKa, p)
        end
        ZKa.ismaximal = 1
        if compute_LLL_basis
          ZKa = lll(ZKa)
        end
        Hecke._set_maximal_order_of_nf(Ka, ZKa)
      end
      c.Kr = Kr
      c.Ka = Ka
      c.mp = (abs2rel, small2abs)
    end
    if cached
      push!(Ac, c)
      Hecke._set_cyclotomic_ext_nf(k, Ac)
    end
    return c
  end

  if gcd(degree(fk), degree(k)) != 1 && !isone(gcd(discriminant(maximal_order(k)), n))
    lf = factor(fk)
    fk = first(keys(lf.fac))
  end

  Kr, Kr_gen = number_field(fk, "z_$n", cached = false, check = false)
  if degree(fk) != 1
    Ka, abs2rel, small2abs = Hecke.absolute_field(Kr, false)

    if compute_maximal_order
      # An equation order defined from a factor of a
      # cyclotomic polynomial is always maximal by Dedekind
      # hence the product basis should be maximal at all primes except
      # for all p dividing the gcd()
      Zk = maximal_order(k)
      b_k = basis(Zk, k)
      B_k = Vector{nf_elem}(undef, degree(Ka))
      for i = 1:length(b_k)
        B_k[i] = small2abs(b_k[i])
      end
      g = abs2rel\(Kr_gen)
      s = length(b_k)+1
      for i = 1:degree(fk)-1
        for j = 1:degree(k)
          B_k[s] = B_k[s-degree(k)]*g
          s += 1
        end
      end
      ZKa = Hecke.NfOrd(B_k)
      if degree(Kr) == euler_phi(n)
        ZKa.disc = (discriminant(Zk)^euler_phi(n))*discriminant(f)^degree(k)
        ZKa.index = root(divexact(numerator(discriminant(Ka)), discriminant(ZKa)), 2)
        ZKa.gen_index = fmpz(ZKa.index)
      else
        ZKa.disc = (discriminant(Zk)^degree(Kr))*numerator(norm(discriminant(fk)))
        ZKa.index = root(divexact(numerator(discriminant(Ka)), discriminant(ZKa)), 2)
        ZKa.gen_index = fmpz(ZKa.index)
      end
      for (p, v) in factor(gcd(discriminant(Zk), fmpz(n)))
        ZKa = pmaximal_overorder(ZKa, p)
      end
      if compute_LLL_basis
        ZKa = lll(ZKa)
      end
      ZKa.ismaximal = 1
      Hecke._set_maximal_order_of_nf(Ka, ZKa)
    end
  else
    Ka = k
    abs2rel = NfToNfRel(Ka, Kr, gen(Ka), -coeff(fk, 0), Kr(gen(Ka)))
    small2abs = id_hom(k)
    if compute_maximal_order
      Ok = maximal_order(k)
      if compute_LLL_basis
        lll(Ok)
      end
    end
  end

  c.Kr = Kr
  c.Ka = Ka
  c.mp = (abs2rel, small2abs)
  if cached
    push!(Ac, c)
    Hecke._set_cyclotomic_ext_nf(k, Ac)
  end
  return c

end

absolute_field(C::CyclotomicExt) = C.Ka
base_field(C::CyclotomicExt) = C.k

################################################################################
#
#  Automorphisms for cyclotomic extensions
#
################################################################################
@doc Markdown.doc"""
    automorphisms(C::CyclotomicExt; gens::Vector{NfToNfMor}) -> Vector{NfToNfMor}
Computes the automorphisms of the absolute field defined by the cyclotomic extension, i.e. of absolute_field(C).
It assumes that the base field is normal. gens must be a set of generators for the automorphism group of the base field of C
"""
function automorphisms(C::CyclotomicExt; gens::Vector{NfToNfMor} = small_generating_set(automorphisms(base_field(C))), copy::Bool = true)

  if degree(absolute_field(C)) == degree(base_field(C))
    return automorphisms(C.Ka, copy = copy)
  end
  genK = C.mp[1](gen(C.Ka))
  gnew = Hecke.NfToNfMor[]
  #First extend the old generators
  for g in gens
    ng = Hecke.extend_to_cyclotomic(C, g)
    na = hom(C.Ka, C.Ka, C.mp[1]\(ng(genK)), check = false)
    push!(gnew, na)
  end
  #Now add the automorphisms of the relative extension
  R = ResidueRing(FlintZZ, C.n, cached = false)
  U, mU = unit_group(R)
  if iscyclic(U)
    k = degree(C.Kr)
    expo = divexact(euler_phi(fmpz(C.n)), k)
    l = hom(C.Kr, C.Kr, gen(C.Kr)^Int(lift(mU(U[1])^expo)), check = false)
    l1 = hom(C.Ka, C.Ka, C.mp[1]\(l(C.mp[1](gen(C.Ka)))), check = false)
    push!(gnew, l1)
  else
    f = C.Kr.pol
    s, ms = sub(U, [x for x in U if iszero(f(gen(C.Kr)^Int(lift(mU(x)))))], false)
    S, mS = snf(s)
    for t = 1:ngens(S)
      l = hom(C.Kr, C.Kr, gen(C.Kr)^Int(lift(mU(ms(mS(S[t]))))), check = false)
      push!(gnew, hom(C.Ka, C.Ka, C.mp[1]\(l(genK)), check = false))
    end
  end
  auts = closure(gnew, degree(C.Ka))
  Hecke._set_automorphisms_nf(C.Ka, auts)
  if copy
    return Base.copy(auts)
  else
    return auts
  end

end


text:
n = rand(1:10)
A = matrixdepot("hilb", n)
B = matrixdepot("lotkin", n)

@test A[2:n, :] == B[2:n, :]
@test B[1,:][:] == ones(n)
println("'lotkin' passed test...")


text:
using PSCOPF

using Test
using Dates

@testset verbose=true "test_ech_generator" begin

    @testset "mode_1" begin
        network = PSCOPF.Networks.Network()
        ts1 = DateTime("2015-01-01T11:00:00")
        TS = PSCOPF.create_target_timepoints(ts1)
        mode = PSCOPF.PSCOPF_MODE_1

        list_ech = PSCOPF.generate_ech(network, TS, mode)

        EXPECTED_ECH = [ts1 - Dates.Hour(4),
                        ts1 - Dates.Hour(1),
                        ts1 - Dates.Minute(30),
                        ts1 - Dates.Minute(15),
                        ts1]
        @test length(list_ech) == length(EXPECTED_ECH)
        @test list_ech == EXPECTED_ECH
    end

    #=
        generate_ech sorts the received TS, and looks at the first target timepoint
    =#
    @testset "mode_1_unsorted_TS" begin
        network = PSCOPF.Networks.Network()
        #unsorted TS
        TS = [DateTime("2015-01-01T11:30:00"), DateTime("2015-01-01T11:45:00"), DateTime("2015-01-01T11:00:00"), DateTime("2015-01-01T11:15:00")]
        mode = PSCOPF.PSCOPF_MODE_1

        list_ech = PSCOPF.generate_ech(network, TS, mode)

        ts1 = DateTime("2015-01-01T11:00:00")
        EXPECTED_ECH = [ts1 - Dates.Hour(4),
                        ts1 - Dates.Hour(1),
                        ts1 - Dates.Minute(30),
                        ts1 - Dates.Minute(15),
                        ts1]
        @test length(list_ech) == length(EXPECTED_ECH)
        @test list_ech == EXPECTED_ECH
    end

    @testset "mode_2" begin
        network = PSCOPF.Networks.Network()
        ts1 = DateTime("2015-01-01T11:00:00")
        TS = PSCOPF.create_target_timepoints(ts1)
        mode = PSCOPF.PSCOPF_MODE_2

        list_ech = PSCOPF.generate_ech(network, TS, mode)

        EXPECTED_ECH = [ts1 - Dates.Hour(4),
                        ts1 - Dates.Hour(1),
                        ts1 - Dates.Minute(30),
                        ts1 - Dates.Minute(15),
                        ts1]
        @test length(list_ech) == length(EXPECTED_ECH)
        @test list_ech == EXPECTED_ECH
    end

    @testset "mode_3" begin
        network = PSCOPF.Networks.Network()
        ts1 = DateTime("2015-01-01T11:00:00")
        TS = PSCOPF.create_target_timepoints(ts1)
        mode = PSCOPF.PSCOPF_MODE_3

        list_ech = PSCOPF.generate_ech(network, TS, mode)

        EXPECTED_ECH = [ts1 - Dates.Hour(4),
                        ts1 - Dates.Hour(1),
                        ts1 - Dates.Minute(30),
                        ts1 - Dates.Minute(15),
                        ts1]
        @test length(list_ech) == length(EXPECTED_ECH)
        @test list_ech == EXPECTED_ECH
    end

end

#=

#TEst
    planningTSO = planningRefTSO
    planningMarket = planningRefMarket

    data = (Grid, ECH, TS, ...)
    planning = (planningTSO, planningMarket)
    OptimResult = run(MarketMode1, planning, data)

    tester les résutlats

#TEst update (démarrage des unités)

    planningTSO = planningRefTSO
    planningMarket = planningRefMarket

    data = (Grid, ECH, TS, ...)
    planning = (planningTSO, planningMarket)

    update_imposition(planning, OptimResult)
    update_limitation(planning, OptimResult)
    update_units(planning, OptimResult)
    update_reserve(planning, OptimResult)

    test sur planning
    #Decision ferme/
    niveau["unit1"][] = 15
    niveau["unit1"] = 15

=#

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
## Running code and tests for Dieter.jl

# %% Load packages
using JuMP
using CPLEX
using Dieter
using ColorSchemes
using Plots

# %% Instantiate data model
projectpath = joinpath(ENV["HOME"],"Documents/Projects/ESM/Dieter.jl/")
datapath = joinpath(projectpath,"testdata/")
rdir = joinpath(projectpath,"results_of_test_case/")

data_instance = Dict{String,Any}()
dtr = InitialiseDieterModel(DieterModel, data_instance)
dtr.settings[:datapath] = datapath
# data_instance["files"] = Dict{String,String}()
# fileDict = dtr.data["files"]

initialise_data_file_dict!(dtr)
check_files_exist(dtr.data["files"])

# Modify the :min_res setting over [0,100] and rerun to see comparison.
dtr.settings[:min_res] = 50
# dtr.settings[:ev] = missing
# dtr.settings[:heat] = missing
# dtr.settings[:h2] = missing

dfDict = parse_data_to_model!(dtr)

# %% Initialise model
# include("model.jl")

# Construct an optimizer factory
# solver = with_optimizer(Clp.Optimizer)
# solver = with_optimizer(Gurobi.Optimizer)
solver = JuMP.with_optimizer(CPLEX.Optimizer)
build_model!(dtr,solver)

# %% Specialised solver settings
# CPLEX:
JuMP.MOI.set(dtr.model, JuMP.MOI.RawParameter("CPX_PARAM_THREADS"), 4)        #  number of threads
JuMP.MOI.set(dtr.model, JuMP.MOI.RawParameter("CPX_PARAM_PARALLELMODE"), 1)   #   -1: Opportunistic parallel, 0: Automatic parallel, 1: Deterministic
JuMP.MOI.set(dtr.model, JuMP.MOI.RawParameter("CPX_PARAM_LPMETHOD"), 4)       #  0: auto, 1: primal simplex, 2: dual, 3: network, 4: barrier
JuMP.MOI.set(dtr.model, JuMP.MOI.RawParameter("CPX_PARAM_BARCROSSALG"), 2)    #  0: automatic, 1: primal, 2: dual
JuMP.MOI.set(dtr.model, JuMP.MOI.RawParameter("CPX_PARAM_SOLUTIONTYPE"), 2)   #  Specifies type of solution (basic or non basic) that CPLEX produces
JuMP.MOI.set(dtr.model, JuMP.MOI.RawParameter("CPX_PARAM_BAREPCOMP"), 1e-8)   # Sets the tolerance on complementarity for convergence; default: 1e-8.

# %% Solve the model and generate results
solve_model!(dtr)
generate_results!(dtr)

# %% Analysis
# include("analysis.jl")
df_summ = summarize_result(dtr,del_zeros=false)

# %% Save results to file
# include("save.jl")

save_results(dtr, rdir)

# %% Merge results with other runs
# include("merge.jl")
post_process_results(rdir)

# %% Plotting
color_dict = Dict()
default_colour = ColorSchemes.leonardo
# ColorScheme([Colors.RGB(0.0, 0.0, 0.0), Colors.RGB(1.0, 1.0, 1.0)],
      # "custom", "twotone, black and white");
      color_dict["Wind Offshore"] = [default_colour[10]];
      color_dict["Wind Onshore"] = [default_colour[20]];
      color_dict["Solar PV"] = [default_colour[30]];
      color_dict["Storages"] = [default_colour[5]];
      color_dict["Hydrogen"] = [default_colour[15]];
      color_dict["Curtailment"] = [default_colour[25]];

marker = (:hexagon, 5, 0.6, :green, stroke(3, 0.2, :black, :dot))

plot_all(rdir,color_dict,sector=:ev,marker=marker,legend=true)


text:
#Design example for a continuous-time lowpass ΔΣ ADC
using RSDeltaSigmaPort
using RSDeltaSigmaPort.EasyPlot #set, cons
import Printf: @sprintf
j=im

#==Baseband modulator (continuous-time implementation)
===============================================================================#
nlev=2
#dsm = RealDSM(order=3, OSR=100, M=nlev-1, f0=0, opt=0, Hinf=1.3, form=:FB)
#dsm = RealDSM(order=3, OSR=32, M=nlev-1, f0=0, opt=0, Hinf=1.5, form=:FB)
dsm = RealDSM(order=3, OSR=32, M=nlev-1, f0=0, opt=2, Hinf=1.5, form=:FB)
umax = 0.83

#Parameters for the continuous-time implementation:
tdac = [0 1] #DAC timing. [0 1] means zero-delay non-return-to-zero

println("\n*** $(dsm.order)th-Order Continuous-Time Lowpass Example")

@info("Performing NTF synthesis..."); flush(stdout); flush(stderr)
#-------------------------------------------------------------------------------
NTF0 = synthesizeNTF(dsm) #Optimized zero placement
#println("p:"); display(NTF0.p)
#println("z:"); display(NTF0.z)
plotcoll_NTF = documentNTF(dsm, NTF0)
println("\tdone.")
saveimage(:png, "dsexample2_NTF.png", plotcoll_NTF, AR=2/1, width=900)
displaygui(plotcoll_NTF)
displaygui(plotNTF(NTF0, dsm.OSR)) #ADDED MALaforge

@info("Performing time-domain simulations..."); flush(stdout); flush(stderr)
#-------------------------------------------------------------------------------
#Example spectrum
#ftest = 0.3 * 0.5/dsm.OSR/3 #~1/3 of the way across the passband
ftest = (0.5/dsm.OSR)/3 #~1/3 of the way across the passband
plot = plotExampleSpectrum(dsm, NTF0, ftest=ftest,N=2^12) #16
#title('Example Spectrum');
saveimage(:png, "dsexample2_PSD.png", plot, AR=2/1, width=900)
displaygui(plot)

#SQNR plot
snrinfo = calcSNRInfo(dsm, NTF=NTF0)
plot = plotSNR(snrinfo, dsm)
	set(plot, xyaxes=set(xmin=-100, ymin=0, ymax=100))
println("\tdone.")
saveimage(:png, "dsexample2_SQNR.png", plot, AR=2/1, width=900)
displaygui(plot)

@info("Mapping to continuous-time..."); flush(stdout); flush(stderr)
#-------------------------------------------------------------------------------
(ABCDc,tdac2) = realizeNTF_ct(NTF0, dsm.form, tdac)
(Ac, Bc, Cc, Dc) = partitionABCD(ABCDc)
sys_c = _ss(Ac,Bc,Cc,Dc)
println("\tdone.")
#Verify that the sampled pulse response of the CT loop filter
#matches the impulse response of the DT prototype
n_imp = 10
y = -impL1(NTF0, n_imp) #Negate impL1 to make the plot prettier
	y_max = maximum(real.(y))
plot = plotLollipop(0:n_imp, y, color=:blue)
	plot.title = "Loop filter pulse/impulse responses (negated)"
yl = floor(0.9*y_max); color=:blue
	ylw = waveform([0.5, 2], [yl, yl])
simglyph = cons(:a, glyph=set(shape=:o, size=1.5, color=color, fillcolor=color))
push!(plot,
	cons(:wfrm, ylw, simglyph, line=set(style=:solid, color=color, width=2), label=""),
	cons(:atext, "   discrete-time", x=2, y=yl, align=:bl),
)
dt = 1/16
t = collect(0:dt:n_imp)
yy = -pulse(sys_c, [0 0;tdac], dt, n_imp*1.0)
	yy = squeeze(yy)
	yyw = waveform(t, yy)
yl = floor(0.7*y_max); color=:green
	ylw = waveform([0.5, 2], [yl, yl])
push!(plot,
	#Add yy when it works
	cons(:wfrm, yyw, line=set(style=:solid, color=color, width=2), label=""),
	cons(:wfrm, ylw, line=set(style=:solid, color=color, width=2), label=""),
	cons(:atext, "   continuous-time", x=2, y=yl, align=:bl),
)
displaygui(plot)

#Map the cts system to its discrete time equivalent and check the NTF
(sys_d, Gp) = mapCtoD(sys_c, t=tdac)
ABCD = [sys_d.A sys_d.B; sys_d.C sys_d.D]
(NTF, G) = calculateTF(ABCD)
NTF = cancelPZ(NTF)
(plot_NTF, plot_fresp) = plotcoll_NTF.plotlist
plotPZ!(plot_NTF, NTF, color=:cyan)
#Also plot the STF
L0 = _zpk(sys_c[1,1])
@warn("FIXME: Again, numerical error in ss->zpk gives 0 gain")
f = range(0, stop=0.5, length=10)
G = evalTFP(L0, NTF, f)
Gw = waveform(f, dbv.(G))
push!(plot_fresp,
	cons(:wfrm, Gw, line=set(style=:solid, color=:magenta, width=2), label=""),
)
plotcoll_NTF.title = "NTF and STF"
displaygui(plotcoll_NTF)

@info("Performing dynamic range scaling..."); flush(stdout); flush(stderr)
#-------------------------------------------------------------------------------
#!!! This code assumes that the scale factors for the DT equivalent apply
#!!! to the CT system. A system with an RZ DAC will have inter-sample peaks
#!!! that exceed the values an the sampling instants.
(ABCDs, umax, S) = scaleABCD(ABCD, nlev=nlev, f=dsm.f0, xlim=1, umax=umax, N0=10^4)
S = S[1:dsm.order,1:dsm.order] #Don't worry about the extra states used in the d-t model
println("\tdone.")
println("\nScaled ABCD matrix (ABCDs):")
@show umax
display(ABCDs)

#Compute ABCDcs:
Sinv = inv(S)
Acs=S*Ac*Sinv; Bcs=S*Bc;  Ccs=Cc*Sinv
ABCDcs = [Acs Bcs; Ccs Dc]
sys_cs = _ss(Acs, Bcs, Ccs, Dc)
#ABCDcs needs to be checked with CT simulations to
# 1. Verify pulse response
# 2. Verify signal swings
println("\nABCDcs matrix:")
display(ABCDcs)

adc = Dict{Symbol, Any}(
	:order => dsm.order,
	:OSR => dsm.OSR,
	:opt => dsm.opt,
	:M => dsm.M,
	:f0 => dsm.f0,
	:NTF => NTF,
	:ABCD => ABCD,
	:umax => umax,
	:peak_snr => snrinfo.peak[2],
	:form => dsm.form,
	:ABCDc => ABCDc,
	:L0  => L0,
	:sys_c => sys_c,
	:ABCDcs => ABCDcs,
	:sys_cs => sys_cs,
)

println()
display(adc)

:END_OF_EXAMPLE


text:
"""
    couple_models(aero::LiftingLine, stru::GEBT)

Create a coupled model using a lifting line aerodynamic model and a
geometrically exact beam theory model.  This model introduces additional
parameters corresponding to the freestream air density ``\\rho_\\infty``,
followed by the external loads ``F_{x,i}, F_{y,i}, F_{z,i}, M_{x,i}, M_{y,i},
M_{z,i}`` or displacements ``u_{x,i}, u_{y,i}, u_{z,i}, \\theta_{x,i},
\\theta_{y,i}, \\theta_{z,i}`` for each node, followed by the constant
distributed loads ``f_{x,i}, f_{y,i}, f_{z,i}, m_{x,i}, m_{y,i}, m_{z,i}``
applied on each beam element (excluding aerodynamic loads), followed by the
body frame linear and angular velocities ``u, v, w, p, q, r``.

**NOTE: When using this model, the local frame for each beam element should be
oriented with the x-axis along the beam's axis, the y-axis forward, and the
z-axis normal to the surface**
"""
couple_models(aero::LiftingLine, stru::GEBT) = (aero, stru)

# --- traits --- #

function number_of_additional_parameters(aero::LiftingLine, stru::GEBT)
    return 1 + 6*length(stru.icol_point) + 6*length(stru.icol_elem) + 6
end

function coupling_inplaceness(::Type{LiftingLine{NA,TA}}, ::Type{<:GEBT}) where {NA,TA}
    return InPlace()
end

function coupling_rate_jacobian_type(::Type{LiftingLine{NA,TA}}, ::Type{<:GEBT}) where {NA,TA}
    model_types = TA.parameters
    if all(isempty.(coupling_rate_jacobian_type.(model_types, Ref(LiftingLineSection))))
        return Empty()
    elseif all(iszero.(coupling_rate_jacobian_type.(model_types, Ref(LiftingLineSection))))
        return Zeros()
    elseif all(isidentity.(coupling_rate_jacobian_type.(model_types, Ref(LiftingLineSection))))
        return Identity()
    elseif all(isinvariant.(coupling_rate_jacobian_type.(model_types, Ref(LiftingLineSection))))
        return Invariant()
    elseif all(isconstant.(coupling_rate_jacobian_type.(model_types, Ref(LiftingLineSection))))
        return Constant()
    elseif all(islinear.(coupling_rate_jacobian_type.(model_types, Ref(LiftingLineSection))))
        return Linear()
    else
        return Nonlinear()
    end
end

function coupling_state_jacobian_type(::Type{LiftingLine{NA,TA}}, ::Type{<:GEBT}) where {NA,TA}
    model_types = TA.parameters
    if all(isempty.(coupling_state_jacobian_type.(model_types, Ref(LiftingLineSection))))
        return Empty()
    elseif all(iszero.(coupling_state_jacobian_type.(model_types, Ref(LiftingLineSection))))
        return Zeros()
    elseif all(isidentity.(coupling_state_jacobian_type.(model_types, Ref(LiftingLineSection))))
        return Identity()
    elseif all(isinvariant.(coupling_state_jacobian_type.(model_types, Ref(LiftingLineSection))))
        return Invariant()
    elseif all(isconstant.(coupling_state_jacobian_type.(model_types, Ref(LiftingLineSection))))
        return Constant()
    elseif all(islinear.(coupling_state_jacobian_type.(model_types, Ref(LiftingLineSection))))
        return Linear()
    else
        return Nonlinear()
    end
end

coupling_parameter_jacobian_type(::Type{<:LiftingLine}, ::Type{<:GEBT}) = Nonlinear()

function coupling_time_gradient_type(::Type{LiftingLine{NA,TA}}, ::Type{<:GEBT}) where {NA,TA}
    model_types = TA.parameters
    if all(isempty.(coupling_time_gradient_type.(model_types, Ref(LiftingLineSection))))
        return Empty()
    elseif all(iszero.(coupling_time_gradient_type.(model_types, Ref(LiftingLineSection))))
        return Zeros()
    elseif all(isinvariant.(coupling_time_gradient_type.(model_types, Ref(LiftingLineSection))))
        return Invariant()
    elseif all(isconstant.(coupling_time_gradient_type.(model_types, Ref(LiftingLineSection))))
        return Constant()
    elseif all(islinear.(coupling_time_gradient_type.(model_types, Ref(LiftingLineSection))))
        return Linear()
    else
        return Nonlinear()
    end
end

# --- methods --- #

function get_coupling_inputs!(y, aero::LiftingLine{NA,TA}, stru::GEBT, dx, x, p, t) where {NA,TA}

    # extract number of state variables, inputs, and parameters
    nua = number_of_states(aero) # number of aerodynamic states
    nus = number_of_states(stru) # number of structural states
    nya = number_of_inputs(aero) # number of aerodynamic inputs
    nys = number_of_inputs(stru) # number of structural inputs
    npa = number_of_parameters(aero) # number of aerodynamic parameters
    nps = number_of_parameters(stru) # number of structural parameters
    npadd = number_of_additional_parameters(aero, stru) # number of additional parameters

    # get indices for state variables, inputs, and parameters
    iua = 1:nua # indices of aerodynamic states
    iya = 1:nya # indices of aerodynamic inputs
    ipa = 1:npa # indices of aerodynamic parameters
    ius = nua + 1 : nua + nus # indices of structural states
    iys = nya + 1 : nya + nys # indices of structural inputs
    ips = npa + 1 : npa + nps # indices of structural parameters
    ipadd = npa + nps + 1 : npa + nps + npadd # indices of additional parameters
    iuas = state_indices(aero.models) # indices of aerodynamic states for each section
    iyas = input_indices(aero.models) # indices of aerodynamic inputs for each section
    ipas = parameter_indices(aero.models) # indices of aerodynamic inputs for each section

    # separate state variables, inputs, and parameters
    dua = view(dx, iua) # aerodynamic rates
    ua = view(x, iua) # aerodynamic states
    ya = view(y, iya) # aerodynamic inputs
    pa = view(p, ipa) # aerodynamic parameters
    dus = view(dx, ius) # structural rates
    us = view(x, ius) # structural states
    ys = view(y, iys) # structural inputs
    ps = view(p, ips) # structural parameters
    padd = view(p, ipadd) # additional parameters for coupled model
    duas = view.(Ref(dua), iuas) # aerodynamic rates for each section
    uas = view.(Ref(ua), iuas) # aerodynamic states for each section
    yas = view.(Ref(ya), iyas) # aerodynamic inputs for each section
    pas = view.(Ref(pa), ipas) # aerodynamic parameters for each section

    # number of points and elements
    npoint = length(stru.icol_point)
    nelem = length(stru.icol_elem)

    # global parameters
    ρ = padd[1]
    ur = padd[1 + 6*npoint + 6*nelem + 1]
    vr = padd[1 + 6*npoint + 6*nelem + 2]
    wr = padd[1 + 6*npoint + 6*nelem + 3]
    pr = padd[1 + 6*npoint + 6*nelem + 4]
    qr = padd[1 + 6*npoint + 6*nelem + 5]
    rr = padd[1 + 6*npoint + 6*nelem + 6]

    # body linear and angular velocity
    V = SVector(ur, vr, wr)
    Ω = SVector(pr, qr, rr)

    # construct assembly from structural parameters
    assembly = gxbeam_assembly(ps, npoint, nelem, stru.start, stru.stop)

    # set prescribed point loads/displacements
    for ip = 1:npoint
        yoff = 6*(ip-1)
        poff = 1 + 6*(ip-1)
        ys[yoff+1] = padd[poff+1]
        ys[yoff+2] = padd[poff+2]
        ys[yoff+3] = padd[poff+3]
        ys[yoff+4] = padd[poff+4]
        ys[yoff+5] = padd[poff+5]
        ys[yoff+6] = padd[poff+6]
    end

    # set aerodynamic inputs and distributed loads
    for i = 1:NA
        # models for this section
        section_aero = aero.models[i]
        section_stru = LiftingLineSection()
        section_models = (section_aero, section_stru)

        # model dimensions for this section
        Nuai = number_of_states(section_aero)
        Nyai = number_of_inputs(section_aero)
        Npai = number_of_parameters(section_aero)
        Nusi = number_of_states(section_stru)
        Nysi = number_of_inputs(section_stru)
        Npsi = number_of_parameters(section_stru)

        # section properties
        icol = stru.icol_elem[i]
        element = assembly.elements[i]
        u_elem = SVector(us[icol], us[icol+1], us[icol+2]) # linear displacement
        θ_elem = SVector(us[icol+3], us[icol+4], us[icol+5]) # angular displacement
        F_elem = SVector(us[icol+6], us[icol+7], us[icol+8]) .* stru.force_scaling # internal force
        M_elem = SVector(us[icol+9], us[icol+10], us[icol+11]) .* stru.force_scaling # internal moment
        P_elem = SVector(us[icol+12], us[icol+13], us[icol+14]) .* stru.mass_scaling # linear momentum
        H_elem = SVector(us[icol+15], us[icol+16], us[icol+17]) .* stru.mass_scaling # angular momentum

        du_elem = SVector(dus[icol], dus[icol+1], dus[icol+2]) # linear displacement
        dθ_elem = SVector(dus[icol+3], dus[icol+4], dus[icol+5]) # angular displacement
        dF_elem = SVector(dus[icol+6], dus[icol+7], dus[icol+8]) .* stru.force_scaling # internal force
        dM_elem = SVector(dus[icol+9], dus[icol+10], dus[icol+11]) .* stru.force_scaling # internal moment
        dP_elem = SVector(dus[icol+12], dus[icol+13], dus[icol+14]) .* stru.mass_scaling # linear momentum
        dH_elem = SVector(dus[icol+15], dus[icol+16], dus[icol+17]) .* stru.mass_scaling # angular momentum

        # convert rotation parameter to Wiener-Milenkovic parameters
        scaling = GXBeam.rotation_parameter_scaling(θ_elem)
        θ_elem *= scaling # angular displacement (Wiener-Milenkovic parameters)
        dθ_elem *= scaling # angular displacement rate (Wiener-Milenkovic parameters)

        # local to body transformation
        CtCab = GXBeam.get_C(θ_elem)'*element.Cab # local to body transformation

        # local structural to aerodynamic transformation
        R = @SMatrix [0 -1 0; 1 0 0; 0 0 1] # local structural to aerodynamic transformation

        # local freestream linear and angular velocities
        vi = -R*(CtCab'*V + GXBeam.element_linear_velocity(element, P_elem, H_elem)) # local linear freestream velocity
        ωi = R*GXBeam.element_angular_velocity(element, P_elem, H_elem) # local angular freestream velocity

        # local freestream linear and angular accelerations
        dvi = -R*GXBeam.element_linear_velocity(element, dP_elem, dH_elem)
        dωi = R*GXBeam.element_angular_velocity(element, dP_elem, dH_elem)

        # NOTE: We assume that the local beam frame y-axis is oriented towards
        # the leading edge and the z-axis is oriented up
        # NOTE: We assume the local aerodynamic frame is oriented with the x-axis
        # in the chordwise direction and the z-axis in the (airfoil) normal
        # direction.

        # section state variable rates
        duai = SVector{Nuai}(duas[i]) # aerodynamic state variables
        dusi = vcat(dvi, dωi) # structural state variables
        dui = vcat(duai, dusi)

        # section state variables
        uai = SVector{Nuai}(uas[i]) # aerodynamic state variables
        usi = vcat(vi, ωi) # structural state variables
        ui = vcat(uai, usi)

        # section parameters
        pai = SVector{Npai}(pas[i]) # aerodynamic parameters
        psi = SVector(ρ) # structural parameters
        pi = vcat(pai, psi)

        # section inputs
        yi = get_coupling_inputs(section_models, dui, ui, pi, t)

        # separate inputs
        yai = view(yi, 1:Nyai)
        ysi = view(yi, Nyai + 1 : Nyai + Nysi)

        # save local inputs
        y[iya[iyas[i]]] .= yai # aerodynamic inputs

        # section aerodynamic loads (in body frame)
        fi = CtCab*R'*SVector(ysi[1], ysi[2], ysi[3])
        mi = CtCab*R'*SVector(ysi[4], ysi[5], ysi[6])

        # add constant distributed loads (in body frame)
        poff = 1 + 6*npoint + 6*(i-1)
        fi += SVector(padd[poff+1], padd[poff+2], padd[poff+3])
        mi += SVector(padd[poff+4], padd[poff+5], padd[poff+6])

        # save distributed loads for this element (in body frame)
        yoff = 6*npoint + 6*(i-1)
        y[iys[yoff+1 : yoff+3]] = fi
        y[iys[yoff+4 : yoff+6]] = mi
    end

    # save body frame linear/angular velocities
    yoff = 6*npoint + 6*nelem
    y[iys[yoff+1]] = ur
    y[iys[yoff+2]] = vr
    y[iys[yoff+3]] = wr
    y[iys[yoff+4]] = pr
    y[iys[yoff+5]] = qr
    y[iys[yoff+6]] = rr

    return y
end

# --- performance overloads --- #

# TODO

# --- convenience methods --- #

function set_additional_parameters!(padd, aero::LiftingLine, stru::GEBT;
    rho, point_conditions, element_loads, u, v, w, p, q, r) where {NA,TA}

    np = length(stru.icol_point)
    ne = length(stru.icol_elem)

    padd[1] = rho

    for ip = 1:np
        padd[1+6*(ip-1)+1] = point_conditions[1,ip]
        padd[1+6*(ip-1)+2] = point_conditions[2,ip]
        padd[1+6*(ip-1)+3] = point_conditions[3,ip]
        padd[1+6*(ip-1)+4] = point_conditions[4,ip]
        padd[1+6*(ip-1)+5] = point_conditions[5,ip]
        padd[1+6*(ip-1)+6] = point_conditions[6,ip]
    end

    for ie = 1:ne
        padd[1+6*np+6*(ie-1)+1] = element_loads[1,ie]
        padd[1+6*np+6*(ie-1)+2] = element_loads[2,ie]
        padd[1+6*np+6*(ie-1)+3] = element_loads[3,ie]
        padd[1+6*np+6*(ie-1)+4] = element_loads[4,ie]
        padd[1+6*np+6*(ie-1)+5] = element_loads[5,ie]
        padd[1+6*np+6*(ie-1)+6] = element_loads[6,ie]
    end

    padd[1+6*np+6*ne+1] = u
    padd[1+6*np+6*ne+2] = v
    padd[1+6*np+6*ne+3] = w
    padd[1+6*np+6*ne+4] = p
    padd[1+6*np+6*ne+5] = q
    padd[1+6*np+6*ne+6] = r

    return padd
end

function separate_additional_parameters(aero::LiftingLine, stru::GEBT, padd)

    np = length(stru.icol_point)
    ne = length(stru.icol_elem)

    rho = padd[1]

    point_conditions = reshape(view(padd, 1 + 1 : 1 + 6*np), 6, np)

    element_loads = reshape(view(padd, 1 + 6*np + 1 : 1 + 6*np + 6*ne), 6, ne)

    u = padd[1 + 6*np + 6*ne + 1]
    v = padd[1 + 6*np + 6*ne + 2]
    w = padd[1 + 6*np + 6*ne + 3]
    p = padd[1 + 6*np + 6*ne + 4]
    q = padd[1 + 6*np + 6*ne + 5]
    r = padd[1 + 6*np + 6*ne + 6]

    return (rho = rho, point_conditions = point_conditions, element_loads = element_loads,
        u = u, v = v, w = w, p = p, q = q, r = r)
end


text:
# blind deconvolution
__precompile__()

module Decon
import JuMIT.DSP
import JuMIT.Misfits
using Optim, LineSearches

type Param
	gf::Array{Float64,2}
	dgf::Array{Float64,2}
	ntgf::Int64
	dobs::Array{Float64,2}
	dcal::Array{Float64,2}
	ddcal::Array{Float64,2}
	wav::Vector{Float64}
	wavmat::Array{Float64,2}
	dwavmat::Array{Float64,2}
	spow2::Array{Complex{Float64}}
	rpow2::Array{Complex{Float64}}
	wpow2::Array{Complex{Float64}}
	nt::Int64
	nr::Int64
	attrib_inv::Symbol
	verbose::Bool
	fftplan
	ifftplan
end

function Param(ntgf, nt, nr, gfobs, wavobs; verbose=false, attrib_inv=:gf,)
	dobs=zeros(nt, nr)
	dcal=zeros(nt, nr)
	ddcal=zeros(nt, nr)

	# initial values are random
	wav=randn(nt)
	wavmat=repeat(wav,outer=(1,nr))
	dwavmat=similar(wavmat)
	gf=randn(ntgf,nr)
	dgf=similar(gf)

	# use maximum threads for fft
	#FFTW.set_num_threads(Sys.CPU_CORES)
	# fft dimension for plan
	np2=nextpow2(maximum([2*nt, 2*ntgf]))
	# create plan
	fftplan = plan_fft!(complex.(zeros(np2,nr)),[1],flags=FFTW.PATIENT)
	ifftplan = plan_ifft!(complex.(zeros(np2,nr)),[1],flags=FFTW.PATIENT)

	# preallocate pow2 vectors
	spow2=complex.(zeros(np2,nr))
	rpow2=complex.(zeros(np2,nr))
	wpow2=complex.(zeros(np2,nr))

	pa = Param(gf, dgf, ntgf, dobs, dcal, ddcal, wav, wavmat, dwavmat,
	     spow2, rpow2, wpow2,
	     nt, nr, attrib_inv, verbose, fftplan, ifftplan)

	F!(pa, mattrib=:gfwav, dattrib=:dobs, gf=gfobs, wav=wavobs)

	return pa
end
function ninv(pa)
	if(pa.attrib_inv == :wav)
		return pa.nt
	else(pa.attrib_inv == :gf)
		return pa.ntgf*pa.nr
	end
end

function model_to_x!(x, pa)
	if(pa.attrib_inv == :wav)
		copy!(x, pa.wav)
	else(pa.attrib_inv == :gf)
		copy!(x, pa.gf)
	end
	return x
end

function x_to_model!(x, pa)
	if(pa.attrib_inv == :wav)
		copy!(pa.wav, x)
	else(pa.attrib_inv == :gf)
		copy!(pa.gf, x)
	end
	return pa
end

function F!(
			pa::Param,
			x::Vector{Float64}=zeros(ninv(pa)),
			last_x::Vector{Float64}=ones(ninv(pa));
			gf=nothing,
			wav=nothing,
			mattrib::Symbol=:x,
			dattrib::Symbol=:dcal,
			   )

	if(x!=last_x)

		if(mattrib == :x)
			x_to_model!(x, pa)
			gf=pa.gf
			wav=pa.wav
		elseif(mattrib == :gf)
			wav=pa.wav
			(gf===nothing) && error("need gf")
		elseif(mattrib == :wav)
			gf=pa.gf
			(wav===nothing) && error("need wav")
		elseif(mattrib == :gfwav)
			(gf===nothing) && error("need gf")
			(wav===nothing) && error("need wav")
		else
			error("invalid mattrib")
		end

		#pa.verbose && println("updating buffer")
		copy!(last_x, x)
		for ir in 1:pa.nr
			pa.wavmat[:,ir]=wav
		end
		DSP.fast_filt!(pa.ddcal, gf, pa.wavmat, :s, nwplags=pa.nt-1,
		 nwnlags=0,nsplags=pa.nt-1,nsnlags=0,nrplags=pa.ntgf-1,nrnlags=0,
		  np2=size(pa.fftplan,1), fftplan=pa.fftplan, ifftplan=pa.ifftplan,
		  spow2=pa.spow2, rpow2=pa.rpow2, wpow2=pa.wpow2,
		  			)
		copy!(getfield(pa, dattrib), pa.ddcal)
		return pa
	end
end

function func_grad!(storage, x::Vector{Float64},
			last_x::Vector{Float64},
			pa)

	#pa.verbose && println("computing gradient...")

	# x to w
	x_to_model!(x, pa)

	F!(pa, x, last_x)

	if(storage === nothing)
		# compute misfit and δdcal
		f = Misfits.fg_cls!(nothing, pa.dcal, pa.dobs)
		return f
	else
		f = Misfits.fg_cls!(pa.ddcal, pa.dcal, pa.dobs)
		Fadj!(pa, storage, pa.ddcal)
		return f
	end
end

"""
Apply Fadj to
"""
function Fadj!(pa, storage, dcal)
	if(pa.attrib_inv == :wav)
		DSP.fast_filt!(dcal, pa.gf, pa.wavmat, :w, nwplags=pa.nt-1,
		 nwnlags=0,nsplags=pa.nt-1,nsnlags=0,nrplags=pa.ntgf-1,nrnlags=0,
		  np2=size(pa.fftplan,1), fftplan=pa.fftplan, ifftplan=pa.ifftplan,
		  spow2=pa.spow2, rpow2=pa.rpow2, wpow2=pa.wpow2,
			)
		copy!(storage, sum(pa.wavmat,2))
	else(pa.attrib_inv == :gf)
		for ir in 1:pa.nr
			pa.wavmat[:,ir]=pa.wav
		end
		DSP.fast_filt!(dcal, pa.dgf, pa.wavmat, :r, nwplags=pa.nt-1,
		 nwnlags=0,nsplags=pa.nt-1,nsnlags=0,nrplags=pa.ntgf-1,nrnlags=0,
		  np2=size(pa.fftplan,1), fftplan=pa.fftplan, ifftplan=pa.ifftplan,
		  spow2=pa.spow2, rpow2=pa.rpow2, wpow2=pa.wpow2,
			)
		copy!(storage, pa.dgf)
	end
	return storage
end

# core algorithm
function update!(pa::Param, x, last_x, df; store_trace::Bool=true, extended_trace::Bool=false,
	     time_limit=Float64=2.0*60.,
	     f_tol::Float64=1e-8, g_tol::Float64=1e-8, x_tol::Float64=1e-8)

	# initial w to x
	model_to_x!(x, pa)

	"""
	Unbounded LBFGS inversion, only for testing
	"""
	res = optimize(df, x,
		       LBFGS(),
		       Optim.Options(g_tol = g_tol, f_tol=f_tol, x_tol=x_tol,
		       iterations = 200, store_trace = store_trace,
		       extended_trace=extended_trace, show_trace = false))
	pa.verbose && println(res)

	x_to_model!(Optim.minimizer(res), pa)

	return res
end

function update_all!(pa)

	# convert initial model to the inversion variable
	pa.attrib_inv=:gf
	xgf = zeros(ninv(pa));
	last_xgf = randn(size(xgf)) # reset last_x

	# create df for optimization
	dfgf = OnceDifferentiable(x -> func_grad!(nothing, x, last_xgf, pa),
	  (storage, x) -> func_grad!(storage, x, last_xgf, pa))

	pa.attrib_inv=:wav
	xwav = zeros(ninv(pa));
	last_xwav = randn(size(xwav)) # reset last_x

	# create df for optimization
	dfwav = OnceDifferentiable(x -> func_grad!(nothing, x, last_xwav, pa),
	  (storage, x) -> func_grad!(storage, x, last_xwav, pa))

	converged_all=false
	itr_all=0
	maxitr_all=10
	while !converged_all && itr_all < maxitr_all
		itr_all += 1
		pa.verbose && (itr_all > 1) && println("\t", "failed to converge.. reintializing round trips")
		# starting models
		randn!(pa.wav)
		randn!(pa.gf)

		fgfmin=Inf
		fgfmax=0.0
		fwavmin=Inf
		fwavmax=0.0

		maxitr=1000
		itr=0
		converged=false
		tol=1e-8
		while !converged && itr < maxitr
			pa.verbose && (itr > 1) && println("\t", "round trip: ", itr)
			itr += 1
			pa.attrib_inv=:gf
			resgf = update!(pa, xgf, last_xgf, dfgf)
			fgf = Optim.minimum(resgf)
			fgfmin = min(fgfmin, fgf)
			fgfmax = max(fgfmax, fgf)
			pa.attrib_inv=:wav
			reswav = update!(pa, xwav, last_xwav, dfwav)
			fwav = Optim.minimum(reswav)
			fwavmin = min(fwavmin, fwav)
			fwavmax = max(fwavmax, fwav)

			converged = ((abs(fwav-fgf)) < tol)
		end
		converged_all = (fgfmin/fgfmax < 1e-3) & (fwavmin/fwavmax < 1e-3)
	end
end

# initialize gf and ssf

# create unknown vector

end # module


text:
function Phi_s(CellData,s,z,Def)
    """
    Solid Potential Transfer Function
    # Add Ins and Outs
        # Cell Data
        # Frequency Vector
        # Discretisation Locations
        # Electrode Definition
    """

 if Def == "Pos"
    Electrode = CellData.Pos #Electrode Length
 else
    Electrode = CellData.Neg #Electrode Length
 end

as = 3*Electrode.ϵ_s/Electrode.Rs # Specific interfacial surf. area
κ_eff = CellData.Const.κ*Electrode.ϵ_e^Electrode.κ_brug #Effective Electrolyte Conductivity
σ_eff = Electrode.σ*Electrode.ϵ_s^Electrode.σ_brug #Effective Electrode Conductivity
comb_cond_eff = κ_eff+σ_eff #combining into single variable

#Defining SOC
θ = CellData.Const.SOC * (Electrode.θ_100-Electrode.θ_0) + Electrode.θ_0

#Beta's
β = @. Electrode.Rs*sqrt(s/Electrode.Ds)

#Prepare for j0
cs0 = Electrode.cs_max * θ

#Current Flux Density
if CellData.Const.CellTyp == "Doyle_94"
   κ = Electrode.k_norm/Electrode.cs_max/CellData.Const.ce0^(1-Electrode.α)
   j0 = κ*(CellData.Const.ce0*(Electrode.cs_max-cs0))^(1-Electrode.α)*cs0^Electrode.α
else
   j0 = Electrode.k_norm*(CellData.Const.ce0*(cs0/Electrode.cs_max*(1-cs0/Electrode.cs_max)))^(1-Electrode.α)
end

#Resistance
Rtot = R*CellData.Const.T/(j0*F^2) + Electrode.RFilm

#∂Uocp_Def
∂Uocp_elc = CellData.Const.∂Uocp(Def,θ)/Electrode.cs_max

ν = @. Electrode.L*sqrt((as/σ_eff+as/κ_eff)/(Rtot+∂Uocp_elc*(Electrode.Rs/(F*Electrode.Ds))*(tanh(β)/(tanh(β)-β)))) #Condensing Variable - eq. 4.13
ν_∞ = @. Electrode.L*sqrt((as/κ_eff+as/σ_eff)/(Rtot))

ϕ_tf = @. (-Electrode.L*(κ_eff*(cosh(ν)-cosh((z-1)*ν)))-Electrode.L*(σ_eff*(1-cosh(z*ν)+z*ν*sinh(ν))))/(CellData.Const.CC_A*σ_eff*(comb_cond_eff)*ν*sinh(ν)) #Transfer Function - eq. 4.19
D = @. (-Electrode.L*(κ_eff*(cosh(ν_∞)-cosh(z-1)*ν_∞))-Electrode.L*(σ_eff*(1-cosh(z*ν_∞)+z*ν_∞*sinh(ν_∞))))/(CellData.Const.CC_A*σ_eff*(comb_cond_eff)*ν_∞*sinh(ν_∞)) # Contribution to D as G->∞
D_term = "@. -$(Electrode.L)*($κ_eff*(cosh($ν_∞)-cosh($z-1)*$ν_∞))/($(CellData.Const.CC_A)*$σ_eff*($comb_cond_eff)*$ν_∞*sinh($ν_∞))-$(Electrode.L)*($σ_eff*(1-cosh($z*$ν_∞)+$z*$ν_∞*sinh($ν_∞)))/($(CellData.Const.CC_A)*$σ_eff*($comb_cond_eff)*$ν_∞*sinh($ν_∞))"
zero_tf = @. Electrode.L*(z-2)*z/(2*CellData.Const.CC_A*σ_eff)
ϕ_tf[:,findall(s.==0)] .= zero_tf[:,findall(s.==0)]
res0 = zeros(length(z))

if Def == "Pos"
   ϕ_tf = -ϕ_tf
   D = -D
   D_term = "@. $(Electrode.L)*($κ_eff*(cosh($ν_∞)-cosh($z-1)*$ν_∞))/($(CellData.Const.CC_A)*$σ_eff*($comb_cond_eff)*$ν_∞*sinh($ν_∞))-$(Electrode.L)*($σ_eff*(1-cosh($z*$ν_∞)+$z*$ν_∞*sinh($ν_∞)))/($(CellData.Const.CC_A)*$σ_eff*($comb_cond_eff)*$ν_∞*sinh($ν_∞))"
   if Debug == 1
      println("D:Phi_s:Pos:",D)
      println("zero_tf:Phi_s:Pos:",zero_tf)
      println("zero_tf:Phi_s:Pos:",ϕ_tf[1:5])
      println("z:Phi_s:Pos:",z)
      println("ν_∞:Phi_s:Pos:",ν_∞)
      println("as:Phi_s:Pos:",as)
      println("L:Phi_s:Pos:",Electrode.L)
      println("ν:Phi_s:Pos:",ν[1:5])
      println("θ:Phi_s:Pos:",θ)
      println("κ:Phi_s:Pos:",κ)
      println("Rtot:Phi_s:Pos:",Rtot)
      println("j0:Phi_s:Pos:",j0)
      println("κ_eff:Phi_s:Pos:",κ_eff)
      println("σ_eff:Phi_s:Pos:",σ_eff)
   end
else
   if Debug == 1
      println("D:Phi_s:Neg:",D)
      println("D_check:Phi_s:Neg:",D_check)
      println("ν_∞:Phi_s:Neg:",ν_∞)
      println("θ:Phi_s:Neg:",θ)
      println("κ:Phi_s:Neg:",κ)
      println("Rtot:Phi_s:Neg:",Rtot)
      println("j0:Phi_s:Neg:",j0)
      println("κ_eff:Phi_s:Neg:",κ_eff)
      println("σ_eff:Phi_s:Neg:",σ_eff)
      println("zero_tf:Phi_s:Neg:",zero_tf)
      #println("z:Phi_s:Neg:",z)
      #println("ν:Phi_s:Neg:",ν)
   end
end
return ϕ_tf, D, res0, D_term

end

max_stars_repo_name: BradyPlanden/RealisationAlgorithms

max_issues_repo_name: BradyPlanden/RealisationAlgorithms

max_forks_repo_name: BradyPlanden/RealisationAlgorithms


text:
import RowEchelon

#=
Homework 3
Section 1.7
Problem 41
=#

import RowEchelon

#=
Homework 3
Section 1.7
Problem 5
=#

# Create the matrix
A = [3 -4 10 7 -4 ;
	 -5 -3 -7 -11 15 ;
	 4 3 5 2 1 ;
	 8 -7 23 4 15]

# Row Reduce the matrix
A_rref = RowEchelon.rref(A)

# Store the rows and columns of A
rows = size(A)[1]
cols = size(A)[2]

# Find the pivot columns
pivotCols = zeros(0)
for i = 1:rows
	for j = 1:cols
		if A_rref[i,j] != 0
			append!(pivotCols,j)
			break
		end
	end
end

# The columns we want are the pivot columns
B = Array{Float64}(undef, 4, 0)
for col in pivotCols
	global B
	B = hcat(B, A[:,Int64(col)])
end

# Row Reduce
B_rref = RowEchelon.rref(B)

# Print results:
println("Original Matrix:")
println(A)
print("Matrix made from columns ")
print(pivotCols)
println(":")
println(B)
println("Row Reduced:")
println(B_rref)