text:
using NKFtool

folder = "C:\\Playlists\\" # 変換元M3Uファイルを置いてあるフォルダパス
list = readdir(folder)
for file in list
    if occursin("m3u", file)
        # M3Uファイルの読込元と出力先を開く
        filein = folder * file
        fileout = folder*"UTF-8\\"*file #
        fw = open(fileout, "w")
        fr = open(filein, "r")

        # fileinの文字列のエンコードをUTF-8に変換
        lines = nkf_convert(fr)
        # フォルダパスをクラウドストレージに合わせる
        lines = replace(lines, "C:\\MUSIC"=>"/Ubuntu/MUSIC")
        # 区切り文字をスラッシュに変換する
        lines = replace(lines, "\\"=>"/")
        # コメント行と"\r"だけの行を除いて書き込み
        for line in split(lines, "\n")
            if !occursin("#", line) && !(line=="\r")
                println(fw, line)
            end
        end
        close(fr)
        close(fw)
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ b787a8ca-7bd4-40a7-8574-334c72eda77d
begin
	using PrettyTables
	using PlutoUI
end

# ╔═╡ 67909d8a-421f-41cb-bc87-c1488be36074
md"""
### Practice Problem: Gas Phase Chemical Reaction Equilibrium Calculation

The gas-phase reaction:

$$A_{(g)}+B_{(g)} ~{\rightleftarrows}~ 2C_{(g)}$$

occurs in a closed well-mixed reactor with constant volume V = 10~L, a constant T = 500$^{\circ}$C, but with an unknown pressure (measured in bar). Initailly, there is 2 mol of A and 1 mol of B and no product C. Let the index of A = 1, B = 2 and C = 3. The Gibbs energies of formation at standard temperature and pressure are given by: $G_{1}^{\circ}$ = -2.0 kJ/mol,  $G_{2}^{\circ}$ = -3.6 kJ/mol and $G_{3}^{\circ}$ = -6.7 kJ/mol, respectively.

The teaching team measured the equilibrium extent of reaction $\epsilon^{\star}_{1}$ = 0.67 mol for these conditions. However, the _mythical mole meter_ used by the teaching team to measure this extent has been acting up (and may be giving the wrong reading).

__Assume__:
* The gas phase is _ideal_ for all components
* The ideal gas law is valid for this system
* The reference pressure $P^{\circ}$ = 1 bar.
* R = 8.314 J mol$^{-1}$ K$^{-1}$ and R = 0.08314 bar L mol$^{-1}$ K$^{-1}$.

__Question__:
* Is the measured equilibrium extent of reaction consistent with the problem setup?

"""

# ╔═╡ 80b572d7-3400-4373-af66-a394fe45296f
begin

	# setup the problem -
	V = 10.0 			# units: L
	T = (500+273.15)	# units: K
	R = 0.008314 		# units: kJ mol^-1 K^-1
	R_bar = 0.08314 	# units: bar L mol^-1 K^-1
	Pₒ = 1.0 			# units: bar

	# G_formation_array -
	G_formation_array = [
		-2.0 ; # 1 G_A units: kJ/mol
		-3.6 ; # 2 G_B units: kJ/mol
		-6.7 ; # 3 G_C units: kJ/mol
	];

	# reaction stoichiometry -
	S = [
		-1.0 ; 	# 1 A
		-1.0 ; 	# 2 B
		2.0 ; 	# 3 C
	];

	# initial mol -
	n_initial = [
		2.0 ; # 1 A units: mol
		1.0 ; # 2 B units: mol
		0.0 ; # 3 C units: mol
	];

	# what is the measured extent?
	ϵ₁ = 0.67 # units: mol

	# show -
	nothing
end

# ╔═╡ 9c7e21c6-dfe3-4f3e-b2e0-23d7ab578619
md"""
#### Solution: Compute the Eqilibrium constant $K$ using two different approaches and compare

To determine if the _mythical mole meter_ is correct, we can compute the equilibrium constant using the extent $\epsilon_{1}$, and then compare that value with the equilibrium constant computed using the $\Delta{G}^{\circ}$. If the values are the same, the _mythical mole meter_ is working, otherwise the _mythical mole meter_ is broken (and giving bad measurements of the extent).

"""

# ╔═╡ 4a63eb66-7090-41d0-80b3-94b56f5edebf
md"""
###### Compute K$_{\Delta{G}}$: equilibrium constant using the ΔG values
"""

# ╔═╡ 4e10a9b5-cb5e-4b4a-87b4-0ce024a1ad6f
begin

	# compute the ΔG -
	ΔG = sum(G_formation_array.*S)
	RT = R*T

	# compute the equlibrium constant using the ΔG -
	K_ΔG = exp(-ΔG/RT)

	# display -
	with_terminal() do
		println("K_ΔG = $(K_ΔG)")
	end
end

# ╔═╡ 9f253794-d773-48bc-a1b5-94ca242ff798
md"""
###### Compute K$_{\epsilon}$: equilibrium constant using the extent values

Let's compute the equilibrium constant from the extent of reaction. Setup the state table using the species mole balance:

$$n_{i} = n^{\circ}_{i} + \sigma_{i1}\epsilon_{1}$$

"""

# ╔═╡ 47e684a6-773b-4009-b7cc-23dc2fbf4a4b
begin

	# compute the final number of mols, mol fraction -
	n_final = n_initial .+ S*ϵ₁
	n_total = sum(n_final)
	y_final = (1/n_total)*n_final

	# setup data array =
	state_table_array = Array{Any,2}(undef, 3,5)

	# A -
	state_table_array[1,1] = "A"
	state_table_array[1,2] = S[1]
	state_table_array[1,3] = n_initial[1]
	state_table_array[1,4] = n_final[1]
	state_table_array[1,5] = y_final[1]

	# B -
	state_table_array[2,1] = "B"
	state_table_array[2,2] = S[2]
	state_table_array[2,3] = n_initial[2]
	state_table_array[2,4] = n_final[2]
	state_table_array[2,5] = y_final[2]

	# C -
	state_table_array[3,1] = "C"
	state_table_array[3,2] = S[3]
	state_table_array[3,3] = n_initial[3]
	state_table_array[3,4] = n_final[3]
	state_table_array[3,5] = y_final[3]

	# setup a pretty table =
	# header row -
	state_table_header_row = (["Species","σ_i1","nᵢ_initial","nᵢ_final","yᵢ_final"],["","","mol", "mol", ""]);

	with_terminal() do
		# write the table -
		pretty_table(state_table_array; header=state_table_header_row)
	end

end

# ╔═╡ 8dca4509-6d6b-41c2-9d22-81c4563a5136
md"""
Now that we have the state table, we can use the gas phase equlibrium expression:

$$\prod_{i=1}^{\mathcal{M}}\left(y_{i}\hat{\phi}_{i}^{v}\right)^{\sigma_{i1}} = \left(\frac{P}{P^{\circ}}\right)^{-\sigma_{1}}K$$

To compute the equlibrium constant using the extent. Becuase we have assumed ideality, all the fugacity coefficient terms $\hat{\phi}^{v}_{i}\rightarrow\forall{i}$. In addition the pressure term is 1 becuase $\sigma_{1} = 0$ (the number of moles is conserved).
Thus, the equilibrium constant is given by:

$$K = \frac{y_{3}^2}{y_{1}y_{2}}$$
"""

# ╔═╡ b7c489c1-513b-4c1c-b22a-9c0738e4a62e
begin

	# compute the overall stoichiometric coefficient -
	σ₁ = sum(S);

	# Compute pressure term -
	P = (1/V)*(R_bar*T)*n_total
	P_term = (P/Pₒ)^(-σ₁)

	# compute the K -
	K_ϵ = (y_final[1]^S[1])*(y_final[2]^S[2])*(y_final[3]^S[3])*(1/P_term)

	with_terminal() do
		println("K_ϵ = $(K_ϵ)")
	end
end

# ╔═╡ 9ce4fb80-a6da-4a3d-88f2-83092be1f340
md"""
### Bummer. K$_{\Delta{G}}$ $\neq$ K$_{\epsilon}$ $\longrightarrow$ the _mythical mole meter_ is broken!
"""

# ╔═╡ ccc364fe-3d4a-11ec-3415-0753629989e9
html"""
<style>
main {
    max-width: 1200px;
    width: 85%;
    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]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"

[compat]
PlutoUI = "~0.7.18"
PrettyTables = "~1.2.3"
"""

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

[[AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "0ec322186e078db08ea3e7da5b8b2885c099b393"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.0"

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

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

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

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

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

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

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

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

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

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

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

[[IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.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"

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

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

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

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

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

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

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

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

[[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 = ["AbstractPlutoDingetjes", "Base64", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "57312c7ecad39566319ccf5aa717a20788eb8c1f"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.18"

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

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

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

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

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

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

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

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

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

[[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:
# ╟─67909d8a-421f-41cb-bc87-c1488be36074
# ╠═80b572d7-3400-4373-af66-a394fe45296f
# ╟─9c7e21c6-dfe3-4f3e-b2e0-23d7ab578619
# ╟─4a63eb66-7090-41d0-80b3-94b56f5edebf
# ╠═4e10a9b5-cb5e-4b4a-87b4-0ce024a1ad6f
# ╟─9f253794-d773-48bc-a1b5-94ca242ff798
# ╠═47e684a6-773b-4009-b7cc-23dc2fbf4a4b
# ╠═8dca4509-6d6b-41c2-9d22-81c4563a5136
# ╠═b7c489c1-513b-4c1c-b22a-9c0738e4a62e
# ╟─9ce4fb80-a6da-4a3d-88f2-83092be1f340
# ╠═b787a8ca-7bd4-40a7-8574-334c72eda77d
# ╟─ccc364fe-3d4a-11ec-3415-0753629989e9
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

ext: jl


text:
struct DTSA <: NeXLAlgorithm end

"""
   mac(zz::Int, ev::Float64, ::Type{DTSA})::Float64

Calculate the elemental MAC using Heinrich's IXCOM 11 formula as implemented by Myklebust in DTSA.
"""
function mac(elm::Element, ev::Float64, ::Type{DTSA})::Float64
    # Ref: Heinrich's formula as implemented by Myklebust translated into Julia
    # This expression only works for x-ray energies below the K-edge and
    # above the K-edge for Z < 50. Energies above the K-edge for elements
    # Z > 49 are completely nuts.
    zz = z(elm)
    if ev <= 10.0
        return 1.0e6
    end
    if (zz < 3) || (zz > 95)
        return 0.001
    end
    ee = collect(
        has(elm, sh) ? energy(AtomicSubShell(zz, sh)) : 0.0 for sh in allsubshells[1:10]
    )
    # formula  in eV units
    nm, cc, az, bias = 0.0, 0.0, 0.0, 0.0
    if ev > ee[1] # ev is above the K edge.
        if zz < 6
            cc = 0.001 * ((1.808599 * zz) - 0.287536)
            az = (((-14.15422 * zz) + 155.6055) * zz) + 24.4545
            bias = (18.2 * zz) - 103
            nm = (((-0.01273815 * zz) + 0.02652873) * zz) + 3.34745
        else
            cc =
                1.0E-5 * (
                    (
                        ((525.3 + (133.257 * zz)) - (7.5937 * zz * zz)) +
                        (0.169357 * zz * zz * zz)
                    ) - (0.0013975 * zz * zz * zz * zz)
                )
            az = ((((-0.152624 * zz) + 6.52) * zz) + 47) * zz
            nm = 3.112 - (0.0121 * zz)
            if (ev > ee[1]) && (zz >= 50)
                az = ((((-0.015 * zz) + 3.52) * zz) + 47) * zz
            end
            if (ev > ee[1]) && (zz >= 57)
                cc = 1.0E-6 * ((200.0 + (100.0 * zz)) - (zz * zz))
            end
        end
    elseif ev > ee[4]
        # ev is below K-edge & above L3-edge
        c =
            0.001 * (
                ((-0.0924 + (0.141478 * zz)) - (0.00524999 * zz * zz)) +
                (9.85296E-5 * zz * zz * zz)
            )
        c = (c - (9.07306E-10 * zz * zz * zz * zz)) + (3.19245E-12 * zz * zz * zz * zz * zz)
        cc = c
        az = ((((((-0.000116286 * zz) + 0.01253775) * zz) + 0.067429) * zz) + 17.8096) * zz
        nm = (((-4.982E-5 * zz) + 0.001889) * zz) + 2.7575
        if (ev < ee[2]) && (ev > ee[3])
            cc = c * 0.858
        end
        if ev < ee[3]
            cc = c * ((0.8933 - (0.00829 * zz)) + (6.38E-5 * zz * zz))
        end
    elseif (ev < ee[4]) && (ev > ee[5])
        # ev is below L3 and above M1
        nm = (((((4.4509E-6 * zz) - 0.00108246) * zz) + 0.084597) * zz) + 0.5385
        if zz < 30
            c =
                (
                    (
                        (((((0.072773258 * zz) - 11.641145) * zz) + 696.02789) * zz) -
                        18517.159
                    ) * zz
                ) + 188975.7
        else
            c =
                (
                    (
                        (((((0.001497763 * zz) - 0.40585911) * zz) + 40.424792) * zz) -
                        1736.63566
                    ) * zz
                ) + 30039
        end
        cc = 1.0E-7 * c
        az =
            (
                (((((-0.00018641019 * zz) + 0.0263199611) * zz) - 0.822863477) * zz) +
                10.2575657
            ) * zz
        if zz < 61
            bias =
                (
                    (((((-0.0001683474 * zz) + 0.018972278) * zz) - 0.536839169) * zz) +
                    5.654
                ) * zz
        else
            bias =
                (
                    (((((0.0031779619 * zz) - 0.699473097) * zz) + 51.114164) * zz) -
                    1232.4022
                ) * zz
        end
    elseif ev >= ee[9]
        az = (4.62 - (0.04 * zz)) * zz
        c = 1.0E-8 * ((((((-0.129086 * zz) + 22.09365) * zz) - 783.544) * zz) + 7770.8)
        c = c * ((((((4.865E-6 * zz) - 0.0006561) * zz) + 0.0162) * zz) + 1.406)
        cc = c * ((((-0.0001285 * zz) + 0.01955) * zz) + 0.584)
        bias = ((((0.000378 * zz) - 0.052) * zz) + 2.51) * ee[8]
        nm = 3 - (0.004 * zz)
        if (ev < ee[5]) && (ev >= ee[6])
            cc = c * ((0.001366 * zz) + 1.082)
        end
        if (ev < ee[7]) && (ev >= ee[8])
            cc = 0.95 * c
        end
        if (ev < ee[8]) && (ev >= ee[9])
            cc = 0.8 * c * ((((0.0005083 * zz) - 0.06) * zz) + 2.0553)
        end
    elseif ev < ee[9]
        cc = 1.08E-7 * ((((((-0.0669827 * zz) + 17.07073) * zz) - 1465.3) * zz) + 43156)
        az = ((((0.00539309 * zz) - 0.61239) * zz) + 19.64) * zz
        bias = (4.5 * zz) - 113.0
        nm = 0.3736 + (0.02401 * zz)
    end
    if ev > ee[10]
        mu = (cc * exp(nm * log(12397.0 / ev)) * zz * zz * zz * zz) / a(elm)
        return mu * (1.0 - exp((bias - ev) / az))
    else
        mu = (cc * exp(nm * log(12397.0 / ev)) * zz * zz * zz * zz) / a(elm)
        return (1.02 * mu * (ev - 10.0)) / (ee[10] - 10.0)
    end
end


text:
function MixtureModel(hmm::HMM)
    sdists = hmm.A[1, :]
    MixtureModel([hmm.B...], sdists)
end

function HMM(m::MixtureModel)
    K = ncomponents(m)
    a = probs(m)
    A = repeat(permutedims(m.prior.p), K, 1)
    B = m.components
    HMM(a, A, B)
end
#
# function PeriodicHMM(vec_mix::Vector{MixtureModel})
#     T = length(vec_mix)
#     K = ncomponents(vec_mix[1])
#     a = probs(vec_mix[1])
#     A = repeat.(permutedims.(probs.(vec_mix)), K, 1)
#     B = components.(vec_mix)
#     PeriodicHMM(a, [A[t][k,l] for k in 1:K, l in 1:K, t in 1:T], [B[t][k] for k in 1:K, t in 1:T])
# end

function copy(mix::MixtureModel)
    MixtureModel(components(mix), probs(mix))
end

function fit_mle(mix::MixtureModel, observations; init = :none, kwargs...)
    mix = copy(mix)

    if init == :kmeans
        kmeans_init!(mix, observations, display = get(kwargs, :display, :none))
    end
    fit_em(mix, observations; kwargs...)
end

function update_a!(
    a::AbstractVector,
    A::AbstractMatrix
)
    a[:] = A[1, :]
end

function update_A!(
    A::AbstractMatrix,
    γ::AbstractMatrix
)
    K = size(A, 1)
    A[:, :] = repeat(mean(γ, dims = 1), K, 1)
end

function posteriors!(
    γ::AbstractMatrix,
    c::AbstractVector,
    A::AbstractMatrix,
    LL::AbstractMatrix
)
    K = size(A, 1)
    a = similar(LL)
    α = log.(A[1, :])
    [a[:, k] = LL[:, k] .+ α[k] for k = 1:K]
    c[:] = logsumexp(a, dims = 2)
    γ[:, :] = exp.(a .- c)
end

function fit_em(mix::MixtureModel, y; display = :none, maxiter = 100, tol = 1e-3, robust = false)

    @argcheck display in [:none, :iter, :final]
    @argcheck maxiter >= 0

    N, K = length(y), ncomponents(mix)
    history = EMHistory(false, 0, [])

    # Allocate memory for in-place updates

    LL = zeros(N, K)
    γ = similar(LL)
    c = zeros(N)
    types = typeof.(components(mix))

    # Initial parameters
    α = copy(probs(mix))
    dists = copy(components(mix))

    # E-step
    # evaluate likelihood for each type k
    [LL[:, k] = log(α[k]) .+ logpdf.(dists[k], y) for k = 1:K]
    robust && replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))
    # get posterior of each category
    c[:] = logsumexp(LL, dims = 2)
    γ[:, :] = exp.(LL .- c)

    # Loglikelihood
    logtot = sum(c)
    (display == :iter) && println("Iteration 0: logtot = $logtot")

    for it = 1:maxiter

        # M-step
        # with γ in hand, maximize (update) the parameters
        α[:] = mean(γ, dims = 1)
        dists[:] = [fit_mle(types[k], y, γ[:, k]) for k = 1:K]

        # E-step
        # evaluate likelihood for each type k
        [LL[:, k] = log(α[k]) .+ logpdf.(dists[k], y) for k = 1:K]
        robust && replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))
        # get posterior of each category
        c[:] = logsumexp(LL, dims = 2)
        γ[:, :] = exp.(LL .- c)

        # Loglikelihood
        logtotp = sum(c)
        (display == :iter) && println("Iteration $it: logtot = $logtotp")

        push!(history.logtots, logtotp)
        history.iterations += 1

        if abs(logtotp - logtot) < tol
            (display in [:iter, :final]) &&
                println("EM converged in $it iterations, logtot = $logtotp")
            history.converged = true
            break
        end

        logtot = logtotp
    end

    if !history.converged
        if display in [:iter, :final]
            println("EM has not converged after $(history.iterations) iterations, logtot = $logtot")
        end
    end

    return MixtureModel(dists, α), history
end

# function fit_em_GE(mix::MixtureModel, y, M1; display = :none, maxiter = 100, tol = 1e-3, robust = false)
#
# 	@argcheck display in [:none, :iter, :final]
#     @argcheck maxiter >= 0
#
# 	N, K = length(y), ncomponents(mix)
# 	history = EMHistory(false, 0, [])
#
# 	# Allocate memory for in-place updates
#
#     LL = zeros(N,K)
#     γ = similar(LL)
# 	c = zeros(N)
# 	types = typeof.(components(mix))
# 	α = copy(probs(mix))
# 	dists = copy(components(mix))
#
# 	# evaluate likelihood for each type k
# 	[LL[:,k] = log(α[k]).+logpdf.(dists[k], y) for k in 1:K]
# 	robust && replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))
#
# 	# get posterior of each type
# 	c[:] = logsumexp(LL, dims = 2)
# 	γ[:,:] =  exp.(LL .- c)
#
# 	logtot = sum(c)
# 	(display == :iter) && println("Iteration 0: logtot = $logtot")
#
# 	for it = 1:maxiter
#
# 		# with p in hand, update
# 		α[:] = mean(γ, dims=1)
# 		θ = [mean(y, Weights(γ[:,k])) for k in 1:K]
# 		alpha = (M1 - θ[2]α[2])/(α[1]θ[1])
# 		dists[1] = Gamma(alpha, θ[1])
# 		dists[2] = Exponential(θ[2])
#
#         # evaluate likelihood for each type k
# 		[LL[:,k] = log(α[k]) .+ logpdf.(dists[k], y) for k in 1:K]
# 		robust && replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))
#
# 		# get posterior of each type
#         c[:] = logsumexp(LL, dims = 2)
#         γ[:,:] =  exp.(LL .- c)
#
#         # Log likelihood
# 		logtotp = sum(c)
# 		(display == :iter) && println("Iteration $it: logtot = $logtotp")
#
# 		push!(history.logtots, logtotp)
#         history.iterations += 1
#
# 		if abs(logtotp - logtot) < tol
#             (display in [:iter, :final]) &&
#                 println("EM converged in $it iterations, logtot = $logtotp")
#             history.converged = true
#             break
#         end
#
#         logtot = logtotp
#     end
#
#     return MixtureModel(dists, α), history
# end

# same as fit_em but a tad slower (howeve)
function fit_mle_mixture!(
    hmm::HMM,
    observations,
    ;
    display = :none,
    maxiter = 100,
    tol = 1e-3,
    robust = false,
    estimator = fit_mle
)
    @argcheck display in [:none, :iter, :final]
    @argcheck maxiter >= 0

    N, K = size(observations, 1), size(hmm, 1)
    history = EMHistory(false, 0, [])

    # Allocate memory for in-place updates

    c = zeros(N)
    γ = zeros(N, K)
    LL = zeros(N, K)

    loglikelihoods!(LL, hmm, observations)
    robust && replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))

    posteriors!(γ, c, hmm.A, LL)

    logtot = sum(c)
    (display == :iter) && println("Iteration 0: logtot = $logtot")

    for it = 1:maxiter

        update_A!(hmm.A, γ)
        update_B!(hmm.B, γ, observations, estimator)

        # Ensure the "connected-ness" of the states,
        # this prevents case where there is no transitions
        # between two extremely likely observations.
        robust && (hmm.A .+= eps())

        @check istransmat(hmm.A)

        loglikelihoods!(LL, hmm, observations)
        robust && replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))

        posteriors!(γ, c, hmm.A, LL)

        logtotp = sum(c)
        (display == :iter) && println("Iteration $it: logtot = $logtotp")

        push!(history.logtots, logtotp)
        history.iterations += 1

        if abs(logtotp - logtot) < tol
            (display in [:iter, :final]) &&
                println("EM converged in $it iterations, logtot = $logtotp")
            history.converged = true
            break
        end

        logtot = logtotp
    end
    update_a!(hmm.a, hmm.A)

    if !history.converged
        if display in [:iter, :final]
            println("EM has not converged after $(history.iterations) iterations, logtot = $logtot")
        end
    end
    history
end

## * Parametric model: Trigonometric function

function model_for_dist(dist::Exponential, θ_Y::AbstractArray, γ::AbstractVector, y, n2t, T::Int; silence = true)
    d = (size(θ_Y, 1) - 1) ÷ 2
    @argcheck length(y) == length(n2t)
    K = size(γ, 2)
    N = length(n2t)
    model = Model(Ipopt.Optimizer)
    silence && set_silent(model)
    f = 2π / T

    cos_nj = [cos(f * j * n2t[n]) for n = 1:N, j = 1:d]
    sin_nj = [sin(f * j * n2t[n]) for n = 1:N, j = 1:d]

    trig = [[1; interleave2(cos_nj[n, :], sin_nj[n, :])] for n = 1:N]

    @variable(model, θ_jump[j = 1:(2d+1)])
    # Polynomial P
    @NLexpression(model, Pθ[n = 1:N], sum(trig[n][j] * θ_jump[j] for j = 1:(2d+1)))

    @NLparameter(model, πk[n = 1:N] == γ[n])
    @NLparameter(model, R[n = 1:N] == y[n])

    @NLobjective(
        model, Max,
        -sum(πk[n] * (Pθ[n] + R[n] * exp(-Pθ[n])) for n = 1:N)
    )

    # I don't know if it is the best but https://discourse.julialang.org/t/jump-updating-nlparameter-of-a-model-in-a-loop/35081/3
    model[:πk] = πk
    model[:R] = R

    return model
end

function model_for_dist(dist::Gamma, θ_Y::AbstractArray, γ::AbstractVector, y, n2t, T::Int; silence = true)
    d = (size(θ_Y, 2) - 1) ÷ 2
    K = size(γ, 2)
    N = length(n2t)
    model = Model(Ipopt.Optimizer)
    set_optimizer_attributes(model, "max_iter" => 20)
    silence && set_silent(model)
    f = 2π / T
    cos_nj = [cos(f * j * n2t[n]) for n = 1:N, j = 1:d]
    sin_nj = [sin(f * j * n2t[n]) for n = 1:N, j = 1:d]

    trig = [[1; interleave2(cos_nj[n, :], sin_nj[n, :])] for n = 1:N]

    @variable(model, θ_jump[j = 1:(2d+1)])
    @variable(model, a_jump[j = 1:(2d+1)])
    # set_lower_bound(θ_jump[1], 0)
    # set_upper_bound.(θ_jump, 6)
    # set_lower_bound(a_jump[1], 0)
    # set_upper_bound(a_jump[1], 3)

    # Polynomial P
    @NLexpression(model, Pθ[n = 1:N], sum(trig[n][j] * θ_jump[j] for j = 1:length(trig[n])))
    @NLexpression(model, Pa[n = 1:N], sum(trig[n][j] * a_jump[j] for j = 1:length(trig[n])))

    @NLparameter(model, πk[n = 1:N] == γ[n])
    @NLparameter(model, R[n = 1:N] == y[n])
    register(model, :loggamma, 1, loggamma; autodiff = true) # cf the JuMP warning that prevent loggamma with more than one argument
    @NLobjective(
        model, Max,
        # -sum(πk[n]*(Pθ[n]+R[n]*exp(-Pθ[n])) for n in 1:N)
        # -sum(πk[n]*(exp(Pa[n])*Pθ[n]+R[n]*exp(-Pθ[n])+log(gamma(exp(Pa[n])))-expm1(Pa[n])*log(R[n])) for n in 1:N)
        -sum(πk[n] * (exp(Pa[n]) * Pθ[n] + R[n] * exp(-Pθ[n]) + loggamma(exp(Pa[n])) - expm1(Pa[n]) * log(R[n])) for n = 1:N)
    )
    # @NLexpression(model, mle,
    # sum(πk[n]*(-exp(Pa[n])*Pθ[n]-R[n]*exp(-Pθ[n])-loggamma(exp(Pa[n]))+expm1(Pa[n])*log(R[n])) for n in N)
    # )
    # I don't know if it is the best but https://discourse.julialang.org/t/jump-updating-nlparameter-of-a-model-in-a-loop/35081/3
    model[:πk] = πk
    model[:R] = R
    return model
end

function model_for_α(θ_α::AbstractMatrix, γ::AbstractMatrix, n2t, T::Int; silence = true)
    N, K = size(γ)
    d = (size(θ_α, 2) - 1) ÷ 2
    model = Model(Ipopt.Optimizer)
    silence && set_silent(model)
    f = 2π / T
    cos_nj = [cos(f * j * n2t[n]) for n = 1:N, j = 1:d]
    sin_nj = [sin(f * j * n2t[n]) for n = 1:N, j = 1:d]

    trig = [[1; interleave2(cos_nj[n, :], sin_nj[n, :])] for n = 1:N]
    @variable(model, θ_jump[k = 1:K-1, j = 1:(2d+1)])
    # Polynomial P
    @NLexpression(model, Pn[n = 1:N, k = 1:K-1], sum(trig[n][j] * θ_jump[k, j] for j = 1:length(trig[n])))

    @NLparameter(model, πk[n = 1:N, k = 1:K] == γ[n, k])
    @NLobjective(
        model,
        Max,
        sum(sum(πk[n, k] * Pn[n, k] for k = 1:K-1) - log1p(sum(exp(Pn[n, k]) for k = 1:K-1)) for n = 1:N) #Σ_k πk[n,k] = 1
    )
    # I don't know if it is the best but https://discourse.julialang.org/t/jump-updating-nlparameter-of-a-model-in-a-loop/35081/3
    model[:πk] = πk
    return model
end

function update_α!(
    α::AbstractArray where {T},
    θ::AbstractArray where {T},
    γ::AbstractMatrix,
    model::Model;
    warm_start = true
)
    N, K = size(γ)
    T = size(α, 2)

    # ξ are the filtering probablies
    θ_jump = model[:θ_jump]
    πk = model[:πk]
    ## Update the smoothing parameters in the JuMP model
    [set_value(πk[n, k], γ[n, k]) for n = 1:N, k = 1:K]
    warm_start && set_start_value.(θ_jump, θ)
    # Optimize the updated model
    set_silent(model)

    optimize!(model)
    # Obtained the new parameters
    θ[:, :] = value.(θ_jump)
    [α[k, t] = exp(polynomial_trigo(t, θ[k, :], T = T)) for k = 1:K-1, t = 1:T]
    [α[K, t] = 1 for t = 1:T] # last colum is 1/normalization (one could do otherwise)
    normalization_polynomial = [1 + sum(α[k, t] for k = 1:K-1) for t = 1:T]
    [α[k, t] /= normalization_polynomial[t] for k = 1:K, t = 1:T]
end

function update_dist(dist::AbstractVector{F} where {F<:Exponential}, θ_Y::AbstractArray, γ::AbstractVector, observations, model; warm_start = true)
    @argcheck size(γ, 1) == size(observations, 1)
    dist = copy(dist)
    θ_Y = copy(θ_Y)

    N = size(γ, 1)
    T = size(dist, 1)
    θ_jump = model[:θ_jump]
    πk = model[:πk]
    R = model[:R]

    ## Update the smoothing parameters in the JuMP model
    for n = 1:N
        set_value(πk[n], γ[n])
        set_value(R[n], observations[n])
    end

    warm_start && set_start_value.(θ_jump, θ_Y[:])
    set_silent(model)

    optimize!(model)
    θ_Y[:] = value.(θ_jump)

    p = [exp(polynomial_trigo(t, θ_Y[:], T = T)) for t = 1:T]
    dist[:] = Exponential.(p)
    return dist, θ_Y
end

function update_dist(dist::AbstractVector{F} where {F<:Gamma}, θ_Y::AbstractArray, γ::AbstractVector, observations, model; warm_start = true)
    @argcheck size(γ, 1) == size(observations, 1)
    dist = copy(dist)
    θ_Y = copy(θ_Y)

    N = size(γ, 1)
    T = size(dist, 1)
    θ_jump = model[:θ_jump]
    a_jump = model[:a_jump]
    πk = model[:πk]
    R = model[:R]

    ## Update the smoothing parameters in the JuMP model
    for n = 1:N
        set_value(πk[n], γ[n])
        set_value(R[n], observations[n])
    end

    warm_start && set_start_value.(a_jump, θ_Y[1, :]) # a parameter of Gamma
    warm_start && set_start_value.(θ_jump, θ_Y[2, :]) # θ parameter of Gamma

    optimize!(model)
    θ_Y[1, :] = value.(a_jump)
    θ_Y[2, :] = value.(θ_jump)

    a = [exp(polynomial_trigo(t, θ_Y[1, :], T = T)) for t = 1:T]
    p = [exp(polynomial_trigo(t, θ_Y[2, :], T = T)) for t = 1:T]
    dist[:] = [Gamma(a[t], p[t]) for t = 1:T]
    return dist, θ_Y

end

# * Improved and distributed compatible version

function model_for_dist_d(dist::Exponential, θ_Y::AbstractArray, T::Int; silence = true)
    d = (size(θ_Y, 1) - 1) ÷ 2

    model = Model(Ipopt.Optimizer)
    silence && set_silent(model)
    f = 2π / T

    cos_nj = [cos(f * j * t) for t = 1:T, j = 1:d]
    sin_nj = [sin(f * j * t) for t = 1:T, j = 1:d]
    trig = [[1; interleave2(cos_nj[t, :], sin_nj[t, :])] for t = 1:T]

    @variable(model, θ_jump[j = 1:(2d+1)])
    # Polynomial P
    @NLexpression(model, Pθ[t = 1:T], sum(trig[t][j] * θ_jump[j] for j = 1:(2d+1)))

    @NLparameter(model, s_π[t = 1:T] == 1)
    @NLparameter(model, s_π_R[t = 1:T] == 1) # same initialization but will be changed later

    """
        -sum(πk[n]*(Pθ[n]+R[n]*exp(-Pθ[n])) for n in 1:N)
    """
    @NLobjective(
        model, Max,
        -sum(s_π[t] * Pθ[t] + s_π_R[t] * exp(-Pθ[t]) for t = 1:T)
    )

    # I don't know if it is the best but https://discourse.julialang.org/t/jump-updating-nlparameter-of-a-model-in-a-loop/35081/3
    model[:s_π] = s_π
    model[:s_π_R] = s_π_R

    return model
end

function update_dist_d(dist::AbstractVector{F} where {F<:Exponential}, θ_Y::AbstractArray, γ::AbstractVector, observations::AbstractArray, n_in_t, model::Model; warm_start = true)
    @argcheck size(γ, 1) == size(observations, 1)
    dist = copy(dist)
    θ_Y = copy(θ_Y)

    N = size(γ, 1)
    T = size(dist, 1)
    θ_jump = model[:θ_jump]
    s_π = model[:s_π]
    s_π_R = model[:s_π_R]

    ## Update the smoothing parameters in the JuMP model
    for t = 1:T
        set_value(s_π[t], sum(γ[n] for n in n_in_t[t]; init = 0))
        set_value(s_π_R[t], sum(γ[n] * observations[n] for n in n_in_t[t]; init = 0))
    end

    warm_start && set_start_value.(θ_jump, θ_Y[:])

    optimize!(model)
    θ_Y[:] = value.(θ_jump)

    p = [exp(polynomial_trigo(t, θ_Y[:], T = T)) for t = 1:T]
    dist[:] = Exponential.(p)
    return dist, θ_Y
end

function model_for_dist_d(dist::Gamma, θ_Y::AbstractArray, T::Int; silence = true)
    d = (size(θ_Y, 2) - 1) ÷ 2
    model = Model(Ipopt.Optimizer)
    set_optimizer_attribute(model, "max_cpu_time", 60.0)
    set_optimizer_attribute(model, "max_iter", 100)
    silence && set_silent(model)
    f = 2π / T
    cos_nj = [cos(f * j * t) for t = 1:T, j = 1:d]
    sin_nj = [sin(f * j * t) for t = 1:T, j = 1:d]
    trig = [[1; interleave2(cos_nj[t, :], sin_nj[t, :])] for t = 1:T]

    @variable(model, θ_jump[j = 1:(2d+1)])
    @variable(model, a_jump[j = 1:(2d+1)])
    # set_lower_bound(θ_jump[1], 0)
    # set_upper_bound.(θ_jump, 6)
    # set_lower_bound(a_jump[1], 0)
    # set_upper_bound(a_jump[1], 3)

    # Polynomial P
    @NLexpression(model, Pθ[t = 1:T], sum(trig[t][j] * θ_jump[j] for j = 1:length(trig[t])))
    @NLexpression(model, Pa[t = 1:T], sum(trig[t][j] * a_jump[j] for j = 1:length(trig[t])))

    @NLparameter(model, s_π[t = 1:T] == 1)
    @NLparameter(model, s_π_R[t = 1:T] == 1)
    @NLparameter(model, s_π_logR[t = 1:T] == 1)

    register(model, :loggamma, 1, loggamma; autodiff = true) # cf the JuMP warning that prevent loggamma with more than one argument
    """
        - sum(πk[n]*(exp(Pa[n])*Pθ[n]+R[n]*exp(-Pθ[n])+loggamma(exp(Pa[n]))-expm1(Pa[n])*log(R[n])) for n in 1:N)
    """
    @NLobjective(
        model, Max,
        -sum(s_π[t] * exp(Pa[t]) * Pθ[t] + s_π_R[t] * exp(-Pθ[t])
             + s_π[t] * loggamma(exp(Pa[t])) - s_π_logR[t] * expm1(Pa[t]) for t = 1:T)
    )

    model[:s_π] = s_π
    model[:s_π_R] = s_π_R
    model[:s_π_logR] = s_π_logR

    return model
end

function update_dist_d(dist::AbstractVector{F} where {F<:Gamma}, θ_Y::AbstractArray, γ::AbstractVector, observations::AbstractArray, n_in_t, model::Model; warm_start = true)
    @argcheck size(γ, 1) == size(observations, 1)
    dist = copy(dist)
    θ_Y = copy(θ_Y)

    N = size(γ, 1)
    T = size(dist, 1)
    θ_jump = model[:θ_jump]
    a_jump = model[:a_jump]
    s_π = model[:s_π]
    s_π_R = model[:s_π_R]
    s_π_logR = model[:s_π_logR]

    ## Update the smoothing parameters in the JuMP model
    for t = 1:T
        set_value(s_π[t], sum(γ[n] for n in n_in_t[t]; init = 0))
        set_value(s_π_R[t], sum(γ[n] * observations[n] for n in n_in_t[t]; init = 0))
        set_value(s_π_logR[t], sum(γ[n] * log(observations[n]) for n in n_in_t[t]; init = 0))
    end

    warm_start && set_start_value.(a_jump, θ_Y[1, :]) # a parameter of Gamma
    warm_start && set_start_value.(θ_jump, θ_Y[2, :]) # θ parameter of Gamma

    optimize!(model)
    θ_Y[1, :] = value.(a_jump)
    θ_Y[2, :] = value.(θ_jump)

    a = [exp(polynomial_trigo(t, θ_Y[1, :], T = T)) for t = 1:T]
    p = [exp(polynomial_trigo(t, θ_Y[2, :], T = T)) for t = 1:T]
    dist[:] = [Gamma(a[t], p[t]) for t = 1:T]
    return dist, θ_Y

end

function s_γα!(s_ξ, ξ, n_in_t)
    T, K = size(s_ξ, 1), size(s_ξ, 2)
    for t = 1:T
        for k = 1:K
            s_ξ[t, k] = sum(ξ[n, k] for n in n_in_t[t]; init = 0)
        end
    end
end

function model_for_α_d(θ_α::AbstractMatrix, s_γ::AbstractMatrix; silence = true)
    T, K = size(s_γ)
    d = (size(θ_α, 2) - 1) ÷ 2
    model = Model(Ipopt.Optimizer)
    set_optimizer_attribute(model, "max_cpu_time", 60.0)
    set_optimizer_attribute(model, "max_iter", 100)
    silence && set_silent(model)
    f = 2π / T
    cos_nj = [cos(f * j * t) for t = 1:T, j = 1:d]
    sin_nj = [sin(f * j * t) for t = 1:T, j = 1:d]
    trig = [[1; interleave2(cos_nj[t, :], sin_nj[t, :])] for t = 1:T]

    @variable(model, θ_jump[k = 1:K-1, j = 1:(2d+1)])
    # Polynomial P
    @NLexpression(model, P[t = 1:T, k = 1:K-1], sum(trig[t][j] * θ_jump[k, j] for j = 1:(2d+1)))

    @NLparameter(model, s_π[t = 1:T, k = 1:K] == s_γ[t, k])
    @NLobjective(
        model,
        Max,
        sum(sum(s_π[t, k] * P[t, k] for k = 1:K-1) - sum(s_π[t, k] for k = 1:K) * log1p(sum(exp(P[t, k]) for k = 1:K-1)) for t = 1:T) #Σ_k πk[n,k] = 1
    )
    # I don't know if it is the best but https://discourse.julialang.org/t/jump-updating-nlparameter-of-a-model-in-a-loop/35081/3
    model[:s_π] = s_π
    return model
end

function update_α_d!(
    α::AbstractArray where {T},
    θ::AbstractArray where {T},
    s_γ::AbstractMatrix,
    model::Model;
    warm_start = true
)
    T, K = size(s_γ)

    # ξ are the filtering probablies
    θ_jump = model[:θ_jump]
    s_π = model[:s_π]
    ## Update the smoothing parameters in the JuMP model
    for t = 1:T, k = 1:K
        set_value(s_π[t, k], s_γ[t, k])
    end
    warm_start && set_start_value.(θ_jump, θ)
    # Optimize the updated model
    optimize!(model)
    # Obtained the new parameters
    θ[:, :] = value.(θ_jump)
    [α[k, t] = exp(polynomial_trigo(t, θ[k, :], T = T)) for k = 1:K-1, t = 1:T]
    [α[K, t] = 1 for t = 1:T] # last colum is 1/normalization (one could do otherwise)
    normalization_polynomial = [1 + sum(α[k, t] for k = 1:K-1) for t = 1:T]
    [α[k, t] /= normalization_polynomial[t] for k = 1:K, t = 1:T]
end

# TODO use the in place version (but it was bugy)

function fit_mle(mix::Vector{F} where {F<:MixtureModel}, θ_α::AbstractArray, θ_Y::AbstractArray, y, n2t;
    display = :none, maxiter = 100, tol = 1e-3, robust = false, silence = true, warm_start = true)
    @argcheck display in [:none, :iter, :final]
    @argcheck maxiter >= 0

    N, K, T = length(y), ncomponents(mix[1]), size(mix, 1)
    # TODO deal with K=1 case
    @argcheck K > 1
    d = (length(θ_α[1, :]) - 1) ÷ 2 #
    history = EMHistory(false, 0, [])
    # Allocate memory for in-place updates

    LL = zeros(N, K)
    γ = similar(LL)
    s_γ = zeros(T, K)

    c = zeros(N)

    # Initial parameters
    θ_α = copy(θ_α)
    θ_Y = copy(θ_Y)
    n_in_t = [findall(n2t .== t) for t = 1:T]

    α = copy(hcat(probs.(mix)...))
    dists = copy(hcat(components.(mix)...))
    types = typeof.(dists[:, 1])
    dists = [convert(Vector{types[k, 1]}, dists[k, :]) for k = 1:K] # I have to separate the array to preserve types

    model_α = model_for_α_d(θ_α, s_γ; silence = silence)
    model_dist = [model_for_dist_d(dists[k][1], θ_Y[k], T; silence = silence) for k = 1:K]

    # E-step
    # evaluate likelihood for each type k
    for k = 1:K, n = 1:N
        LL[n, k] = log(α[k, n2t[n]]) + logpdf(dists[k][n2t[n]], y[n])
    end
    robust && replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))
    # get posterior of each category
    c[:] = logsumexp(LL, dims = 2)
    γ[:, :] = exp.(LL .- c)

    # Loglikelihood
    logtot = sum(c)
    (display == :iter) && println("Iteration 0: logtot = $logtot")

    for it = 1:maxiter

        # M-step
        # with γ in hand, maximize (update) the parameters
        s_γα!(s_γ, γ, n_in_t)
        update_α_d!(α, θ_α, s_γ, model_α; warm_start = warm_start)
        for k = 1:K
            dists[k][:], θ_Y[k] = update_dist_d(dists[k], θ_Y[k], γ[:, k], y, n_in_t, model_dist[k], warm_start = warm_start)
        end
        # E-step
        # evaluate likelihood for each type k
        for k = 1:K, n = 1:N
            LL[n, k] = log(α[k, n2t[n]]) + logpdf(dists[k][n2t[n]], y[n])
        end
        robust && replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))
        # get posterior of each category
        c[:] = logsumexp(LL, dims = 2)
        γ[:, :] = exp.(LL .- c)

        # Loglikelihood
        logtotp = sum(c)
        (display == :iter) && println("Iteration $it: logtot = $logtotp")

        push!(history.logtots, logtotp)
        history.iterations += 1

        if abs(logtotp - logtot) < tol
            (display in [:iter, :final]) &&
                println("EM converged in $it iterations, logtot = $logtotp")
            history.converged = true
            break
        end

        logtot = logtotp
    end

    if !history.converged
        if display in [:iter, :final]
            println("EM has not converged after $(history.iterations) iterations, logtot = $logtot")
        end
    end

    return [MixtureModel([dists[k][t] for k = 1:K], [α[k, t] for k = 1:K]) for t = 1:T], θ_α, θ_Y, history
end

function fit_em_GE(mix::Vector{F} where {F<:MixtureModel}, θ_α::AbstractArray, θ_Y::AbstractArray, y, n2t;
    display = :none, maxiter = 100, tol = 1e-3, robust = false, silence = true, warm_start = true)
    @argcheck display in [:none, :iter, :final]
    @argcheck maxiter >= 0

    N, K, T = length(y), ncomponents(mix[1]), size(mix, 1)
    # TODO deal with K=1 case
    @argcheck K > 1
    d = (length(θ_α[1, :]) - 1) ÷ 2 #
    history = EMHistory(false, 0, [])
    # Allocate memory for in-place updates

    LL = zeros(N, K)
    γ = similar(LL)
    c = zeros(N)

    # Initial parameters
    θ_α = copy(θ_α)
    θ_Y = copy(θ_Y)

    α = copy(hcat(probs.(mix)...))
    dists = copy(hcat(components.(mix)...))
    types = typeof.(dists[:, 1])
    dists = [convert(Vector{types[k, 1]}, dists[k, :]) for k = 1:K] # I have to separate the array to preserve types

    model_α = model_for_α(θ_α, γ, n2t, T; silence = silence)
    model_dist = [model_for_dist(dists[k][1], θ_Y[k], γ[:, k], y, n2t, T; silence = silence) for k = 1:K]

    # E-step
    # evaluate likelihood for each type k
    for k = 1:K, n = 1:N
        LL[n, k] = log(α[k, n2t[n]]) + logpdf(dists[k][n2t[n]], y[n])
    end
    robust && replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))
    # get posterior of each category
    c[:] = logsumexp(LL, dims = 2)
    γ[:, :] = exp.(LL .- c)

    # Loglikelihood
    logtot = sum(c)
    (display == :iter) && println("Iteration 0: logtot = $logtot")

    for it = 1:maxiter

        # M-step
        # with γ in hand, maximize (update) the parameters
        update_α!(α, θ_α, γ, model_α; warm_start = warm_start)
        for k = 1:K
            dists[k][:], θ_Y[k] = update_dist(dists[k], θ_Y[k], γ[:, k], y, model_dist[k], warm_start = warm_start)
            # update_dist!(dists[k], θ_Y[k], γ[:,k], y, model_dist[k], warm_start = warm_start)
        end
        # E-step
        # evaluate likelihood for each type k
        for k = 1:K, n = 1:N
            LL[n, k] = log(α[k, n2t[n]]) + logpdf(dists[k][n2t[n]], y[n])
        end
        robust && replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))
        # get posterior of each category
        c[:] = logsumexp(LL, dims = 2)
        γ[:, :] = exp.(LL .- c)

        # Loglikelihood
        logtotp = sum(c)
        (display == :iter) && println("Iteration $it: logtot = $logtotp")

        push!(history.logtots, logtotp)
        history.iterations += 1

        if abs(logtotp - logtot) < tol
            (display in [:iter, :final]) &&
                println("EM converged in $it iterations, logtot = $logtotp")
            history.converged = true
            break
        end

        logtot = logtotp
    end

    if !history.converged
        if display in [:iter, :final]
            println("EM has not converged after $(history.iterations) iterations, logtot = $logtot")
        end
    end

    return [MixtureModel([dists[k][t] for k = 1:K], [α[k, t] for k = 1:K]) for t = 1:T], θ_α, θ_Y, history
end


text:
"""
Provides a rendering function, [`render`](@ref), for writing each supported
[`Formats.Format`](@ref) to file.

Note that currently `Formats.Markdown` is the **only** supported format.

"""
module Writers

import ..Documenter:
    Anchors,
    Builder,
    Documents,
    Expanders,
    Formats,
    Documenter,
    Utilities

using Compat

"""
A parametric type that allows us to use multiple dispatch to pick the appropriate
writer for each output format.

The parameter `f` should be an instance of the [`Formats.Format`](@ref) enumeration.
"""
immutable Writer{f} end
Writer(f::Formats.Format) = Writer{f}()

"""
Writes a [`Documents.Document`](@ref) object to `.user.build` directory in
the format specified in `.user.format`.

The method should be overloaded in each writer as

    render(::Writer{format}, doc)

where `format` is one of the values of the [`Formats.Format`](@ref) enumeration.
"""
render(doc::Documents.Document) = render(Writer(doc.user.format), doc)

include("MarkdownWriter.jl")
include("HTMLWriter.jl")
include("LaTeXWriter.jl")

end


text:
using CUDArt
using Base.Test

##########################
# Low-level memory tests #
##########################
devlist = CUDArt.devices(dev->true)
for dev in devlist
    CUDArt.device(dev)
    p = CUDArt.malloc(UInt8, 1)
    p2 = CUDArt.malloc(UInt16, 100)
    CUDArt.free(p)
    CUDArt.free(p2)
    a = rand(5,3)
    g = CUDArt.CudaArray(a)
    gp = CUDArt.CudaPitchedArray(a)
    CUDArt.free(g)
    CUDArt.free(gp)
    # Also test finalizer calls
    g = CUDArt.CudaArray(a)
    g = CUDArt.CudaPitchedArray(a)
    CUDArt.device_reset(dev)
end

gc()  # check for finalizer errors

# Test free when the allocated device is not the current device
if length(devlist) > 1
    CUDArt.device(devlist[1])
    p1a = CUDArt.malloc(UInt8, 1)
    p1 = CUDArt.malloc(UInt16, 100)
    CUDArt.device(devlist[2])
    p2 = CUDArt.malloc(UInt16, 100)
    CUDArt.free(p1)
    CUDArt.free(p2)
    CUDArt.free(p1a)
    CUDArt.device_reset(devlist[1])
    CUDArt.device_reset(devlist[2])
end

gc()

# Check whether gc/finalizers lead to bad behavior across
# device resets (the pointer address is reproducible)
# BROKEN: see JuliaGPU/CUDArt.jl#66
if false
    CUDArt.device(0)
    a = rand(5,3)
    g = CUDArt.CudaArray(a)
    CUDArt.device_reset(0)
    @test isempty(CUDArt.cuda_ptrs)
    CUDArt.device(0)
    g = CUDArt.CudaArray(a)
    # In case the next gc() test yields an error, the next two lines let
    # us do some archaeology
    dictcopy = deepcopy(CUDArt.cuda_ptrs)
    gptrcopy = copy(pointer(g))
    gc()   # Check that this doesn't delete the new g
    @test !isempty(CUDArt.cuda_ptrs)
    h_g = CUDArt.to_host(g)
    CUDArt.free(g)
    @test isempty(CUDArt.cuda_ptrs)
    @test h_g == a
    CUDArt.device_reset(0)
end

gc()


text:
"""
The full list of keyword arguments that can be passed to [`Pluto.run`](@ref) (or [`Pluto.Configuration.from_flat_kwargs`](@ref)) is divided into four categories. Take a look at the documentation for:

- [`Pluto.Configuration.CompilerOptions`](@ref) defines the command line arguments for notebook `julia` processes.
- [`Pluto.Configuration.ServerOptions`](@ref) configures the HTTP server.
- [`Pluto.Configuration.SecurityOptions`](@ref) configures the authentication options for Pluto's HTTP server. Change with caution.
- [`Pluto.Configuration.EvaluationOptions`](@ref) is used internally during Pluto's testing.

Note that Pluto is designed to be _zero-configuration_, and most users should not (have to) change these settings. Most 'customization' can be achieved using Julia's wide range of packages! That being said, the available settings are useful if you are using Pluto in a special environment, such as docker, mybinder, etc.
"""
module Configuration

using Configurations # https://github.com/Roger-luo/Configurations.jl

import ..Pluto: tamepath

# Using a ref to avoid fixing the pwd() output during the compilation phase. We don't want this value to be baked into the sysimage, because it depends on the `pwd()`. We do want to cache it, because the pwd might change while Pluto is running.
const pwd_ref = Ref{Union{Nothing,String}}()
function notebook_path_suggestion()
    pwd_val = something(pwd_ref[], pwd())
    preferred_dir = startswith(Sys.BINDIR, pwd_val) ? homedir() : pwd_val
    # so that it ends with / or \
    string(joinpath(preferred_dir, ""))
end

function __init__()
    pwd_ref[] = pwd()
end

const ROOT_URL_DEFAULT = nothing
const HOST_DEFAULT = "127.0.0.1"
const PORT_DEFAULT = nothing
const LAUNCH_BROWSER_DEFAULT = true
const DISMISS_UPDATE_NOTIFICATION_DEFAULT = false
const SHOW_FILE_SYSTEM_DEFAULT = true
const DISABLE_WRITING_NOTEBOOK_FILES_DEFAULT = false
const AUTO_RELOAD_FROM_FILE_DEFAULT = false
const AUTO_RELOAD_FROM_FILE_COOLDOWN_DEFAULT = 0.4
const AUTO_RELOAD_FROM_FILE_IGNORE_PKG_DEFAULT = false
const NOTEBOOK_DEFAULT = nothing
const INIT_WITH_FILE_VIEWER_DEFAULT = false
const SIMULATED_LAG_DEFAULT = 0.0
const SIMULATED_PKG_LAG_DEFAULT = 0.0
const INJECTED_JAVASCRIPT_DATA_URL_DEFAULT = "data:text/javascript;base64,"
const ON_EVENT_DEFAULT = function(a) #= @info "$(typeof(a))" =# end

"""
    ServerOptions([; kwargs...])

The HTTP server options. See [`SecurityOptions`](@ref) for additional settings.

# Arguments

- `root_url::Union{Nothing,String} = $ROOT_URL_DEFAULT`
- `host::String = "$HOST_DEFAULT"`
- `port::Union{Nothing,Integer} = $PORT_DEFAULT`
- `launch_browser::Bool = $LAUNCH_BROWSER_DEFAULT`
- `dismiss_update_notification::Bool = $DISMISS_UPDATE_NOTIFICATION_DEFAULT`
- `show_file_system::Bool = $SHOW_FILE_SYSTEM_DEFAULT`
- `notebook_path_suggestion::String = notebook_path_suggestion()`
- `disable_writing_notebook_files::Bool = $DISABLE_WRITING_NOTEBOOK_FILES_DEFAULT`
- `auto_reload_from_file::Bool = $AUTO_RELOAD_FROM_FILE_DEFAULT` Watch notebook files for outside changes and update running notebook state automatically
- `auto_reload_from_file_cooldown::Real = $AUTO_RELOAD_FROM_FILE_COOLDOWN_DEFAULT` Experimental, will be removed
- `auto_reload_from_file_ignore_pkg::Bool = $AUTO_RELOAD_FROM_FILE_IGNORE_PKG_DEFAULT` Experimental flag, will be removed
- `notebook::Union{Nothing,String} = $NOTEBOOK_DEFAULT` Optional path of notebook to launch at start
- `init_with_file_viewer::Bool = $INIT_WITH_FILE_VIEWER_DEFAULT`
- `simulated_lag::Real=$SIMULATED_LAG_DEFAULT` (internal) Extra lag to add to our server responses. Will be multiplied by `0.5 + rand()`.
- `simulated_pkg_lag::Real=$SIMULATED_PKG_LAG_DEFAULT` (internal) Extra lag to add to operations done by Pluto's package manager. Will be multiplied by `0.5 + rand()`.
- `injected_javascript_data_url::String = "$INJECTED_JAVASCRIPT_DATA_URL_DEFAULT"` (internal) Optional javascript injectables to the front-end. Can be used to customize the editor, but this API is not meant for general use yet.
- `on_event::Function = $ON_EVENT_DEFAULT`
"""
@option mutable struct ServerOptions
    root_url::Union{Nothing,String} = ROOT_URL_DEFAULT
    host::String = HOST_DEFAULT
    port::Union{Nothing,Integer} = PORT_DEFAULT
    launch_browser::Bool = LAUNCH_BROWSER_DEFAULT
    dismiss_update_notification::Bool = DISMISS_UPDATE_NOTIFICATION_DEFAULT
    show_file_system::Bool = SHOW_FILE_SYSTEM_DEFAULT
    notebook_path_suggestion::String = notebook_path_suggestion()
    disable_writing_notebook_files::Bool = DISABLE_WRITING_NOTEBOOK_FILES_DEFAULT
    auto_reload_from_file::Bool = AUTO_RELOAD_FROM_FILE_DEFAULT
    auto_reload_from_file_cooldown::Real = AUTO_RELOAD_FROM_FILE_COOLDOWN_DEFAULT
    auto_reload_from_file_ignore_pkg::Bool = AUTO_RELOAD_FROM_FILE_IGNORE_PKG_DEFAULT
    notebook::Union{Nothing,String,Vector{<:String}} = NOTEBOOK_DEFAULT
    init_with_file_viewer::Bool = INIT_WITH_FILE_VIEWER_DEFAULT
    simulated_lag::Real = SIMULATED_LAG_DEFAULT
    simulated_pkg_lag::Real = SIMULATED_PKG_LAG_DEFAULT
    injected_javascript_data_url::String = INJECTED_JAVASCRIPT_DATA_URL_DEFAULT
    on_event::Function = ON_EVENT_DEFAULT
end

const REQUIRE_SECRET_FOR_OPEN_LINKS_DEFAULT = true
const REQUIRE_SECRET_FOR_ACESS_DEFAULT = true

"""
    SecurityOptions([; kwargs...])

Security settings for the HTTP server.

# Arguments

- `require_secret_for_open_links::Bool = $REQUIRE_SECRET_FOR_OPEN_LINKS_DEFAULT`

    Whether the links `http://localhost:1234/open?path=/a/b/c.jl`  and `http://localhost:1234/open?url=http://www.a.b/c.jl` should be protected.

    Use `true` for almost every setup. Only use `false` if Pluto is running in a safe container (like mybinder.org), where arbitrary code execution is not a problem.

- `require_secret_for_access::Bool = $REQUIRE_SECRET_FOR_ACESS_DEFAULT`

    If false, you do not need to use a `secret` in the URL to access Pluto: you will be authenticated by visiting `http://localhost:1234/` in your browser. An authentication cookie is still used for access (to prevent XSS and deceptive links or an img src to `http://localhost:1234/open?url=badpeople.org/script.jl`), and is set automatically, but this request to `/` is protected by cross-origin policy.

    Use `true` on a computer used by multiple people simultaneously. Only use `false` if necessary.

**Leave these options on `true` for the most secure setup.**

Note that Pluto is quickly evolving software, maintained by designers, educators and enthusiasts — not security experts. If security is a serious concern for your application, then we recommend running Pluto inside a container and verifying the relevant security aspects of Pluto yourself.
"""
@option mutable struct SecurityOptions
    require_secret_for_open_links::Bool = REQUIRE_SECRET_FOR_OPEN_LINKS_DEFAULT
    require_secret_for_access::Bool = REQUIRE_SECRET_FOR_ACESS_DEFAULT
end

const RUN_NOTEBOOK_ON_LOAD_DEFAULT = true
const WORKSPACE_USE_DISTRIBUTED_DEFAULT = true
const LAZY_WORKSPACE_CREATION_DEFAULT = false
const CAPTURE_STDOUT_DEFAULT = true

"""
    EvaluationOptions([; kwargs...])

Options to change Pluto's evaluation behaviour during internal testing. These options are not intended to be changed during normal use.

- `run_notebook_on_load::Bool = $RUN_NOTEBOOK_ON_LOAD_DEFAULT` Whether to evaluate a notebook on load.
- `workspace_use_distributed::Bool = $WORKSPACE_USE_DISTRIBUTED_DEFAULT` Whether to start notebooks in a separate process.
- `lazy_workspace_creation::Bool = $LAZY_WORKSPACE_CREATION_DEFAULT`
- `capture_stdout::Bool = $CAPTURE_STDOUT_DEFAULT`
"""
@option mutable struct EvaluationOptions
    run_notebook_on_load::Bool = RUN_NOTEBOOK_ON_LOAD_DEFAULT
    workspace_use_distributed::Bool = WORKSPACE_USE_DISTRIBUTED_DEFAULT
    lazy_workspace_creation::Bool = LAZY_WORKSPACE_CREATION_DEFAULT
    capture_stdout::Bool = CAPTURE_STDOUT_DEFAULT
end

const COMPILE_DEFAULT = nothing
const SYSIMAGE_DEFAULT = nothing
const BANNER_DEFAULT = nothing
const OPTIMIZE_DEFAULT = nothing
const MATH_MODE_DEFAULT = nothing
const STARTUP_FILE_DEFAULT = "no"
const HISTORY_FILE_DEFAULT = "no"

function roughly_the_number_of_physical_cpu_cores()
    # https://gist.github.com/fonsp/738fe244719cae820245aa479e7b4a8d
    if Sys.CPU_THREADS == 1
        1
    elseif Sys.CPU_THREADS == 2 || Sys.CPU_THREADS == 3
        2
    else
        Sys.CPU_THREADS ÷ 2
    end
end

function default_number_of_threads()
    env_value = get(ENV, "JULIA_NUM_THREADS", "")

    all(isspace, env_value) ?
        roughly_the_number_of_physical_cpu_cores() :
        env_value
end

"""
    CompilerOptions([; kwargs...])

These options will be passed as command line argument to newly launched processes. See [the Julia documentation on command-line options](https://docs.julialang.org/en/v1/manual/command-line-options/).

# Arguments
- `compile::Union{Nothing,String} = $COMPILE_DEFAULT`
- `sysimage::Union{Nothing,String} = $SYSIMAGE_DEFAULT`
- `banner::Union{Nothing,String} = $BANNER_DEFAULT`
- `optimize::Union{Nothing,Int} = $OPTIMIZE_DEFAULT`
- `math_mode::Union{Nothing,String} = $MATH_MODE_DEFAULT`
- `startup_file::Union{Nothing,String} = "$STARTUP_FILE_DEFAULT"` By default, the startup file isn't loaded in notebooks.
- `history_file::Union{Nothing,String} = "$HISTORY_FILE_DEFAULT"` By default, the history isn't loaded in notebooks.
- `threads::Union{Nothing,String,Int} = default_number_of_threads()`
"""
@option mutable struct CompilerOptions
    compile::Union{Nothing,String} = COMPILE_DEFAULT
    sysimage::Union{Nothing,String} = SYSIMAGE_DEFAULT
    banner::Union{Nothing,String} = BANNER_DEFAULT
    optimize::Union{Nothing,Int} = OPTIMIZE_DEFAULT
    math_mode::Union{Nothing,String} = MATH_MODE_DEFAULT

    # notebook specified options
    # the followings are different from
    # the default julia compiler options
    startup_file::Union{Nothing,String} = STARTUP_FILE_DEFAULT
    history_file::Union{Nothing,String} = HISTORY_FILE_DEFAULT

    threads::Union{Nothing,String,Int} = default_number_of_threads()
end

"""
Collection of all settings that configure a Pluto session.

`ServerSession` contains a `Configuration`.
"""
@option struct Options
    server::ServerOptions = ServerOptions()
    security::SecurityOptions = SecurityOptions()
    evaluation::EvaluationOptions = EvaluationOptions()
    compiler::CompilerOptions = CompilerOptions()
end

function from_flat_kwargs(;
        root_url::Union{Nothing,String} = ROOT_URL_DEFAULT,
        host::String = HOST_DEFAULT,
        port::Union{Nothing,Integer} = PORT_DEFAULT,
        launch_browser::Bool = LAUNCH_BROWSER_DEFAULT,
        dismiss_update_notification::Bool = DISMISS_UPDATE_NOTIFICATION_DEFAULT,
        show_file_system::Bool = SHOW_FILE_SYSTEM_DEFAULT,
        notebook_path_suggestion::String = notebook_path_suggestion(),
        disable_writing_notebook_files::Bool = DISABLE_WRITING_NOTEBOOK_FILES_DEFAULT,
        auto_reload_from_file::Bool = AUTO_RELOAD_FROM_FILE_DEFAULT,
        auto_reload_from_file_cooldown::Real = AUTO_RELOAD_FROM_FILE_COOLDOWN_DEFAULT,
        auto_reload_from_file_ignore_pkg::Bool = AUTO_RELOAD_FROM_FILE_IGNORE_PKG_DEFAULT,
        notebook::Union{Nothing,String,Vector{<:String}} = NOTEBOOK_DEFAULT,
        init_with_file_viewer::Bool = INIT_WITH_FILE_VIEWER_DEFAULT,
        simulated_lag::Real = SIMULATED_LAG_DEFAULT,
        simulated_pkg_lag::Real = SIMULATED_PKG_LAG_DEFAULT,
        injected_javascript_data_url::String = INJECTED_JAVASCRIPT_DATA_URL_DEFAULT,
        on_event::Function = ON_EVENT_DEFAULT,
        require_secret_for_open_links::Bool = REQUIRE_SECRET_FOR_OPEN_LINKS_DEFAULT,
        require_secret_for_access::Bool = REQUIRE_SECRET_FOR_ACESS_DEFAULT,
        run_notebook_on_load::Bool = RUN_NOTEBOOK_ON_LOAD_DEFAULT,
        workspace_use_distributed::Bool = WORKSPACE_USE_DISTRIBUTED_DEFAULT,
        lazy_workspace_creation::Bool = LAZY_WORKSPACE_CREATION_DEFAULT,
        capture_stdout::Bool = CAPTURE_STDOUT_DEFAULT,
        compile::Union{Nothing,String} = COMPILE_DEFAULT,
        sysimage::Union{Nothing,String} = SYSIMAGE_DEFAULT,
        banner::Union{Nothing,String} = BANNER_DEFAULT,
        optimize::Union{Nothing,Int} = OPTIMIZE_DEFAULT,
        math_mode::Union{Nothing,String} = MATH_MODE_DEFAULT,
        startup_file::Union{Nothing,String} = STARTUP_FILE_DEFAULT,
        history_file::Union{Nothing,String} = HISTORY_FILE_DEFAULT,
        threads::Union{Nothing,String,Int} = default_number_of_threads(),
    )
    server = ServerOptions(;
        root_url,
        host,
        port,
        launch_browser,
        dismiss_update_notification,
        show_file_system,
        notebook_path_suggestion,
        disable_writing_notebook_files,
        auto_reload_from_file,
        auto_reload_from_file_cooldown,
        auto_reload_from_file_ignore_pkg,
        notebook,
        init_with_file_viewer,
        simulated_lag,
        simulated_pkg_lag,
        injected_javascript_data_url,
        on_event,
    )
    security = SecurityOptions(;
        require_secret_for_open_links,
        require_secret_for_access,
    )
    evaluation = EvaluationOptions(;
        run_notebook_on_load,
        workspace_use_distributed,
        lazy_workspace_creation,
        capture_stdout,
    )
    compiler = CompilerOptions(;
        compile,
        sysimage,
        banner,
        optimize,
        math_mode,
        startup_file,
        history_file,
        threads,
    )
    return Options(; server, security, evaluation, compiler)
end

function _merge_notebook_compiler_options(notebook, options::CompilerOptions)::CompilerOptions
    if notebook.compiler_options === nothing
        return options
    end

    kwargs = Dict{Symbol,Any}()
    for each in fieldnames(CompilerOptions)
        # 1. not specified by notebook options
        # 2. general notebook specified options
        if getfield(notebook.compiler_options, each) === nothing
            kwargs[each] = getfield(options, each)
        else
            kwargs[each] = getfield(notebook.compiler_options, each)
        end
    end
    return CompilerOptions(; kwargs...)
end

function _convert_to_flags(options::CompilerOptions)::Vector{String}
    option_list = String[]

    for name in fieldnames(CompilerOptions)
        flagname = string("--", replace(String(name), "_" => "-"))
        value = getfield(options, name)
        if value !== nothing
            push!(option_list, string(flagname, "=", value))
        end
    end

    return option_list
end

end


text:
using Documenter
using Models
using Models.TestUtils

makedocs(;
    modules=[Models],
    format=Documenter.HTML(;
        prettyurls=false,
        assets=[
            "assets/invenia.css",
        ],
    ),
    pages=[
        "Index" => "index.md",
        "Design" => "design.md",
        "API" => "api.md",
        "TestUtils" => "testutils.md",
    ],
    repo="https://github.com/invenia/Models.jl/blob/{commit}{path}#L{line}",
    sitename="Models.jl",
    authors="Invenia Technical Computing Corporation",
    strict=true,
    checkdocs=:exports,
)

deploydocs(;
    repo="github.com/invenia/Models.jl",
    push_preview=true,
)


text:
# This file was generated, do not modify it. # hide
(WtoV ∘ VtoW)(V) ≈ V  &&  (VtoW ∘ WtoV)(W) ≈ W

max_stars_repo_name: cscherrer/InformativePrior

max_issues_repo_name: cscherrer/InformativePrior

max_forks_repo_name: cscherrer/InformativePrior


text:
# code for testing our package
using Hydraulics2, Test

@testset "moody" begin
    # test case moody
	# from https://neutrium.net/fluid-flow/pressure-loss-from-fittings-3k-method/
	# Reynolds number 306900, e = 0.068 mm, D = 102.3
	eD = 0.068 / 102.3
	Reynolds = 306900.0
    	@test (abs(Hydraulics2.calcMoodyF(Reynolds,eD) - 0.018) < 0.001)
end


text:
mutable struct Gaussian{T, RT <: Real, O} <: AbstractBlock
    μ::T
    σ::RT
    grad_μ::T
    grad_σ::RT

    input::O
    output::O

    Gaussian{T, RT, O}(mu, sigma) where {T, RT, O} =
        new{T, RT, O}(mu, sigma, zero(T), zero(RT))
end

Gaussian(::Type{O}, mu::T, sigma::RT) where {T, RT, O} =
    Gaussian{T, real(T), O}(mu, sigma)

Gaussian(mu::T, sigma::RT; nbatch=1) where {T, RT} =
    Gaussian(nbatch > 1 ? Matrix{T} : Vector{T}, mu, sigma)

# TODO: check if generated function will be faster (probably not)
function forward(op::Gaussian{T, RT, O}, input::O) where {T, RT, O}
    op.input = input
    factor = 1 / (sqrt(2 * pi) * op.σ)
    op.output = factor .* exp.(- abs2.(input .- op.μ) ./ (2 * op.σ^2))
    return op.output
end

function backward(op::Gaussian{T, RT, O}, grad::O) where {T <: Complex, RT, O}
    op.grad_μ = sum(real(op.input .- op.μ) .* op.output ./ op.σ^2 .* grad) # Trust your compiler
    op.grad_σ = sum(((abs2(op.input .- op.μ)) .- op.σ^2) ./ op.σ^3 .* op.output .* grad)
    return -(conj(op.input .- op.μ) .* op.output) ./ op.σ^2 .* grad
end


text:
include("../src/PBCM_IO.jl")

cell_general = cell()
cell_design_op = cylindrical_cell_designer(cell_general)
cost = cost_default()
cell_general, cost = convert_all(cell_general, cost, mult)

using PyPlot
using DiffEqSensitivity

per = 10

poro = 0.35
poro_intrvl = [poro*(1 - (per/100)) , poro*(1 + (per/100))]

pos_thic = 150
pos_thic_intrvl = [pos_thic*(1 - (per/100)) , pos_thic*(1 + (per/100))]

neg_thic = 150
neg_thic_intrvl = [neg_thic*(1 - (per/100)) , neg_thic*(1 + (per/100))]

no_cells_yr = 10
no_cells_yr_intrvl = [no_cells_yr*(1 - (per/100)) , no_cells_yr*(1 + (per/100))]

cost_pos_AM = 25
cost_pos_AM_intrvl = [cost_pos_AM*(1 - (per/100)) , cost_pos_AM*(1 + (per/100))]

function gsa_cost(arr)

    cell_general.cathode.por        = arr[1]
    cell_general.cathode.th         = converter(arr[2], mult.pos_th)
    cell_general.anode.th           = converter(arr[3], mult.neg_th)
    cost.general_costs.no_units_mfg = converter(arr[4], mult.units_mfg)
    cost.cell_costs.cathode.AM[1]   = arr[5]

    dollars_kWh = cost_calc(cell_general, cost, system="Cell", cost_verbosity=0)[1]

    return dollars_kWh
end

# gsa_cost([0.25, 150.0, 130.0, 25.0, 19.0])

gsa_result = gsa(gsa_cost, Sobol(), [poro_intrvl, pos_thic_intrvl, neg_thic_intrvl, no_cells_yr_intrvl, cost_pos_AM_intrvl], N = 10000)

print("\n")
S1 = gsa_result.S1
print(S1)
clf()
x = ["Pos Porosity", "Pos Th", "Neg Th", "Cells MFG", "Pos AM Cost"]
bar(x, S1, color="#0f87bf", align="center", alpha=0.8)
ylabel("Sobol Indices")

title(string("GSA: 1st Order : ", per, " %"))
figure(2)


text:
# conv test

function prepare_input(input_size)
	return (rand(input_size),rand(input_size))
end

function test_function(input_data)
	return conv(input_data[1], input_data[2])
end


text:
#=
  Test set 3: Compare the results from test sets 1 and 2.
=#

println("\nTEST SET 3")

outfile = joinpath(outdir, "linkage_comparison.csv")
compare_linkage_runs(outdir1b, outdir2, outfile)
result  = table_to_set_of_dicts(outfile)

@test length(result) == 14
@test in(Dict(:tablename => "LINKAGE RUNS", :status1 => outdir1b, :status2 => outdir2, :nrecords => missing), result)
@test in(Dict(:tablename => "spine", :status1 => "existent",    :status2 => "existent",    :nrecords => 3), result)
@test in(Dict(:tablename => "spine", :status1 => "existent",    :status2 => "nonexistent", :nrecords => 1), result)
@test in(Dict(:tablename => "spine", :status1 => "nonexistent", :status2 => "existent",    :nrecords => 3), result)
@test in(Dict(:tablename => "hospital_admissions", :status1 => "linked with criteria ID 3", :status2 => "linked with criteria ID 1", :nrecords => 3), result)
@test in(Dict(:tablename => "hospital_admissions", :status1 => "unlinked",                  :status2 => "linked with criteria ID 1", :nrecords => 2), result)
@test in(Dict(:tablename => "emergency_presentations", :status1 => "linked with criteria ID 1", :status2 => "linked with criteria ID 2", :nrecords => 2), result)
@test in(Dict(:tablename => "emergency_presentations", :status1 => "linked with criteria ID 2", :status2 => "linked with criteria ID 3", :nrecords => 1), result)
@test in(Dict(:tablename => "emergency_presentations", :status1 => "unlinked", :status2 => "linked with criteria ID 2", :nrecords => 2), result)
@test in(Dict(:tablename => "influenza_cases", :status1 => "linked with criteria ID 4", :status2 => "linked with criteria ID 4", :nrecords => 3), result)
@test in(Dict(:tablename => "influenza_cases", :status1 => "linked with criteria ID 4", :status2 => "linked with criteria ID 7", :nrecords => 1), result)
@test in(Dict(:tablename => "influenza_cases", :status1 => "linked with criteria ID 5", :status2 => "linked with criteria ID 5", :nrecords => 2), result)
@test in(Dict(:tablename => "influenza_cases", :status1 => "linked with criteria ID 6", :status2 => "linked with criteria ID 6", :nrecords => 1), result)
@test in(Dict(:tablename => "influenza_cases", :status1 => "linked with criteria ID 7", :status2 => "linked with criteria ID 7", :nrecords => 1), result)


text:
#demo_snp_rw.jl: sNp (Touchstone) file tests
#-------------------------------------------------------------------------------

using CMDimCircuits
CMDimCircuits.@using_CData()

#Get a demo display:
include(CMDimCircuits.demoplotcfgscript)
#Normally use something like:
#CMDimData.@includepkg EasyPlotInspect

#==Constants
===============================================================================#
LBLAX_MAGDB = "Magnitude (dB)"
LBLAX_FREQ = "Frequency (Hz)"
color1 = cons(:a, line = set(color=:red, width=2))
color2 = cons(:a, line = set(color=:blue, width=2))

#==Input data
===============================================================================#
Z0 = 50 #Ohms (reference impedance)
ℓ = 5e-3 #Meters
Zc = 200 #Ohms (line impedance)
μ0 = 4pi*1e-7 #F/m
ϵ0 = 8.854e-12 #H/m
f = collect(1:.1:100)*1e9

#==Computations
===============================================================================#
len_mm = ℓ/1e-3
filepath_S = "lineS_$(len_mm)m.s2p"
filepath_Z = "lineZ_$(len_mm)m.s2p"

ω = 2pi*f
α = 0
β = ω*sqrt(μ0*ϵ0)
γ = α .+ im*β

#Convert γ type to DataF1 (function of 1 argument):
γ = DataF1(f, γ)

#ABCD matrix:
A = cosh(γ*ℓ)      ;B = sinh(γ*ℓ)*Zc;
C = sinh(γ*ℓ)/Zc   ;D = cosh(γ*ℓ);
T = Network(:ABCD, [A B; C D])

#Convert to S-paramters:
data = Network(:S, T)

#Save data to .s2p file:
EDAData.write_snp(filepath_S, data)
EDAData.write_snp(filepath_Z, Network(:Z, data))

#Re-load data:
data = EDAData.read_snp(filepath_S, numports=2)
(s11, s12, s21, s22) = mx2elem(data)
@show Symbol(NPType(data))

#==Generate plot
===============================================================================#
i𝛤 = 1; i𝛵 = 2; #Strip indices

plot = cons(:plot, nstrips=2,
	ystrip1 = set(axislabel=LBLAX_MAGDB, striplabel="Reflection Coefficients (𝛤)"),
	ystrip2 = set(axislabel=LBLAX_MAGDB, striplabel="Transmission Coefficients (𝛵)"),
	xaxis = set(label=LBLAX_FREQ),
)
push!(plot,
	cons(:wfrm, dB20(s11), color1, label="s11", strip=i𝛤),
	cons(:wfrm, dB20(s22), color2, label="s22", strip=i𝛤),
	cons(:wfrm, dB20(s12), color1, label="s12", strip=i𝛵),
	cons(:wfrm, dB20(s21), color2, label="s21", strip=i𝛵),
)

#Overlay result of reading/writing Z-parameters:
data = EDAData.read_snp(filepath_Z, numports=2)
@show s = Symbol(NPType(data))
(s11, s12, s21, s22) = mx2elem(Network(:S, data))
push!(plot,
	cons(:wfrm, dB20(s11), color1, label="s11 ($(s)P)", strip=i𝛤),
	cons(:wfrm, dB20(s22), color2, label="s22 ($(s)P)", strip=i𝛤),
	cons(:wfrm, dB20(s12), color1, label="s12 ($(s)P)", strip=i𝛵),
	cons(:wfrm, dB20(s21), color2, label="s21 ($(s)P)", strip=i𝛵),
)

pcoll = push!(cons(:plotcoll, title="EDAData Tests: sNp (Touchstone) Format"), plot)
	pcoll.ncolumns = 1

#==Display results as a plot
===============================================================================#
EasyPlot.displaygui(pcoll)

#==Return pcoll to user (call evalfile(...))
===============================================================================#
pcoll


text:
using FastFilter, Test, BenchmarkTest
import Base: Fix1

@testset "FastFilter.jl" begin
    ary = map(x -> rand(Bool) ? nothing : x, 1:10_000_000)

    pred(a) = !isnothing(a)
    tran(a) = pred(a) ? a : 0

    @test all(filter(pred, ary) .== fastfilter(pred, ary))
    @test all(replace(tran, ary) .== fastreplace(tran, ary))

    @btime filter($pred, $ary);
    @btime fastfilter($pred, $ary);
    @btime replace($tran, $ary);
    @btime fastreplace($tran, $ary);

    summer(f, ary) = sum(f(ary))

    @test summer(Fix1(filter, pred), ary) == summer(Fix1(fastfilter, pred), ary)
    @test summer(Fix1(replace, pred), ary) == summer(Fix1(fastreplace, pred), ary)

    @btime summer(Fix1(filter, pred), $(ary))
    @btime summer(Fix1(fastfilter, pred), $(ary))
    @btime summer(Fix1(replace, tran), $(ary))
    @btime summer(Fix1(fastreplace, tran), $(ary))

    # type inference failed
    @btime filter($ary) do a
        string(a) != "nothing"
    end
    @btime fastfilter($ary) do a
        string(a) != "nothing"
    end
end


text:
module Cards
# This module generates cards. Its contents will have to be changed to account
# for variations in the card template.

export mk_card

using Printf
using Dates
using Colors
using ColorTypes

# This function returns a function that converts linear indices to Cartesian (x
# and y) indices for the given array.
# Usage: map(ci_mapfunc(xs), collect(enumerate(xs)))
ci_mapfunc(xs) = ((li, x),) -> (CartesianIndices(xs)[li], x)

mk_gridcell((i, n)) = if i.I === (3,3) ""
  else @sprintf("%u %u gat moveto (%02u) show\n", i.I[2]-1, i.I[1]-1, n) end

mk_grid(ns) = map(mk_gridcell ∘ ci_mapfunc(ns), enumerate(ns))

handle_date(date::Nothing)        = ""
handle_date(date::Date)           = handle_date(Dates.format(date,
  dateformat"YYYY-mm-dd"))
handle_date(date::DateTime)       = handle_date(Date(date)) *
  handle_time(Time(date))
handle_date(date::AbstractString) = """
  3 -0.53 ga moveto
  (Date) show
  4 -0.53 ga moveto
  ($(date)) show
  """

handle_time(time::Nothing)        = ""
handle_time(time::Time)           = handle_time(Dates.format(time,
  dateformat"II:MM pp"))
handle_time(time::AbstractString) = """
  3 -0.64 ga moveto
  (Starting At) show
  4 -0.64 ga moveto
  ($(time)) show
  """

# The interface by which card details are passed to this function is subject
# to change.
# (I might make it a struct.)
mk_card(ns, colorant::String; kwargs...) =
  mk_card(ns,
          (convert(RGB, parse(Colorant, colorant)), colorant);
          kwargs...)

mk_card(ns, (col, col_name); num=0, round=0, date=nothing, draft=false) = """
  %!PS
  % Bingo card template
  % Written by Tucker R. Twomey
  /in {72 mul} def

  % Hardware margin
  /hwmargin  0.1389 in            def

  % Get paper dimensions
  currentpagedevice /PageSize get
  dup
  /width    exch 0 get def
  /height   exch 1 get def

  /DeviceCMYK setcolorspace
  % C M   Y   K
  0.0 0.0 0.0 1.0 setcolor

  % Function to draw box
  /rect
  {
    % Convert coord pairs to arrays to make stack ops easier
    4 array astore /coord exch def % Shame!
    % Stack is currently [ x2 y2 x1 y1 ] (top)
    % X1 and Y1
    %aload 4 2 roll moveto
    %dup aload pop moveto
    % X1          Y1
    coord 2 get   coord 3 get moveto
    % X1          Y2
    coord 2 get   coord 1 get lineto
    % X2          Y2
    coord 0 get   coord 1 get lineto
    % X2          Y1
    coord 0 get   coord 3 get lineto
    closepath
  } def

  /strike
  {
    % Convert coord pairs to arrays to make stack ops easier
    4 array astore /coord exch def % Shame!
    % Stack is currently [ x2 y2 x1 y1 ] (top)
    % X1 and Y1
    %aload 4 2 roll moveto
    %dup aload pop moveto
    % X1          Y1
    coord 2 get   coord 3 get moveto
    % X2          Y2
    coord 0 get   coord 1 get lineto stroke
    % X1          Y2
    coord 2 get   coord 1 get moveto
    % X2          Y1
    coord 0 get   coord 3 get lineto stroke
  } def

  %%% Grid drawing
  %%  Usage: [ <x1> <y1> <x2> <y2> <box-length> <box-height> ] grid-draw
  /grid
  {
    /args exch def
    %% Vertical lines
    % Start               Inc            End
    args 0 get            args 4 get     args 2 get
    {
      dup
      % Y
      args 1 get moveto
      args 3 get  lineto
      stroke
    } for

    %% Horizontal lines
    % Start               Inc            End
    args 1 get            args 5 get     args 3 get
    {
      dup
      % X
      args 0 get exch moveto
      args 2 get exch lineto
      stroke
    } for
  } def

  /gridaddr
  {
    /args exch def
    args 5 get mul args 1 get add
    exch
    args 4 get mul args 0 get add
    exch
  } def

  /gridx1 1 in def
  /gridy1 height 1 in sub 1.3 in 5.6 mul sub def
  /gridx2 7.5 in def
  /gridy2 height 1 in sub 1.3 in 0.6 mul sub def
  /gridbw 1.3 in def
  /gridbh 1.3 in def
  /ga  { [ gridx1 gridy1 gridx2 gridy2 gridbw gridbh ] gridaddr } def
  /gat { ga 25 add exch 05 add exch } def

  /DeviceRGB setcolorspace
  $(float(col.r)) $(float(col.g)) $(float(col.b)) setcolor

  /URWGothic-Book 20 selectfont
  0 -0.25 ga moveto
  (Supporting) show
  /URWGothic-Demi 45 selectfont
  0 -0.65 ga moveto
  (El Cerrito) show
  0 -1.05 ga moveto
  (Forensics) show

  $(
  if draft
    """
    gsave
    0.8 setgray % Faint gray
    /Helvetica 175 selectfont
    2 2 ga moveto
    58 rotate
    (Draft) show
    grestore
    gsave
    0.0 setgray % Black
    20 setlinewidth
    0 0 ga 2.25 -1.05 ga strike
    grestore
    """
  else
    ""
  end
  )

  %1 in 1 in width 1 in sub height 1 in sub rect stroke

  $(float(col.r)) $(float(col.g)) $(float(col.b)) setcolor
  /Helvetica-Bold 12 selectfont
  3 -0.20 ga moveto
  (Name) show
  %3.50 -0.42 ga 5 -0.33 ga rect stroke
  3.50 -0.20 ga moveto 5 -0.20 ga lineto stroke

  3 -0.31 ga moveto
  (Card) show
  4 -0.31 ga moveto
  ($(@sprintf("%05u", num))) show

  3 -0.42 ga moveto
  (Round) show
  4 -0.42 ga moveto
  ($(@sprintf("%02u (%s)", round, col_name))) show

  $(handle_date(date))

  gridx1 gridy1 gridx2 gridy2 rect stroke
  [ gridx1 gridy1 gridx2 gridy2 gridbw gridbh ] grid

  /URWGothic-Book 46 selectfont
  2 2 ga moveto (Free) show

  /URWGothic-Book 72 selectfont
  0.2 5 ga moveto
  (BINGO)
  [1.2 gridbw mul
   0.8 gridbw mul
   0.9 gridbw mul
   1.0 gridbw mul
   1.0 gridbw mul]
  xshow

  /URWGothic-Book 72 selectfont
  $(mk_grid(ns)...)

  showpage
  """

end

max_stars_repo_name: tuckertwo/ECbingo

max_issues_repo_name: tuckertwo/ECbingo

max_forks_repo_name: tuckertwo/ECbingo


text:
# ----------
# Dye Module, Experimental feature!
# ----------

"""
    mutable struct Dye

Represents the clock of a problem.
$(TYPEDFIELDS)
"""
mutable struct Dye{Bool,Function,Array,Tmp}
    "Dye Flag"
     dyeflag :: Bool
    "Dye Function"
stepforward! :: Function
    "Dye Array in complexed Form"
           ρ :: Array
    "Dye Sol and cache"
         tmp :: Tmp
end

show(io::IO, dye::Dye) =
     print(io, "Dye\n",
               "  ├────── dye flag: ", dye.dyeflag, '\n',
               "  ├────── function: ", dye.stepforward!, '\n',
               "  └─────────── dye: ρ - size :", size(dye.ρ));

mutable struct Tmp{Atrans}
  sol₀ :: Atrans
  sol₁ :: Atrans
  RHS₁ :: Atrans
  RHS₂ :: Atrans
  RHS₃ :: Atrans
  RHS₄ :: Atrans
end

function DyeContructer(::Dev, DyeFlag::Bool, grid::AbstractGrid) where Dev
  T = eltype(grid)

  if (DyeFlag == true )
  	@devzeros Dev T (grid.nx, grid.ny, grid.nz) ρ
  	@devzeros Dev Complex{T} (grid.nkr, grid.nl, grid.nm) sol₀ sol₁ RHS₁ RHS₂ RHS₃ RHS₄
  	tmp = Tmp(sol₀, sol₁, RHS₁, RHS₂, RHS₃, RHS₄);
  	dye = Dye(true, Dyestepforward!, ρ, tmp);
  else
  	ρ  = zeros(T,1);
  	tmp = Tmp(zeros(1),zeros(1),zeros(1),zeros(1),zeros(1),zeros(1));
  	dye = Dye(false, Dyestepforward!, ρ, tmp);
  end
  return dye
end

function Dyestepforward!(prob)
  dye, clock, vars, params, grid, clock = prob.dye, prob.clock, prob.vars, prob.params, prob.grid, prob.clock;
  ts = dye.tmp;
  RK4substeps!(ts.sol₀, clock, dye.tmp, vars, params, grid, clock.t, clock.dt)
  RK4update!(dye.tmp.sol₀, ts.RHS₁, ts.RHS₂, ts.RHS₃, ts.RHS₄, clock.dt)
  FourierFlows.dealias!(ts.sol₀,grid);
  DyeVarsUpdates!(ts.sol₀, dye.ρ, prob.grid);
end

function RK4update!(sol, RHS₁, RHS₂, RHS₃, RHS₄, dt)
  @. sol += dt*(RHS₁ / 6 + RHS₂ / 3  + RHS₃ / 3 + RHS₄ / 6)

  return nothing
end

function RK4substeps!(sol, clock, ts,  vars, params, grid, t, dt)
  # Substep 1
  DyeEqn!(ts.RHS₁, sol, t, clock, vars, params, grid)

  # Substep 2
  substepsol!(ts.sol₁, sol, ts.RHS₁, dt/2)
  DyeEqn!(ts.RHS₂, ts.sol₁, t+dt/2, clock, vars, params, grid)

  # Substep 3
  substepsol!(ts.sol₁, sol, ts.RHS₂, dt/2)
  DyeEqn!(ts.RHS₃, ts.sol₁, t+dt/2, clock, vars, params, grid)

  # Substep 4
  substepsol!(ts.sol₁, sol, ts.RHS₃, dt)
  DyeEqn!(ts.RHS₄, ts.sol₁, t+dt, clock, vars, params, grid)

  return nothing
end

function substepsol!(newsol, sol, RHS, dt)
  @. newsol = sol + dt*RHS

  return nothing
end

function DyeVarsUpdates!(sol,ρ, grid)
	ldiv!(ρ, grid.rfftplan, deepcopy(sol))

	return nothing
end

function DyeEqn!(N, sol, t, clock, vars, params, grid)
    # δ function
    δ(a::Int,b::Int) = ( a == b ? 1 : 0 );

    # To emulate the dye tracing,  We define dye density as ρ, and conserve the continuity equation
    # Such that, ∂ρ/∂t = ∇ ⋅ (ρ \vec{v})
    # So that, in spectral space,
    # ∂ρₖ/∂t = ∑ᵢ -im*kᵢ*(ρₖ vₖᵢ)

    # Initialization of rho
    @. N*=0;

    for (uᵢ,kᵢ) ∈ zip([vars.ux,vars.uy,vars.uz],[grid.kr,grid.l,grid.m])
         # Initialization
        @. vars.nonlin1 *= 0;
        ρuᵢ  = vars.nonlin1;
        ρuᵢh = vars.nonlinh1;

        # get back the updated sol in real space using fft
        ldiv!(ρuᵢ, grid.rfftplan, deepcopy(sol))

        # Pre-Calculation in Real Space
        @. ρuᵢ *= uᵢ;

        # Fourier transform
        mul!(ρuᵢh, grid.rfftplan, ρuᵢ);

        # Perform the actual calculation
        @. N += -im*kᵢ*ρuᵢh;
    end
    return nothing
end

max_stars_repo_name: MHDFlows/MHDFlows

max_issues_repo_name: MHDFlows/MHDFlows

max_forks_repo_name: MHDFlows/MHDFlows


text:
struct Xpack{T} end
Xpack(T) = Xpack{T}

function make_url(::Type{Xpack{:_sql}}, info::Esinfo)
	"$(info.url)/_xpack/sql/"
end

function make_url(::Type{Xpack{:_translate}}, info::Esinfo)
	"$(info.url)/_xpack/sql/translate"
end

function make_url(::Type{Xpack{:_sql_close}}, info::Esinfo)
	"$(info.url)/_xpack/sql/close"
end

@genfunction "POST" xpack_sql Xpack{:_sql} 1
@genfunction "POST" xpack_sql_close Xpack{:_sql_close} 1
@genfunction "POST" xpack_translate Xpack{:_translate} 1

max_stars_repo_name: frankill/ES

max_issues_repo_name: frankill/ES

max_forks_repo_name: frankill/ES


text:
using Pkg
Pkg.add("JuMP")
Pkg.add("StructArrays")
Pkg.add("DataFrames")
Pkg.add("XLSX")
Pkg.add("StatsBase")
Pkg.add("Query")
Pkg.add("Gadfly")
Pkg.add("Cairo")
Pkg.add("Fontconfig")
Pkg.add("CSV")
Pkg.add("JLD")

# these packages need their software installed at the following locations
ENV["CPLEX_STUDIO_BINARIES"] = "/opt/ibm/ILOG/CPLEX_Studio1210/cplex/bin/x86-64_linux/"
Pkg.add("CPLEX")
ENV["GUROBI_HOME"] = "/opt/gurobi911/linux64/"
Pkg.add("Gurobi")

max_stars_repo_name: Dibillilia/JuliaOptHeuristics

max_issues_repo_name: Dibillilia/JuliaOptHeuristics

max_forks_repo_name: Dibillilia/JuliaOptHeuristics


text:
struct Object
    surface::Surface
    frontMaterial::Union{Nothing, Material}
    backMaterial::Union{Nothing, Material}

    function Object(surface::Surface, material::Material)
        new(surface, material, nothing)
    end
end

function hit(obj::Object, ray::Ray, range::NamedTuple{(:min, :max), Tuple{T, T}} where T <: AbstractFloat)::Tuple{Bool, AbstractFloat, Union{Nothing, HitPoint}, Union{Nothing, Material}}
    succeeded::Bool, t::AbstractFloat, hp::Union{Nothing, HitPoint} = hit(obj.surface, ray, range)
    if !succeeded
        return false, NaN, nothing, nothing
    end
    local m::Union{Nothing, Material}
    if hp.front
        m = obj.frontMaterial
    else
        m = obj.backMaterial
    end
    m != nothing, t, hp, m
end

max_stars_repo_name: goooxu/JuliaRayTracing

max_issues_repo_name: goooxu/JuliaRayTracing

max_forks_repo_name: goooxu/JuliaRayTracing


text:
@testset "colons, mapslices etc" begin

    M = randn(3,4)

    @cast C1[i][j] := M[i,j]
    @cast C2[i][j] |= M[i,j]
    @cast C3[i] := M[i,:]
    @cast C4[i] |= M[i,:4] assert

    @test first(C1) isa SubArray
    @test first(C2) isa Array
    @test first(C3) isa SubArray
    @test first(C4) isa Array

    @cast R1[j]{i} := M[i,j]     # size from M
    @cast R2[j]{i:3} := M[i,j]   # size inside
    @cast R3[j]{i} := M[i,j] i:3 # size outside
    @cast R4[j] := M{:3,j}       # size on the colon
    @cast R5[j] := M{:,j}        # size from M

    @test first(R1) isa SArray
    @test first(R2) isa SArray
    @test first(R3) isa SArray
    @test first(R4) isa SArray
    @test first(R5) isa SArray

    γ = 3
    @cast R6[j]{i:γ} |= M[i,j] # test both |= and γ
    @cast R7[j] |= M{:γ,j} assert

    @test R5 isa Base.ReinterpretArray
    @test R6 isa Array
    @test R7 isa Array

    S1 = vec(sum(M, dims=1))
    @cast S2[j] := sum(M[:,j])
    @cast S3[j] := sum(M{:3,j})
    @reduce S4[j] := sum(i) M[i,j]
    @test S1 == S2 == S3 == S4

    @cast M1[i,j] := copy(M[:,j])[i]
    @cast M2[i,j] := copy(M{:3,j})[i]
    @cast M3[i,j] := copy(M{:,j})[i]  # size from M
    @cast M4[i,j] := copy(M{:,j}){i}  # static glue
    @test M == M1 == M2 == M3 == M4

    fun(x) = (sum=sum(x), same=x, one=1)

    @cast M5[i,j] := fun(M[:,j]).same[i]
    @cast M6[i,j] := fun(M{:3,j}).same{i}
    @cast M7[i,j] := fun(M[:,j] * 1).same[i]
    @test M == M5 == M6 == M7

    @cast S5[j] := fun(M[:,j]).sum
    @test S1 == S5

    @cast D1[i,j] := 2(C1[i][j])
    @cast D2[i,j] := copy(C1[i] + 1 * C1[i])[j]
    @test 2M == D1 == D2

end
@testset "fields & sizes" begin

    B = [ (scal=i, vect=[1,2,3,4]) for i=1:3 ]
    @cast A[j,i] := B[i].vect[j]  i:3, j:4
    @test A[4,1] == 4

    @cast C[i] := B[i].scal
    @test C == 1:3

end
@testset "int/bool indexing" begin

    ind = [1,2,1,2]
    M = rand(1:99, 3,4)

    N1 = M[1,ind]
    @cast N2[i] := M[1,ind[i]]
    @test N1 == N2

    N3 = M[ind,ind]
    @cast N4[i,j] := M[ind[i],ind[j]]
    @test N3 == N4

    bool = ind .== 1
    B1 = M[:,bool]
    @cast B2[i,j] := M[i,bool[j]]
    @test B1 == B2

    V1 = M[CartesianIndex.(ind, ind)]
    # @cast V2[i] := M[ind[i],ind[i]]

end
@testset "indexing with a range" begin

    M = rand(1:99, 3,4)

    # ranges are treated first like constants, then like colons:
    @cast C1[i,j] := identity(M[i, 1:2])[j]
    R1 = 1:2
    @cast C2[i,j] := M[i, R1[j]] # explicit inner indexing
    @test C1 == C2 == M[:, 1:2]

    β = 2
    @cast C3[i,j] := identity(M[i, 1:β])[j]
    @test C3 == C1

end
@testset "in-place += and *=" begin

    A = ones(2,3)
    B = rand(3,2)

    @cast A[i,j] += B[j,i]
    @test A == 1 .+ B'

    C = 2 .* ones(3)
    D = rand(3)

    @cast C[i] *= D[i] + 1
    @test C == 2 .* (D .+ 1)

    M = ones(3,3)

    @cast M[i,i] -= D[i] + 5
    @test diag(M) == 1 .- (D .+ 5)

end
@testset "fixing inner indices" begin

    imgs = [ rand(28,28,1,4) for i=1:2 ];
    @cast Z[i,j,k,l] := imgs[l][i,j,_,k];

    @test size(Z) == (28,28,4,2)
    @test Z[2,3,4,1] == imgs[1][2,3,1,4]

    @cast out[l][i,j,_,k] := Z[i,j,k,l];
    @test out == imgs

    uno = 1 # try again with $constant
    @cast Z[i,j,k,l] := imgs[l][i,j,$uno,k];

    @test size(Z) == (28,28,4,2)
    @test Z[2,3,4,1] == imgs[1][2,3,1,4]

    @cast out[l][i,j,1,k] := Z[i,j,k,l];
    @test out == imgs

end
@testset "scalars & 0-arrays" begin

    B = rand(2,3)
    @reduce S := sum(i,j) sqrt(B[j,i])
    S′ = @reduce sum(i,j) sqrt(B[j,i])
    @test S isa Number
    @test S ≈ S′ ≈ sum(sqrt, B)

    @reduce T[] := sum(i,j) sqrt(B[j,i])
    @test T isa Array{Float64,0}
    @test T[] ≈ sum(sqrt, B)

    @reduce T1[_] := sum(i,j) sqrt(B[j,i])
    @test T1 isa Array{Float64,1}
    @test T1[1] ≈ sum(sqrt, B)

    V = rand(4)
    @cast TV[i] := V[i] + T[]
    @test TV ≈ V .+ S

    U =  similar(T)
    @reduce U[] = sum(i,j) B[j,i]/2 # inplace
    @test U[] ≈ sum(B)/2

    @cast V[] := 2 * T[]
    @test_broken V isa Array{Float64,0}

end
@testset "tupple broadcasting" begin

    V = 1:4
    @cast T[i,j] := (V[i], V[j]^2 + 10^3, 99)
    @test T[3,4] == (3,1016,99)

    @cast N[i,j] := (a = V[i], b = V[j]^2 + 10^3, c = 99)
    @test N[3,4] == (a=3, b=1016, c=99)

end
@testset "arrays of functions" begin

    funs = [sin, cos, tan]
    vals = 1:3
    @cast M[i,j] := (funs[i])(vals[j])
    @cast V[i] := (funs[i])(vals[i])
    @test M[1,2] == sin(2)
    @test V[2] == cos(2)

    op = [+, -, *, /]
    x = 1:5
    y = 1:6
    @cast T[i,j,k] := op[i](x[j], y[k]) + 10
    @test T[2,3,4] == 3-4 + 10
    @test T[4,5,6] == 5/6 + 10

    @cast A[i] := op[i](π,ℯ)
    @test A[end] == π / ℯ

end
@testset "A * B * C" begin

    A = rand(2,2)
    B = rand(2,2)
    C = rand(2,2)

    @matmul Z[i,l] := sum(j,k) A[i,j] * B[j,k] * C[k,l]
    @test_broken Z == A * B * C
    @test vec(Z) == vec(A * B * C)

end
@testset "string macros" begin

    A = rand(2,10); B = rand(10,10);

    reduce" C_ii := sum_k A_1k * log(A_2i * B_ik) "
    @test C isa Diagonal

end
@testset "from todo list" begin

    list = [ i .* ones(2,2,1) for i=1:8 ];
    @cast mat[x\i, y\j] := Int(list[i\j][x,y,1])  i:2
    @cast mat2[x\i, y\j] := Int(list[i\j][x,y,1])  i:2, lazy # crazy type
    @test mat[3,5] == 6
    @test mat == mat2

    @cast C1[i,i',k] := (1:4)[i⊗i′⊗k] + im  (i:2, i′:2)  # two tensor signs
    @cast C2[i,i',k] := (1:4)[i⊗i'⊗k] + im  (i:2, i′:2)  # more primes
    @cast C3[i,i',k] := (1:4)[i⊗i'⊗k] + im  (i:2, i':2)
    @test C1[1,2,1] == C2[1,2,1] == C3[1,2,1] == 3+1im

    @cast C4[i,i'⊗k] := (1:4)[i⊗i′⊗k] + im  (i:2, i':2)
    @test size(C4) == (2,2)
    @test C4[2,1] == 2+im

    ∇λ = ones(3); topd = rand(3); logind=3; λ=rand(3); d=2
    @cast ∇λ[c] = ∇λ[c] - λ[$d] * exp(topd[c]) / logind # all const
    @test ∇λ[1] == 1 - λ[2] * exp(topd[1]) / 3

    W = rand(3); X = [rand(3) for _=1:3, _=1:2];
    @cast Z[i,j] := W[i] * exp(X[1,1][i] - X[2,2][j]) # all const
    @test Z == @. W * exp(X[1,1] - X[2,2]')

    mat = rand(3,4)
    outer(v::AbstractVector) = v * v'
    @cast out[i⊗i',j] := outer(mat[:,j])[i,i′] i:3, i':3
    @test size(out) == (9,4)

    yy = rand(3,4)
    t = 4
    xx = zeros(9)
    @cast xx[μ⊗ν] = yy[μ,$t] * yy[ν,$t]
    @test xx ≈ vec(yy[:,4] * yy[:,4]')

    A = [collect((1:4) .+ 10i) for i=1:5]
    @cast B[i,j] := getindex(A[i], (2:3)[j])
    @test B == getindex.(A, (2:3)')

    @cast C[i,j] |= mat[j,i] # now identity.() not collect()
    C[2,1] = 99
    @test mat[1,2] != 99

end
@static if VERSION >= v"1.1"
@testset "@avx" begin

    using LoopVectorization

    A = rand(4,5)
    @test exp.(A) ≈ @cast B[i,j] := exp(A[i,j]) avx

    @test exp.(A') ≈ @cast B[i,j] := exp(A[j,i]) avx
    @test exp.(A.+1) ≈ @cast B[i,j] := exp(A[i,j]+1) avx

end
end
@testset "OffsetArrays" begin

    # https://github.com/mcabbott/TensorCast.jl/issues/11
    A = OffsetArray(rand(3, 3, 5, 5), 1:3, 1:3, -2:2, -2:2);
    # A = OffsetArray(rand(3, 3, 5, 5), 1:3, 11:13, -2:2, -2:2);
    B = OffsetArray(rand(3,    5)   , 1:3,        -2:2)
    δ1 = Matrix(I,3,3)
    # δ1 = OffsetArray(Matrix(I,3,3), 11:13, 11:13)
    δ2 = OffsetArray(Matrix(I,5,5), -2:2, -2:2)
    Ap, Bp, δ1p, δ2p = parent(A), parent(B), δ1, parent(δ2);
    # Ap, Bp, δ1p, δ2p = parent(A), parent(B), parent(δ1), parent(δ2);

    @cast Cp[a, b, c, d] := Ap[a, b, c, d] + Bp[a, c] * δ1p[a, b] * δ2p[c, d]
    @cast C[ a, b, c, d] := A[ a, b, c, d] + B[ a, c] * δ1[ a, b] * δ2[ c, d]
    @test C.parent == Cp

end
@testset "Strided" begin

    bcde = rand(2,3,4,5);
    bcde2 = similar(bcde);

    @cast oo[d,e,b,c] |= bcde[b,c,d,e] strided;
    @test size(oo) == (4,5,2,3)
    @test oo[1,3,1,3] == bcde[1,3,1,3]

    oo[1] = 99
    @test bcde[1] != 99 # copy not a view

    @cast g[(y,e),(b,c),x] := bcde[b,c,(x,y),e]  x:1, strided;
    @test size(g) == (20, 6, 1)
    @cast bcde2[b,c,(x,y),e] = g[(y,e),(b,c),x]  strided;
    @test all(bcde2 .== bcde) # used to say "fails when using Strided"

    # https://github.com/mcabbott/TensorCast.jl/issues/2
    A = rand(3,3); B = rand(3,5);
    @reduce C[i, j] := sum(l) A[i, l] * B[l, j] strided
    @test C ≈ A * B
    @reduce D[i, j] |= sum(l) A[i, l] * B[l, j] strided
    @test D isa Array

end
@testset "parse-time errors" begin

    using TensorCast: MacroError, _macro, CallInfo

    @test_throws MacroError _macro(:(  A[i,j,i] := B[i,j]  )) # repeated
    @test_throws MacroError _macro(:(  A[i,j] := B[i,j,-i]  ))
    @test_throws MacroError _macro(:(  A[k] := sum(k)  ),:(  B[k]  ), call=CallInfo(:reduce))

    @test_throws MacroError _macro(:(  A[i,j] := B[k]  )) # "can't find index k on the left"
    @test_throws MacroError _macro(:(  A[i] := sum(k)  ),:(  (B[j]+ C[j])[i]  ), call=CallInfo(:reduce))

    @test_throws MacroError _macro(:(  A[i] := (B[i], c=3)  ))

    @test_throws MacroError _macro(:( A[h, w\n] := B[n][2,-h,w] )) # can't reverse inner
    @test_throws MacroError _macro(:( A[i,j] := B[i][~j] )) # can't shuffle inner
    @test_throws MacroError _macro(:( A[i,-j] := B[i,j] )) # can't reverse left

end
@testset "run-time errors" begin

    B = [ (scal=i, vect=[1,2,3,4]) for i=1:3 ]
    @test_throws DimensionMismatch @cast A[j,i] := B[i].vect[j]  i:99, j:4 # wrong size
    @test_throws DimensionMismatch @cast A[j,i] := B[i].vect[j]  i:3, j:99

    M = randn(3,4)
    fun(x) = (sum=sum(x), same=x, one=1)
    @test_throws DimensionMismatch  @cast M5[i,j] := fun(M[:,j]).same[i]  i:3, j:99
    @test_throws DimensionMismatch  @cast M5[i,j] := fun(M[:99,j]).same[i]  j:4
    # @test_throws DimensionMismatch  @cast M5[i,j] := fun(M[:,j]).same[i]  i:99, j:4 # TODO make this check canonical length?

end


text:
using Base.Test
import Base: REPL, LineEdit
import Q.KdbMode: install_kdb_mode

mutable struct FakeTerminal <: Base.Terminals.UnixTerminal
    in_stream::Base.IO
    out_stream::Base.IO
    err_stream::Base.IO
    hascolor::Bool
    raw::Bool
    size::Tuple{Int,Int}
    FakeTerminal(stdin, stdout, stderr, hascolor=false) =
        new(stdin, stdout, stderr, hascolor, false, (24, 80))
end

Base.Terminals.hascolor(t::FakeTerminal) = t.hascolor
Base.Terminals.raw!(t::FakeTerminal, raw::Bool) = t.raw = raw
Base.Terminals.size(t::FakeTerminal) = t.size

function fake_repl()
    p = [Pipe() for _ in 1:3]
    map(Base.link_pipe, p)
    term = FakeTerminal(p[1].out, p[2].in, p[3].in)
    repl = Base.REPL.LineEditREPL(term)
    repl.history_file = false
    p[1].in, p[2].out, p[3].out, repl
end

stdin, stdout, stderr, repl = fake_repl()

repltask = @async begin
    Base.REPL.run_repl(repl)
end
@testset "REPL" begin
    write(stdin, "42\r")
    @test (out = readuntil(stdout, "\n\n"); contains(out, "42\n"))
    @test install_kdb_mode(repl) &&
        (contains(readline(stderr), "q)") ||
         contains(readline(stderr), "q)"))
    write(stdin, "\\")
    readuntil(stdout, "q)")
    write(stdin, "til 10\r")
    @test (out = readline(stdout); contains(out, "til 10"))
    @test (out = readline(stdout); contains(out, "0 1 2 3"))
end

# Close REPL ^D
write(stdin, '\x04')
wait(repltask)


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

using .Markdown
import .Markdown: MD, Paragraph, Header, Italic, Bold, plain, term, html, Table, Code
import Base: writemime

# Basics
# Equality is checked by making sure the HTML output is
# the same – the structure itself may be different.

@test md"foo" == MD(Paragraph("foo"))
@test md"foo *bar* baz" == MD(Paragraph(["foo ", Italic("bar"), " baz"]))

@test md"#no title" == MD(Paragraph(["#no title"]))
@test md"# title" == MD(Header{1}("title"))
@test md"""
  #
  empty
  """ == MD(Header{1}(""), Paragraph("empty"))
@test md"## section" == MD(Header{2}("section"))
@test md"# title *foo* `bar` **baz**" ==
    MD(Header{1}(["title ", Italic("foo")," ",Code("bar")," ",Bold("baz")]))
@test md"""
h1
===""" == md"# h1"
@test md"""
h2
   ---""" == md"## h2"

@test md"**foo *bar* baz**" == MD(Paragraph(Bold(["foo ", Italic("bar"), " baz"])))
@test md"*foo **bar** baz*" == MD(Paragraph(Italic(["foo ", Bold("bar"), " baz"])))

@test md"""```julia
foo
```
""" == MD(Code("julia", "foo"))
@test md"``code```more code``" == MD(Any[Paragraph(Any[Code("","code```more code")])])
@test md"``code``````more code``" == MD(Any[Paragraph(Any[Code("","code``````more code")])])

@test md"""
* one
* two

1. pirate
2. ninja
3. zombie""" == Markdown.MD([Markdown.List(["one", "two"]),
                             Markdown.List(["pirate", "ninja", "zombie"], true)])

@test md"Foo [bar]" == MD(Paragraph("Foo [bar]"))
@test md"Foo [bar](baz)" != MD(Paragraph("Foo [bar](baz)"))
@test md"Foo \[bar](baz)" == MD(Paragraph("Foo [bar](baz)"))

# Basic plain (markdown) output

@test md"foo" |> plain == "foo\n"
@test md"foo *bar* baz" |> plain == "foo *bar* baz\n"
@test md"# title" |> plain == "# title\n"
@test md"## section" |> plain == "## section\n"
@test md"## section `foo`" |> plain == "## section `foo`\n"
@test md"""Hello

---
World""" |> plain == "Hello\n\n–––\n\nWorld\n"

# HTML output

@test md"foo *bar* baz" |> html == "
foo bar baz
\n"
@test md"something ***" |> html == "
something ***
\n"
@test md"# h1## " |> html == "
h1##
\n"
@test md"## h2 ### " |> html == "
h2
\n"
@test md"###### h6" |> html == "
h6
\n"
@test md"####### h7" |> html == "
####### h7
\n"
@test md"   >" |> html == "
\n
\n"
@test md"1. Hello" |> html == "
\n
Hello
\n
\n"
@test md"* World" |> html == "
\n
World
\n
\n"
@test md"# title *blah*" |> html == "
title blah
\n"
@test md"## title *blah*" |> html == "
title blah
\n"
@test md"""Hello

---
World""" |> html == "
Hello
\n
\n
World
\n"
@test md"`escape</code>`" |> html == "
escape</code>
\n"

@test md"""
    code1

    code2
""" |> html == "
code1\n\ncode2
\n" # single code block

# @test md"""
# - Foo
#  ---
# - Bar""" |> html == "
\n
Foo
\n
\n
\n
\n
Bar
\n
\n"
@test md"""
h1
===
h2
---
not
== =""" |> html == "
h1
\n
h2
\n
not == =
\n"

# Latex output
book = md"""
# Title

Some discussion

> A quote

## Section *important*

Some **bolded**

- list1
- list2
"""
@test latex(book) == "\\section{Title}\nSome discussion\n\\begin{quote}\nA quote\n\\end{quote}\n\\subsection{Section \\emph{important}}\nSome \\textbf{bolded}\n\\begin{itemize}\n\\item list1\n\\item list2\n\\end{itemize}\n"

# Interpolation / Custom types

type Reference
    ref
end

ref(x) = Reference(x)

# if Base.USE_GPL_LIBS

ref(fft)

writemime(io::IO, m::MIME"text/plain", r::Reference) =
    print(io, "$(r.ref) (see Julia docs)")

fft_ref = md"Behaves like $(ref(fft))"
@test plain(fft_ref) == "Behaves like fft (see Julia docs)\n"
@test html(fft_ref) == "
Behaves like fft (see Julia docs)
\n"

writemime(io::IO, m::MIME"text/html", r::Reference) =
    Markdown.withtag(io, :a, Markdown.Pair(:href, "test")) do
        Markdown.htmlesc(io, Markdown.plaininline(r))
    end
@test html(fft_ref) == "
Behaves like fft (see Julia docs)
\n"

# end # USE_GPL_LIBS

@test md"""
````julia
foo()
````""" == md"""
```julia
foo()
```"""

# GH tables
@test md"""
    a  | b
    ---|---
    1  | 2""" == MD(Table(Any[["a","b"],
                              ["1","2"]], [:r, :r]))

@test md"""
    | a  |  b | c |
    | :-- | --: | --- |
    | d`gh`hg | hgh**jhj**ge | f |""" == MD(Table(Any[["a","b","c"],
                                                      Any[["d",Code("gh"),"hg"],
                                                          ["hgh",Bold("jhj"),"ge"],
                                                          "f"]],
                                                  [:l, :r, :r]))
@test md"""
no|table
no error
""" == MD([Paragraph(Any["no|table no error"])])

t = """a   |   b
:-- | --:
1   |   2
"""
@test plain(Markdown.parse(t)) == t