text:
## Euclid 1

"""
    $TYPEDEF

Continuous Euclidean variable of dimension `N` representing a Position in cartesian space.
"""
struct Position{N} <: InferenceVariable end

Position(N::Int) = Position{N}()

# not sure if these overloads are necessary since DFG 775?
DFG.getManifold(::InstanceType{Position{N}}) where N = TranslationGroup(N)
DFG.getDimension(val::InstanceType{Position{N}}) where N = manifold_dimension(getManifold(val))
DFG.getPointType(::Type{Position{N}}) where N = Vector{Float64}
DFG.getPointIdentity(M_::Type{Position{N}}) where N = zeros(N) # identity_element(getManifold(M_), zeros(N))

Base.convert(::Type{<:ManifoldsBase.AbstractManifold}, ::InstanceType{Position{N}}) where N = TranslationGroup(N)

#

"""
$(TYPEDEF)

Most basic continuous scalar variable in a `::DFG.AbstractDFG` object.

Alias of `Position{1}`
"""
const ContinuousScalar = Position{1}
const ContinuousEuclid{N} = Position{N}

const Position1 = Position{1}
const Position2 = Position{2}
const Position3 = Position{3}
const Position4 = Position{4}

## Circular

"""
$(TYPEDEF)

Circular is a `Manifolds.Circle{ℝ}` mechanization of one rotation, with `theta in [-pi,pi)`.
"""
@defVariable Circular RealCircleGroup() [0.0;]

#


text:
using PackageCompiler
using Pkg

Pkg.activate(".")
create_sysimage(:MasterEquationSolve,sysimage_path="mesolve_sysimage.so",precompile_execution_file="precompile.jl")

ext: jl


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

using Markdown
using InteractiveUtils

# ╔═╡ b70ee5f2-3c05-11ec-37db-1704f9c112b5
begin
	using ValidatedNumerics
	using StaticArrays
	using DelimitedFiles
	using IntervalArithmetic
	using IntervalRootFinding
	using LinearAlgebra
	using PyCall
	np = pyimport("numpy")
	cel = pyimport("celluloid")
	using PyPlot
end

# ╔═╡ f7a8506d-2e53-4e72-aaa8-b4133a6520cf
md"# Interval arithmetic"

# ╔═╡ 5397e0bd-378c-42fe-a23e-d5eca0b5a010
md"## Self positive feedback loop"

# ╔═╡ 17757d95-5205-4620-b0ad-c6dd844c13cc
begin
	function fself(x,c1,c2)
	    fself1 = c1*x^2/(c2+x^2) - x
	    return fself1
	end

	Sself = [];
	Sselfx = [];
	Sselfy = [];

	for c1 in 0.1:1:20.0
	    for c2 in 0.1:1:20.0
	        X0 = (0.0..c1)
	        rts = roots(xx -> fself(xx, c1, c2), X0, Krawczyk, 0.000001)
	        nsol = 0
	        for rt in rts
	            if rt.status == :unique
	                nsol += 1
	            end
	        end
	        if nsol >= 2
	            append!(Sself,[[c1,c2]])
	            append!(Sselfx,c1)
	            append!(Sselfy,c2)
	        end

	    end
	end
end

# ╔═╡ c92b4d43-33a1-4ca1-8403-e2900ca6d2be
begin
	tself = np.arange(0.0, 10.0, 0.2)
	uself = tself.^2/4
end;

# ╔═╡ 957a1119-3ffa-4f29-aae6-8d4991e0ebfc
begin
	plt.figure(figsize=(7,7))
	plt.scatter(Sselfx,Sselfy,label = "Bistable configuration",color="#D24D3E",s=25)
	plt.plot(tself,uself,label = "Critical curve",color="#781C81")
	plt.legend()
	plt.gca().set_aspect("equal")
	plt.xlabel("c (a.u.)")
	plt.ylabel("K (a.u.)")
	plt.xlim(0,20)
	plt.ylim(0,20)
	plt.show()
end

# ╔═╡ 64334315-5d05-4032-9ad3-8130d66fe008
md"## Double positive feedback loop"

# ╔═╡ 0b88f3f9-8e5f-479a-9c06-ee0aebfd200e
begin
	function fpos(x,c1,c2)
	    fpos1 = c1*x[2]^2/(1+x[2]^2) - x[1]
	    fpos2 = c2*x[1]^2/(1+x[1]^2) - x[2]
	    return SVector(fpos1, fpos2)
	end

	Spos = [];
	Sposx = [];
	Sposy = [];

	for c1 in 0.0:0.4:10.0
	    for c2 in 0.0:0.4:10.0

	        X0 = IntervalBox(0.0..c1, 0.0..c2)
	        rts = roots(xx -> fpos(xx, c1, c2), X0, Krawczyk, 0.000001)
	        nsol = 0
	        for rt in rts
	            if rt.status == :unique
	                nsol += 1
	            end
	        end
	        if nsol >= 2
	            append!(Spos,[[c1,c2]])
	            append!(Sposx,c1)
	            append!(Sposy,c2)
	        end

	    end
	end
end

# ╔═╡ 20a1adeb-9286-45d7-8f27-b5a41d6052d9
begin
	tpos = np.arange(0.0, 1.6, 0.01)
	apos = (4 .*tpos) ./(3 .- tpos.^2)
	bpos = (((4 ./(1 .+tpos.^2)).-1).+0im).^ (1/2) .* (1 .+ 1 ./tpos.^2)
end;

# ╔═╡ c0d8791e-190e-4799-9ec9-224fcd7bfcfb
begin
	plt.figure(figsize=(7,7))
	plt.scatter(Sposx,Sposy,label = "Bistable configuration",color="#D24D3E",s=20)
	plt.scatter(0,0,color="#3D52A1",s=50)
	plt.plot(bpos,apos,label = "Critical curve",color="#781C81")
	plt.legend()
	plt.gca().set_aspect("equal")
	plt.xlabel("c (a.u.)")
	plt.ylabel("K (a.u.)")
	plt.xlim(-0.5,10)
	plt.ylim(-0.5,10)
	plt.show()
end

# ╔═╡ 6ee9dd70-2e76-4fb3-8972-14edbe3f4103
md"## Double negative feedback loop"

# ╔═╡ 6218e2ed-902b-4a3c-9f19-8fced206df1c
begin
	function fneg(x,c1,c2)
	    fneg1 = c1/(1+x[2]^2) - x[1]
	    fneg2 = c2/(1+x[1]^2) - x[2]
	    return SVector(fneg1, fneg2)
	end

	Sneg = [];
	Snegx = [];
	Snegy = [];

	for c1 in 0.1:0.4:10.0
	    for c2 in 0.1:0.4:10.0

	        X0 = IntervalBox(0.0..c1, 0.0..c2)
	        #rts = roots(xx -> f(xx, c1, c2), X0,0.000001)
	        rts = roots(xx -> fneg(xx, c1, c2), X0, Krawczyk, 0.000001)
	        nsol = 0
	        for rt in rts
	            if rt.status == :unique
	                nsol += 1
	            end
	        end
	        if nsol >= 2
	            append!(Sneg,[[c1,c2]])
	            append!(Snegx,c1)
	            append!(Snegy,c2)
	        end

	    end
	end
end

# ╔═╡ f5f2d5b5-effb-46ac-ac23-a10f88f00b46
begin
	tneg = np.arange(0.6, 5.0, 0.01)
	aneg = (4 .*tneg.^3) ./(3 .*tneg.^2 .- 1)
	bneg = ((1 .+ tneg.^2)./(3 .*tneg.^2 .- 1).+0im).^ (1/2) .* (1 .+ tneg.^2)
end;

# ╔═╡ 954f9963-1249-43b4-800e-6d46604a27b0
begin
	plt.figure(figsize=(7,7))
	plt.scatter(Snegx,Snegy,label = "Bistable configuration",color="#D24D3E",s=20)
	plt.scatter(0,0,color="#3D52A1",s=50)
	plt.plot(aneg,bneg,label = "Critical curve",color="#781C81")
	plt.legend()
	plt.gca().set_aspect("equal")
	plt.xlabel("c (a.u.)")
	plt.ylabel("K (a.u.)")
	plt.xlim(-0.5,10)
	plt.ylim(-0.5,10)
	plt.show()
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
IntervalRootFinding = "d2bf35a9-74e0-55ec-b149-d360ff49b807"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
ValidatedNumerics = "d621b6e3-7715-5857-9c6f-c67000ef6083"

[compat]
IntervalArithmetic = "~0.17.8"
IntervalRootFinding = "~0.5.7"
PyCall = "~1.92.5"
PyPlot = "~2.10.0"
StaticArrays = "~1.2.13"
ValidatedNumerics = "~0.11.0"
"""

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

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

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

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

[[CRlibm]]
deps = ["Libdl"]
git-tree-sha1 = "9d1c22cff9c04207f336b8e64840d0bd40d86e0e"
uuid = "96374032-68de-5a5b-8d9e-752f78720389"
version = "0.8.0"

[[ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "3533f5a691e60601fe60c90d8bc47a27aa2907ec"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.11.0"

[[ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

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

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

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

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

[[Conda]]
deps = ["Downloads", "JSON", "VersionParsing"]
git-tree-sha1 = "6cdc8832ba11c7695f494c9d9a1c31e90959ce0f"
uuid = "8f4d0f93-b110-5947-807f-2305c1781a2d"
version = "1.6.0"

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

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

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

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

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

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

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

[[ErrorfreeArithmetic]]
git-tree-sha1 = "d6863c556f1142a061532e79f611aa46be201686"
uuid = "90fa49ef-747e-5e6f-a989-263ba693cf1a"
version = "0.5.2"

[[ExprTools]]
git-tree-sha1 = "b7e3d17636b348f005f11040025ae8c6f645fe92"
uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04"
version = "0.1.6"

[[FastRounding]]
deps = ["ErrorfreeArithmetic", "Test"]
git-tree-sha1 = "224175e213fd4fe112db3eea05d66b308dc2bf6b"
uuid = "fa42c844-2597-5d31-933b-ebd51ab2693f"
version = "0.2.0"

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

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

[[InlineStrings]]
deps = ["Parsers"]
git-tree-sha1 = "19cb49649f8c41de7fea32d089d37de917b553da"
uuid = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48"
version = "1.0.1"

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

[[IntervalArithmetic]]
deps = ["CRlibm", "FastRounding", "LinearAlgebra", "Markdown", "Random", "RecipesBase", "RoundingEmulator", "SetRounding", "StaticArrays"]
git-tree-sha1 = "00cce14aeb4b256f2f57caf3f3b9354c27d93259"
uuid = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
version = "0.17.8"

[[IntervalConstraintProgramming]]
deps = ["IntervalArithmetic", "IntervalContractors", "IntervalRootFinding", "MacroTools", "Test"]
git-tree-sha1 = "e5234b7d91eaaa80305cd8e74a5f2231d4d0e05e"
uuid = "138f1668-1576-5ad7-91b9-7425abbf3153"
version = "0.9.1"

[[IntervalContractors]]
deps = ["IntervalArithmetic"]
git-tree-sha1 = "02e049c761e7fd911566c6b0eff9deffe7de876b"
uuid = "15111844-de3b-5229-b4ba-526f2f385dc9"
version = "0.4.5"

[[IntervalOptimisation]]
deps = ["IntervalArithmetic"]
git-tree-sha1 = "7928f5018b5f54adb6780a43b975880633ead850"
uuid = "c7c68f13-a4a2-5b9a-b424-07d005f8d9d2"
version = "0.4.5"

[[IntervalRootFinding]]
deps = ["ForwardDiff", "IntervalArithmetic", "LinearAlgebra", "Polynomials", "StaticArrays"]
git-tree-sha1 = "2b3c1cbe892ceb2f3936fd734c1c98ea97dfb18e"
uuid = "d2bf35a9-74e0-55ec-b149-d360ff49b807"
version = "0.5.7"

[[Intervals]]
deps = ["Dates", "Printf", "RecipesBase", "Serialization", "TimeZones"]
git-tree-sha1 = "323a38ed1952d30586d0fe03412cde9399d3618b"
uuid = "d8418881-c3e1-53bb-8760-2df7ec849ed5"
version = "1.5.0"

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

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

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

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

[[LaTeXStrings]]
git-tree-sha1 = "f2355693d6778a178ade15952b7ac47a4ff97996"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.3.0"

[[LazyArtifacts]]
deps = ["Artifacts", "Pkg"]
uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"

[[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", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

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

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

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

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

[[Mocking]]
deps = ["Compat", "ExprTools"]
git-tree-sha1 = "29714d0a7a8083bba8427a4fbfb00a540c681ce7"
uuid = "78c3b35d-d492-501b-9361-3d52fe80e533"
version = "0.7.3"

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

[[MutableArithmetics]]
deps = ["LinearAlgebra", "SparseArrays", "Test"]
git-tree-sha1 = "8d9496b2339095901106961f44718920732616bb"
uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
version = "0.2.22"

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

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

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

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

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

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

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

[[Polynomials]]
deps = ["Intervals", "LinearAlgebra", "MutableArithmetics", "RecipesBase"]
git-tree-sha1 = "7499556d31417baeabaa55d266a449ffe4ec5a3e"
uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
version = "2.0.17"

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

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

[[PyCall]]
deps = ["Conda", "Dates", "Libdl", "LinearAlgebra", "MacroTools", "Serialization", "VersionParsing"]
git-tree-sha1 = "4ba3651d33ef76e24fef6a598b63ffd1c5e1cd17"
uuid = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
version = "1.92.5"

[[PyPlot]]
deps = ["Colors", "LaTeXStrings", "PyCall", "Sockets", "Test", "VersionParsing"]
git-tree-sha1 = "14c1b795b9d764e1784713941e787e1384268103"
uuid = "d330b81b-6aea-500a-939a-2ce795aea3ee"
version = "2.10.0"

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

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

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

[[Reexport]]
deps = ["Pkg"]
git-tree-sha1 = "7b1d07f411bc8ddb7977ec7f377b97b158514fe0"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "0.2.0"

[[RoundingEmulator]]
git-tree-sha1 = "40b9edad2e5287e05bd413a38f61a8ff55b9557b"
uuid = "5eaf0fd0-dfba-4ccb-bf02-d820a40db705"
version = "0.2.1"

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

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

[[SetRounding]]
git-tree-sha1 = "d7a25e439d07a17b7cdf97eecee504c50fedf5f6"
uuid = "3cc68bcd-71a2-5612-b932-767ffbe40ab0"
version = "0.2.1"

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

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

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

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

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

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

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

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

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

[[TimeZones]]
deps = ["Dates", "Downloads", "InlineStrings", "LazyArtifacts", "Mocking", "Pkg", "Printf", "RecipesBase", "Serialization", "Unicode"]
git-tree-sha1 = "8de32288505b7db196f36d27d7236464ef50dba1"
uuid = "f269a46b-ccf7-5d73-abea-4c690281aa53"
version = "1.6.2"

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

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

[[ValidatedNumerics]]
deps = ["IntervalArithmetic", "IntervalConstraintProgramming", "IntervalContractors", "IntervalOptimisation", "IntervalRootFinding", "Reexport", "Test"]
git-tree-sha1 = "3e913fa00b3943675ac8ddd4121adeb1754e99f9"
uuid = "d621b6e3-7715-5857-9c6f-c67000ef6083"
version = "0.11.0"

[[VersionParsing]]
git-tree-sha1 = "e575cf85535c7c3292b4d89d89cc29e8c3098e47"
uuid = "81def892-9a0e-5fdd-b105-ffc91e053289"
version = "1.2.1"

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

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

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

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

# ╔═╡ Cell order:
# ╟─f7a8506d-2e53-4e72-aaa8-b4133a6520cf
# ╠═b70ee5f2-3c05-11ec-37db-1704f9c112b5
# ╟─5397e0bd-378c-42fe-a23e-d5eca0b5a010
# ╠═17757d95-5205-4620-b0ad-c6dd844c13cc
# ╠═c92b4d43-33a1-4ca1-8403-e2900ca6d2be
# ╠═957a1119-3ffa-4f29-aae6-8d4991e0ebfc
# ╟─64334315-5d05-4032-9ad3-8130d66fe008
# ╠═0b88f3f9-8e5f-479a-9c06-ee0aebfd200e
# ╠═20a1adeb-9286-45d7-8f27-b5a41d6052d9
# ╠═c0d8791e-190e-4799-9ec9-224fcd7bfcfb
# ╟─6ee9dd70-2e76-4fb3-8972-14edbe3f4103
# ╠═6218e2ed-902b-4a3c-9f19-8fced206df1c
# ╠═f5f2d5b5-effb-46ac-ac23-a10f88f00b46
# ╠═954f9963-1249-43b4-800e-6d46604a27b0
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: Flamanjaune/PHD-Thesis

max_issues_repo_name: Flamanjaune/PHD-Thesis

max_forks_repo_name: Flamanjaune/PHD-Thesis


text:
using DataStructures

# this is the main heurisctic used to find the test fragments
function generate_basis(eq, x; homotopy=true)
    # if homotopy return generate_homotopy(eq, x) end
    eq = expand(eq)
    S = 0 #Set{Any}()
    for t in terms(eq)
        q = equivalent(t, x)
        f = kernel(q)

        C₁ = closure(f, x)
        p = q * inverse(f)

        if homotopy
            if isdependent(p, x)
                C₂ = generate_homotopy(p, x)
            else
                C₂ = 1
            end
        else
            C₂ = find_candidates(p, x)
        end

        S += sum(c₁*c₂ for c₁ in C₁ for c₂ in C₂)
    end
    return unique([one(x); [equivalent(t,x) for t in terms(S)]])
end

function expand_basis(basis, x)
    b = sum(basis)
    eq = (1 + x) * (b + Differential(x)(b))
    eq = expand(expand_derivatives(eq))
    S = Set{Any}()
    enqueue_expr!(S, eq, x)
    return [one(x); [s for s in S]]
end

function closure(eq, x; max_terms=50)
    if !isdependent(eq, x) return [one(x)] end
    D = Differential(x)
    S = Set{Any}()
    q = Queue{Any}()
    enqueue_expr!(S, q, eq, x)

    while !isempty(q) && length(S) < max_terms
        y = dequeue!(q)
        enqueue_expr!(S, q, expand_derivatives(D(y)), x)
    end
    unique([one(x); [s for s in S]; [s*x for s in S]])
end

function find_candidates(eq, x)
    eq = apply_d_rules(eq)
    D = Differential(x)

    S = Set{Any}()
    q = Queue{Any}()
    enqueue_expr!(S, q, eq, x)

    for y in q
        ∂y = expand_derivatives(D(y))
        enqueue_expr!(S, ∂y, x)
    end

    return unique([one(x); [s for s in S]])
end

function candidate_pow_minus(p, k; abstol=1e-8)
    if isnan(poly_deg(p))
        return [p^k, p^(k+1), log(p)]
    end

    x = var(p)
    d = poly_deg(p)

    for j = 1:10  # will try 10 times to find the roots
        r, s = find_roots(p, x)
        if length(r) + length(s) >= d break end
    end
    s = s[1:2:end]
    r = nice_parameter.(r)
    s = nice_parameter.(s)

    # ∫ 1 / ((x-z₁)(x-z₂)) dx = ... + c₁ * log(x-z₁) + c₂ * log(x-z₂)
    q = Any[log(x - u) for u in r]
    for i in eachindex(s)
        β = s[i]
        if abs(imag(β)) > abstol
            push!(q, atan((x - real(β))/imag(β)))
            push!(q, (1+x)*log(x^2 - 2*real(β)*x + abs2(β)))
        else
            push!(q, log(x - real(β)))
        end
    end
    q = unique(q)

    # return [[p^k, p^(k+1)]; candidates(q₁, x)]
    if k ≈ -1
        return [[p^k]; q]
    else
        return [[p^k, p^(k+1)]; q]
    end
end

function candidate_sqrt(p, k)
    x = var(p)

    h = Any[p^k, p^(k+1)]

    if poly_deg(p) == 2
        r, s = find_roots(p, x)
        l = leading(p, x)
        if length(r) == 2
            if sum(r) ≈ 0
                r₁ = abs(r[1])
                if l > 0
                    push!(h, acosh(x/r₁))
                else
                    push!(h, asin(x/r₁))
                end
            end
        elseif real(s[1]) ≈ 0
            push!(h, asinh(x/imag.(s[1])))
        end
    end

    Δ = expand_derivatives(Differential(x)(p))
    push!(h, log(0.5*Δ + sqrt(p)))
    return h
end

###############################################################################

function enqueue_expr!(S, q, eq::SymbolicUtils.Add, x)
    for t in arguments(eq)
        enqueue_expr!(S, q, t, x)
    end
end

function enqueue_expr!(S, q, eq, x)
    y = eq / coef(eq, x)
    if y ∉ S && isdependent(y, x)
        enqueue!(q, y)
        push!(S, y)
    end
end

function enqueue_expr!(S, eq::SymbolicUtils.Add, x)
    for t in arguments(eq)
        enqueue_expr!(S, t, x)
    end
end

function enqueue_expr!(S, eq, x)
    y = eq / coef(eq, x)
    if y ∉ S && isdependent(y, x)
        push!(S, y)
    end
end


text:
"""
exrss
"""

function exrss(Q, r; maxiter=5000, steps=20, rtol=1.0e-6) where {Tv<:Number}
    dot(sprob(Q, maxiter=maxiter, steps=steps, rtol=rtol), r)
end

function exrss(m::CTMCModel{Tv}; reward, maxiter=5000, steps=20, rtol=1.0e-6) where Tv
    exrss(m.Q, m.reward[reward], maxiter=maxiter, steps=steps, rtol=rtol)
end

"""
exrt, cexrt
"""

function exrt(ts::Tt, Q, x0, r; ufact=1.01, eps=1.0e-8, rmax=500) where {Tt<:Number}
    dot(tprob(ts, Q, x0, ufact=ufact, eps=eps, rmax=rmax), r)
end

function exrt(ts::Tt, m::CTMCModel{Tv}; reward, ufact=1.01, eps=1.0e-8, rmax=500) where {Tv,Tt<:Number}
    exrt(ts, m.Q, m.initv, m.reward[reward], ufact=ufact, eps=eps, rmax=rmax)
end

function cexrt(ts::Tt, Q, x0, r; ufact=1.01, eps=1.0e-8, rmax=500) where {Tt<:Number}
    dot(ctprob(ts, Q, x0, ufact=ufact, eps=eps, rmax=rmax), r)
end

function cexrt(ts::Tt, m::CTMCModel{Tv}; reward, ufact=1.01, eps=1.0e-8, rmax=500) where {Tv,Tt<:Number}
    cexrt(ts, m.Q, m.initv, m.reward[reward], ufact=ufact, eps=eps, rmax=rmax)
end

###

function exrt(ts::AbstractVector{Tt}, Q, x0, r; forward=:T, ufact=1.01, eps=1.0e-8, rmax=500) where {Tt<:Number}
    _exrt(ts, Q, x0, r, forward, ufact, eps, rmax)
end

function cexrt(ts::AbstractVector{Tt}, Q, x0, r; forward=:T, ufact=1.01, eps=1.0e-8, rmax=500) where {Tt<:Number}
    _cexrt(ts, Q, x0, r, forward, ufact, eps, rmax)
end

function exrt(ts::AbstractVector{Tt}, m::CTMCModel{Tv}; reward, forward=:T, ufact=1.01, eps=1.0e-8, rmax=500) where {Tv,Tt<:Number}
    exrt(ts, m.Q, m.initv, m.reward[reward], forward=forward, ufact=ufact, eps=eps, rmax=rmax)
end

function cexrt(ts::AbstractVector{Tt}, m::CTMCModel{Tv}; reward, forward=:T, ufact=1.01, eps=1.0e-8, rmax=500) where {Tv,Tt<:Number}
    cexrt(ts, m.Q, m.initv, m.reward[reward], forward=forward, ufact=ufact, eps=eps, rmax=rmax)
end

"""
exrt, cexrt for ts
"""

function _exrt(ts::AbstractVector{Tt}, Q::AbstractMatrix{Tv}, x::Vector{Tv}, r::Vector{Tv}, forward, ufact, eps, rmax) where {Tv<:Number,Tt<:Number}
    result, _, = tran(Q, x, r, ts, forward=forward, ufact=ufact, eps=eps, rmax=rmax)
    result
end

function _cexrt(ts::AbstractVector{Tt}, Q::AbstractMatrix{Tv}, x::Vector{Tv}, r::Vector{Tv}, forward, ufact, eps, rmax) where {Tv<:Number,Tt<:Number}
    _, cresult, = tran(Q, x, r, ts, forward=forward, ufact=ufact, eps=eps, rmax=rmax)
    cresult
end

"""
symbolic
"""

mutable struct SymbolicCTMCTranExpression{Tv} <: AbstractVectorSymbolic{Tv}
    params::Set{Symbol}
    op::Symbol
    x0::AbstractVectorSymbolic{Tv}
    Q::AbstractMatrixSymbolic{Tv}
    r::AbstractVectorSymbolic{Tv}
    ts::AbstractVector{<:Number}
    options::Dict{Symbol,Any}
    dim::Int
end

function _toexpr(x::SymbolicCTMCTranExpression)
    Expr(:call, x.op, _toexpr(x.x0), _toexpr(x.Q), _toexpr(x.r))
end

function Base.show(io::IO, x::SymbolicCTMCTranExpression{Tv}) where Tv
    Base.show(io, "SymbolicCTMCTranExpression $(objectid(x))")
end

"""
exrt, cexrt
"""

function _exrt(ts::AbstractVector{Tt}, Q::AbstractMatrix{<:AbstractSymbolic{Tv}}, x::Vector{<:AbstractSymbolic{Tv}}, r::Vector{<:AbstractSymbolic{Tv}}, forward, ufact, eps, rmax) where {Tv<:Number,Tt<:Number}
    _exrt(ts, convert(AbstractMatrixSymbolic{Tv}, Q), convert(AbstractVectorSymbolic{Tv}, x), convert(AbstractVectorSymbolic{Tv}, r), forward, ufact, eps, rmax)
end

function _cexrt(ts::AbstractVector{Tt}, Q::AbstractMatrix{<:AbstractSymbolic{Tv}}, x::Vector{<:AbstractSymbolic{Tv}}, r::Vector{<:AbstractSymbolic{Tv}}, forward, ufact, eps, rmax) where {Tv<:Number,Tt<:Number}
    _cexrt(ts, convert(AbstractMatrixSymbolic{Tv}, Q), convert(AbstractVectorSymbolic{Tv}, x), convert(AbstractVectorSymbolic{Tv}, r), forward, ufact, eps, rmax)
end

function _exrt(ts::AbstractVector{Tt}, Q::AbstractMatrixSymbolic{Tv}, x::AbstractVectorSymbolic{Tv}, r::AbstractVectorSymbolic{Tv}, forward, ufact, eps, rmax) where {Tv<:Number,Tt<:Number}
    s = union(x.params, Q.params, r.params)
    SymbolicCTMCTranExpression{Tv}(s, :exrt, x, Q, r, ts, Dict(:forward=>forward, :ufact=>ufact, :eps=>eps, :rmax=>rmax), length(ts))
end

function _cexrt(ts::AbstractVector{Tt}, Q::AbstractMatrixSymbolic{Tv}, x::AbstractVectorSymbolic{Tv}, r::AbstractVectorSymbolic{Tv}, forward, ufact, eps, rmax) where {Tv<:Number,Tt<:Number}
    s = union(x.params, Q.params, r.params)
    SymbolicCTMCTranExpression{Tv}(s, :cexrt, x, Q, r, ts, Dict(:forward=>forward, :ufact=>ufact, :eps=>eps, :rmax=>rmax), length(ts))
end

"""
seval(f, env, cache)
Return the value for expr f
"""

function _eval(::Val{:exrt}, f::SymbolicCTMCTranExpression{Tv}, env::SymbolicEnv, cache::SymbolicCache) where Tv
    Q = seval(f.Q, env, cache)
    x = seval(f.x0, env, cache)
    r = seval(f.r, env, cache)
    ts = f.ts
    result, cresult, = tran(Q, x, r, ts, forward=f.options[:forward], ufact=f.options[:ufact], eps=f.options[:eps], rmax=f.options[:rmax])
    result
end

function _eval(::Val{:exrt}, f::SymbolicCTMCTranExpression{Tv}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache) where Tv
    Q = seval(f.Q, env, cache)
    x = seval(f.x0, env, cache)
    r = seval(f.r, env, cache)
    dQ = seval(f.Q, dvar, env, cache)
    dx = seval(f.x0, dvar, env, cache)
    dr = seval(f.r, dvar, env, cache)
    ts = f.ts

    n = length(x)

    QQ = SparseCSC(BlockCOO(2, 2, [(1,1,Q), (1,2,dQ), (2,2,Q)]))
    xx = [x..., dx...]
    rr = [dr..., r...]
    result, cresult, = tran(QQ, xx, rr, ts, forward=f.options[:forward], ufact=f.options[:ufact], eps=f.options[:eps], rmax=f.options[:rmax])
    result
end

function _eval(::Val{:exrt}, f::SymbolicCTMCTranExpression{Tv}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache) where Tv
    Q = seval(f.Q, env, cache)
    x = seval(f.x0, env, cache)
    r = seval(f.r, env, cache)

    dQ_a = seval(f.Q, dvar[1], env, cache)
    dx_a = seval(f.x0, dvar[1], env, cache)
    dr_a = seval(f.r, dvar[1], env, cache)

    dQ_b = seval(f.Q, dvar[2], env, cache)
    dx_b = seval(f.x0, dvar[2], env, cache)
    dr_b = seval(f.r, dvar[2], env, cache)

    dQ_ab = seval(f.Q, dvar, env, cache)
    dx_ab = seval(f.x0, dvar, env, cache)
    dr_ab = seval(f.r, dvar, env, cache)
    ts = f.ts

    n = length(x)
    QQ = SparseCSC(BlockCOO(4, 4, [(1,1,Q), (2,2,Q), (3,3,Q), (4,4,Q), (1,2,dQ_a), (1,3,dQ_b), (1,4,dQ_ab), (2,4,dQ_b), (3,4,dQ_a)]))
    xx = [x..., dx_a..., dx_b..., dx_ab...]
    rr = [dr_ab..., dr_b..., dr_a..., r...]
    result, cresult, = tran(QQ, xx, rr, ts, forward=f.options[:forward], ufact=f.options[:ufact], eps=f.options[:eps], rmax=f.options[:rmax])
    result
end

"""
seval for cexrt
"""

function _eval(::Val{:cexrt}, f::SymbolicCTMCTranExpression{Tv}, env::SymbolicEnv, cache::SymbolicCache) where Tv
    Q = seval(f.Q, env, cache)
    x = seval(f.x0, env, cache)
    r = seval(f.r, env, cache)
    ts = f.ts
    result, cresult, = tran(Q, x, r, ts, forward=f.options[:forward], ufact=f.options[:ufact], eps=f.options[:eps], rmax=f.options[:rmax])
    cresult
end

function _eval(::Val{:cexrt}, f::SymbolicCTMCTranExpression{Tv}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache) where Tv
    Q = seval(f.Q, env, cache)
    x = seval(f.x0, env, cache)
    r = seval(f.r, env, cache)
    dQ = seval(f.Q, dvar, env, cache)
    dx = seval(f.x0, dvar, env, cache)
    dr = seval(f.r, dvar, env, cache)
    ts = f.ts

    n = length(x)

    QQ = SparseCSC(BlockCOO(2, 2, [(1,1,Q), (1,2,dQ), (2,2,Q)]))
    xx = [x..., dx...]
    rr = [dr..., r...]
    result, cresult, = tran(QQ, xx, rr, ts, forward=f.options[:forward], ufact=f.options[:ufact], eps=f.options[:eps], rmax=f.options[:rmax])
    cresult
end

function _eval(::Val{:cexrt}, f::SymbolicCTMCTranExpression{Tv}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache) where Tv
    Q = seval(f.Q, env, cache)
    x = seval(f.x0, env, cache)
    r = seval(f.r, env, cache)

    dQ_a = seval(f.Q, dvar[1], env, cache)
    dx_a = seval(f.x0, dvar[1], env, cache)
    dr_a = seval(f.r, dvar[1], env, cache)

    dQ_b = seval(f.Q, dvar[2], env, cache)
    dx_b = seval(f.x0, dvar[2], env, cache)
    dr_b = seval(f.r, dvar[2], env, cache)

    dQ_ab = seval(f.Q, dvar, env, cache)
    dx_ab = seval(f.x0, dvar, env, cache)
    dr_ab = seval(f.r, dvar, env, cache)
    ts = f.ts

    n = length(x)
    QQ = SparseCSC(BlockCOO(4, 4, [(1,1,Q), (2,2,Q), (3,3,Q), (4,4,Q), (1,2,dQ_a), (1,3,dQ_b), (1,4,dQ_ab), (2,4,dQ_b), (3,4,dQ_a)]))
    xx = [x..., dx_a..., dx_b..., dx_ab...]
    rr = [dr_ab..., dr_b..., dr_a..., r...]
    result, cresult, = tran(QQ, xx, rr, ts, forward=f.options[:forward], ufact=f.options[:ufact], eps=f.options[:eps], rmax=f.options[:rmax])
    cresult
end


text:
##########################################################
"""
$(TYPEDEF)

Abstract type for finite volume system structure.
"""
abstract type AbstractSystem{Tv<:Number, Ti <:Integer, Tm <:Integer} end

default_check_allocs()= haskey(ENV,"VORONOIFVM_CHECK_ALLOCS") ? parse(Bool,ENV["VORONOIFVM_CHECK_ALLOCS"]) :  false

##################################################################
"""
````
function System(grid,physics;
                unknown_storage=:dense;
                matrixindextype=Int32)
````

Create Finite Volume System.

- `grid`: 1D/2D/3D discretization grid
- `physics`: Physics struct containing node and edge callbacks
- `unknown_storage`: string or symbol:
     - `:dense` :  solution vector is an  `nspecies` x `nnodes`  dense matrix
     - `:sparse` :  solution vector is an `nspecies` x `nnodes`  sparse matrix
- `matrixindextype` : Index type for sparse matrices created in the system
"""
function System(grid,physics::Physics; unknown_storage=:dense, matrixindextype=Int64, check_allocs=default_check_allocs())
    if Symbol(unknown_storage)==:dense
        sys=DenseSystem(grid,physics, matrixindextype=matrixindextype)
    elseif Symbol(unknown_storage)==:sparse
        sys=SparseSystem(grid,physics, matrixindextype=matrixindextype)
    else
        throw("specify either unknown_storage=:dense  or unknown_storage=:sparse")
    end
    check_allocs!(sys,check_allocs)
end

# Test if allocated is zero and if allocation check is enabled
# However we need to be aware that @allocated reports allocations
# during the compilation phase. So we need to wait for at least
# one run of the system in order enact the checking.
function _check_allocs(system, allocated)
       if system.allocs>=0 # we had a couple of runs before to bridge the compilation phase
           system.allocs=allocated
           return system.allocs==0
       elseif system.allocs > -100 # probably still in compiling phase
           system.allocs=system.allocs+1
           return true
       else
           # otherwise, checking has been switched off.
           return true
       end
end

"""
```
check_allocs!(system,true_or_false)
```

Enable/disable checking for time-consuming allocations in the assembly loop.
By default,  this check  is switched off.  By setting  the environment
variable `ENV["VORONOIFVM_CHECK_ALLOCS"]="true"`, this  default can be
changed.

Unless  the   matrix  pattern  changes,  there   shouldn't  occur  any
allocations in this loop. The check  method is aware of matrix pattern
changes. As a consequence, allocations in the assembly loop are mostly
due to type instabilities in physics callbacks, see the the discussion
[here](../runexamples/#Performance-with-closures).  Type instabilities
can be debugged via the `@time`  macro applied to expressions in a
physics callback.

The following  cases provide some ideas  where to look for  reasons of
the problem and possible remedies:

Case 1: a parameter changes its value, and Julia is not sure about the type.
```julia
eps=1.0

flux(f,_u,edge)
    u=unkowns(edge,_u)
    f[1]=eps*(u[1,1]-[1,2])
end
... solve etc ...
eps=2.0
```
Remedy: use a type annotation `eps::Float64=...` to signalize your intent to Julia.
This behaviour is explained in the [Julia documentation](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-captured).

Case 2: variables in the closure have the same name as a variable
introduced in a callback.
```julia
flux(f,_u,edge)
    u=unkowns(edge,_u)
    f[1]=(u[1,1]-[1,2])
end

... create etc ...

u=solve(...)
```
Remedy: rename e.g. `u=solve()` to `sol=solve()`

"""
function check_allocs!(system::AbstractSystem,chk::Bool)
    if chk
        system.allocs=-2
    else
        system.allocs=-1000
    end
    system
end

##################################################################

# Constant to be used as boundary condition factor
# to mark Dirichlet boundary conditons.
const Dirichlet=1.0e30

##################################################################
"""
````
is_boundary_species(AbstractSystem, ispec) -> Bool
````

Check if species number corresponds to a boundary species.
"""
function is_boundary_species(system::AbstractSystem, ispec)
    isbspec=false
    for ibreg=1:num_bfaceregions(system.grid)
        if system.bregion_species[ispec,ibreg]>0
            isbspec=true
        end
    end
    return isbspec
end

##################################################################
"""
````
is_bulk_species(AbstractSystem, ispec) -> Bool
````

Check if species number corresponds to a bulk species.
"""
function is_bulk_species(system::AbstractSystem, ispec)
    isrspec=false
    for ixreg=1:num_cellregions(system.grid)
        if system.region_species[ispec,ixreg]>0
            isrspec=true
        end
    end
    return isrspec
end

##################################################################
"""
````
enable_species!(system,ispec,regions)
````

Add species `ispec` to a list of bulk regions. Species numbers for
bulk and boundary species have to be distinct.
"""
function enable_species!(system::AbstractSystem,ispec::Integer, regions::AbstractVector)
    if ispec>num_species(system)
        throw(DomainError(ispec,"Number of species exceeded"))
    end
    if is_boundary_species(system,ispec)
        throw(DomainError(ispec,"Species is already boundary species"))
    end
    _cellregions=cellregions(system.grid)
    _cellnodes=cellnodes(system.grid)
    for i in eachindex(regions)
        ireg=regions[i]
        system.region_species[ispec,ireg]=ispec
        for icell=1:num_cells(system.grid)
            if _cellregions[icell]==ireg
                for iloc=1:num_targets(_cellnodes,icell)
                    iglob=_cellnodes[iloc,icell]
                    system.node_dof[ispec,iglob]=ispec
                end
            end
        end
    end
end

##################################################################
"""
````
enable_boundary_species!(system,ispec,regions)
````

Add species `ispec` to a list of boundary regions. Species numbers for
bulk and boundary species have to be distinct.

"""
function enable_boundary_species!(system::AbstractSystem, ispec::Integer, bregions::AbstractVector)
    if ispec>num_species(system)
        throw(DomainError(ispec,"Number of species exceeded"))
    end
    if is_bulk_species(system,ispec)
        throw(DomainError(ispec,"Species is already bulk species"))
    end
    _bfaceregions=bfaceregions(system.grid)
    _bfacenodes=bfacenodes(system.grid)

    for i in eachindex(bregions)
        ireg=bregions[i]
        system.bregion_species[ispec,ireg]=ispec
        for ibface=1:num_bfaces(system.grid)
            if _bfaceregions[ibface]==ireg
                for iloc=1:size(_bfacenodes,1)
                    iglob=_bfacenodes[iloc,ibface]
                    system.node_dof[ispec,iglob]=ispec
                end
            end
        end
    end
end

# Create matrix in system and figure out if species
# distribution is homgeneous
function _complete!(system::AbstractSystem{Tv,Ti, Tm};create_newtonvectors=false) where {Tv,Ti, Tm}

    if isdefined(system,:matrix)
        return
    end
    system.matrix=ExtendableSparseMatrix{Tv,Tm}(num_dof(system), num_dof(system))
    system.species_homogeneous=true
    species_added=false
    for inode=1:size(system.node_dof,2)
        for ispec=1:size(system.node_dof,1)
            if system.node_dof[ispec,inode]==ispec
                species_added=true
            else
                system.species_homogeneous=false
            end
        end
    end

    if (!species_added)
        error("No species enabled.\n Call enable_species(system,species_number, list_of_regions) at least once.")
    end
    if create_newtonvectors
        system.residual=unknowns(system)
        system.update=unknowns(system)
    end
    update_grid!(system)
    if has_generic_operator(system)
        if has_generic_operator_sparsity(system)
            system.generic_matrix=system.physics.generic_operator_sparsity(system)
        else
            generic_operator(f,u)=system.physics.generic_operator(f,u,system)
            input=rand(num_dof(system))
            output=similar(input)
            tdetect=@elapsed begin
                sparsity_pattern = jacobian_sparsity(generic_operator,output,input)
                system.generic_matrix = Float64.(sparse(sparsity_pattern))
            end
            println("sparsity detection for generic operator: $(tdetect) s")
            if nnz(system.generic_matrix)==0
                error("Sparsity detection failed: no pattern found")
            end
        end
        tdetect=@elapsed begin
            system.generic_matrix_colors = matrix_colors(system.generic_matrix)
        end
        println("matrix coloring for generic operator: $(tdetect) s")
    end
end

"""
$(SIGNATURES)
Set generic operator sparsity, in the case where a generic operator has been
defined in physics.
"""
function generic_operator_sparsity!(system::AbstractSystem, sparsematrix::SparseMatrixCSC)
    system.generic_matrix=sparsematrix
end

"""
````
update_grid!(system; grid=system.grid)
````

Update grid (e.g. after rescaling of coordinates).
"""
function update_grid!(system::AbstractSystem{Tv,Ti}; grid=system.grid) where {Tv,Ti}

    geom        = grid[CellGeometries][1]
    csys        = grid[CoordinateSystem]
    coord       = grid[Coordinates]
    cellnodes   = grid[CellNodes]
    cellregions = grid[CellRegions]
    bgeom       = grid[BFaceGeometries][1]
    bfacenodes  = grid[BFaceNodes]
    nbfaces     = num_bfaces(grid)
    ncells      = num_cells(grid)

    system.cellnodefactors  = zeros(Tv, num_nodes(geom), ncells)
    system.celledgefactors  = zeros(Tv, num_edges(geom), ncells)
    system.bfacenodefactors = zeros(Tv, num_nodes(bgeom), nbfaces)
    system.bfaceedgefactors = zeros(Tv, num_edges(bgeom), nbfaces)

    for icell=1:ncells
        @views cellfactors!(geom, csys, coord, cellnodes, icell,
                            system.cellnodefactors[:,icell],system.celledgefactors[:,icell])
    end

    for ibface=1:nbfaces
        @views bfacefactors!(bgeom,csys,coord,bfacenodes,ibface,
                             system.bfacenodefactors[:,ibface], system.bfaceedgefactors[:, ibface])
    end
end

##################################################################
"""
$(SIGNATURES)

Check if degree of freedom is defined.
"""
isdof(system::AbstractSystem,ispec,inode)= system.node_dof[ispec,inode]==ispec

##################################################################
"""
$(SIGNATURES)

Number of species in system
"""
num_species(system::AbstractSystem) = system.physics.num_species

##################################################################
"""
$(SIGNATURES)

Retrieve user data record.
"""
data(system::AbstractSystem) = system.physics.data

##################################################################
"""
$(SIGNATURES)

Set Dirichlet boundary condition for species ispec at boundary ibc:

``u_{ispec}=v`` on ``\\Gamma_{ibc}``
"""
function boundary_dirichlet!(system::AbstractSystem, ispec, ibc, v)
    system.boundary_factors[ispec,ibc]=Dirichlet
    system.boundary_values[ispec,ibc]=v
end

##################################################################
"""
$(SIGNATURES)

Set Neumann boundary condition for species ispec at boundary ibc:

``\\mathrm{flux}_{ispec}\\cdot \\vec n=v`` on ``\\Gamma_{ibc}``
"""
function boundary_neumann!(system::AbstractSystem, ispec, ibc, v)
    system.boundary_factors[ispec,ibc]=0.0
    system.boundary_values[ispec,ibc]=v
end

##################################################################
"""
$(SIGNATURES)

Set Robin boundary condition for species ispec at boundary ibc:

``\\mathrm{flux}_{ispec}\\cdot \\vec n + \\alpha u_{ispec}=v`` on ``\\Gamma_{ibc}``

"""
function boundary_robin!(system::AbstractSystem, ispec, ibc, α, v)
    system.boundary_factors[ispec,ibc]=α
    system.boundary_values[ispec,ibc]=v
end

##################################################################
"""
$(SIGNATURES)

Number of species (size of first dimension) of solution array.
"""
num_species(a::AbstractArray)=size(a,1)

#
# Initialize Dirichlet BC
#
function _initialize_dirichlet!(U::AbstractMatrix,system::AbstractSystem)
    _bfaceregions=bfaceregions(system.grid)
    _bfacenodes=bfacenodes(system.grid)
    for ibface=1:num_bfaces(system.grid)
        ibreg=_bfaceregions[ibface]
        for ispec=1:num_species(system)
            if system.boundary_factors[ispec,ibreg]≈ Dirichlet
                for inode=1:dim_grid(system.grid)
                    U[ispec,_bfacenodes[inode,ibface]]=system.boundary_values[ispec,ibreg]
                end
            end
        end
    end
end

function _initialize!(U::AbstractMatrix,system::AbstractSystem)
    _initialize_dirichlet!(U,system)
    _initialize_inactive_dof!(U,system)
end

function Base.show(io::IO,sys::AbstractSystem)
    str=@sprintf("%s(num_species=%d)",typeof(sys),sys.physics.num_species)
    println(io,str)
end

#####################################################
has_generic_operator(sys::AbstractSystem) = sys.physics.generic_operator!=nofunc_generic
has_generic_operator_sparsity(sys::AbstractSystem) =  sys.physics.generic_operator_sparsity!=nofunc_generic_sparsity


text:
function fva(S, b, lb, ub, ridxs = eachindex(lb);
        check_obj::Union{Nothing, Int} = nothing,
        check_obj_atol::Real = 1e-4,
        verbose = true,
        batchlen::Int = 50,
        on_non_optimal_sol::Function = (idx, lp_model) -> error("FBA failed, non OPTIMAL solution returned!!!"),
        nths = nthreads(),
        solver = Clp.Optimizer
    )

    ridxs = collect(ridxs)

    # bounds
    fvalb = lb[ridxs]
    fvaub = ub[ridxs]

    verbose && (prog = Progress(length(ridxs); desc = "Doing FVA (-t$nths)  "))
    ii_ch = Channel{Int}(nths) do ch_
        for ii in eachindex(ridxs)
            put!(ch_, ii)
            verbose && next!(prog)
        end
    end

    @threads for _ in 1:nths

        lp_model = build_lp_model(S, b, lb, ub, solver)
        x = _get_vars(lp_model)

        for ii in ii_ch
            ridx = ridxs[ii]

            # optimize max
            JuMP.@objective(lp_model, MIN_SENSE, x[ridx])
            JuMP.optimize!(lp_model)
            status = JuMP.termination_status(lp_model)
            if status == JuMP.MOI.OPTIMAL
                fvalb[ii] = JuMP.objective_value(lp_model)
                else; on_non_optimal_sol(ridx, lp_model)
            end

            # optimize min
            JuMP.@objective(lp_model, MAX_SENSE, x[ridx])
            JuMP.optimize!(lp_model)
            status = JuMP.termination_status(lp_model)
            if status == JuMP.MOI.OPTIMAL
                fvaub[ii] = max(JuMP.objective_value(lp_model), fvalb[ii])
                else; on_non_optimal_sol(ridx, lp_model)
            end
        end
    end
    verbose && finish!(prog)

    # Check bounds
    if !isnothing(check_obj)
        return MetLP.check_newbounds(S, b, lb, ub, fvalb, fvaub,
            check_obj, ridxs; check_obj_atol, verbose, batchlen
        )
    end

    return fvalb, fvaub
end

function fva(model::MetNet, iders = eachindex(model.lb);
        check_obj = nothing, kwargs...
    )
    obj_idx = isnothing(check_obj) ? nothing : rxnindex(model, check_obj)
    idxs = [rxnindex(model, idx) for idx in iders]
    model_fields = _extract_dense(model, [:S, :b, :lb, :ub])
    return fva(model_fields..., idxs; check_obj = obj_idx, kwargs...);
end

fva(model::MetNet, ider::IDER_TYPE; kwargs...) = fva(model, [ider]; kwargs...)


text:
@testset verbose = true "Tests for exceptions" begin

    using JavaCall.JNI
    using JavaCall: JavaCodeGeneration

    @testset "Integer Parsing Exception" begin
        exceptions_mod = Module(:Exceptions)
        Core.eval(exceptions_mod, :(using JavaCall))
        Core.eval(exceptions_mod, JavaCodeGeneration.loadclass(
            exceptions_mod,
            Symbol("java.lang.Integer")
        ))
        Core.eval(exceptions_mod, JavaCodeGeneration.loadclass(
            exceptions_mod,
            Symbol("java.lang.String")
        ))

        invalidnumber = exceptions_mod.JString(Char['1', '!', '3'])

        @test_throws exceptions_mod.JNumberFormatExceptionJuliaImpl exceptions_mod.j_parse_int(
            exceptions_mod.JInteger,
            invalidnumber
        )
        @test_throws exceptions_mod.JNumberFormatException exceptions_mod.j_parse_int(
            exceptions_mod.JInteger,
            invalidnumber
        )
        @test_throws exceptions_mod.JRuntimeException exceptions_mod.j_parse_int(
            exceptions_mod.JInteger,
            invalidnumber
        )

        try
            exceptions_mod.j_parse_int(exceptions_mod.JInteger, invalidnumber)
        catch e
            @test isa(e, exceptions_mod.JNumberFormatException)
            @test isa(e, exceptions_mod.JThrowable)
            @test exceptions_mod.j_length(exceptions_mod.j_get_message(e)) > 0
        end
    end
end


text:
function logsumexp(logx::Array{Float64,1})
    maxlog = maximum(logx)
    maxlog + log(sum(exp.(logx - maxlog)))
end

function logsumexp(logx::Array{Float64,2}, dim::Int)
    maxlog = maximum(logx, dim)
    maxlog + log.(sum(exp.(logx .- maxlog), dim))
end


text:
using PackageCompiler
@time build_executable("scripts/executable.jl", "insarlos")


text:
module Banana

struct Node{T}
    data::T
    next::Int
end

struct LinkedList{T}
    nodes::Vector{Node{T}}
    head::Int
end

export LinkedList

function LinkedList(els::Vector)
    nodes = [Node(el,i+1) for (i,el) in enumerate(els)]
    nodes[end] = Node(last(els),0)
    LinkedList(nodes,1)
end

start(ll::LinkedList) = (0, ll.head)
done(ll::LinkedList, state) = state[2] < 1
next(ll::LinkedList, (prev, next)) = ll.nodes[next].data, (next, ll.nodes[next].next)

function Base.iterate(ll::LinkedList, state=start(ll))
    done(ll, state) ? nothing : next(ll, state)
end

export start, next, done

"""
This is how you insert elements.
"""
function Base.insert!(ll::LinkedList, pos, item)
    prev, next = pos
    push!(ll.nodes, Node(item, next))
    prev > 0 && (ll.nodes[prev] = Node(ll.nodes[prev].data, lastindex(ll.nodes)))
    prev == 0  && (ll.head = length(ll.nodes))
    return ll
end

end # module


text:
# TO RUN:
# include("main.jl")
#
# AUTHOR
# Dr. Umberto Saetti, Assistant Professor, Department of Aerospace Engineering,
# Auburn university
# saetti@auburn.edu
#
# LAST UPDATED
# 1/3/2022
#
# DESCRIPTION
# Complex helicopter model based on the General Helicopter (GenHel) flight
# dynamics simulation model (Ref. 1). This model is representative of a utility
# helicopter similar to a UH-60. The model contains a 6-DoF rigid-body dynamic
# model of the fuselage, nonlinear aerodynamic lookup tables for the fuselage,
# rotor blades, and empennage, rigid flap and lead-lag rotor blade dynamics,
# a three-state Pitt-Peters inflow model, and a one-state tail rotor model.
# The trim, linearization routines, and overall architecture of the simulation
# follows the teachings of Dr. Joe Horn at Penn State. If using the code,
# please cite Ref. 2.
#
# REFERENCES
# 1) Howlett, J. J., “UH-60A Black Hawk Engineering Simulation Program.
#    Volume 1: Mathematical Model,” Tech. rep. NASA-CR-166309, 1980.
# 2) Saetti, U., and Horn, J. F., "Flight Simulation and Control using the Julia
#    Language", AIAA Scitech Forum, San Diego, CA, Jan 3-7, 2022.
#    DOI: https://arc.aiaa.org/doi/10.2514/6.2022-2354.
#
# STATES             | NINDICES    | DESCRIPTION
# ______________________________________________________________________________
# u v w              |  1  2  3    | body velocities [ft/s]
# p q r              |  4  5  6    | angula rates [rad]
# phi theta psi      |  7  8  9    | Euler angles [rad]
# x  y  z            | 10 11 12    | position [ft]
# b0 b1c b1s b0d     | 13 14 15 16 | flapping angles [rad]
# db0 db1c db1s db0d | 17 18 19 20 | flapping angles derivatives [rad/s]
# z0 z1c z1s z0d     | 21 22 23 24 | lead-lag angles [rad]
# dz0 dz1c dz1s dz0d | 25 26 27 28 | lead-lag angles derivative [rad/s]
# lam0 lam1c lam1s   | 29 30 31    | inflow angles [rad]
# psi                | 32          | rotor azimuth [rad]
# dynamic twist      | 33          | dynamic twist force? [?]
# lam0T              | 34          | tail rotor inflow [rad]
#
# CONTROL INPUTS     | INDICED     | DESCRIPTION
# ______________________________________________________________________________
# delta_lat          | 1           | lateral stick
# delta_lon          | 2           | longitudinal stick
# delta_col          | 3           | collective stick
# delta_ped          | 4           | pedals
#
# ------------------------------------------------------------------------------

# include packages
#using DifferentialEquations
using LinearAlgebra
using SparseArrays
using ControlSystems
using MAT
using Interpolations
using ApproxFun
using Printf
#ENV["MPLBACKEND"]="tkagg"
using Plots
using JLD
#using PyPlot
#using GR

@printf("\n                  J-GenHel           \n")
@printf("\n          Author: Dr. Umberto Saetti\n")
@printf("\n-----------------------------------------------\n")

# include functions
include("control_mixing.jl")
include("table_lookup.jl")
include("stabSched.jl")
include("atan2.jl")
include("eqnmot.jl")
include("interpp1.jl")
include("linearize.jl")
include("interpp2.jl")
include("wrapper.jl")
include("trimmer.jl")
include("ModRed.jl")
include("ModRed8.jl")
include("ModRed10.jl")
include("GenHel.jl")
include("GenHel_DI.jl")
include("rk4.jl")
include("finp.jl")
include("simulate.jl")
include("simulate_DI.jl")
# load H60 constants
include("H60_constants.jl")

# ------------------------ TRIM FLIGHT DYNAMICS MODEL --------------------------

# forward, lateral, and velocities in trim [ft/s]
VXTRIM=80*1.688
VYTRIM=0.0
VZTRIM=.0
# turn rate in trim [rad/s]
PSIDTRIM=0.0
# initial heading [rad]
PSITRIM = 0.0*D2R
# initial position [ft]
ALTTRIM = 0.0
XNTRIM = 0.0
YETRIM = 0.0
# initialize trim variables
include("trimInit.jl")
# trim aircraft
x0, u0 = trimmer(GenHel!,x0,u0,xdot_targ,t)
# linearized aircraft dynamics
A, B = linearize(H60!,x0,u0,xdot0,t)

# ----------------------------- SPECTRAL ANALYSIS ------------------------------

# eigenvalues
eigs=eigvals(A)
# reduce to 8-state model
sysRed = ModRed8(A,B)
A8=sysRed[1]
B8=sysRed[2]
eigs8=eigvals(A8)
# reduce to 10-state model
sysRed = ModRed10(A,B)
A10=sysRed[1]
B10=sysRed[2]
eigs10=eigvals(A10)
# store coefficient matrices and eigenvalues
#save("AB_J-GenHel.jld","A",A,"B",B,"A8",A8,"B8",B8)
#save("eigs_J-GenHel.jld","eigs",eigs,"eigs8",eigs8)

# plot eigenvalues
gr(size=(800,600))
plt_eigs=plot(real(eigs),imag(eigs),seriestype=:scatter,
      label="Full-Order", legendfontsize=12, legend=:topleft,
      marker = (:circle, :royalblue, 6), markerstrokecolor = :royalblue)
plot!(real(eigs10),imag(eigs10),seriestype=:scatter,
      label="Reduced-Order (10-state)", legendfontsize=12,
      marker = (:utriangle, :brown3, 6), markerstrokecolor = :brown3,
      reuse = false, color=:purple)
plot!(real(eigs8),imag(eigs8),seriestype=:scatter,
      label="Reduced-Order (8-state)", legendfontsize=12,
      marker = (:star5, :forestgreen, 6), markerstrokecolor = :forestgreen)
display(plt_eigs)
xaxis!("Real",xguidefontsize=14)
yaxis!("Imag",yguidefontsize=14)

# ---------------------------- FREQUENCY RESPONSES -----------------------------

# transfer functions (full-order model)
sys=ss(A,B,Matrix{Float64}(I, 37, 37),zeros(37,4))
pdlat=tf(sys[4,1])
# transfer functions (10-state model)
sys10=ss(A10,B10,Matrix{Float64}(I, 10, 10),zeros(10,4))
pdlat10=tf(sys10[4,1])
# transfer functions (8-state model)
sys8=ss(A8,B8,Matrix{Float64}(I, 8, 8),zeros(8,4))
pdlat8=tf(sys8[4,1])
# freq response (pdlat)
mag, phase, w = bode(pdlat, 10.0.^(range(-1,stop=2,length=1000)))
mag10, phase10, w10 = bode(pdlat10, 10.0.^(range(-1,stop=2,length=1000)))
mag8, phase8, w8 = bode(pdlat8, 10.0.^(range(-1,stop=2,length=1000)))

# plot (pdlat)
gr(size=(800,600))
p1_mag=plot(vec(w), 20*log10.(vec(mag)/perc2in[1,1]),
      label="Full Order", reuse = false, line=(3, :royalblue, :solid),
      legendfontsize=12)
plot!(vec(w10), 20*log10.(vec(mag10)/perc2in[1,1]), label="Reduced Order (10-state)",
                          reuse = false, line=(3, :brown3, :dash))
plot!(vec(w8), 20*log10.(vec(mag8)/perc2in[1,1]), label="Reduced Order (8-state)",
                          reuse = false, line=(3, :forestgreen, :dashdot))
xaxis!(:log10)
yaxis!("Magnitude [dB]",yguidefontsize=14)
xlims!((0.1,30))
ylims!((-50,0))
p1_phase=plot(vec(w), wrapper(vec(phase)).-360, label="", line=(3, :royalblue, :solid))
plot!(vec(w10), wrapper(vec(phase10)).-360, label="", line=(3, :brown3, :dash))
plot!(vec(w8), wrapper(vec(phase8)).-360, label="", line=(3, :forestgreen, :dashdot))
yaxis!("Phase [deg]",yguidefontsize=14)
xaxis!(:log10)
xlims!((0.1,30))
p1=plot(p1_mag, p1_phase, layout=(2,1), leftmargin=3Plots.mm)
display(p1)

# ------------------------------ TIME SIMULATION -------------------------------

# time step [s]
dt=0.01
# length of simulation [s]
Tsim=10
# run open-loop simulation
state_OL, time_OL = simulate(GenHel!,finp,Tsim,dt,x0,xdot0,u0)
# number of dynamic inverse states
NDISTATES=9
# initial state vector of closed-loop simulation
x0_DI=[x0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; u0[1]; u0[2]; u0[4]]
# initial state derivative vector of closed-loop simulation
xdot0_DI=[xdot0; zeros(NDISTATES,1)]
# run closed-loop simulation
state_CL, time_CL = simulate_DI(GenHel_DI!,finp,Tsim,dt,x0,xdot0,u0)
# store open- and closed-loop simulations
save("resp_GenHel.jld","time_OL",time_OL,"time_CL",time_CL,
    "state_OL",state_OL,"state_CL",state_CL)

# plot attitude
gr(size=(800,600))
i=7
p1=plot(time_OL,state_OL[i,:]*180/pi,label="Open Loop",legend=:topright,
    legendfontsize=12,line=(3, :royalblue, :solid))
plot!(time_OL,state_CL[i,:]*180/pi,label="Closed Loop",legend=:topright,
    legendfontsize=12,line=(3, :brown3, :dash))
yaxis!("φ [deg]",yguidefontsize=14)
xlims!((0,10))
i=8
p2=plot(time_OL,state_OL[i,:]*180/pi,label="",line=(3, :royalblue, :solid))
plot!(time_OL,state_CL[i,:]*180/pi,label="",line=(3, :brown3, :dash))
yaxis!("θ [deg]",yguidefontsize=14)
xlims!((0,10))
i=9
p3=plot(time_OL,state_OL[i,:]*180/pi,label="",line=(3, :royalblue, :solid))
plot!(time_OL,state_CL[i,:]*180/pi,label="",line=(3, :brown3, :dash))
xaxis!("Time [s]",xguidefontsize=14)
yaxis!("ψ [deg]",yguidefontsize=14)
xlims!((0,10))
p4=plot(p1, p2, p3, layout=(3,1), leftmargin=3Plots.mm)
display(p4)
#savefig(".\\Plots\\att_J-GenHel.svg")

# plot angular rates
gr(size=(800,600))
i=4
p5=plot(time_OL,state_OL[i,:],label="Open Loop",legend=:bottomleft,
    legendfontsize=12,line=(3, :royalblue, :solid))
plot!(time_OL,state_CL[i,:],label="Closed Loop",legend=:bottomleft,
    legendfontsize=12,line=(3, :brown3, :dash))
yaxis!("p [rad/s]",yguidefontsize=14)
xlims!((0,10))
i=5
p6=plot(time_OL,state_OL[i,:],label="",line=(3, :royalblue, :solid))
plot!(time_OL,state_CL[i,:],label="",line=(3, :brown3, :dash))
yaxis!("q [rad/s]",yguidefontsize=14)
xlims!((0,10))
i=6
p7=plot(time_OL,state_OL[i,:],label="",line=(3, :royalblue, :solid))
plot!(time_OL,state_CL[i,:],label="",line=(3, :brown3, :dash))
xaxis!("Time [s]",xguidefontsize=14)
yaxis!("r [rad/s]",yguidefontsize=14)
xlims!((0,10))
p8=plot(p5, p6, p7, layout=(3,1), leftmargin=3Plots.mm)
display(p8)
#savefig(".\\Plots\\ang_J-GenHel.svg")

max_stars_repo_name: umbertosaetti/J-GenHel

max_issues_repo_name: umbertosaetti/J-GenHel

max_forks_repo_name: umbertosaetti/J-GenHel


text:
import Base.length, Base.size

export getSF, SVD, tikhonovLU, setlambda

function converttoreal(S::AbstractArray{Complex{T}},f) where T
  N = prod(calibSize(f))
  M = div(length(S),N)
  S = reshape(S,N,M)
  S = reshape(reinterpret(T,vec(S)),(2*N,M))
  for l=1:M
    tmp = S[:,l]
    S[1:N,l] = tmp[1:2:end]
    S[N+1:end,l] = tmp[2:2:end]
  end
  return S
end

setlambda(S::AbstractMatrix, λ) = nothing

function getSF(bSF, frequencies, sparseTrafo, solver; kargs...)
  SF, grid = getSF(bSF, frequencies, sparseTrafo; kargs...)
  if solver == "kaczmarz"
    return transpose(SF), grid
  elseif solver == "pseudoinverse"
    return SVD(svd(transpose(SF))...), grid
  elseif solver == "cgnr" || solver == "lsqr" || solver == "fusedlasso"
    return copy(transpose(SF)), grid
  elseif solver == "direct"
    return RegularizedLeastSquares.tikhonovLU(copy(transpose(SF))), grid
  else
    return SF, grid
  end
end

getSF(bSF::Union{T,Vector{T}}, frequencies, sparseTrafo::Nothing; kargs...) where {T<:MPIFile} =
   getSF(bSF, frequencies; kargs...)

# TK: The following is a hack since colored and sparse trafo are currently not supported
getSF(bSF::Vector{T}, frequencies, sparseTrafo::AbstractString; kargs...) where {T<:MPIFile} =
  getSF(bSF[1], frequencies, sparseTrafo; kargs...)

getSF(bSFs::MultiContrastFile, frequencies; kargs...) =
   getSF(bSFs.files, frequencies; kargs...)

function getSF(bSFs::Vector{T}, frequencies; kargs...) where {T<:MPIFile}
  maxFov = [0.0,0.0,0.0]
  maxSize = [0,0,0]
  for l=1:length(bSFs)
    maxFov = max.(maxFov, calibFov(bSFs[l]))
    maxSize = max.(maxSize, collect(calibSize(bSFs[l])) )
  end
  data = [getSF(bSF, frequencies; gridsize=maxSize, fov=maxFov, kargs...) for bSF in bSFs]
  return cat([d[1] for d in data]..., dims=1), data[1][2] # grid of the first SF
end

getSF(f::MPIFile; recChannels=1:numReceivers(f), kargs...) = getSF(f, filterFrequencies(f, recChannels=recChannels); kargs...)

function repairDeadPixels(S, shape, deadPixels)
  shapeT = tuple(shape...)
  for k=1:size(S,2)
    for dp in deadPixels
      ix,iy,iz = ind2sub(shapeT,dp)
      if 1<ix<shape[1] && 1<iy<shape[2] && 1<iz<shape[3]
         lval = S[ sub2ind(shapeT,ix,iy+1,iz)  ,k]
         rval = S[ sub2ind(shapeT,ix,iy-1,iz)  ,k]
         S[dp,k] = 0.5*(lval+rval)
         fval = S[ sub2ind(shapeT,ix+1,iy,iz)  ,k]
         bval = S[ sub2ind(shapeT,ix-1,iy,iz)  ,k]
         S[dp,k] = 0.5*(fval+bval)
         uval = S[ sub2ind(shapeT,ix,iy,iz+1)  ,k]
         dval = S[ sub2ind(shapeT,ix,iy,iz-1)  ,k]
         S[dp,k] = 0.5*(uval+dval)
      end
    end
  end
end

function getSF(bSF::MPIFile, frequencies; returnasmatrix = true, procno::Integer=1,
               bgcorrection=false, bgCorrection=bgcorrection, loadasreal=false,
	       gridsize=collect(calibSize(bSF)),
	       fov=calibFov(bSF), center=[0.0,0.0,0.0], deadPixels=Int[], kargs...)

  nFreq = rxNumFrequencies(bSF)

  S = getSystemMatrix(bSF, frequencies, bgCorrection=bgCorrection; kargs...)

  if !isempty(deadPixels)
    repairDeadPixels(S,gridsize,deadPixels)
  end

  if collect(gridsize) != collect(calibSize(bSF)) ||
    center != [0.0,0.0,0.0] ||
    fov != calibFov(bSF)
    @debug "Perform SF Interpolation..."

    origin = RegularGridPositions(calibSize(bSF),calibFov(bSF),[0.0,0.0,0.0])
    target = RegularGridPositions(gridsize,fov,center)

    SInterp = zeros(eltype(S),prod(gridsize),length(frequencies)*acqNumPeriodsPerFrame(bSF))
    for k=1:length(frequencies)*acqNumPeriodsPerFrame(bSF)
      A = MPIFiles.interpolate(reshape(S[:,k],calibSize(bSF)...), origin, target)
      SInterp[:,k] = vec(A)
    end
    S = SInterp
    grid = target
  else
    grid = RegularGridPositions(calibSize(bSF),calibFov(bSF),[0.0,0.0,0.0])
  end

  if loadasreal
    S = converttoreal(S,bSF)
    resSize = [gridsize..., 2*length(frequencies)*acqNumPeriodsPerFrame(bSF)]
  else
    resSize = [gridsize...,length(frequencies)*acqNumPeriodsPerFrame(bSF)]
  end

  if returnasmatrix
    return reshape(S,prod(resSize[1:end-1]),resSize[end]), grid
  else
    return squeeze(reshape(S,resSize...)), grid
  end
end

"""
This function reads the component belonging to the given mixing-factors und channel from the given MPIfile.\n
In the case of a 2D systemfunction, the third dimension is added.
"""
function getSF(bSF::MPIFile, mx::Int, my::Int, mz::Int, recChan::Int)
  maxFreq  = rxNumFrequencies(bSF)
  freq = mixFactorToFreq(bSF, mx, my, mz)
  freq = clamp(freq,0,maxFreq-1)
  k = freq + (recChan-1)*maxFreq
  @debug "Frequency = $k"
  SF, grid = getSF(bSF,[k+1],returnasmatrix = true, bgcorrection=true)
  N = calibSize(bSF)
  SF = reshape(SF,N[1],N[2],N[3])

  return SF
end

include("Sparse.jl")


text:
#= This line is a comment containing operators like &, - and %
 A comment with the German word "Übergrößengeschäft" (store for oversizes) =#
# This line contains some special unicode characters: €, α, γ, w², Δₓ, ∩
# A comment with some numbers: 424, 1.23, 0.2E-5, -9.9e+9
function Style_3rd_Test(x, y)
    myver = v"2.00"
    mystr = "A string with \"Übergrößengeschäft\", α, π, ∪ and the + operator."
    myset = Set( [2, 9, 12, 25, 33])
    x_in_myset = x ∈ myset
    myset⁽²⁾ = myset ∪ Set( [4, 8, 12, 33])
    z1vec = rand(Int8, 3)
    z₂vec = Array{Int8}(undef, 3)
    z₂vec[1:2] = [x % y, y \ x]
    t = x % 2 == 0 ? x : x + 1
    t = ~(t & x | y) ⊻ y
    α = √0.3
    βᵅ = 3.2e+5^α
    myβvar = ∛0.12E-2 * βᵅ
    z₂vec[3] = y^2 + 3.4x*y - (α + myβvar) * t/2
    z₂vec = (z₂vec + z1vec).^2
    if !(0.1 ≤ norm(z1vec') < norm(z₂vec') + e ÷ pi + γ + φ)
        mystr = String( mystr, " signed ")
        println( mystr)
        return true;
    elseif 3.2 ≥ norm(z₂vec - z1vec) > 2.69
        if norm(z₂vec - z1vec) ≠ 3.0
            println( String( "Error in ", myver, "!"))
        end
        return false;
    end
end


text:
module Chronometer
    using Dates
    export message_with_time

    """Функция, которая возвращает строку, которая сформирована
    из служебных символов, и времени вызова.

    message: сообщение, которое требуется обогатить датой и временем
    """
    function message_with_time(message)::String
        @info "Вызвана функция хронометра"

        @debug "Параметры вызова функции:"
        @debug "> сообщение, которое требуется передать" message

        return string("<", Dates.now(), "> ", message)
    end
end


text:
module CenterIndexedArrays

using Interpolations, OffsetArrays
using OffsetArrays: IdentityUnitRange

export CenterIndexedArray

include("symrange.jl")

"""
A `CenterIndexedArray` is one for which the array center has indexes
`0,0,...`. Along each coordinate, allowed indexes range from `-n:n`.

CenterIndexedArray(A) "converts" `A` into a CenterIndexedArray. All
the sizes of `A` must be odd.
"""
struct CenterIndexedArray{T,N,A<:AbstractArray} <: AbstractArray{T,N}
    data::A
    halfsize::NTuple{N,Int}

    function CenterIndexedArray{T,N,A}(data::A) where {T,N,A<:AbstractArray}
        new{T,N,A}(data, _halfsize(data))
    end
end

CenterIndexedArray(A::AbstractArray{T,N}) where {T,N} = CenterIndexedArray{T,N,typeof(A)}(A)
CenterIndexedArray{T,N}(::UndefInitializer, sz::Vararg{<:Integer,N}) where {T,N} =
    CenterIndexedArray(Array{T,N}(undef, sz...))
CenterIndexedArray{T,N}(::UndefInitializer, sz::NTuple{N,Integer}) where {T,N} =
    CenterIndexedArray(Array{T,N}(undef, sz))
CenterIndexedArray{T}(::UndefInitializer, sz::Vararg{<:Integer,N}) where {T,N} =
    CenterIndexedArray{T,N}(undef, sz...)
CenterIndexedArray{T}(::UndefInitializer, sz::NTuple{N,Integer}) where {T,N} =
    CenterIndexedArray{T,N}(undef, sz)

# This is the AbstractArray default, but do this just to be sure
Base.IndexStyle(::Type{A}) where {A<:CenterIndexedArray} = IndexCartesian()

Base.size(A::CenterIndexedArray) = size(A.data)
Base.axes(A::CenterIndexedArray) = map(SymRange, A.halfsize)

const SymAx = Union{SymRange, Base.Slice{SymRange}}
Base.axes(r::Base.Slice{SymRange}) = (r.indices,)

function Base.similar(A::CenterIndexedArray, ::Type{T}, inds::Tuple{SymAx,Vararg{SymAx}}) where T
    data = Array{T}(undef, map(length, inds))
    CenterIndexedArray(data)
end
function Base.similar(::Type{T}, inds::Tuple{SymAx, Vararg{SymAx}}) where T<:AbstractArray
    data = Array{eltype(T)}(undef, map(length, inds))
    CenterIndexedArray(data)
end

# This is incomplete: ideally we wouldn't need SymAx in the first slot
# as long as there was at least one SymAx.
function Base.similar(A::CenterIndexedArray, ::Type{T}, inds::Tuple{SymAx,Vararg{Union{Int,<:IdentityUnitRange,SymAx}}}) where T
    torange(n) = isa(n, Int) ? Base.OneTo(n) : n
    return OffsetArray{T}(undef, map(torange, inds))
end

function _halfsize(A::AbstractArray)
    all(isodd, size(A)) || error("Must have all-odd sizes")
    map(n->n>>UInt(1), size(A))
end

@inline function Base.getindex(A::CenterIndexedArray{T,N}, i::Vararg{Int,N}) where {T,N}
    @boundscheck checkbounds(A, i...)
    @inbounds val = A.data[map(offset, A.halfsize, i)...]
    val
end

Base.@propagate_inbounds Base.getindex(A::CenterIndexedArray{T,N,I}, i::Vararg{Int,N}) where {T,N,I<:AbstractInterpolation} =
    _getindex(A, i...)
Base.@propagate_inbounds Base.getindex(A::CenterIndexedArray{T,N,I}, i::Vararg{Number,N}) where {T,N,I<:AbstractInterpolation} =
    _getindex(A, i...)

@inline function _getindex(A::CenterIndexedArray{T,N,I}, i::Vararg{Number,N}) where {T,N,I<:AbstractInterpolation}
    @boundscheck checkbounds(A, i...)
    @inbounds val = A.data(map(offset, A.halfsize, i)...)
    val
end
Base.throw_boundserror(A::CenterIndexedArray, I) = (Base.@_noinline_meta; throw(BoundsError(A, I)))

offset(off, i) = off+i+1

@inline function Base.setindex!(A::CenterIndexedArray{T,N}, v, i::Vararg{Int,N}) where {T,N}
    @boundscheck checkbounds(A, i...)
    @inbounds A.data[map(offset, A.halfsize, i)...] = v
    v
end

Base.BroadcastStyle(::Type{<:CenterIndexedArray}) = Broadcast.ArrayStyle{CenterIndexedArray}()
function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{CenterIndexedArray}}, ::Type{ElType}) where ElType
    similar(Array{ElType}, axes(bc))
end

Base.parent(A::CenterIndexedArray) = A.data

function Base.showarg(io::IO, A::CenterIndexedArray, toplevel)
    print(io, "CenterIndexedArray(")
    Base.showarg(io, parent(A), false)
    print(io, ')')
    toplevel && print(io, " with eltype ", eltype(A))
end

include("deprecated.jl")

end  # module


text:
############################################################################################
# methods for AbstractSem
############################################################################################

function objective!(model::AbstractSemSingle, parameters)
    objective!(imply(model), parameters, model)
    return objective!(loss(model), parameters, model)
end

function gradient!(gradient, model::AbstractSemSingle, parameters)
    fill!(gradient, zero(eltype(gradient)))
    gradient!(imply(model), parameters, model)
    gradient!(gradient, loss(model), parameters, model)
end

function hessian!(hessian, model::AbstractSemSingle, parameters)
    fill!(hessian, zero(eltype(hessian)))
    hessian!(imply(model), parameters, model)
    hessian!(hessian, loss(model), parameters, model)
end

function objective_gradient!(gradient, model::AbstractSemSingle, parameters)
    fill!(gradient, zero(eltype(gradient)))
    objective_gradient!(imply(model), parameters, model)
    objective_gradient!(gradient, loss(model), parameters, model)
end

function objective_hessian!(hessian, model::AbstractSemSingle, parameters)
    fill!(hessian, zero(eltype(hessian)))
    objective_hessian!(imply(model), parameters, model)
    objective_hessian!(hessian, loss(model), parameters, model)
end

function gradient_hessian!(gradient, hessian, model::AbstractSemSingle, parameters)
    fill!(gradient, zero(eltype(gradient)))
    fill!(hessian, zero(eltype(hessian)))
    gradient_hessian!(imply(model), parameters, model)
    gradient_hessian!(gradient, hessian, loss(model), parameters, model)
end

function objective_gradient_hessian!(gradient, hessian, model::AbstractSemSingle, parameters)
    fill!(gradient, zero(eltype(gradient)))
    fill!(hessian, zero(eltype(hessian)))
    objective_gradient_hessian!(imply(model), parameters, model)
    return objective_gradient_hessian!(gradient, hessian, loss(model), parameters, model)
end

############################################################################################
# methods for SemFiniteDiff
############################################################################################

gradient!(gradient, model::SemFiniteDiff, par) =
    FiniteDiff.finite_difference_gradient!(gradient, x -> objective!(model, x), par)

hessian!(hessian, model::SemFiniteDiff, par) =
    FiniteDiff.finite_difference_hessian!(hessian, x -> objective!(model, x), par)

function objective_gradient!(
        gradient,
        model::SemFiniteDiff,
        parameters)
    gradient!(gradient, model, parameters)
    return objective!(model, parameters)
end

# other methods
function gradient_hessian!(
        gradient,
        hessian,
        model::SemFiniteDiff,
        parameters)
    gradient!(gradient, model, parameters)
    hessian!(hessian, model, parameters)
end

function objective_hessian!(hessian, model::SemFiniteDiff, parameters)
    hessian!(hessian, model, parameters)
    return objective!(model, parameters)
end

function objective_gradient_hessian!(
        gradient,
        hessian,
        model::SemFiniteDiff,
        parameters)
    hessian!(hessian, model, parameters)
    return objective_gradient!(gradient, model, parameters)
end

############################################################################################
# methods for SemLoss
############################################################################################

function objective!(loss::SemLoss, par, model)
    return mapreduce(
        fun_weight -> fun_weight[2]*objective!(fun_weight[1], par, model),
        +,
        zip(loss.functions, loss.weights)
        )
end

function gradient!(gradient, loss::SemLoss, par, model)
    for (lossfun, w) in zip(loss.functions, loss.weights)
        new_gradient = gradient!(lossfun, par, model)
        gradient .+= w*new_gradient
    end
end

function hessian!(hessian, loss::SemLoss, par, model)
    for (lossfun, w) in zip(loss.functions, loss.weights)
        hessian .+= w*hessian!(lossfun, par, model)
    end
end

function objective_gradient!(gradient, loss::SemLoss, par, model)
    return mapreduce(
        fun_weight -> objective_gradient_wrap_(gradient, fun_weight[1], par, model, fun_weight[2]),
        +,
        zip(loss.functions, loss.weights)
        )
end

function objective_hessian!(hessian, loss::SemLoss, par, model)
    return mapreduce(
        fun_weight -> objective_hessian_wrap_(hessian, fun_weight[1], par, model, fun_weight[2]),
        +,
        zip(loss.functions, loss.weights)
        )
end

function gradient_hessian!(gradient, hessian, loss::SemLoss, par, model)
    for (lossfun, w) in zip(loss.functions, loss.weights)
        new_gradient, new_hessian = gradient_hessian!(lossfun, par, model)
        gradient .+= w*new_gradient
        hessian .+= w*new_hessian
    end
end

function objective_gradient_hessian!(gradient, hessian, loss::SemLoss, par, model)
    return mapreduce(
        fun_weight -> objective_gradient_hessian_wrap_(gradient, hessian, fun_weight[1], par, model, fun_weight[2]),
        +,
        zip(loss.functions, loss.weights)
        )
end

# wrapper to update gradient/hessian and return objective value
function objective_gradient_wrap_(gradient, lossfun, par, model, w)
    new_objective, new_gradient = objective_gradient!(lossfun, par, model)
    gradient .+= w*new_gradient
    return w*new_objective
end

function objective_hessian_wrap_(hessian, lossfun, par, model, w)
    new_objective, new_hessian = objective_hessian!(lossfun, par, model)
    hessian .+= w*new_hessian
    return w*new_objective
end

function objective_gradient_hessian_wrap_(gradient, hessian, lossfun, par, model, w)
    new_objective, new_gradient, new_hessian = objective_gradient_hessian!(lossfun, par, model)
    gradient .+= w*new_gradient
    hessian .+= w*new_hessian
    return w*new_objective
end

############################################################################################
# methods for SemEnsemble
############################################################################################

function objective!(ensemble::SemEnsemble, par)
    return mapreduce(
        model_weight -> model_weight[2]*objective!(model_weight[1], par),
        +,
        zip(ensemble.sems, ensemble.weights)
        )
end

function gradient!(gradient, ensemble::SemEnsemble, par)
    fill!(gradient, zero(eltype(gradient)))
    for (model, w) in zip(ensemble.sems, ensemble.weights)
        gradient_new = similar(gradient)
        gradient!(gradient_new, model, par)
        gradient .+= w*gradient_new
    end
end

function hessian!(hessian, ensemble::SemEnsemble, par)
    fill!(hessian, zero(eltype(hessian)))
    for (model, w) in zip(ensemble.sems, ensemble.weights)
        hessian_new = similar(hessian)
        hessian!(hessian_new, model, par)
        hessian .+= w*hessian_new
    end
end

function objective_gradient!(gradient, ensemble::SemEnsemble, par)
    fill!(gradient, zero(eltype(gradient)))
    return mapreduce(
        model_weight -> objective_gradient_wrap_(gradient, model_weight[1], par, model_weight[2]),
        +,
        zip(ensemble.sems, ensemble.weights)
        )
end

function objective_hessian!(hessian, ensemble::SemEnsemble, par)
    fill!(hessian, zero(eltype(hessian)))
    return mapreduce(
        model_weight -> objective_hessian_wrap_(hessian, model_weight[1], par, model_weight[2]),
        +,
        zip(ensemble.sems, ensemble.weights)
        )
end

function gradient_hessian!(gradient, hessian, ensemble::SemEnsemble, par)
    fill!(gradient, zero(eltype(gradient)))
    fill!(hessian, zero(eltype(hessian)))
    for (model, w) in zip(ensemble.sems, ensemble.weights)

        new_gradient = similar(gradient)
        new_hessian = similar(hessian)

        gradient_hessian!(new_gradient, new_hessian, model, par)

        gradient .+= w*new_gradient
        hessian .+= w*new_hessian

    end
end

function objective_gradient_hessian!(gradient, hessian, ensemble::SemEnsemble, par)
    fill!(gradient, zero(eltype(gradient)))
    fill!(hessian, zero(eltype(hessian)))
    return mapreduce(
        model_weight -> objective_gradient_hessian_wrap_(gradient, hessian, model_weight[1], par, model, model_weight[2]),
        +,
        zip(ensemble.sems, ensemble.weights)
        )
end

# wrapper to update gradient/hessian and return objective value
function objective_gradient_wrap_(gradient, model::AbstractSemSingle, par, w)
    gradient_pre = similar(gradient)
    new_objective = objective_gradient!(gradient_pre, model, par)
    gradient .+= w*gradient_pre
    return w*new_objective
end

function objective_hessian_wrap_(hessian, model::AbstractSemSingle, par, w)
    hessian_pre = similar(hessian)
    new_objective = objective_hessian!(hessian_pre, model, par)
    hessian .+= w*new_hessian
    return w*new_objective
end

function objective_gradient_hessian_wrap_(gradient, hessian, model::AbstractSemSingle, par, w)
    gradient_pre = similar(gradient)
    hessian_pre = similar(hessian)
    new_objective = objective_gradient_hessian!(gradient_pre, hessian_pre, model, par)
    gradient .+= w*new_gradient
    hessian .+= w*new_hessian
    return w*new_objective
end

############################################################################################
# generic methods for loss functions
############################################################################################

function objective_gradient!(lossfun::SemLossFunction, par, model)
    objective = objective!(lossfun::SemLossFunction, par, model)
    gradient = gradient!(lossfun::SemLossFunction, par, model)
    return objective, gradient
end

function objective_hessian!(lossfun::SemLossFunction, par, model)
    objective = objective!(lossfun::SemLossFunction, par, model)
    hessian = hessian!(lossfun::SemLossFunction, par, model)
    return objective, hessian
end

function gradient_hessian!(lossfun::SemLossFunction, par, model)
    gradient = gradient!(lossfun::SemLossFunction, par, model)
    hessian = hessian!(lossfun::SemLossFunction, par, model)
    return gradient, hessian
end

function objective_gradient_hessian!(lossfun::SemLossFunction, par, model)
    objective = objective!(lossfun::SemLossFunction, par, model)
    gradient = gradient!(lossfun::SemLossFunction, par, model)
    hessian = hessian!(lossfun::SemLossFunction, par, model)
    return objective, gradient, hessian
end

# throw an error by default if gradient! and hessian! are not implemented

#= gradient!(lossfun::SemLossFunction, par, model) =
    throw(ArgumentError("gradient for $(typeof(lossfun).name.wrapper) is not available"))

hessian!(lossfun::SemLossFunction, par, model) =
    throw(ArgumentError("hessian for $(typeof(lossfun).name.wrapper) is not available")) =#

############################################################################################
# generic methods for imply
############################################################################################

function objective_gradient!(semimp::SemImply, par, model)
    objective!(semimp::SemImply, par, model)
    gradient!(semimp::SemImply, par, model)
    return nothing
end

function objective_hessian!(semimp::SemImply, par, model)
    objective!(semimp::SemImply, par, model)
    hessian!(semimp::SemImply, par, model)
    return nothing
end

function gradient_hessian!(semimp::SemImply, par, model)
    gradient!(semimp::SemImply, par, model)
    hessian!(semimp::SemImply, par, model)
    return nothing
end

function objective_gradient_hessian!(semimp::SemImply, par, model)
    objective!(semimp::SemImply, par, model)
    gradient!(semimp::SemImply, par, model)
    hessian!(semimp::SemImply, par, model)
    return nothing
end

############################################################################################
# Documentation
############################################################################################
"""
    objective!(model::AbstractSem, parameters)

Returns the objective value at `parameters`.
The model object can be modified.

# Implementation
To implement a new `SemImply` or `SemLossFunction` subtype, you need to add a method for
    objective!(newtype::MyNewType, parameters, model::AbstractSemSingle)

To implement a new `AbstractSem` subtype, you need to add a method for
    objective!(model::MyNewType, parameters)
"""
function objective! end

"""
    gradient!(gradient, model::AbstractSem, parameters)

Writes the gradient value at `parameters` to `gradient`.

# Implementation
To implement a new `SemImply` or `SemLossFunction` type, you can add a method for
    gradient!(newtype::MyNewType, parameters, model::AbstractSemSingle)

To implement a new `AbstractSem` subtype, you can add a method for
    gradient!(gradient, model::MyNewType, parameters)
"""
function gradient! end

"""
    hessian!(hessian, model::AbstractSem, parameters)

Writes the hessian value at `parameters` to `hessian`.

# Implementation
To implement a new `SemImply` or `SemLossFunction` type, you can add a method for
    hessian!(newtype::MyNewType, parameters, model::AbstractSemSingle)

To implement a new `AbstractSem` subtype, you can add a method for
    hessian!(hessian, model::MyNewType, parameters)
"""
function hessian! end


text:
module TontiTest
  using AlgebraicDynamics.Tonti
  using CombinatorialSpaces
  using GeometryBasics
  using Random
  using Test
  td = TontiDiagram(2, [:C =>:IP,  :dC =>:TP,  :ΔC=>:IL,
                         :∑ϕ=>:TS2, :ϕ=>:TL2]);

  s = EmbeddedDeltaSet2D{Bool, Point{3,Float64}}()
  points = [(0,0,0),(0,0,1),(0,1,0),(0,1,1)]

  add_vertices!(s, 4, point=points)
  glue_sorted_triangle!(s, 1,2,3)
  glue_sorted_triangle!(s, 2,3,4)

  sd = EmbeddedDeltaDualComplex2D{Bool, Float64, Point{3,Float64}}(s)
  subdivide_duals!(sd, Barycenter())
  star0 = ⋆(0,sd)

  @test gen_form(s, x->1) == [1,1,1,1]

  addSpace!(td, s)
  addTime!(td)
  const k = 0.003

  diffusion(diff) = ΔC -> diff*ΔC
  mass_conv() = ∑ϕ->∑ϕ

  addTransform!(td, s, [:ΔC], diffusion(k), [:ϕ])
  addTransform!(td, s, [:∑ϕ], mass_conv(), [:dC]);

  sim, data_syms = vectorfield(td, s)

  @test data_syms[1][1] == :C
  @test data_syms[1][2] == 4

  Random.seed!(42)
  u  = gen_form(s,x->rand())
  du = zeros(Float64,nv(s))

  tempu = copy(u)
  dt = 0.001
  for i in 10000
    sim(du, tempu, [],0)
    tempu .+= du
    # Check that mass is conserved
  end
  @test all(tempu .!= u)
  @test all(sum(star0*tempu)-sum(star0*u) < 1e-6)
end


text:
module ToggledLH

export @toggled, @tog_assert, toggle

assert_toggle() = true

"""
    @toggled(cond, expr)
    @toggled(expr)

Execute `expr` if `toggle` is ON and `cond` is true.
`cond` may be omitted.

`expr` can be several expressions as in
```
@toggled (1 == 1.0) (println("Start"); global x = 2; println("End"))
@toggled (1 == 1.0) begin
    println("Start");
    global x = 3;
    println("End");
end
```
"""
macro toggled(cond, expr)
    quote
        if assert_toggle()
            if $cond
                $(esc(expr))
            else
                nothing
            end
        else
            nothing
        end
        # :(assert_toggle() ?  ($cond  &&  $expr)  : nothing)
    end
end

macro toggled(expr)
    esc(:( @toggled(true, $expr) ))
end

"""
    @tog_assert(cond, text = nothing)

Toggled assert.
This is verbatim from `ToggleableAsserts.jl`
"""
macro tog_assert(cond, text = nothing)
    if text==nothing
        assert_stmt = esc(:(@assert $cond))
    else
        assert_stmt = esc(:(@assert $cond $text))
    end
    :(assert_toggle() ? $assert_stmt  : nothing)
end

const toggle_lock = ReentrantLock()

function toggle(enable::Bool)
    lock(toggle_lock) do
        @eval ToggledLH assert_toggle() = $enable
        on_or_off = enable ? "on." : "off."
        @info "ToggledLH turned "*on_or_off
    end
end

end # module


text:
using EquivalentCircuits, CSV

"""
Circuitidentification is done with the circuitevolution function. This function can take several keyword arguments,
allowing users to tune the circuitidentification if necessary. These can be found in the function documentation.
The measurements and frequencies should be provided either as a CSV file, or as arrays.

"""
#Measurements input as CSV file with as columns: real impedances, imaginary impedances and frequencies.
circuitevolution("example_measurements.csv")

#Example of input as arrays.
data = readdlm(raw"C:\Users\Dell\Documents\EquivalentCircuits.jl\example_measurements.csv",',')
measurements = data[:,1] .+ data[:,2]im
frequencies = data[:,3]
circuitevolution(measurements,frequencies)

"""
When the equivalent circuit to be used is known, its parameters can be fit using the parameteroptimisation function.

"""

parameteroptimisation("[[C1,R2]-R3,P4]", measurements, frequencies)


text:
using Unitful,Datax
using Test

@testset "Datax.jl" begin
    io = IOBuffer();

    # Basic printing
    Datax.printdata(io,"String","%.4g")
    @test String(take!(io)) == "String"
    Datax.printdata(io,4.1,"%.4g")
    @test String(take!(io)) == "\\num{4.1}"
    Datax.printdata(io,[1.0,2.4,3.0],"%.4g")
    @test String(take!(io)) == "\\numlist{1;2.4;3;}"

    # Print units
    Datax.printdata(io,(1.0,"\\meter"),"%.4g")
    @test String(take!(io)) == "\\SI{1}{\\meter}"
    Datax.printdata(io,([1,2,3],"\\kilogram"),"%.4g")
    @test String(take!(io)) == "\\SIlist{1;2;3;}{\\kilogram}"

    # Format strings
    Datax.printdata(io,pi,"%.3g")
    @test String(take!(io)) == "\\num{3.14}"
    Datax.printdata(io,("%.2g",pi),"%.4g")
    @test String(take!(io)) == "\\num{3.1}"
    Datax.printdata(io,("%.5g",pi,"\\liter"),"%.4g")
    @test String(take!(io)) == "\\SI{3.1416}{\\liter}"

    # Unitful
    Datax.printdata(io,2u"m","%.4g")
    @test String(take!(io)) == "\\SI{2}{\\meter}"
    Datax.printdata(io,[1.4,2.8]u"g/J","%.4g")
    @test String(take!(io)) == "\\SIlist{1.4;2.8;}{\\gram\\per\\joule}"

end


text:
# Test function

function testit(name, result, check::Function)
    if check(result)
        @info "$name successful, with result $result"
    else
        @warn "$name failed, with result $result"
    end
end

# Part 2

# spaces:
#    1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 # hallway
#    12, 13, 14, 15 # Room for As
#    16, 17, 18, 19 # Room for Bs
#    20, 21, 22, 23 # Room for Cs
#    24, 25, 26, 27 # Room for Ds

const AMPHIPODS = Dict(
    'A' => (:A1, :A2, :A3, :A4),
    'B' => (:B1, :B2, :B3, :B4),
    'C' => (:C1, :C2, :C3, :C4),
    'D' => (:D1, :D2, :D3, :D4)
)

const SIBLINGS = reduce(
    merge,
    [Dict(a => Set(filter(≠(a), AMPHIPODS[s])) for a in AMPHIPODS[s]) for s in 'A':'D']
)

const MULTIPLIER = reduce(
    merge,
    Dict(k => 10^(j-1) for k in AMPHIPODS[c]) for (j, c) in enumerate('A':'D')
)

function entry(k)
    for (i, c) in enumerate('A':'D')
        if k in AMPHIPODS[c]
            return (1 + 2 * i, 12 + 4 * (i-1))
        end
    end
end

exit(n) = div(n - 12, 4) |> m -> (3 + 2 * m, 12 + 4 * m)

# initial configuration for the test

#############
#...........#
###B#C#B#D###
  #D#C#B#A#
  #D#B#A#C#
  #A#D#C#A#
  #########

initial_positions_test = Dict(
    :A1 => 15,
    :A2 => 22,
    :A3 => 25,
    :A4 => 27,
    :B1 => 12,
    :B2 => 18,
    :B3 => 20,
    :B4 => 21,
    :C1 => 16,
    :C2 => 17,
    :C3 => 23,
    :C4 => 26,
    :D1 => 13,
    :D2 => 14,
    :D3 => 19,
    :D4 => 24
)

# An easy one for testing 4 x 1 + 4 x 10 + 2 x 1 = 46 points

#############
#...........#
###B#A#C#D###
  #A#B#C#D#
  #A#B#C#D#
  #A#B#C#D#
  #########

initial_positions_test_almostgameoverAB = Dict(
    :A1 => 16,
    :A2 => 13,
    :A3 => 14,
    :A4 => 15,
    :B1 => 12,
    :B2 => 17,
    :B3 => 18,
    :B4 => 19,
    :C1 => 20,
    :C2 => 21,
    :C3 => 22,
    :C4 => 23,
    :D1 => 24,
    :D2 => 25,
    :D3 => 26,
    :D4 => 27
)

initial_positions_test_almostgameoverCD = Dict(
    :A1 => 12,
    :A2 => 13,
    :A3 => 14,
    :A4 => 15,
    :B1 => 16,
    :B2 => 17,
    :B3 => 18,
    :B4 => 19,
    :C1 => 24,
    :C2 => 21,
    :C3 => 22,
    :C4 => 23,
    :D1 => 20,
    :D2 => 25,
    :D3 => 26,
    :D4 => 27
)

initial_positions_test_almostgameoverAADD = Dict(
    :A1 => 24,
    :A2 => 25,
    :A3 => 14,
    :A4 => 15,
    :B1 => 16,
    :B2 => 17,
    :B3 => 18,
    :B4 => 19,
    :C1 => 20,
    :C2 => 21,
    :C3 => 22,
    :C4 => 23,
    :D1 => 12,
    :D2 => 13,
    :D3 => 26,
    :D4 => 27
)

# initial configuration for the challenge

#############
#...........#
###B#D#C#A###
  #D#C#B#A#
  #D#B#A#C#
  #C#D#B#A#
  #########

initial_positions_challenge = Dict(
    :A1 => 22,
    :A2 => 24,
    :A3 => 25,
    :A4 => 27,
    :B1 => 12,
    :B2 => 18,
    :B3 => 21,
    :B4 => 23,
    :C1 => 15,
    :C2 => 17,
    :C3 => 20,
    :C4 => 26,
    :D1 => 13,
    :D2 => 14,
    :D3 => 16,
    :D4 => 19
)

# Part 2

gameover(positions) = all(
    u -> Set(positions[a] for a in AMPHIPODS[last(u)]) == Set((8 + 4*first(u) + j for j in 0:3)),
    enumerate('A':'D')
)

function allowedmoves(positions, k)
    pk = positions[k]
    mvs = []
    enh, enr = entry(k) # entry point in hallway and in room
    exh, exr = exit(positions[k]) # exit point in hallway and in room
    if pk ≤ 11
        for i in 3:-1:0
            if isdisjoint(values(positions), union(enh:pk-1, pk+1:enh, enr:enr+i)) &&
                    enr+i+1:enr+3 ⊆ (positions[sk] for sk in SIBLINGS[k])
                push!(mvs, (enr+i, length(enh:pk-1) + length(pk+1:enh) + i + 1))
                return mvs
            end
        end
    elseif pk ∉ enr:enr+3 || pk+1:enr+3 ⊈ (positions[sk] for sk in SIBLINGS[k])
        for i in 3:-1:0
            if isdisjoint(values(positions), union(exr:pk-1, exh:enh, enh:exh, enr:enr+i)) &&
                enr+i+1:enr+3 ⊆ (positions[sk] for sk in SIBLINGS[k])
                push!(mvs, (enr+i, length(exr:pk-1) + length(exh:enh) + length(enh:exh) + i + 1))
                return mvs
            end
        end
        for p in (1, 2, 4, 6, 8, 10, 11)
            if isdisjoint(values(positions), union(exr:pk-1, exh:p, p:exh))
                push!(mvs, (p, length(exr:pk-1) + length(exh:p) + length(p:exh)))
            end
        end
    end
    return mvs
end

function move(positions, points = 0, minpoints = typemax(Int),
    moves = Vector{Tuple{Symbol, Pair{Int, Int}, Int}}(),
    minmoves = Vector{Tuple{Symbol, Pair{Int, Int}, Int}}()
)
    if points < minpoints && gameover(positions)
        minpoints = points
        minmoves = moves
    else
        for k in keys(positions)
            for (pos, steps) in allowedmoves(positions, k)
                step_pts = MULTIPLIER[k] * steps
                newpoints = points + step_pts
                if newpoints < minpoints
                    newpositions = copy(positions)
                    newpositions[k] = pos
                    newmoves = [moves; (k, positions[k] => pos, step_pts)]
                    minmoves, minpoints = move(newpositions, newpoints, minpoints, newmoves, minmoves)
                end
            end
        end
    end
    return minmoves, minpoints
end

function part2(initial_positions)
    move(initial_positions)
end

testit("Test Part 2 AB", part2(initial_positions_test_almostgameoverAB)[2], ==(46))
testit("Test Part 2 CD", part2(initial_positions_test_almostgameoverCD)[2], ==(4600))
testit("Test Part 2 AADD", part2(initial_positions_test_almostgameoverAADD)[2], ==(18022))
testit("Test Part 2", part2(initial_positions_test)[2], ==(44169))
testit("Challenge Part 2", part2(initial_positions_challenge)[2], ==(15322))

# Visualization

# A better approach is to make `Board` a subtype of AbstractVector, with length 27, and
# just set directly the element at a given position on the board, instead of rebuilding
# the entire board after each move.

struct Board
    mat::Matrix{Char}
end

function Board(positions::Dict = Dict())
    b = Board(fill('#', 7, 13))
    b.mat[4:end, 1:2] .= ' '
    b.mat[4:end, 12:13] .= ' '
    b.mat[2, 2:12] .= '.'
    b.mat[3:6, [4, 6, 8, 10]] .= '.'
    for (k, v) in positions
        if 1 ≤ v ≤ 11
            b.mat[2, 1 + v] = first(string(k))
        else
            d, r = divrem(v - 12, 4)
            b.mat[3 + r, 4 + 2*d] = first(string(k))
        end
    end
    return b
end

function Base.show(io::IO, b::Board)
    for bi in eachrow(b.mat)
        println(join(bi))
    end
end

function showsolution(positions, moves)
    b = Board(positions)
    println(b)
    sleep(5)
    for p in moves
        positions[p[1]] = last(p[2])
        b = Board(positions)
        println(b)
        sleep(1)
    end
    sleep(5)
end

# positions = copy(initial_positions_test)
# minmoves, minpoints = part2(positions)
# showsolution(positions, minmoves)


text:
module EulerCharacteristic

import AnnealingAPI

include("calc.jl")
include("tracker.jl")

export euler_characteristic,
    EulerTracker

end # module


text:
# Deps
using CpelTdm

# Set dir of test data
pack_dir = dirname(dirname(pathof(CpelTdm)))
data_dir = "$(pack_dir)/test/"

# Set paths
fa = "$(data_dir)/fasta/ref.fa"
bed = "$(data_dir)/bed/roi.bed"
bams1 = "$(data_dir)/bam/g1_s" .* string.(collect(1:8)) .* ".bam"
bams2 = "$(data_dir)/bam/g2_s" .* string.(collect(1:8)) .* ".bam"

# Output
outprefix = "cpeltdm"
outdir = "$(data_dir)/cpeltdm/"

# Run analysis
rm(outdir;force=true,recursive=true)
cpel_tdm(bams1,bams2,bed,fa,outdir,outprefix;matched=false)


text:
"""
    MeasurementOperator{M,G}

General measurement operator. `M` is the number of wires affected by the measurement operator, `G` is the actual operator type (gate or sparse matrix).
```
julia> m = MeasurementOperator(X, (1,))
MeasurementOperator{1,XGate}(XGate(), (1,))

julia> m = MeasurementOperator([0 1; 1 0], (1,))
MeasurementOperator{1,SparseArrays.SparseMatrixCSC{Complex{Float64},Int64}}(
  [2, 1]  =  1.0+0.0im
  [1, 2]  =  1.0+0.0im, (1,))
```
"""
struct MeasurementOperator{M,G}
    operator::G
    iwire::NTuple{M,Int}
    function MeasurementOperator(g::G, iwire::NTuple{M,Integer}) where {M,G <: AbstractGate}
        M == num_wires(g) || error("CircuitGate affecting $(num_wires(g)) wires given, `iwire` of length $(length(iwire)) provided")
        g == adjoint(g) || error("Measurement operator must be Hermitian.")
        new{M,G}(g, iwire)
    end

    function MeasurementOperator(m::G, iwire::NTuple{M,Integer}) where {M,G <: AbstractMatrix}
        d = 2
        size(m) == (d^M, d^M) || error("Measurement operator must be a $(d^M) × $(d^M) matrix.")
        m ≈ adjoint(m) || error("Measurement operator must be Hermitian.")
        new{M,SparseMatrixCSC{Complex{Float64},Int}}(sparse(m), iwire)
    end

    function MeasurementOperator(m, iwire::Integer...)
        MeasurementOperator(m, iwire)
    end
end

mop(operator, iwire) = MeasurementOperator(operator, iwire)

@memoize function num_wires(m::MeasurementOperator)
    length(m.iwire)
end
@memoize function Base.size(m::MeasurementOperator)
    length(m.iwire)
end

@memoize function check_commute(ops::Vector{<:MeasurementOperator})
    for (i, operator) in enumerate(ops)
        for operator2 in ops[i+1:end]
            if !iscommuting(operator, operator2)
                return false
            end
        end
    end
    true
end

"""
    sparse_matrix(m::MeasurementOperator{M,G}, N::Integer=0) where {M,G}

returns matrix representation of a `MeasurementOperator{M,G}` object that can applied to a state vector of `N` qubits.
"""
function sparse_matrix(m::MeasurementOperator{M,G}, N::Integer=0) where {M,G}
    # convert to array
    iwire = collect(m.iwire)

    if N == 0
        N = maximum([iwire' M])
    else
        N >= maximum(iwire) || error("CircuitGate applied to iwires, $iwire. Input circuit size `N` is $N")
    end

    gmat = sparse_matrix(m.operator)

    distribute_to_wires(gmat, iwire, N, M)
end

"""
    req_wires(m::MeasurementOperator{M,G})

Minimum number of required qubit wires in a circuit to host the measurement operator `m`.
```jldoctest
julia> meas = mop(X, (4,));
julia> req_wires(meas)
4
```
"""
req_wires(m::MeasurementOperator{M,G}) where {M,G} = maximum(m.iwire)