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

using Markdown
using InteractiveUtils

# ╔═╡ da5ad68e-1757-11eb-37ee-5ff74db1f6ca
using LinearAlgebra, Images, Plots

# ╔═╡ 6437c7fa-178e-11eb-3966-31c2f7ea3408
md"# Deconstructing PCA using a DR perspective

Results from the paper by [Narasimhan and Bhatt (2015)](https://www.sciencedirect.com/science/article/abs/pii/S0098135415000873)
"

# ╔═╡ 3a7ed398-1790-11eb-0662-07879f0beaf1
begin
	img = load("Flow_process.png")
	plot(img, axis=false)
end

# ╔═╡ 5661af4e-1791-11eb-0ded-3541c01717b3
md"### Constraint matrix"

# ╔═╡ 73ed6bc6-1759-11eb-024c-83cfae5b1b81
A = [1 1 -1 0 0 0;
	0 0 1 -1 0 0;
	0 0 0 1 -1 -1;
	0 -1 0 0 0 1]

# ╔═╡ 2fdcdcc0-175c-11eb-2a4e-5d913ceb7f52
md"## Measurements and Fluctuations
- Both, $F_1$ and $F_2$ are measured with fluctuations
- Base value for both is $10$
- Fluctuation are $\mathcal{N}(0,1)$ and $\mathcal{N}(0,2^2)$, respectively
- A total of $1000$ data points are available for the steady state flow balance equations
"

# ╔═╡ 1fb5fc10-175b-11eb-139c-8126aaca3a00
begin
	N = 1000
	F1_base, F1_fluct = 10, 1
	F2_base, F2_fluct = 10, 2
	F1ₜ = F1_base .+ (F1_fluct) * randn(N)
	F2ₜ = F2_base .+ (F2_fluct) * randn(N)
end;

# ╔═╡ 7225218e-1791-11eb-2c02-617ad622b901
md"Breaking the constraint matrix into the measured $A_x$ and the unmeasured $A_u$ matrices."

# ╔═╡ d478db56-175e-11eb-31b6-d32eadc31c46
begin
	Aₓ = A[:,1:2]
	Aᵤ = A[:,3:end]
end;

# ╔═╡ bca0ca10-1791-11eb-3de6-b3fc37014809
md"For deriving the true values of the other variables, we need to substitute the true values of $F_1$ and $F_2$ (all $1000$ values) and obtain the corresponding values of the derived/dependent variables, $F_3$, $F_4$, $F_5$ and $F_6$.

$A_x [F_1, F_2]^T + A_u [F_3, F_4, F_5, F_6]^T = 0$
or

$A_u [F_3, F_4, F_5, F_6]^T = - A_x [F_1, F_2]^T$

"

# ╔═╡ ac79bc0a-175f-11eb-1205-df0aab8c88d6
b = -Aₓ*[F1ₜ F2ₜ]'

# ╔═╡ ea3993f8-175f-11eb-308a-11c6e6d08303
dependent_vars = Aᵤ\b

# ╔═╡ 7b461812-1760-11eb-1ca4-b9d591c18345
begin
	F3ₜ = dependent_vars[1,:]
	F4ₜ = dependent_vars[2, :]
	F5ₜ = dependent_vars[3, :]
	F6ₜ = dependent_vars[4, :]
end;

# ╔═╡ c06ad3ec-1760-11eb-178e-397e40a42dc4
size(F3ₜ), size(F4ₜ), size(F5ₜ), size(F6ₜ)

# ╔═╡ 78a1293e-1761-11eb-1d7c-67574000a6a4
md"#### Verify the true values"

# ╔═╡ 635e0f74-1761-11eb-0fc9-7b1c6bfc24df
X = [F1ₜ, F2ₜ, F3ₜ, F4ₜ, F5ₜ, F6ₜ]

# ╔═╡ 0590fbea-1761-11eb-0b0a-ad1cfe527582
A * X

# ╔═╡ c26af9fa-1813-11eb-34e7-c1b7bbdc230b
all(x -> x < 1e-12, reduce(hcat, A*X))

# ╔═╡ 1f53b0c2-175c-11eb-0ad4-3b99740221e2
md"### SDE"

# ╔═╡ 089c4dde-175d-11eb-3702-3f3ad6fbc824
begin
	sde = [0.1, 0.08, 0.15, 0.2, 0.18, 0.1]
	Σ = diagm(sde.^2)
end

# ╔═╡ 9e3ed1e0-1762-11eb-3b40-0de3fe9e6c1f
md"#### Simulating measurements around the true values"

# ╔═╡ ac38bda6-1762-11eb-3c95-51a83eaa5e3d
Y = X + [sd * randn(N) for sd in sde]

# ╔═╡ 288f1b4e-175d-11eb-3c8b-1554fbd84e14
md"## Applying DR"

# ╔═╡ 45f72cbc-175d-11eb-0055-01c9815689c0
md"### LA route

- Find the $W$ matrix
- Reconciled estimates, $X_e = WY$
- Updated covariance, $Σ_e = WΣW^T$
"

# ╔═╡ 5002007e-175d-11eb-2391-8fd69f92794a
W = I -Σ * A'* inv(A * Σ * A') * A

# ╔═╡ b44187ea-175e-11eb-2921-25e3c2345ca1
size(W)

# ╔═╡ c1e879ce-175e-11eb-0944-61396ca26a54
md"#### Estimates"

# ╔═╡ be181a9a-1761-11eb-1989-c1f15900ffdf
Xₑ = W * Y

# ╔═╡ 6ef7277e-176d-11eb-00ca-c95bc64fd8ed
md"##### Error covariance matrix"

# ╔═╡ 6a1c0092-1765-11eb-0ca7-7d6d5ee76b66
Σₑ = W * Σ * W'

# ╔═╡ 7d0befd4-176d-11eb-2881-69eea63cb4e9
diag(Σₑ)

# ╔═╡ 918e9662-1765-11eb-30b1-030cd8a03991
diag(Σₑ) .< diag(Σ)

# ╔═╡ 256ad438-1764-11eb-0783-874a2b0a8486
md"#### Verifying whether the estimates satisfy the constraints"

# ╔═╡ 19fb961c-1764-11eb-1eee-b9300189ab65
A * Xₑ

# ╔═╡ f09e3672-1764-11eb-3b4e-bff973f423b5
begin
	residual = reduce(hcat, A*Xₑ);
	size(residual)
end

# ╔═╡ 7ceeb166-1764-11eb-11ca-851e9d7fbdf6
all(x -> isapprox(x,0,atol=1e-12), residual)

# ╔═╡ 626ffc8a-1766-11eb-0830-4f2533596b99
md"#### Effect of DR"

# ╔═╡ 72f4527c-1766-11eb-3f89-03fd973eeaf3
begin
	rmse_before = [a.^2 for a in Y - X]
	rmse_before = [√sum(a) for a in rmse_before]

	rmse_after = [a.^2 for a in Xₑ - X]
	rmse_after = [√sum(a) for a in rmse_after]
end;

# ╔═╡ ba9ccd86-1767-11eb-0924-6154f763fa2e
rmse_before, rmse_after

# ╔═╡ a504e824-1768-11eb-0479-8b0fd218b00f
rmse_before .> rmse_after

# ╔═╡ 395f44de-1781-11eb-2f57-9326eb43c0c4
md"### PCA route

**Wrong approach**
"

# ╔═╡ 42e79e2a-1781-11eb-07c3-5595346dc7da
begin
	Y_mat = reduce(hcat, Y)'./ √N;
	size(Y_mat)
end

# ╔═╡ bb6cd70a-1814-11eb-2beb-19d04c8aef8b
md"The error in the above implementation is the violation of the assumption in Section 4.
- The errors must be i.i.d.
- So, we may express the error-covariance as: $Σ = σ^2I$
Here, the identically distributed assumption is violated since the simulated measurement data is generated using different variances, as given by the error-covariance matrix $Σ$ above, which is a diagonal matrix with diagonal:

$[s for s in sde]
"

# ╔═╡ fc20aefa-1816-11eb-0700-efaa98fa8fdb
md"
Therefore, we need to turn to Section 4.1 and 4.2, where we include the effect of a known error-covariance matrix, with different diagonal entries. The trick is to transform (scale) the measurement matrix ($Y$) by the *cholesky* decomposed lower triangular matrix of Σ.
"

# ╔═╡ 4fa77514-1816-11eb-39bf-fb7c4ec67469
begin
	C = cholesky(Σ)
	L = C.L
end

# ╔═╡ 1df550d0-1815-11eb-2d14-c709b0d64eac
begin
	Y_mat2 = reduce(hcat, Y)'./ √N;
	Y_mat_scaled = inv(L) * Y_mat2
	size(Y_mat)
end

# ╔═╡ b7db36f8-1784-11eb-18e7-9df8eff09198
U, S, Vt = svd(Y_mat_scaled)

# ╔═╡ 35266c1c-1786-11eb-1dfc-89f1fdfc9f60
plot(S, marker="0", label="Singular values")

# ╔═╡ 946dd5a8-1785-11eb-3661-9734f2aaa832
threshold_var = 0.01;

# ╔═╡ 13c16910-1785-11eb-07b3-014ff0eb15bc
begin
	rel_var = 100*cumsum(reverse(S.^2))/sum(S.^2)
	plot(rel_var)
end

# ╔═╡ cb1d5724-1792-11eb-05ba-9bf4030ef2c0
rel_var

# ╔═╡ 66d7c00a-1787-11eb-2657-8f9e938a2346
m = findlast( rel_var .< threshold_var )

# ╔═╡ 2bc8d784-1788-11eb-3e11-99ee9641592b
begin
	n = size(A,2)
	p = n - m
end

# ╔═╡ c8986b5c-1787-11eb-2840-1d97394cdd34
md"### Denoised estimates
$X_e = √N (U_1  S_1  V_1 ^T)$
"

# ╔═╡ 96736630-1788-11eb-0528-c165152d4b6a
begin
	X_denoised_s = √N * U[:,1:p] * diagm(S[1:p]) * Vt[:,1:p]'
	X_denoised = L * X_denoised_s
end

# ╔═╡ 35913aa4-1819-11eb-1316-83e9fca7344c
begin
	rmse_PCA = (X_denoised - reduce(hcat, X)').^2
	rmse_PCA = sqrt.(sum(rmse_PCA, dims=2))
	rmse_PCA = reduce(vcat,rmse_PCA)
end

# ╔═╡ 023cd9fc-181d-11eb-3344-03a6ac092538
rmse_before, rmse_after

# ╔═╡ 38078eb0-181d-11eb-163e-7f72c7963b17
md"As can be seen, the RMSE in the estimates w.r.t. the true values has reduced by employing the PCA approach. This is even closely comparable to the DR estimates' RMSE.
"

# ╔═╡ 15dec576-178a-11eb-1063-0fd1e64d75db
md"### Recovered Model

$A_m = {U_2}^T$
"

# ╔═╡ b68f2b3a-1789-11eb-18ed-c7a7668db14c
begin
	U₂s = U[:,p+1:end]
	Aₘ = U₂s' * inv(L)
end

# ╔═╡ 88691d56-178a-11eb-178d-4b85c730b107
md"*Ofcourse, this estimated model is rotated w.r.t. to the true constraint matrix*

$A_m = QA$
But it is orthogonal to the measured estimates
"

# ╔═╡ 64ec59c2-1789-11eb-1b29-f7fe37cc910f
out_denoised = Aₘ*X_denoised

# ╔═╡ e1dc9cd4-178a-11eb-00c5-314ca4d343f7
all(x -> isapprox(x, 0, atol=1e-12), out_denoised)

# ╔═╡ d4299626-178c-11eb-2a26-d35ce8d672c9
md"**But does this satisfy the estimates obtained from LA based DR?**"

# ╔═╡ e1607cf6-178c-11eb-30c3-c50e55e2f919
out_der_model = Aₘ * Xₑ

# ╔═╡ ea6b8804-178c-11eb-310a-27b5f7374079
out_model = reduce(hcat, out_der_model)'

# ╔═╡ 0b6db504-178d-11eb-33cd-e13fdb6aa7dd
maximum(abs.(out_model)), minimum(abs.(out_model))

# ╔═╡ 5cd79d60-178d-11eb-0793-3329c24fcecb
md"*As is observed, the derived model from data alone, does not gel well with traditional DR estimates, but the residuals are still quite small*"

# ╔═╡ Cell order:
# ╟─6437c7fa-178e-11eb-3966-31c2f7ea3408
# ╠═da5ad68e-1757-11eb-37ee-5ff74db1f6ca
# ╟─3a7ed398-1790-11eb-0662-07879f0beaf1
# ╟─5661af4e-1791-11eb-0ded-3541c01717b3
# ╠═73ed6bc6-1759-11eb-024c-83cfae5b1b81
# ╟─2fdcdcc0-175c-11eb-2a4e-5d913ceb7f52
# ╠═1fb5fc10-175b-11eb-139c-8126aaca3a00
# ╟─7225218e-1791-11eb-2c02-617ad622b901
# ╠═d478db56-175e-11eb-31b6-d32eadc31c46
# ╟─bca0ca10-1791-11eb-3de6-b3fc37014809
# ╠═ac79bc0a-175f-11eb-1205-df0aab8c88d6
# ╠═ea3993f8-175f-11eb-308a-11c6e6d08303
# ╠═7b461812-1760-11eb-1ca4-b9d591c18345
# ╠═c06ad3ec-1760-11eb-178e-397e40a42dc4
# ╟─78a1293e-1761-11eb-1d7c-67574000a6a4
# ╠═635e0f74-1761-11eb-0fc9-7b1c6bfc24df
# ╠═0590fbea-1761-11eb-0b0a-ad1cfe527582
# ╠═c26af9fa-1813-11eb-34e7-c1b7bbdc230b
# ╟─1f53b0c2-175c-11eb-0ad4-3b99740221e2
# ╠═089c4dde-175d-11eb-3702-3f3ad6fbc824
# ╟─9e3ed1e0-1762-11eb-3b40-0de3fe9e6c1f
# ╠═ac38bda6-1762-11eb-3c95-51a83eaa5e3d
# ╟─288f1b4e-175d-11eb-3c8b-1554fbd84e14
# ╟─45f72cbc-175d-11eb-0055-01c9815689c0
# ╠═5002007e-175d-11eb-2391-8fd69f92794a
# ╠═b44187ea-175e-11eb-2921-25e3c2345ca1
# ╟─c1e879ce-175e-11eb-0944-61396ca26a54
# ╠═be181a9a-1761-11eb-1989-c1f15900ffdf
# ╟─6ef7277e-176d-11eb-00ca-c95bc64fd8ed
# ╠═6a1c0092-1765-11eb-0ca7-7d6d5ee76b66
# ╠═7d0befd4-176d-11eb-2881-69eea63cb4e9
# ╠═918e9662-1765-11eb-30b1-030cd8a03991
# ╟─256ad438-1764-11eb-0783-874a2b0a8486
# ╠═19fb961c-1764-11eb-1eee-b9300189ab65
# ╠═f09e3672-1764-11eb-3b4e-bff973f423b5
# ╠═7ceeb166-1764-11eb-11ca-851e9d7fbdf6
# ╟─626ffc8a-1766-11eb-0830-4f2533596b99
# ╠═72f4527c-1766-11eb-3f89-03fd973eeaf3
# ╠═ba9ccd86-1767-11eb-0924-6154f763fa2e
# ╠═a504e824-1768-11eb-0479-8b0fd218b00f
# ╟─395f44de-1781-11eb-2f57-9326eb43c0c4
# ╠═42e79e2a-1781-11eb-07c3-5595346dc7da
# ╟─bb6cd70a-1814-11eb-2beb-19d04c8aef8b
# ╟─fc20aefa-1816-11eb-0700-efaa98fa8fdb
# ╠═4fa77514-1816-11eb-39bf-fb7c4ec67469
# ╠═1df550d0-1815-11eb-2d14-c709b0d64eac
# ╠═b7db36f8-1784-11eb-18e7-9df8eff09198
# ╠═35266c1c-1786-11eb-1dfc-89f1fdfc9f60
# ╠═946dd5a8-1785-11eb-3661-9734f2aaa832
# ╠═13c16910-1785-11eb-07b3-014ff0eb15bc
# ╠═cb1d5724-1792-11eb-05ba-9bf4030ef2c0
# ╠═66d7c00a-1787-11eb-2657-8f9e938a2346
# ╠═2bc8d784-1788-11eb-3e11-99ee9641592b
# ╟─c8986b5c-1787-11eb-2840-1d97394cdd34
# ╠═96736630-1788-11eb-0528-c165152d4b6a
# ╠═35913aa4-1819-11eb-1316-83e9fca7344c
# ╠═023cd9fc-181d-11eb-3344-03a6ac092538
# ╟─38078eb0-181d-11eb-163e-7f72c7963b17
# ╟─15dec576-178a-11eb-1063-0fd1e64d75db
# ╠═b68f2b3a-1789-11eb-18ed-c7a7668db14c
# ╟─88691d56-178a-11eb-178d-4b85c730b107
# ╠═64ec59c2-1789-11eb-1b29-f7fe37cc910f
# ╠═e1dc9cd4-178a-11eb-00c5-314ca4d343f7
# ╟─d4299626-178c-11eb-2a26-d35ce8d672c9
# ╠═e1607cf6-178c-11eb-30c3-c50e55e2f919
# ╠═ea6b8804-178c-11eb-310a-27b5f7374079
# ╠═0b6db504-178d-11eb-33cd-e13fdb6aa7dd
# ╟─5cd79d60-178d-11eb-0793-3329c24fcecb

max_stars_repo_name: samiit/DR-PCA

max_issues_repo_name: samiit/DR-PCA

max_forks_repo_name: samiit/DR-PCA


text:
using DiffEqSensitivity,OrdinaryDiffEq, ParameterizedFunctions,
      RecursiveArrayTools, DiffEqBase, ForwardDiff, Calculus, QuadGK,
      LinearAlgebra
using Test

fb = @ode_def_bare begin
  dx = a*x - b*x*y
  dy = -c*y + x*y
end a b c
f = @ode_def begin
  dx = a*x - b*x*y
  dy = -c*y + x*y
end a b c

p = [1.5,1.0,3.0]
prob = ODEProblem(f,[1.0;1.0],(0.0,10.0),p)
sol = solve(prob,Vern9(),abstol=1e-14,reltol=1e-14)
probb = ODEProblem(fb,[1.0;1.0],(0.0,10.0),p)
solb = solve(probb,Vern9(),abstol=1e-14,reltol=1e-14)

# Do a discrete adjoint problem
println("Calculate discrete adjoint sensitivities")
t = 0.0:0.5:10.0 # TODO: Add end point handling for callback
# g(t,u,i) = (1-u)^2/2, L2 away from 1
function dg(out,u,p,t,i)
  (out.=2.0.-u)
end

easy_res = adjoint_sensitivities(sol,Vern9(),dg,t,abstol=1e-14,
                                 reltol=1e-14,iabstol=1e-14,ireltol=1e-12)
easy_res2 = adjoint_sensitivities(solb,Vern9(),dg,t,abstol=1e-14,
                                 reltol=1e-14,iabstol=1e-14,ireltol=1e-12,sensealg=SensitivityAlg(quad=true,backsolve=false))
easy_res3 = adjoint_sensitivities(solb,Vern9(),dg,t,abstol=1e-14,
                                  reltol=1e-14,iabstol=1e-14,ireltol=1e-12,sensealg=SensitivityAlg(quad=false,backsolve=false))
easy_res4 = adjoint_sensitivities(solb,Vern9(),dg,t,abstol=1e-14,
                                  reltol=1e-14,iabstol=1e-14,ireltol=1e-12,sensealg=SensitivityAlg(backsolve=true))
easy_res5 = adjoint_sensitivities(sol,Kvaerno5(nlsolve=NLAnderson(), smooth_est=false),dg,t,abstol=1e-14,
                                 reltol=1e-14,iabstol=1e-14,ireltol=1e-12,sensealg=SensitivityAlg(backsolve=true))
easy_res6 = adjoint_sensitivities(solb,Vern9(),dg,t,abstol=1e-14,
                                  reltol=1e-14,iabstol=1e-14,ireltol=1e-12,
                                  sensealg=SensitivityAlg(checkpointing=true,quad=true),
                                  checkpoints=sol.t[1:5:end])
easy_res7 = adjoint_sensitivities(solb,Vern9(),dg,t,abstol=1e-14,
                                  reltol=1e-14,iabstol=1e-14,ireltol=1e-12,
                                  sensealg=SensitivityAlg(checkpointing=true,quad=false),
                                  checkpoints=sol.t[1:5:end])

adj_prob = ODEAdjointProblem(sol,dg,t)
adj_sol = solve(adj_prob,Vern9(),abstol=1e-14,reltol=1e-14)
integrand = AdjointSensitivityIntegrand(sol,adj_sol)
res,err = quadgk(integrand,0.0,10.0,atol=1e-14,rtol=1e-12)

@test isapprox(res, easy_res, rtol = 1e-10)
@test isapprox(res, easy_res2, rtol = 1e-10)
@test isapprox(res, easy_res3, rtol = 1e-10)
@test isapprox(res, easy_res4, rtol = 1e-10)
@test isapprox(res, easy_res5, rtol = 1e-9)
@test isapprox(res, easy_res6, rtol = 1e-9)
@test isapprox(res, easy_res7, rtol = 1e-9)

println("Calculate adjoint sensitivities from autodiff & numerical diff")
function G(p)
  tmp_prob = remake(prob,u0=convert.(eltype(p),prob.u0),p=p)
  sol = solve(tmp_prob,Vern9(),abstol=1e-14,reltol=1e-14,saveat=t)
  A = convert(Array,sol)
  sum(((2 .- A).^2)./2)
end
G([1.5,1.0,3.0])
res2 = ForwardDiff.gradient(G,[1.5,1.0,3.0])
res3 = Calculus.gradient(G,[1.5,1.0,3.0])

using Flux
res4 = Flux.Tracker.gradient(G,[1.5,1.0,3.0])[1]

using ReverseDiff
res5 = ReverseDiff.gradient(G,[1.5,1.0,3.0])

@test norm(res' .- res2) < 1e-8
@test norm(res' .- res3) < 1e-6
@test norm(res' .- res4) < 1e-6
@test norm(res' .- res5) < 1e-6

# Do a continuous adjoint problem

# Energy calculation
g(u,p,t) = (sum(u).^2) ./ 2
# Gradient of (u1 + u2)^2 / 2
function dg(out,u,p,t)
  out[1]= u[1] + u[2]
  out[2]= u[1] + u[2]
end

adj_prob = ODEAdjointProblem(sol,g,nothing,dg)
adj_sol = solve(adj_prob,Tsit5(),abstol=1e-14,reltol=1e-10)
integrand = AdjointSensitivityIntegrand(sol,adj_sol)
res,err = quadgk(integrand,0.0,10.0,atol=1e-14,rtol=1e-10)

println("Test the `adjoint_sensitivities` utility function")
easy_res = adjoint_sensitivities(sol,Tsit5(),g,nothing,dg,abstol=1e-14,
                                 reltol=1e-14,iabstol=1e-14,ireltol=1e-12)
easy_res2 = adjoint_sensitivities(sol,Tsit5(),g,nothing,dg,abstol=1e-14,
                                  reltol=1e-14,iabstol=1e-14,ireltol=1e-12,
                                  sensealg=SensitivityAlg(quad=false))
easy_res3 = adjoint_sensitivities(sol,Tsit5(),g,nothing,abstol=1e-14,
                                 reltol=1e-14,iabstol=1e-14,ireltol=1e-12)
easy_res4 = adjoint_sensitivities(sol,Tsit5(),g,nothing,abstol=1e-14,
                                  reltol=1e-14,iabstol=1e-14,ireltol=1e-12,
                                  sensealg=SensitivityAlg(autodiff=false))
easy_res5 = adjoint_sensitivities(sol,Tsit5(),g,nothing,abstol=1e-14,
                                  reltol=1e-14,iabstol=1e-14,ireltol=1e-12,
                                  sensealg=SensitivityAlg(checkpointing=true))
@test norm(easy_res .- res) < 1e-8
@test norm(easy_res2 .- res) < 1e-8
@test norm(easy_res3 .- res) < 1e-8
@test norm(easy_res4 .- res) < 1e-8
@test norm(easy_res5 .- res) < 1e-8

println("Calculate adjoint sensitivities from autodiff & numerical diff")
function G(p)
  tmp_prob = remake(prob,u0=eltype(p).(prob.u0),p=p,
                    tspan=eltype(p).(prob.tspan))
  sol = solve(tmp_prob,Vern9(),abstol=1e-14,reltol=1e-14)
  res,err = quadgk((t)-> (sum(sol(t)).^2)./2,0.0,10.0,atol=1e-14,rtol=1e-10)
  res
end
res2 = ForwardDiff.gradient(G,[1.5,1.0,3.0])
res3 = Calculus.gradient(G,[1.5,1.0,3.0])

@test norm(res' .- res2) < 1e-8
@test norm(res' .- res3) < 1e-6

# Buffer length test
f = (du, u, p, t) -> du .= 0
p = zeros(3); u = zeros(50)
prob = ODEProblem(f,u,(0.0,10.0),p)
sol = solve(prob,Vern9(),abstol=1e-14,reltol=1e-14)
@test_nowarn res = adjoint_sensitivities(sol,Vern9(),dg,t,abstol=1e-14,
                                 reltol=1e-14,iabstol=1e-14,ireltol=1e-12)

using DiffEqSensitivity: adjoint_sensitivities_u0

function dg(out,u,p,t,i)
  out .= 1 .- u
end

ū0 = adjoint_sensitivities_u0(sol,Vern9(),dg,t,abstol=1e-14,
                         reltol=1e-14,iabstol=1e-14,ireltol=1e-12)[1]

ū0 ≈ ForwardDiff.gradient(prob.u0) do u0
  tmp_prob = remake(prob,u0=u0)
  sol = solve(tmp_prob,Vern9(),abstol=1e-14,reltol=1e-14,saveat=t)
  A = convert(Array,sol)
  sum(((1 .- A).^2)./2)
end

ext: jl


text:
module Util

s3sync(from, to) = run(`aws s3 sync $(from) $(to)`)
s3sync(from, to, tags) = run(`aws s3 sync $(from) $(to) $(tags)`)

function redirect_all(f::Function, fout::String, ferr::String)
  open(fout, "w") do ioout
    open(ferr, "w") do ioerr
      redirect_stdout(ioout) do
        redirect_stderr(ioerr) do
          f()
        end
      end
    end
  end
end

"""
example usage:

```julia
redirect_all("path/to/out") do
  dosomething()
end
```
"""
function redirect_all(f::Function, path::String)
  redirect_all("$(path)-out.txt", "$(path)-err.txt") do
    f()
  end
end

end


text:
""" Categories of C-sets and attributed C-sets.
"""
module CSets
export ACSetTransformation, CSetTransformation, components, force, is_natural,
  homomorphism, homomorphisms, isomorphism, isomorphisms, migrate!,
  generate_json_acset, parse_json_acset, read_json_acset, write_json_acset

using AutoHashEquals
using JSON
using Reexport
using StaticArrays: SVector

@reexport using ...CSetDataStructures
using ...GAT, ..FreeDiagrams, ..Limits, ..Sets, ..FinSets
import ..Limits: limit, colimit, universal
import ..FinSets: FinSet, FinFunction, FinDomFunction, force
using ...Theories: Category, CatDesc, AttrDesc, ob, hom, attr, adom, acodom
import ...Theories: dom, codom, compose, ⋅, id

# FinSets interop
#################

""" Create `FinSet` for part of attributed C-set.
"""
FinSet(X::ACSet, type::Symbol) = FinSet{Int,Int}(nparts(X, type))

""" Create `FinFunction` for part or subpart of attributed C-set.

The subpart must be of kind `Ob`. For indexed subparts, the index is included.
"""
FinFunction(X::ACSet, name::Symbol) = fin_function(X, Val{name})

@generated function fin_function(X::ACSet{CD,AD,Ts,Idxed},
    ::Type{Val{name}}) where {CD,AD,Ts,Idxed,name}
  if name ∈ ob(CD)
    quote
      FinFunction(identity, FinSet(X, $(QuoteNode(name))))
    end
  elseif name ∈ hom(CD)
    quote
      FinFunction(subpart(X, $(QuoteNode(name))),
                  FinSet(X, $(QuoteNode(codom(CD, name)))),
                  index=$(name ∈ Idxed ? :(X.indices.$name) : false))
    end
  else
    throw(ArgumentError("$(repr(name)) not in $(ob(CD)) or $(hom(CD))"))
  end
end

""" Create `FinDomFunction` for part or subpart of attributed C-set.

The codomain is always of type `TypeSet`, regardless of whether the subpart is
of kind `Ob` or `Data`. For indexed subparts, the index is included.
"""
FinDomFunction(X::ACSet, name::Symbol) = fin_dom_function(X, Val{name})

@generated function fin_dom_function(X::ACSet{CD,AD,Ts,Idxed},
    ::Type{Val{name}}) where {CD,AD,Ts,Idxed,name}
  if name ∈ ob(CD)
    quote
      n = nparts(X, $(QuoteNode(name)))
      FinDomFunction(1:n, FinSet(n), TypeSet{Int}())
    end
  elseif name ∈ hom(CD) || name ∈ attr(AD)
    quote
      FinDomFunction(subpart(X, $(QuoteNode(name))),
                     index=$(name ∈ Idxed ? :(X.indices.$name) : false))
    end
  else
    throw(ArgumentError(
      "$(repr(name)) not in $(ob(CD)), $(hom(CD)), or $(attr(AD))"))
  end
end

# C-set transformations
#######################

""" Transformation between attributed C-sets.

A homomorphism of C-sets is a natural transformation: a transformation between
functors C → Set satisfying the naturality axiom for all morphisms in C. This
struct records the data of a transformation; it does not enforce naturality, but
see [`is_natural`](@ref).

A C-set transformation has a component for every object in C. When C-sets have
attributes, the data types are assumed to be fixed. Thus, the naturality axiom
for data attributes is a commutative triangle, rather than a commutative square.
"""
@auto_hash_equals struct ACSetTransformation{CD <: CatDesc, AD <: AttrDesc{CD},
    Comp <: NamedTuple, Dom <: AbstractACSet{CD,AD}, Codom <: AbstractACSet{CD,AD}}
  components::Comp
  dom::Dom
  codom::Codom

  function ACSetTransformation{CD,AD}(components::NamedTuple, X::Dom, Y::Codom) where
      {Ob, CD <: CatDesc{Ob}, AD <: AttrDesc{CD},
       Dom <: AbstractACSet{CD,AD}, Codom <: AbstractACSet{CD,AD}}
    @assert keys(components) ⊆ Ob
    coerced_components = NamedTuple{Ob}(
      coerce_component(ob, get(components, ob) do; Int[] end, X, Y)
      for ob in Ob)
    new{CD,AD,typeof(coerced_components),Dom,Codom}(coerced_components, X, Y)
  end
end

function coerce_component(ob::Symbol, f::FinFunction{Int,Int}, X, Y)
  @assert length(dom(f)) == nparts(X,ob) "Domain error in component $ob"
  @assert length(codom(f)) == nparts(Y,ob) "Codomain error in component $ob"
  return f
end
function coerce_component(ob::Symbol, f, X, Y)::FinFunction{Int,Int}
  FinFunction(f, nparts(X,ob), nparts(Y,ob))
end

ACSetTransformation(components, X::Dom, Y::Codom) where
    {CD, AD, Dom <: AbstractACSet{CD,AD}, Codom <: AbstractACSet{CD,AD}} =
  ACSetTransformation{CD,AD}(components, X, Y)
ACSetTransformation(X::Dom, Y::Codom; components...) where
    {CD, AD, Dom <: AbstractACSet{CD,AD}, Codom <: AbstractACSet{CD,AD}} =
  ACSetTransformation{CD,AD}((; components...), X, Y)

const CSetTransformation{CD, Comp,
                         Dom <: AbstractCSet{CD}, Codom <: AbstractCSet{CD}} =
  ACSetTransformation{CD,AttrDesc{CD,(),(),(),()},Comp,Dom,Codom}

CSetTransformation(components, X::AbstractCSet, Y::AbstractCSet) =
  ACSetTransformation(components, X, Y)
CSetTransformation(X::AbstractCSet, Y::AbstractCSet; components...) =
  ACSetTransformation(X, Y; components...)

components(α::ACSetTransformation) = α.components
Base.getindex(α::ACSetTransformation, ob) = α.components[ob]

""" Is the transformation between C-sets a natural transformation?

Uses the fact that to check whether a transformation is natural, it suffices to
check the naturality equation on a generating set of morphisms.
"""
function is_natural(α::ACSetTransformation{CD,AD}) where {CD,AD}
  X, Y = dom(α), codom(α)
  for (f, c, d) in zip(hom(CD), dom(CD), codom(CD))
    Xf, Yf, α_c, α_d = subpart(X,f), subpart(Y,f), α[c], α[d]
    all(Yf[α_c(i)] == α_d(Xf[i]) for i in eachindex(Xf)) || return false
  end
  for (f, c) in zip(attr(AD), adom(AD))
    Xf, Yf, α_c = subpart(X,f), subpart(Y,f), α[c]
    all(Yf[α_c(i)] == Xf[i] for i in eachindex(Xf)) || return false
  end
  return true
end

force(α::ACSetTransformation) =
  ACSetTransformation(map(force, components(α)), dom(α), codom(α))

# Finding C-set transformations
###############################

""" Find a homomorphism between two attributed ``C``-sets.

Returns `nothing` if no homomorphism exists. For many categories ``C``, the
``C``-set homomorphism problem is NP-complete and thus this procedure generally
runs in exponential time. It works best when the domain object is small.

This procedure uses the classic backtracking search algorithm for a
combinatorial constraint satisfaction problem (CSP). As is well known, the
homomorphism problem for relational databases is equivalent to CSP. Since the
C-set homomorphism problem is "the same" as the database homomorphism problem
(insofar as attributed C-sets are "the same" as relational databases), it is
also equivalent to CSP. Backtracking search for CSP is described in many
computer science textbooks, such as (Russell & Norvig 2010, *Artificial
Intelligence*, Third Ed., Chapter 6: Constraint satisfaction problems, esp.
Algorithm 6.5). In our implementation, the search tree is ordered using the
popular heuristic of "minimum remaining values" (MRV), also known as "most
constrained variable."

Set the keyword argument `monic=true` to restrict to monomorphisms:
homomorphisms whose components are all injective functions.

See also: [`homomorphisms`](@ref), [`isomorphism`](@ref).
"""
function homomorphism(X::AbstractACSet, Y::AbstractACSet; monic::Bool=false)
  result = backtracking_search(X, Y, findall=false, monic=monic)
  isnothing(result) ? nothing : ACSetTransformation(result, X, Y)
end

""" Find all homomorphisms between two attributed ``C``-sets.

This function is at least as expensive as [`homomorphism`](@ref) and when no
homomorphisms exist, it is exactly as expensive.
"""
function homomorphisms(X::AbstractACSet, Y::AbstractACSet; monic::Bool=false)
  results = backtracking_search(X, Y, findall=true, monic=monic)
  map(components -> ACSetTransformation(components, X, Y), results)
end

""" Find an isomorphism between two attributed ``C``-sets, if one exists.

See [`homomorphism`](@ref) for more information about the algorithms involved.
"""
function isomorphism(X::AbstractACSet, Y::AbstractACSet)
  result = backtracking_search(X, Y, findall=false, iso=true)
  isnothing(result) ? nothing : ACSetTransformation(result, X, Y)
end

""" Find all isomorphisms between two attributed ``C``-sets.

This function is at least as expensive as [`isomorphism`](@ref) and when no
homomorphisms exist, it is exactly as expensive.
"""
function isomorphisms(X::AbstractACSet, Y::AbstractACSet)
  results = backtracking_search(X, Y, findall=true, iso=true)
  map(components -> ACSetTransformation(components, X, Y), results)
end

""" Internal state for backtracking search for ACSet homomorphisms.
"""
struct BacktrackingState{CD <: CatDesc, AD <: AttrDesc{CD},
    Assign <: NamedTuple, Dom <: AbstractACSet{CD,AD}, Codom <: AbstractACSet{CD,AD}}
  """ The current assignment, a partially-defined homomorphism of ACSets. """
  assignment::Assign
  """ Depth in search tree at which assignments were made. """
  assignment_depth::Assign
  """ Inverse assignment if finding a monomorphism, otherwise `nothing`. """
  inv_assignment::Union{Nothing,Assign}
  """ Domain ACSet: the "variables" in the CSP. """
  dom::Dom
  """ Codomain ACSet: the "values" in the CSP. """
  codom::Codom
  """ List of completed assignments if finding all homomorphisms;
  always `nothing` if finding one homomorphism.
  """
  results::Union{Nothing,Vector{Assign}}
end

function backtracking_search(X::AbstractACSet{CD}, Y::AbstractACSet{CD};
                             monic::Bool=false, iso::Bool=false,
                             findall::Bool=false) where {Ob, CD<:CatDesc{Ob}}
  # Fail early if no monic/iso exists on cardinality grounds.
  Assign = NamedTuple{Ob,Tuple{(Vector{Int} for c in Ob)...}}
  results = findall ? Assign[] : nothing
  if iso
    all(nparts(X,c) == nparts(Y,c) for c in Ob) || return results
    # Injections between finite sets are bijections, so reduce to that case.
    monic = true
  elseif monic
    all(nparts(X,c) <= nparts(Y,c) for c in Ob) || return results
  end

  assignment = NamedTuple{Ob}(zeros(Int, nparts(X, c)) for c in Ob)
  assignment_depth = map(copy, assignment)
  inv_assignment = monic ?
    NamedTuple{Ob}(zeros(Int, nparts(Y, c)) for c in Ob) : nothing
  state = BacktrackingState(assignment, assignment_depth, inv_assignment,
                            X, Y, results)
  backtracking_search(state, 0)
end

function backtracking_search(state::BacktrackingState, depth::Int)
  # Choose the next unassigned element.
  mrv, mrv_elem = find_mrv_elem(state, depth)
  if isnothing(mrv_elem)
    # No unassigned elements remain, so we have a complete assignment.
    return isnothing(state.results) ? state.assignment :
      push!(state.results, map(copy, state.assignment))
  elseif mrv == 0
    # An element has no allowable assignment, so we must backtrack.
    return state.results
  end
  c, x = mrv_elem

  # Attempt all assignments of the chosen element.
  Y = state.codom
  for y in parts(Y, c)
    if assign_elem!(state, depth, c, x, y)
      result = backtracking_search(state, depth + 1)
      if isnothing(state.results) && !isnothing(result)
        return result
      end
    end
    unassign_elem!(state, depth, c, x)
  end
  state.results
end

""" Find an unassigned element having the minimum remaining values (MRV).
"""
function find_mrv_elem(state::BacktrackingState{CD}, depth) where CD
  mrv, mrv_elem = Inf, nothing
  Y = state.codom
  for c in ob(CD), (x, y) in enumerate(state.assignment[c])
    y == 0 || continue
    n = count(can_assign_elem(state, depth, c, x, y) for y in parts(Y, c))
    if n < mrv
      mrv, mrv_elem = n, (c, x)
    end
  end
  (mrv, mrv_elem)
end

""" Check whether element (c,x) can be assigned to (c,y) in current assignment.
"""
function can_assign_elem(state::BacktrackingState, depth, c, x, y)
  # Although this method is nonmutating overall, we must temporarily mutate the
  # backtracking state, for several reasons. First, an assignment can be a
  # consistent at each individual subpart but not consistent for all subparts
  # simultaneously (consider trying to assign a self-loop to an edge with
  # distinct vertices). Moreover, in schemas with non-trivial endomorphisms, we
  # must keep track of which elements we have visited to avoid looping forever.
  ok = assign_elem!(state, depth, c, x, y)
  unassign_elem!(state, depth, c, x)
  return ok
end

""" Attempt to assign element (c,x) to (c,y) in the current assignment.

Returns whether the assignment succeeded. Note that the backtracking state can
be mutated even when the assignment fails.
"""
function assign_elem!(state::BacktrackingState{CD,AD}, depth, c, x, y) where {CD, AD}
  y′ = state.assignment[c][x]
  y′ == y && return true  # If x is already assigned to y, return immediately.
  y′ == 0 || return false # Otherwise, x must be unassigned.
  if !isnothing(state.inv_assignment) && state.inv_assignment[c][y] != 0
    # Also, y must unassigned in the inverse assignment.
    return false
  end

  # Check attributes first to fail as quickly as possible.
  X, Y = state.dom, state.codom
  for f in out_attr(AD, Val{c})
    subpart(X,x,f) == subpart(Y,y,f) || return false
  end

  # Make the assignment and recursively assign subparts.
  state.assignment[c][x] = y
  state.assignment_depth[c][x] = depth
  if !isnothing(state.inv_assignment)
    state.inv_assignment[c][y] = x
  end
  for (f, d) in out_hom(CD, Val{c})
    assign_elem!(state, depth, d,
                 subpart(X,x,f), subpart(Y,y,f)) || return false
  end
  return true
end

""" Unassign the element (c,x) in the current assignment.
"""
function unassign_elem!(state::BacktrackingState{CD}, depth, c, x) where CD
  state.assignment[c][x] == 0 && return
  assign_depth = state.assignment_depth[c][x]
  @assert assign_depth <= depth
  if assign_depth == depth
    X = state.dom
    if !isnothing(state.inv_assignment)
      y = state.assignment[c][x]
      state.inv_assignment[c][y] = 0
    end
    state.assignment[c][x] = 0
    state.assignment_depth[c][x] = 0
    for (f, d) in out_hom(CD, Val{c})
      unassign_elem!(state, depth, d, subpart(X,x,f))
    end
  end
end

@generated function out_hom(::Type{CD}, ::Type{Val{c}}) where {CD<:CatDesc, c}
  Expr(:tuple, (:($(QuoteNode(f)) => $(QuoteNode(codom(CD, f))))
                for f in hom(CD) if dom(CD, f) == c)...)
end
@generated function out_attr(::Type{AD}, ::Type{Val{c}}) where {AD<:AttrDesc, c}
  Expr(:tuple, (QuoteNode(f) for f in attr(AD) if dom(AD, f) == c)...)
end

# Category of C-sets
####################

@instance Category{ACSet, ACSetTransformation} begin
  dom(α::ACSetTransformation) = α.dom
  codom(α::ACSetTransformation) = α.codom

  id(X::ACSet) = ACSetTransformation(map(id, fin_sets(X)), X, X)

  function compose(α::ACSetTransformation, β::ACSetTransformation)
    # Question: Should we incur cost of checking that codom(β) == dom(α)?
    ACSetTransformation(map(compose, components(α), components(β)),
                        dom(α), codom(β))
  end
end

fin_sets(X::ACSet) = map(table -> FinSet(length(table)), tables(X))

# Limits and colimits
#####################

""" Limit of C-sets that stores the pointwise limits in FinSet.
"""
struct CSetLimit{Ob <: AbstractCSet, Diagram, Cone <: Multispan{Ob},
                 Limits <: NamedTuple} <: AbstractLimit{Ob,Diagram}
  diagram::Diagram
  cone::Cone
  limits::Limits
end

""" Colimit of attributed C-sets that stores the pointwise colimits in FinSet.
"""
struct ACSetColimit{Ob <: AbstractACSet, Diagram, Cocone <: Multicospan{Ob},
                    Colimits <: NamedTuple} <: AbstractColimit{Ob,Diagram}
  diagram::Diagram
  cocone::Cocone
  colimits::Colimits
end

# Compute limits and colimits of C-sets by reducing to those in FinSet using the
# "pointwise" formula for (co)limits in functor categories.

function limit(diagram::AbstractFreeDiagram{ACS}) where
    {CD <: CatDesc, ACS <: AbstractCSet{CD}}
  limits = map(limit, unpack_diagram(diagram))
  Xs = cone_objects(diagram)
  Y = ACS()
  for (c, lim) in pairs(limits)
    add_parts!(Y, c, length(ob(lim)))
  end
  for (f, c, d) in zip(hom(CD), dom(CD), codom(CD))
    Yfs = map(legs(limits[c]), Xs) do π, X
      compose(π, FinFunction(subpart(X, f), nparts(X, d)))
    end
    Yf = universal(limits[d], Multispan(ob(limits[c]), Yfs))
    set_subpart!(Y, f, collect(Yf))
  end
  πs = pack_components(map(legs, limits), map(X -> Y, Xs), Xs)
  CSetLimit(diagram, Multispan(Y, πs), limits)
end

function universal(lim::CSetLimit, cone::Multispan)
  components = map(universal, lim.limits, unpack_diagram(cone))
  CSetTransformation(components, apex(cone), ob(lim))
end

function colimit(diagram::AbstractFreeDiagram{ACS}) where
    {CD <: CatDesc, AD <: AttrDesc{CD}, Ts, ACS <: AbstractACSet{CD,AD,Ts}}
  # Colimit of C-set without attributes.
  colimits = map(colimit, unpack_diagram(diagram))
  Xs = cocone_objects(diagram)
  Y = ACS()
  for (c, colim) in pairs(colimits)
    add_parts!(Y, c, length(ob(colim)))
  end
  for (f, c, d) in zip(hom(CD), dom(CD), codom(CD))
    Yfs = map(legs(colimits[d]), Xs) do ι, X
      compose(FinFunction(subpart(X, f), nparts(X, d)), ι)
    end
    Yf = universal(colimits[c], Multicospan(ob(colimits[d]), Yfs))
    set_subpart!(Y, f, collect(Yf))
  end
  ιs = pack_components(map(legs, colimits), Xs, map(X -> Y, Xs))

  # Set data attributes by canonical inclusion from attributes in diagram.
  for (attr, c, d) in zip(attr(AD), adom(AD), acodom(AD))
    T = Ts.parameters[d]
    data = Vector{Union{Some{T},Nothing}}(nothing, nparts(Y, c))
    for (ι, X) in zip(ιs, Xs)
      for i in parts(X, c)
        j = ι[c](i)
        if isnothing(data[j])
          data[j] = Some(subpart(X, i, attr))
        else
          val1, val2 = subpart(X, i, attr), something(data[j])
          val1 == val2 || error(
            "ACSet colimit does not exist: $attr attributes $val1 != $val2")
        end
      end
    end
    set_subpart!(Y, attr, map(something, data))
  end

  ACSetColimit(diagram, Multicospan(Y, ιs), colimits)
end

function universal(colim::ACSetColimit, cocone::Multicospan)
  components = map(universal, colim.colimits, unpack_diagram(cocone))
  ACSetTransformation(components, ob(colim), apex(cocone))
end

""" Diagram in C-Set → named tuple of diagrams in FinSet
"""
unpack_diagram(discrete::DiscreteDiagram{<:AbstractACSet}) =
  map(DiscreteDiagram, unpack_sets(ob(discrete)))
unpack_diagram(span::Multispan{<:AbstractACSet}) =
  map(Multispan, fin_sets(apex(span)), unpack_components(legs(span)))
unpack_diagram(cospan::Multicospan{<:AbstractACSet}) =
  map(Multicospan, fin_sets(apex(cospan)), unpack_components(legs(cospan)))
unpack_diagram(para::ParallelMorphisms{<:AbstractACSet}) =
  map(ParallelMorphisms, unpack_components(hom(para)))

function unpack_diagram(d::Union{FreeDiagram{ACS},BipartiteFreeDiagram{ACS}}) where
    {Ob, CD <: CatDesc{Ob}, ACS <: AbstractACSet{CD}}
  NamedTuple{Ob}([ map(d, Ob=X -> FinSet(X, ob), Hom=α -> α[ob]) for ob in Ob ])
end

""" Vector of C-sets → named tuple of vectors of FinSets
"""
unpack_sets(Xs::AbstractVector{<:AbstractACSet{CD}}) where
    {Ob, CD <: CatDesc{Ob}} =
  NamedTuple{Ob}([ map(X -> FinSet(X, ob), Xs) for ob in Ob ])

""" Vector of C-set transformations → named tuple of vectors of FinFunctions
"""
unpack_components(αs::AbstractVector{<:ACSetTransformation{CD}}) where
    {Ob, CD <: CatDesc{Ob}} =
  NamedTuple{Ob}([ map(α -> α[ob], αs) for ob in Ob ])

""" Named tuple of vectors of FinFunctions → vector of C-set transformations
"""
function pack_components(fs::NamedTuple{Ob}, doms, codoms) where Ob
  components = map((x...) -> NamedTuple{Ob}(x), fs...) # XXX: Is there a better way?
  map(ACSetTransformation, components, doms, codoms)
end

""" Objects in diagram that will have explicit legs in limit cone.

Encodes common conventions such as, when taking a pullback of a cospan, not
explicitly including a cone leg for the cospan apex since it can be computed
from the other legs.

FIXME: Should this function be part of the official limits interface?
"""
cone_objects(diagram) = ob(diagram)
cone_objects(diagram::BipartiteFreeDiagram) = ob₁(diagram)
cone_objects(cospan::Multicospan) = feet(cospan)
cone_objects(para::ParallelMorphisms) = SVector(dom(para))

""" Objects in diagram that will have explicit legs in colimit cocone.
"""
cocone_objects(diagram) = ob(diagram)
cocone_objects(diagram::BipartiteFreeDiagram) = ob₂(diagram)
cocone_objects(span::Multispan) = feet(span)
cocone_objects(para::ParallelMorphisms) = SVector(codom(para))

# Functorial data migration
###########################

""" Pullback functorial data migration from one ACSet to another.

Note that this operation is contravariant: the data is transferred from `X` to
`Y` but the functor, represented by two dictionaries, maps the schema for `Y`
to the schema for `X`.

When the functor is the identity, this function is equivalent to
[`copy_parts!`](@ref).
"""
function migrate!(Y::ACSet{CD, AD}, X::ACSet,
                  FOb::AbstractDict, FHom::AbstractDict) where {CD, AD}
  ob(CD) ⊆ keys(FOb)     || error("Every object in $CD must be a key in $FOb")
  hom(CD) ⊆ keys(FHom)   || error("Every morphism in $CD must be a key in $FHom")
  attr(AD) ⊆ keys(FHom)  || error("Every attribute in $AD must be a key in $FHom")

  partsY = NamedTuple{ob(CD)}(map(ob(CD)) do obY
    add_parts!(Y, obY, nparts(X, FOb[obY]))
  end)
  for homY in hom(CD)
    domY, codomY = dom(CD, homY), codom(CD, homY)
    set_subpart!(Y, partsY[domY], homY, partsY[codomY][subpart(X, FHom[homY])])
  end
  for attrY in attr(AD)
    domY = dom(AD, attrY)
    set_subpart!(Y, partsY[domY], attrY, subpart(X, FHom[attrY]))
  end
  return Y
end

function (::Type{T})(X::ACSet, FOb::AbstractDict,
                     FHom::AbstractDict) where T <: AbstractACSet
  Y = T()
  migrate!(Y, X, FOb, FHom)
end

""" Serialize an ACSet object to a JSON string
"""
function generate_json_acset(x::T) where T <: AbstractACSet
  JSON.json(x.tables)
end

""" Deserialize a dictionary from a parsed JSON string to an object of the given ACSet type
"""
function parse_json_acset(type::Type{T}, input::Dict) where T <: AbstractACSet
  out = type()
  for (k,v) ∈ input
    add_parts!(out, Symbol(k), length(v))
  end
  for l ∈ values(input)
    for (i, j) ∈ enumerate(l)
      for (k,v) ∈ j
        vtype = eltype(out[Symbol(k)])
        if !(typeof(v) <: vtype)
          v = vtype(v)
        end
        set_subpart!(out, i, Symbol(k), v)
      end
    end
  end
  out
end

""" Deserialize a JSON string to an object of the given ACSet type
"""
function parse_json_acset(type::Type{T}, input::String) where T <: AbstractACSet
  parse_json_acset(type, JSON.parse(input))
end

""" Read a JSON file to an object of the given ACSet type
"""
function read_json_acset(type::Type{T}, input::String) where T <: AbstractACSet
  parse_json_acset(type, open(f->read(f, String), input))
end

""" Serialize an ACSet object to a JSON file
"""
function write_json_acset(x::T, fname::AbstractString) where T <: AbstractACSet
  open(string(fname), "w") do f
    write(f, generate_json_acset(x))
  end
end

end


text:
module FoldableTest

using Foldable
using Base.Test

function main()
    # Array
    let
        xs = [1,2,4,8,16]
        ys = [1,2,3,4,5]
        @test freduce(+, 0, xs) == 31
        @test freduce((r,x,y) -> r+x*y, 0, xs, ys) == 129
        x2s = [1 2 4; 8 16 32]
        y2s = [1 2 3; 4 5 6]
        @test freduce(+, 0, x2s) == 63
        @test freduce((r,x,y) -> r+x*y, 0, x2s, y2s) == 321
    end

    # Set
    let
        xs = Set([1,2,4,8,16])
        @test freduce(+, 0, xs) == 31
    end

    # Tuple
    let
        xs = (1,2,4,8,16)
        ys = (1,2,3,4,5)
        @test freduce(+, 0, xs) == 31
        @test freduce((x,y,z) -> x+y*z, 0, xs, ys) == 129
    end
end

end


text:
using Starlight
using Starlight.Log
using Test

# the namespace in front of the method name is important, apparently
Starlight.update!(r::Root, Δ::AbstractFloat) = r.updated = true

mutable struct TestEntity <: Entity end

Starlight.update!(t::TestEntity, Δ::AbstractFloat) = t.updated = true

@testset "Starlight" begin
  # load test file
  a = App()

  # no systems should be running
  @test off(a)

  # kick off
  awake!(a)

  # all systems should be running
  @test on(a)

  # test logging
  Log.test()

  # close
  shutdown!(a)

  # systems should no longer be running
  @test off(a)

  # ok back on
  awake!(a)

  # visual testing
  # NOTE in SDL, (0,0) is top-left and y goes down
  shapes = [
    ColorRect(200, 200; color=colorant"blue", pos=XYZ(200, 200)),
    Sprite(artifact"test/test/sprites/charmap-cellphone_black.png"; color=RGBA(1.0, 0, 0, 0.5), pos=XYZ(100, 100)),
  ]

  # test manipulations on root
  root = getEntityById(0)

  root.updated = false

  @test !root.updated

  # takes just a little longer than a second
  # for the first update to propagate when
  # JULIA_DEBUG=Starlight, this fixes it
  sleep(1.5)

  @test root.updated

  # manipulations on test entity
  tst = TestEntity()
  instantiate!(tst, props=Dict(
    :updated=>false
  ))

  @test !tst.updated

  sleep(1)

  @test tst.updated

  destroy!(tst)

  @test length(ecs()) == 3

  destroy!(shapes...)

  @test length(ecs()) == 1

  shutdown!(a)

  # root should remain and be unmodified
  @test root.updated

  @test off(a)

end


text:
# ### Image viewer using MiniFB
# This example displays the PNG files in a directory, rotating through the available images on each mouse click.

# This code needs Images.jl and ImageTransformations.jl. Add those packages to your environment before proceeding
using Images
using ImageTransformations
using MiniFB

# Set the size of the window
global const WIDTH = 600
global const HEIGHT = 400

# Global state
global ni = 1
global num_images = 1

# A function that will be called when a mouse button is clicked. It simply rotates the current image by storing a global integer.
function onclick(window, button, mod, isPressed)::Cvoid
    global ni, num_images
    if Bool(isPressed)
        if ni < num_images
            ni = ni+1
        else
            ni = 1
        end
    end
    return nothing
end

# The callback needs to be passed to the C library. The @cfunction macro creates a
onclick_c = @cfunction(onclick, Cvoid, (Ptr{Cvoid}, mfb_mouse_button, mfb_key_mod, UInt8) )

# Populate a MiniFB buffer from the image data.
# First, resize the image to the window dimensions. Calculate it's position so that it is centered in the window.
# Finally, inside the loop, convert each pixel to the 32 bit MiniFB buffer format.
function populate_buffer!(buffer, img)
    global WIDTH, HEIGHT
    buffer[:] = zeros(UInt32, HEIGHT*WIDTH)
    h, w = size(img)
    ratio = min(1, HEIGHT/h, WIDTH/w)
    nimg = imresize(img, ratio=ratio)
    nh, nw = size(nimg)
    oh = (HEIGHT-nh) ÷ 2
    ow = (WIDTH-nw) ÷ 2
    for i in 1:nh
        for j in 1:nw
            buffer[(oh+i-1)*WIDTH + ow+j] = mfb_rgb(nimg[i,j])
        end
    end
end

# This is the main function. It takes a directory, calls `loadfiles` to get the image data.
# Creates the windows, and sets up the callback. Then, inside the `while` loop, renders the buffer into the window.
# The `mfb_update` not only renders the buffer, but also flushes the input event queue.
# The buffer is changed only when the state changes. The `mfb_wait_sync` method enforces the required frame rate.
function imageview(dir::String=".")
    global ni, num_images
    images = loadfiles(dir)
    num_images = length(images)
    buffer = zeros(UInt32, HEIGHT*WIDTH)
    mfb_set_target_fps(10)
    window = mfb_open_ex("Image Viewer", WIDTH, HEIGHT, MiniFB.WF_RESIZABLE);
    mfb_set_mouse_button_callback(window, onclick_c);
    old_ni=0
    while mfb_wait_sync(window)
        if ni != old_ni
            populate_buffer!(buffer, images[ni])
            old_ni = ni
        end
        state = mfb_update(window, buffer);
        if state != MiniFB.STATE_OK
            break;
        end
    end
    mfb_close(window)
end

# Load the files from disk. Select all png files, and then use ImageIO to load them into memory.
function loadfiles(dir::String=".")
    files = readdir(dir)
    filter!(x->occursin(r"\.png"i, x), files)
    [load(joinpath(dir,x)) for x in files]
end

# Finally, call the main method to display the application
imageview(joinpath(dirname(pathof(MiniFB)),"..", "example"))

# ![](../assets/imageviewer.png)


text:
abstract type MicroCoverageException <: Exception
end

struct AlwaysAssertionError <: MicroCoverageException
    msg::String
end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/EAGO.jl
# The main file for EAGO.
#############################################################################

__precompile__()

module EAGO

    import MathOptInterface

    using Reexport, Cassette, IntervalArithmetic, NumericIO, DocStringExtensions
    using FastRounding, SpecialFunctions

    import JuMP
    import JuMP._Derivatives: operators, NodeData, VARIABLE, VALUE, SUBEXPRESSION,
                              Linearity, classify_linearity, UserOperatorRegistry,
                              register_multivariate_operator!
    using JuMP._Derivatives: univariate_operators,
                             univariate_operator_to_id

    import Ipopt, Hypatia, GLPK, Pavito

    using DataStructures: OrderedDict, BinaryMinMaxHeap, popmin!, popmax!, top
    using SparseArrays: SparseMatrixCSC, spzeros, rowvals, nzrange, nonzeros, sparse, findnz
    using LinearAlgebra: eigmin, norm
    using Printf: @printf

    import IntervalArithmetic: mid

    @reexport using McCormick
    @reexport using SpecialFunctions

    using IntervalContractors

    const MOI = MathOptInterface
    const MOIU = MathOptInterface.Utilities
    const MOIB = MathOptInterface.Bridges

    const SV = MOI.SingleVariable
    const SAF = MOI.ScalarAffineFunction
    const SQF = MOI.ScalarQuadraticFunction
    const VECVAR = MOI.VectorOfVariables

    const SQT = MOI.ScalarQuadraticTerm
    const SAT = MOI.ScalarAffineTerm

    const LT = MOI.LessThan
    const GT = MOI.GreaterThan
    const ET = MOI.EqualTo
    const IT = MOI.Interval
    const ZO = MOI.ZeroOne
    const SOC_CONE = MOI.SecondOrderCone
    const POW_CONE = MOI.PowerCone
    const EXP_CONE = MOI.ExponentialCone
    const PSD_CONE = MOI.PositiveSemidefiniteConeTriangle

    const VI = MOI.VariableIndex
    const CI = MOI.ConstraintIndex
    const ListCI = MOI.ListOfConstraintIndices

    const SCoefC = MOI.ScalarCoefficientChange
    const SConsC = MOI.ScalarConstantChange
    const MOIU = MOI.Utilities

    # Add storage types for EAGO optimizers
    export NodeBB, get_history, get_lower_bound, get_upper_bound, get_lower_time,
           get_upper_time, get_preprocess_time, get_postprocess_time, get_lower_bound, get_solution_time,
           get_iteration_number, get_node_count, get_absolute_gap, get_relative_gap

    export register_eago_operators!

    # map/reduce non-allocating no bounds checking map-reduce like utilities
    include("eago_optimizer/utilities/unsafe_utilities.jl")
    include("eago_optimizer/utilities/register_special.jl")

    # defines structure used to store information at each iteration of global optimize
    include("eago_optimizer/logging/log.jl")

    # defines the optimizer structures
    include("eago_optimizer/types/types.jl")

    # defines routines to add saf, sqf, and nlp block constraints
    include("eago_optimizer/math_opt_interface.jl")

    # functions which print information to console
    include("eago_optimizer/display.jl")

    #
    include("eago_optimizer/evaluate/evaluate.jl")

    #
    include("eago_optimizer/domain_reduction.jl")

    #
    include("eago_optimizer/parse_classify/parse_classify.jl")

    #
    include("eago_optimizer/logging/log_iteration.jl")

    #
    include("eago_optimizer/optimize/optimize.jl")

    # import the script solving utilities
    include("eago_script/script.jl")

    # routines for solving SIPs
    export SIPResult, SIPProblem, SIPCallback, SIPSubResult,
           sip_solve, SIPRes, SIPResRev, SIPHybrid,
           build_model, set_tolerance_inner!, set_tolerance!, get_disc_set,
           sip_llp!, sip_bnd!, sip_res!, get_sip_optimizer, check_convergence,
           LowerLevel1, LowerLevel2, LowerLevel3, LowerProblem, UpperProblem,
           ResProblem, AbstractSIPAlgo, AbstractSubproblemType
    include("eago_semiinfinite/semiinfinite.jl")
end


text:
base_path_default = "/gpfs/group/ebf11/default/acc"
base_path_default = "/home/eford/Data/HARPSN"
#manifest_subdir_default = "."
#data_subdir_default = "e2ds"


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

using Markdown
using InteractiveUtils

# ╔═╡ 9cd5e7cc-5c6c-11eb-3485-011143087d8d
begin
	import Pkg
	Pkg.activate(Base.current_project())
	using ReinforcementLearning
	using Flux
	using Statistics
	using Plots, Measures
	using SparseArrays
end

# ╔═╡ 41c57942-5c6c-11eb-3c99-9345b6668a1a
md"""
# Chapter 9 On-policy Prediction with Approximation

In this notebook, we'll focus on the linear approximation methods.
"""

# ╔═╡ be5a557c-5c6c-11eb-035e-39e94df41852
md"""
## Figure 9.1

We've discussed the `RandomWalk1D` environment before. In previous example, the state space is relatively small (`1:7`). Here we expand it into `1:1000` and see how the `LinearVApproximator` will work here.
"""

# ╔═╡ 6a0881f0-5c6d-11eb-143e-0196833abc05
ACTIONS = collect(Iterators.flatten((-100:-1, 1:100)))

# ╔═╡ 7ee0867c-5c6d-11eb-11b4-a7858177564f
NA = length(ACTIONS)

# ╔═╡ 7aae4986-5c6d-11eb-09b0-fd883165bc72
NS = 1002

# ╔═╡ 8fc27a60-5c6d-11eb-05ff-dbdcd106b853
md"""
First, let's roll out a large experiment to calculate the *true* state values of each state:
"""

# ╔═╡ c4d3a538-5c6d-11eb-3369-9bd67cc10bcd
TRUE_STATE_VALUES = begin
    env = RandomWalk1D(N=NS, actions=ACTIONS)
    agent = Agent(
        policy=VBasedPolicy(
            learner=TDLearner(
                approximator=TabularVApproximator(;n_state=NS,opt=Descent(0.01)),
                method=:SRS,
                ),
            mapping=(env,V) -> rand(action_space(env))
            ),
        trajectory=VectorSARTTrajectory()
    )
    run(agent, env, StopAfterEpisode(10^5))
    agent.policy.learner.approximator.table
end

# ╔═╡ df244c4e-5c6d-11eb-23e2-bf718f61180e
plot(TRUE_STATE_VALUES[2:end-1])

# ╔═╡ e0f292ce-5c6d-11eb-0fc8-7fd0cb86d6b5
md"""
Next, we define a preprocessor to map adjacent states into groups.
"""

# ╔═╡ 2e40e8dc-5c6e-11eb-264a-f99d9f5d62e0
N_GROUPS = 10

# ╔═╡ 248d425e-5c6e-11eb-2b34-fff856ca079c
begin
	Base.@kwdef struct GroupMapping
		n::Int
		n_groups::Int = N_GROUPS
		n_per_group::Int=div(n, N_GROUPS)
	end
	function (p::GroupMapping)(x::Int)
		if x == 1
			res = 1
		elseif x == p.n
			res = p.n_groups + 2
		else
			res = div(x - 2, p.n_per_group) + 2
		end
		res
	end
end

# ╔═╡ 484bb1a8-5c6e-11eb-3aaa-ebb84d7c8785
plot([GroupMapping(;n=NS)(i) for i in 1:NS], legend=nothing)

# ╔═╡ be5e114c-5c6e-11eb-007d-8b22a431ee4f
md"""
To count the frequency of each state, we need to write a hook.
"""

# ╔═╡ ed416b12-5c6e-11eb-1d57-b3d730907456
begin
	struct CountStates <: AbstractHook
		counts::Vector{Int}
		CountStates(n) = new(zeros(Int, n))
	end
	(f::CountStates)(::PreActStage, agent, env, action) = f.counts[state(env.env)] += 1
end

# ╔═╡ 091cad38-5c6f-11eb-3f4e-893c28972111
md"""
Now let's kickoff our experiment:
"""

# ╔═╡ 11cf3a2c-5c6f-11eb-231a-99660c29621e
agent_1 = Agent(
    policy=VBasedPolicy(
        learner=MonteCarloLearner(
            approximator=TabularVApproximator(n_state=N_GROUPS+2,opt=Descent(2e-5)),
            kind=EVERY_VISIT,  # this is very important!
            ),
        mapping=(env,V) -> rand(action_space(env))
        ),
    trajectory=VectorSARTTrajectory()
)

# ╔═╡ 198327e2-5c6f-11eb-26c4-19272b3374ae
env_1 = StateTransformedEnv(
    RandomWalk1D(N=NS, actions=ACTIONS),
    state_mapping=GroupMapping(n=NS)
)

# ╔═╡ 248e79a2-5c6f-11eb-30f5-9d417d286927
hook=CountStates(NS)

# ╔═╡ 2dcedb7e-5c6f-11eb-2026-abb9b378ee89
run(agent_1, env_1, StopAfterEpisode(10^5),hook)

# ╔═╡ 3fa64b5a-5c6f-11eb-2421-799d630e5460
begin
	fig_9_1 = plot(legend=:topleft, ylabel="Value scale", xlabel="State", right_margin = 1.5cm)
	fig_9_1_right = twinx(fig_9_1)
	plot!(fig_9_1, hook.counts./sum(hook.counts), color=:gray, label="state distribution")
	plot!(fig_9_1_right, agent_1.policy.learner.approximator.(env_1.state_mapping(s) for s in 2:NS-1), label="MC Learner", legend=:bottomright)
	plot!(fig_9_1_right, TRUE_STATE_VALUES[2:end-1], label="true values",legend=:bottomright, ylabel="Distribution scale")
end

# ╔═╡ 23060d86-5c70-11eb-2faa-a3851e3b5d2f
md"""
## Figure 9.2
"""

# ╔═╡ b19224a0-5c71-11eb-0582-2337e78a5ea9
agent_2 = Agent(
    policy=VBasedPolicy(
        learner=TDLearner(
            approximator=TabularVApproximator(n_state=N_GROUPS+2,opt=Descent(2e-4)),
            method=:SRS,
            ),
        mapping=(env,V) -> rand(action_space(env))
        ),
    trajectory=VectorSARTTrajectory()
)

# ╔═╡ e4ac979e-5c71-11eb-2ef3-8b5d88e38091
run(agent_2, env_1, StopAfterEpisode(10^5))

# ╔═╡ f12833c0-5c71-11eb-0c30-dd7f5c670a66
begin
	fig_9_2_left = plot(legend=:bottomright, xlabel="State")
	plot!(fig_9_2_left, agent_2.policy.learner.approximator.(env_1.state_mapping(s) for s in 2:NS-1), label="TD Learner", legend=:bottomright)
	plot!(fig_9_2_left, TRUE_STATE_VALUES[2:end-1], label="true values",legend=:bottomright)
	fig_9_2_left
end

# ╔═╡ c0afcbc6-5c72-11eb-3d31-69434487eda7
md"""
### Figure 9.2 right
"""

# ╔═╡ b797372c-5c72-11eb-3f71-5d57efd29d38
begin
	struct RecordRMS <: AbstractHook
		rms::Vector{Float64}
		RecordRMS() = new([])
	end
	function (f::RecordRMS)(::PostEpisodeStage, agent, env)
		push!(f.rms, sqrt(mean((agent.policy.learner.approximator.(env.state_mapping.(2:(NS-1))) - TRUE_STATE_VALUES[2:end-1]).^2)))
	end
end

# ╔═╡ 42af660e-5c73-11eb-3b38-bf3bb6bc6800
n_groups = 20

# ╔═╡ 804c0536-5c72-11eb-2294-93bb2512ed7a
function run_once(n, α)
    env = StateTransformedEnv(
		RandomWalk1D(N=NS, actions=ACTIONS),
		state_mapping=GroupMapping(n=NS)
	)
    agent = Agent(
        policy=VBasedPolicy(
            learner=TDLearner(
                approximator=TabularVApproximator(;
					n_state=n_groups+2,
					opt=Descent(α)
					),
                method=:SRS,
                n=n
                ),
            mapping=(env,V) -> rand(action_space(env))
            ),
        trajectory=VectorSARTTrajectory()
    )

    hook = RecordRMS()
    run(agent, env, StopAfterEpisode(10),hook)
    mean(hook.rms)
end

# ╔═╡ 54a7601e-5c73-11eb-2372-513196c58e8f
begin
	A = [0., 0.03, 0.06, 0.1:0.1:1...]
	fig_9_2_right = plot(legend=:bottomright, ylim=[0.25,0.55])
	for n in [2^i for i in 0:9]
		plot!(
			fig_9_2_right,
			A,
			mean(
				[run_once(n, α) for α in A]
				for _ in 1:100
				),
			label="n = $n")
	end
	fig_9_2_right
end

# ╔═╡ be968ab2-5c74-11eb-321f-0186685a329a
md"""
## Figure 9.5
"""

# ╔═╡ dfcef90a-5c74-11eb-04e1-f3320eba9bf8
begin
	struct FourierPreprocessor
		order::Int
	end
	(fp::FourierPreprocessor)(s::Number) = [cos(i * π * s) for i = 0:fp.order]
end

# ╔═╡ 7901bf16-5c75-11eb-13e5-0dd84f259c0a
begin
	struct PolynomialPreprocessor
		order::Int
	end
	(pp::PolynomialPreprocessor)(s::Number) = [s^i for i = 0:pp.order]
end

# ╔═╡ 08d133a0-5c77-11eb-1fbb-ed6b8da42d9f
md"""
## Figure 9.10

Implementing the tile encoding in Julia is quite easy！😀
"""

# ╔═╡ 2ef2aa46-5c77-11eb-1eec-13ad13061214
begin
	struct Tiling{N,Tr<:AbstractRange}
		ranges::NTuple{N,Tr}
		inds::LinearIndices{N,NTuple{N,Base.OneTo{Int}}}
	end

	Tiling(ranges...) =Tiling(
		ranges,
		LinearIndices(Tuple(length(r) - 1 for r in ranges))
	)

	Base.length(t::Tiling) = reduce(*, (length(r) - 1 for r in t.ranges))
end

# ╔═╡ 87c528bc-5c75-11eb-2f2f-adf254afda01
function run_once_MC(preprocessor, order, α)
    env = StateTransformedEnv(
        RandomWalk1D(N=NS, actions=ACTIONS),
        state_mapping=preprocessor
    )
    agent = Agent(
        policy=VBasedPolicy(
            learner=MonteCarloLearner(
                approximator=RLZoo.LinearVApproximator(;n=order+1,opt=Descent(α)),
                kind=EVERY_VISIT,
                ),
            mapping=(env,V) -> rand(1:NA)
            ),
        trajectory=VectorSARTTrajectory(;state=Vector{Float64})
    )

    hook=RecordRMS()
    run(agent, env, StopAfterEpisode(5000;is_show_progress=false),hook)
    hook.rms
end

# ╔═╡ c52bcb44-5c74-11eb-0e2b-fbb72e8edad8
begin

	fig_9_5 = plot(legend=:topright)

	for order in [5, 10, 20]
	plot!(
		fig_9_5,
		mean(
			run_once_MC(
				x -> FourierPreprocessor(order)(x/NS),
				order,
				0.00005
			)
			for _ in 1:5
		),
		label="Fourier $order",
		linestyle=:dash
	)

	plot!(
		fig_9_5,
		mean(
			run_once_MC(
				x -> PolynomialPreprocessor(order)(x/NS),
				order,
				0.0001
			)
			for _ in 1:5
		),
		label="Polynomial $order",
		linestyle=:solid
	)
	end

	fig_9_5
end

# ╔═╡ 592ac4a0-5c78-11eb-3d28-f7b178f4b94f
encode(range::AbstractRange, x) = floor(Int, div(x - range[1], step(range)) + 1)

# ╔═╡ 5c0304ee-5c78-11eb-2394-8fc17938918c
encode(t::Tiling, xs) = t.inds[CartesianIndex(Tuple(map(encode,  t.ranges, xs)))]

# ╔═╡ 3c773ea6-5c78-11eb-1a09-0f1fc560386d
t = Tiling(range(1, step=200, length=7))

# ╔═╡ 925141b4-5c78-11eb-208c-13e289d11f66
tt = [Tiling((range(1-4*(i-1), step=200, length=7))) for i in 1:50]

# ╔═╡ 7ab5686e-5c78-11eb-1067-a3127da36994
function run_once_MC_tiling(preprocessor, α, n)
    env = StateTransformedEnv(
        RandomWalk1D(N=NS, actions=ACTIONS),
        state_mapping=preprocessor
    )
    agent = Agent(
        policy=VBasedPolicy(
            learner=MonteCarloLearner(
                approximator=RLZoo.LinearVApproximator(;n=n,opt=Descent(α)),
                kind=EVERY_VISIT,
                ),
            mapping=(env,V) -> rand(1:NA)
            ),
        trajectory=VectorSARTTrajectory(;state=Vector{Float64})
    )

    hook=RecordRMS()
    run(agent, env, StopAfterEpisode(10000;is_show_progress=true),hook)
    hook.rms
end

# ╔═╡ 8b90595a-5c78-11eb-356b-df8d22646ed1
begin
	fig_9_10 = plot(xlabel="Episodes", ylabel="RMS error")

	plot!(
		fig_9_10,
		run_once_MC_tiling(
			x -> sparse([encode(t, x) for t in tt], 1:50, ones(50), 7, 50) |> vec,
			1e-4/50,
			7*50
		),
		label="50 tilings"
	)

	plot!(
		fig_9_10,
		run_once_MC_tiling(
			x -> Flux.onehot(encode(t, x), 1:7),
			1e-4,
			7
		),
		label = "one tiling"
	)

	fig_9_10
end

# ╔═╡ 248e1648-5c7a-11eb-0a7f-2767d27c80b6
md"""
Feel free to make a PR if you can improve the speed of generating this figure. ❤
"""

# ╔═╡ Cell order:
# ╟─41c57942-5c6c-11eb-3c99-9345b6668a1a
# ╠═9cd5e7cc-5c6c-11eb-3485-011143087d8d
# ╟─be5a557c-5c6c-11eb-035e-39e94df41852
# ╠═6a0881f0-5c6d-11eb-143e-0196833abc05
# ╠═7ee0867c-5c6d-11eb-11b4-a7858177564f
# ╠═7aae4986-5c6d-11eb-09b0-fd883165bc72
# ╟─8fc27a60-5c6d-11eb-05ff-dbdcd106b853
# ╠═c4d3a538-5c6d-11eb-3369-9bd67cc10bcd
# ╠═df244c4e-5c6d-11eb-23e2-bf718f61180e
# ╟─e0f292ce-5c6d-11eb-0fc8-7fd0cb86d6b5
# ╠═2e40e8dc-5c6e-11eb-264a-f99d9f5d62e0
# ╠═248d425e-5c6e-11eb-2b34-fff856ca079c
# ╠═484bb1a8-5c6e-11eb-3aaa-ebb84d7c8785
# ╟─be5e114c-5c6e-11eb-007d-8b22a431ee4f
# ╠═ed416b12-5c6e-11eb-1d57-b3d730907456
# ╟─091cad38-5c6f-11eb-3f4e-893c28972111
# ╠═11cf3a2c-5c6f-11eb-231a-99660c29621e
# ╠═198327e2-5c6f-11eb-26c4-19272b3374ae
# ╠═248e79a2-5c6f-11eb-30f5-9d417d286927
# ╠═2dcedb7e-5c6f-11eb-2026-abb9b378ee89
# ╠═3fa64b5a-5c6f-11eb-2421-799d630e5460
# ╟─23060d86-5c70-11eb-2faa-a3851e3b5d2f
# ╠═b19224a0-5c71-11eb-0582-2337e78a5ea9
# ╠═e4ac979e-5c71-11eb-2ef3-8b5d88e38091
# ╠═f12833c0-5c71-11eb-0c30-dd7f5c670a66
# ╟─c0afcbc6-5c72-11eb-3d31-69434487eda7
# ╠═b797372c-5c72-11eb-3f71-5d57efd29d38
# ╠═42af660e-5c73-11eb-3b38-bf3bb6bc6800
# ╠═804c0536-5c72-11eb-2294-93bb2512ed7a
# ╠═54a7601e-5c73-11eb-2372-513196c58e8f
# ╟─be968ab2-5c74-11eb-321f-0186685a329a
# ╠═dfcef90a-5c74-11eb-04e1-f3320eba9bf8
# ╠═7901bf16-5c75-11eb-13e5-0dd84f259c0a
# ╠═87c528bc-5c75-11eb-2f2f-adf254afda01
# ╠═c52bcb44-5c74-11eb-0e2b-fbb72e8edad8
# ╟─08d133a0-5c77-11eb-1fbb-ed6b8da42d9f
# ╠═2ef2aa46-5c77-11eb-1eec-13ad13061214
# ╠═592ac4a0-5c78-11eb-3d28-f7b178f4b94f
# ╠═5c0304ee-5c78-11eb-2394-8fc17938918c
# ╠═3c773ea6-5c78-11eb-1a09-0f1fc560386d
# ╠═925141b4-5c78-11eb-208c-13e289d11f66
# ╠═7ab5686e-5c78-11eb-1067-a3127da36994
# ╠═8b90595a-5c78-11eb-356b-df8d22646ed1
# ╟─248e1648-5c7a-11eb-0a7f-2767d27c80b6


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""CronJobSpec describes how the job execution will look like and when it will actually run.

    IoK8sApiBatchV1beta1CronJobSpec(;
        concurrencyPolicy=nothing,
        failedJobsHistoryLimit=nothing,
        jobTemplate=nothing,
        schedule=nothing,
        startingDeadlineSeconds=nothing,
        successfulJobsHistoryLimit=nothing,
        suspend=nothing,
    )

    - concurrencyPolicy::String : Specifies how to treat concurrent executions of a Job. Valid values are: - \"Allow\" (default): allows CronJobs to run concurrently; - \"Forbid\": forbids concurrent runs, skipping next run if previous run hasn't finished yet; - \"Replace\": cancels currently running job and replaces it with a new one
    - failedJobsHistoryLimit::Int32 : The number of failed finished jobs to retain. This is a pointer to distinguish between explicit zero and not specified. Defaults to 1.
    - jobTemplate::IoK8sApiBatchV1beta1JobTemplateSpec : Specifies the job that will be created when executing a CronJob.
    - schedule::String : The schedule in Cron format, see https://en.wikipedia.org/wiki/Cron.
    - startingDeadlineSeconds::Int64 : Optional deadline in seconds for starting the job if it misses scheduled time for any reason.  Missed jobs executions will be counted as failed ones.
    - successfulJobsHistoryLimit::Int32 : The number of successful finished jobs to retain. This is a pointer to distinguish between explicit zero and not specified. Defaults to 3.
    - suspend::Bool : This flag tells the controller to suspend subsequent executions, it does not apply to already started executions.  Defaults to false.
"""
mutable struct IoK8sApiBatchV1beta1CronJobSpec <: SwaggerModel
    concurrencyPolicy::Any # spec type: Union{ Nothing, String } # spec name: concurrencyPolicy
    failedJobsHistoryLimit::Any # spec type: Union{ Nothing, Int32 } # spec name: failedJobsHistoryLimit
    jobTemplate::Any # spec type: Union{ Nothing, IoK8sApiBatchV1beta1JobTemplateSpec } # spec name: jobTemplate
    schedule::Any # spec type: Union{ Nothing, String } # spec name: schedule
    startingDeadlineSeconds::Any # spec type: Union{ Nothing, Int64 } # spec name: startingDeadlineSeconds
    successfulJobsHistoryLimit::Any # spec type: Union{ Nothing, Int32 } # spec name: successfulJobsHistoryLimit
    suspend::Any # spec type: Union{ Nothing, Bool } # spec name: suspend

    function IoK8sApiBatchV1beta1CronJobSpec(;concurrencyPolicy=nothing, failedJobsHistoryLimit=nothing, jobTemplate=nothing, schedule=nothing, startingDeadlineSeconds=nothing, successfulJobsHistoryLimit=nothing, suspend=nothing)
        o = new()
        validate_property(IoK8sApiBatchV1beta1CronJobSpec, Symbol("concurrencyPolicy"), concurrencyPolicy)
        setfield!(o, Symbol("concurrencyPolicy"), concurrencyPolicy)
        validate_property(IoK8sApiBatchV1beta1CronJobSpec, Symbol("failedJobsHistoryLimit"), failedJobsHistoryLimit)
        setfield!(o, Symbol("failedJobsHistoryLimit"), failedJobsHistoryLimit)
        validate_property(IoK8sApiBatchV1beta1CronJobSpec, Symbol("jobTemplate"), jobTemplate)
        setfield!(o, Symbol("jobTemplate"), jobTemplate)
        validate_property(IoK8sApiBatchV1beta1CronJobSpec, Symbol("schedule"), schedule)
        setfield!(o, Symbol("schedule"), schedule)
        validate_property(IoK8sApiBatchV1beta1CronJobSpec, Symbol("startingDeadlineSeconds"), startingDeadlineSeconds)
        setfield!(o, Symbol("startingDeadlineSeconds"), startingDeadlineSeconds)
        validate_property(IoK8sApiBatchV1beta1CronJobSpec, Symbol("successfulJobsHistoryLimit"), successfulJobsHistoryLimit)
        setfield!(o, Symbol("successfulJobsHistoryLimit"), successfulJobsHistoryLimit)
        validate_property(IoK8sApiBatchV1beta1CronJobSpec, Symbol("suspend"), suspend)
        setfield!(o, Symbol("suspend"), suspend)
        o
    end
end # type IoK8sApiBatchV1beta1CronJobSpec

const _property_map_IoK8sApiBatchV1beta1CronJobSpec = Dict{Symbol,Symbol}(Symbol("concurrencyPolicy")=>Symbol("concurrencyPolicy"), Symbol("failedJobsHistoryLimit")=>Symbol("failedJobsHistoryLimit"), Symbol("jobTemplate")=>Symbol("jobTemplate"), Symbol("schedule")=>Symbol("schedule"), Symbol("startingDeadlineSeconds")=>Symbol("startingDeadlineSeconds"), Symbol("successfulJobsHistoryLimit")=>Symbol("successfulJobsHistoryLimit"), Symbol("suspend")=>Symbol("suspend"))
const _property_types_IoK8sApiBatchV1beta1CronJobSpec = Dict{Symbol,String}(Symbol("concurrencyPolicy")=>"String", Symbol("failedJobsHistoryLimit")=>"Int32", Symbol("jobTemplate")=>"IoK8sApiBatchV1beta1JobTemplateSpec", Symbol("schedule")=>"String", Symbol("startingDeadlineSeconds")=>"Int64", Symbol("successfulJobsHistoryLimit")=>"Int32", Symbol("suspend")=>"Bool")
Base.propertynames(::Type{ IoK8sApiBatchV1beta1CronJobSpec }) = collect(keys(_property_map_IoK8sApiBatchV1beta1CronJobSpec))
Swagger.property_type(::Type{ IoK8sApiBatchV1beta1CronJobSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiBatchV1beta1CronJobSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiBatchV1beta1CronJobSpec }, property_name::Symbol) =  _property_map_IoK8sApiBatchV1beta1CronJobSpec[property_name]

function check_required(o::IoK8sApiBatchV1beta1CronJobSpec)
    (getproperty(o, Symbol("jobTemplate")) === nothing) && (return false)
    (getproperty(o, Symbol("schedule")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiBatchV1beta1CronJobSpec }, name::Symbol, val)
end


text:
importall Base
using Knet,AutoGrad
using Knet: gpuinfo, gpusync
gpu(true)

AutoGrad._dbg{T}(x::KnetArray{T})="K$(join([id2(x),size(x)...,8*sizeof(T)],'_'))"

include(Pkg.dir("AutoGrad/test/profile.jl"))
gpuinfo("before d0kn,w2kn")
d0kn = map(i->map(KnetArray,i),d0)
w2kn = map(KnetArray,w2)
gpuinfo("after d0kn,w2kn")
# w1kn = map(KnetArray,w1)
# timeall(w2kn,d0kn,10)
# loop(fun[9],w2kn,d0kn,10)
# (x1kn,y1kn)=first(d0kn)
# fun[10](w2kn,x1kn,y1kn)

function timeall_kn(w=w2kn,d=d0kn,t=10)
    for i=1:length(fun)
        print(i); printfun(fun[i])
        for j=1:4
            gc_enable(false)
            @time (loop_kn(fun[i],w,d,t); gpusync())
            # @time loop_kn(fun[i],w,d,t)
            gpuinfo("before gc")
            gc_enable(true)
            gc()
            sleep(2)
            gpuinfo("after gc ")
        end
    end
end

function loop_kn(f,w,d,t)
    for i in 1:t
        for (x,y) in d
            f(w,x,y)
        end
    end
end

# gc_enable(false); @time loop_kn(fun[9],w2kn,d0kn,10); gc_enable(true)
# gc_enable(false); @time loop_kn(fun[9],w2kn,d0kn,10); gc_enable(true)
# gc_enable(false); @time loop_kn(fun[9],w2kn,d0kn,10); gc_enable(true)
# gc_enable(false); @time loop_kn(fun[9],w2kn,d0kn,10); gc_enable(true)
# gc_enable(false); @time loop_kn(fun[9],w2kn,d0kn,10); gc_enable(true)
# gc_enable(false); @profile loop_kn(fun[9],w2kn,d0kn,10); gc_enable(true)
# Profile.print()
# Profile.print(format=:flat)

timeall_kn()


text:
module classification

using Markdown
using InteractiveUtils
using CSV
using DataFrames
export get_img_paths, load_classes, labelFromPattern, examplelabeller, load_train_imgs
# MAIN FUNCTIONS

"""
Get all images in folder
"""
function get_img_paths(data_dir,extra_format = "")
    img_path = []
    for i in readdir(data_dir)
        if split(i, ".")[end] == "jpg" || split(i, ".")[end] == "png" || split(i, ".")[end] == extra_format
            push!(img_path, joinpath(data_dir, i))
        end
    end
    img_path
end

"""
- Identify the number of classes in the data. Return an array with all image paths and class labels
- Does not actually load the images. Just returns paths
- Assuming heirarchy of classes to be:
    - Main Folder
        - class1
            -img1
            -img2...
        - class2
            -img1
            -img2...
"""
function from_folder(data_dir)
    classes = readdir(main_path)
    @info length(classes), classes
    total_files = sum([length(readdir(main_path*classes[x])) for x in 1:length(classes)] )
    @info "Total images: $total_files"
    labels, paths = Array{String}(undef, total_files), Array{String}(undef, total_files)
    current_index = 1
    for path in classes
        current_path = readdir(main_path*path)
        current_length = length(current_path)
        labels[current_index:current_length+current_index-1] .= path
        # @info paths[current_index:current_length]
        paths[current_index:current_length+current_index-1] = main_path*path.*current_path
        current_index += current_length
    end
    labels, paths
end

"""
- Read paths from a csv file along with labels given path column name and label column name
- Further can be passed into labellers and encoders
- Just a convinence function really
"""
function from_csv(fname,pathcol, namecol, header = true,delim = ",", ignore = [])
        df = CSV.File(fname, delim = delim, header = header) |> DataFrame
        return df.namecol , df.pathcol
end

"""
- Takes paths generated from load_classes as input
- Takes a function defining a custom label
- eg: labeller(x) = split(x, "/")[-1]
- Returns labels based on function specified
"""
function labelFromPattern(paths, labeller)
    labeller.(paths)
end

"""
- A function to show a custom labeller function as an example
- This splits the path based on "/" and returns the lower case of the last element
"""
examplelabeller(x) = lowercase(split(x,"/")[end-1])

"""
- Input path -> Read an image -> Convert to array -> Resize to given dimensions -> Convert to float
"""
loadim(path, dims = (128, 128)) = float.(permutedims(channelview(imresize(load(path), dims...)), (2, 3, 1)))

"""
- Return a label encoded vector from 1..number of unique vals
- Also returns a copy of the labels for later use, use as encoded, _ if you do not care about re encoding it again
"""
function encoder(arr)
        uniquevals = unique(arr)
        hot = Dict(uniquevals[x] =>x for x in 1:length(uniquevals))
        @info hot
        return replace(arr, hot...) , reverse.(hot)
end

"""
- Load all images
- Resize them to given size
- one hot encode labels
"""
function load_train_imgs(labels, train_img_paths, dims = (128,128))
    train_imgs = zeros()
    imgs_loaded = 0
    Threads.@spawn for i in train_img_paths
        push!(train_imgs, loadim(i, dims...))
        imgs_loaded += 1
        if imgs_loaded % 1000 == 0
            @info "$imgs_loaded Images have been loaded"
        end
    end
end

end

# EXAMPLE WORKFLOW

# load filepaths, images
#main_path = "/home/subhaditya/Desktop/Datasets/asl-downloaded/"

#labels, paths = from_folder(main_path)

# label it in whatever way you want
#@info labelFromPattern(paths, examplelabeller)[1:3]

#one hot encoded labels (for classification)
#labels,d = encoder(labels)


text:
"""
    function stoichiometry_matrix(m::SBML.Model)

Extract the vector of species (aka metabolite) identifiers, vector of reaction
identifiers, and a sparse stoichiometry matrix (of type `SparseMatrixCSC` from
`SparseArrays` package) from an existing `SBML.Model`. Returns a 3-tuple with
these values.
"""
function stoichiometry_matrix(m::SBML.Model)
    rows = collect(keys(m.species))
    cols = collect(keys(m.reactions))
    row_idx = Dict(k => i for (i, k) in enumerate(rows))
    col_idx = Dict(k => i for (i, k) in enumerate(cols))

    nnz = 0
    for (_, r) in m.reactions
        for _ in r.reactants
            nnz += 1
        end
        for _ in r.products
            nnz += 1
        end
    end

    SI = Int[]
    RI = Int[]
    SV = Float64[]
    sizehint!(SI, nnz)
    sizehint!(RI, nnz)
    sizehint!(SV, nnz)

    for (rid, r) in m.reactions
        ridx = col_idx[rid]
        for (sid, stoi) in r.reactants
            push!(SI, row_idx[sid])
            push!(RI, ridx)
            push!(SV, -stoi)
        end
        for (sid, stoi) in r.products
            push!(SI, row_idx[sid])
            push!(RI, ridx)
            push!(SV, stoi)
        end
    end
    return rows, cols, SparseArrays.sparse(SI, RI, SV, length(rows), length(cols))
end

"""
    flux_bounds(m::SBML.Model)::NTuple{2, Vector{Tuple{Float64,String}}}

Extract the vectors of lower and upper bounds of reaction rates from the model,
in the same order as `keys(m.reactions)`.  All bounds are accompanied with the
unit of the corresponding value (the behavior is based on SBML specification).
Missing bounds are represented by negative/positive infinite values with
empty-string unit.
"""
function flux_bounds(m::SBML.Model)::NTuple{2,Vector{Tuple{Float64,String}}}
    # Now this is tricky. There are multiple ways in SBML to specify a
    # lower/upper bound. There are the "global" model bounds that we completely
    # ignore now because no one uses them. In reaction, you can specify the
    # bounds using "LOWER_BOUND" and "UPPER_BOUND" parameters, but also there
    # may be a FBC plugged-in objective name that refers to the parameters.
    # We extract these, using the units from the parameters. For unbounded
    # reactions this gives -Inf or Inf as a default.

    function get_bound(rxn, fld, param, default)
        p = mayfirst(getfield(rxn, fld), param)
        get(rxn.kinetic_parameters, p, get(m.parameters, p, default))
    end

    (
        get_bound.(values(m.reactions), :lower_bound, "LOWER_BOUND", Ref((-Inf, ""))),
        get_bound.(values(m.reactions), :upper_bound, "UPPER_BOUND", Ref((Inf, ""))),
    )
end

"""
    flux_objective(m::SBML.Model)::Vector{Float64}

Extract the vector of objective coefficients of each reaction, in the same
order as `keys(m.reactions)`.
"""
function flux_objective(m::SBML.Model)::Vector{Float64}
    # As with bounds, this sometimes needs to be gathered from 2 places (maybe
    # even more). FBC-specified OC gets a priority.
    function get_oc(rid::String)
        mayfirst(
            get(m.objective, rid, nothing),
            maylift(
                first,
                get(m.reactions[rid].kinetic_parameters, "OBJECTIVE_COEFFICIENT", nothing),
            ),
            0.0,
        )
    end
    get_oc.(keys(m.reactions))
end

"""
    mayfirst(args::Maybe{T}...)::Maybe{T} where T

Helper to get the first non-`nothing` value from the arguments.
"""
function mayfirst(args...)
    for i in args
        if !isnothing(i)
            return i
        end
    end
    nothing
end

"""
    maylift(f, args::Maybe...)

Helper to lift a function to work on [`Maybe`](@ref), returning `nothing`
whenever there's a `nothing` in args.
"""
maylift(f, args::Maybe...) = any(isnothing, args) ? nothing : f(args...)

"""
    get_compartment_size(m::SBML.Model, compartment; default = nothing)

A helper for easily getting out a defaulted compartment size.
"""
get_compartment_size(m::SBML.Model, compartment; default = nothing) =
    mayfirst(maylift(x -> x.size, get(m.compartments, compartment, nothing)), default)

"""
    initial_amounts(
        m::SBML.Model;
        convert_concentrations = false,
        compartment_size = comp -> get_compartment_size(m, comp),
    )

Return initial amounts for each species as a generator of pairs
`species_name => initial_amount`; the amount is set to `nothing` if not
available. If `convert_concentrations` is true and there is information about
initial concentration available together with compartment size, the result is
computed from the species' initial concentration.

In the current version, units of the measurements are completely ignored.

# Example
```
# get the initial amounts as dictionary
Dict(SBML.initial_amounts(model, convert_concentrations = true))

# suppose the compartment size is 10.0 if unspecified
collect(SBML.initial_amounts(
    model,
    convert_concentrations = true,
    compartment_size = comp -> SBML.get_compartment_size(model, comp, 10.0),
))

# remove the empty entries
Dict(k => v for (k,v) in SBML.initial_amounts(model) if !isnothing(v))
```
"""
initial_amounts(
    m::SBML.Model;
    convert_concentrations = false,
    compartment_size = comp -> get_compartment_size(m, comp),
) = (
    k => mayfirst(
        maylift(first, s.initial_amount),
        if convert_concentrations
            maylift(
                (ic, s) -> first(ic) * s,
                s.initial_concentration,
                compartment_size(s.compartment),
            )
        end,
    ) for (k, s) in m.species
)

"""
    initial_concentrations(
        m::SBML.Model;
        convert_amounts = false,
        compartment_size = comp -> get_compartment_size(m, comp),
    )

Return initial concentrations of the species in the model. Refer to work-alike
[`initial_amounts`](@ref) for details.
"""
initial_concentrations(
    m::SBML.Model;
    convert_amounts = false,
    compartment_size = comp -> get_compartment_size(m, comp),
) = (
    k => mayfirst(
        maylift(first, s.initial_concentration),
        if convert_amounts
            maylift(
                (ia, s) -> first(ia) / s,
                s.initial_amount,
                compartment_size(s.compartment),
            )
        end,
    ) for (k, s) in m.species
)

"""
    function extensive_kinetic_math(
        m::SBML.Model,
        formula::SBML.Math;
        handle_empty_compartment_size = (id::String) -> throw(
            DomainError(
                "Non-substance-only-unit reference to species `\$id' in an unsized compartment `\$(m.species[id].compartment)'",
            ),
        ),
    )

Convert a SBML math `formula` to "extensive" kinetic laws, where the references
to species that are marked as not having only substance units are converted
from amounts to concentrations.

If the data is missing, you can supply a function that adds them. A common way
to handle errors is to assume that unsized compartments have volume 1.0 (of
whatever units), you can specify that behavior by supplying
`handle_empty_compartment_size = _ -> 1.0`.

Handling of units in the conversion process is ignored in this version.
"""
function extensive_kinetic_math(
    m::SBML.Model,
    formula::SBML.Math;
    handle_empty_compartment_size = (id::String) -> throw(
        DomainError(
            "Non-substance-only-unit reference to species `$id' in an unsized compartment `$(m.species[id].compartment)'",
        ),
    ),
)
    conv(x::SBML.MathIdent) = begin
        haskey(m.species, x.id) || return x
        sp = m.species[x.id]
        sp.only_substance_units && return x
        sz = m.compartments[sp.compartment].size
        isnothing(sz) && (sz = handle_empty_compartment_size(x.id))
        SBML.MathApply("/", [x, SBML.MathVal(sz)])
    end
    conv(x::SBML.MathApply) = SBML.MathApply(x.fn, conv.(x.args))
    conv(x::SBML.Math) = x

    conv(formula)
end

"""
    get_error_messages(doc::VPtr, error::Exception, report_severities)

Show the error messages reported by SBML in the `doc` document and throw the
`error` if they are more than 1.

`report_severities` switches the reporting of certain error types defined by
libsbml; you can choose from `["Fatal", "Error", "Warning", "Informational"]`.
"""
function get_error_messages(doc::VPtr, error::Exception, report_severities)
    n_errs = ccall(sbml(:SBMLDocument_getNumErrors), Cuint, (VPtr,), doc)
    do_throw = false
    for i = 1:n_errs
        err = ccall(sbml(:SBMLDocument_getError), VPtr, (VPtr, Cuint), doc, i - 1)
        msg = string(strip(get_string(err, :XMLError_getMessage)))
        sev = string(strip(get_string(err, :XMLError_getSeverityAsString)))
        # keywords from `libsbml/src/sbml/xml/XMLError.cpp` xmlSeverityStringTable:
        if sev == "Fatal"
            sev in report_severities && @error "SBML reported fatal error: $(msg)"
            do_throw = true
        elseif sev == "Error"
            sev in report_severities && @error "SBML reported error: $(msg)"
            do_throw = true
        elseif sev == "Warning"
            sev in report_severities && @warn "SBML reported warning: $(msg)"
        elseif sev == "Informational"
            sev in report_severities && @info "SBML reported: $(msg)"
        end
    end
    do_throw && throw(error)
    nothing
end

"""
    check_errors(
        success::Integer,
        doc::VPtr,
        error::Exception,
        report_severities = ["Fatal", "Error"],
    )

If success is a 0-valued `Integer` (a logical `false`), then call
[`get_error_messages`](@ref) to show the error messages reported by SBML in the
`doc` document and throw the `error` if they are more than 1.  `success` is
typically the value returned by an SBML C function operating on `doc` which
returns a boolean flag to signal a successful operation.
"""
check_errors(
    success::Integer,
    doc::VPtr,
    error::Exception,
    report_severities = ["Fatal", "Error"],
) = Bool(success) || get_error_messages(doc, error, report_severities)

function Base.show(io::IO, ::MIME"text/plain", m::SBML.Model)
    print(
        io,
        repr(typeof(m)),
        " with ",
        length(m.reactions),
        " reactions, ",
        length(m.species),
        " species, and ",
        length(m.parameters),
        " parameters.",
    )
end

# Define equality `==` methods for some of our custom types.
Base.:(==)(a::MathApply, b::MathApply) = a.fn == b.fn && a.args == b.args
Base.:(==)(a::AlgebraicRule, b::AlgebraicRule) = a.math == b.math
Base.:(==)(a::T, b::T) where {T<:Union{AssignmentRule,RateRule}} =
    a.id == b.id && a.math == b.math
Base.:(==)(a::Constraint, b::Constraint) = a.math == b.math && a.message == b.message


text:
using NFFT, DataFrames, LinearAlgebra, LaTeXStrings, DelimitedFiles
using BenchmarkTools
using Plots, StatsPlots, CategoricalArrays
pgfplotsx()
#gr()

include("../Wrappers/NFFT3.jl")
include("../Wrappers/FINUFFT.jl")

const threads = [1,2,4,8]
const precomp = [NFFT.LINEAR, NFFT.POLYNOMIAL, NFFT.TENSOR, NFFT.LINEAR, NFFT.LINEAR, NFFT.TENSOR]
const packagesCtor = [NFFTPlan, NFFTPlan, NFFTPlan, FINUFFTPlan, NFFT3Plan, NFFT3Plan]
const packagesStr = ["NFFT.jl/LINEAR", "NFFT.jl/POLY", "NFFT.jl/TENSOR", "FINUFFT", "NFFT3/LINEAR", "NFFT3/TENSOR"]
const benchmarkTime = [2, 5, 5]

NFFT.FFTW.set_num_threads(Threads.nthreads())
ccall(("omp_set_num_threads",NFFT3.lib_path_nfft),Nothing,(Int64,),convert(Int64,Threads.nthreads()))
@info ccall(("nfft_get_num_threads",NFFT3.lib_path_nfft),Int64,())
NFFT._use_threads[] = (Threads.nthreads() > 1)

function nfft_performance_comparison(m = 4, σ = 2.0)
  println("\n\n ##### nfft_performance ##### \n\n")

  df = DataFrame(Package=String[], Threads=Int[], D=Int[], M=Int[], N=Int[],
                   Undersampled=Bool[], m = Int[], σ=Float64[],
                   TimePre=Float64[], TimeTrafo=Float64[], TimeAdjoint=Float64[] )

  N = [collect(4096* (4 .^(0:3))),collect(128* (2 .^ (0:3))),[32,48,64,72]]
  fftflags = NFFT.FFTW.MEASURE

  for D = 2:2
    for U = 4:4
      NN = ntuple(d->N[D][U], D)
      M = prod(NN) #÷ 8

        @info D, NN, M

        T = Float64

        x = T.(rand(T,D,M) .- 0.5)
        x .= sortslices(x, dims=2) # sort nodes to gain cache locality

        fHat = randn(Complex{T}, M)
        f = randn(Complex{T}, NN)

        for pl = 1:length(packagesStr)

          planner = packagesCtor[pl]
          BenchmarkTools.DEFAULT_PARAMETERS.seconds = benchmarkTime[1]
          b = @benchmark $planner($x, $NN; m=$m, σ=$σ, window=:kaiser_bessel,
                                 precompute=$(precomp[pl]), sortNodes=false, fftflags=$fftflags)
          tpre = minimum(b).time / 1e9
          p = planner(x, NN; m=m, σ=σ, window=:kaiser_bessel,
                      precompute=(precomp[pl]), sortNodes=false, fftflags=fftflags)
          BenchmarkTools.DEFAULT_PARAMETERS.seconds = benchmarkTime[2]
          b = @benchmark mul!($f, $(adjoint(p)), $fHat)
          tadjoint = minimum(b).time / 1e9
          BenchmarkTools.DEFAULT_PARAMETERS.seconds = benchmarkTime[3]
          b = @benchmark mul!($fHat, $p, $f)
          ttrafo = minimum(b).time / 1e9

          push!(df, (packagesStr[pl], Threads.nthreads(), D, M, N[D][U], false, m, σ,
                   tpre, ttrafo, tadjoint))

      end
    end
  end
  return df
end

function plot_performance(df; D=2, N=1024, M=N*N)

  Plots.scalefontsizes()
  Plots.scalefontsizes(1.5)

  labelsA = packagesStr

  tpre = zeros(length(threads),length(labelsA))
  ttrafo = zeros(length(threads), length(labelsA))
  tadjoint = zeros(length(threads), length(labelsA))
  for (i,p) in enumerate(packagesStr)
    for (j,th) in enumerate(threads) #.&& df.Pre.==precomp[i]
      df_ = df[df.Threads .== th .&& df.Package.==p .&& df.D .== D .&& df.M .== M .&& df.N .==N ,:]
      tpre[j,i] = df_[1,:TimePre]
      ttrafo[j,i] = df_[1,:TimeTrafo]
      tadjoint[j,i] = df_[1,:TimeAdjoint]
    end
  end

  labelsB = "t = " .* string.(threads)

  ctg = CategoricalArray(repeat(labelsA, inner = length(threads)))
  levels!(ctg, labelsA)
  name = CategoricalArray(repeat(labelsB, outer = length(labelsA)))
  levels!(name, labelsB)

  maxtime = max(maximum(tpre), maximum(ttrafo), maximum(tadjoint))

  p1 = groupedbar(name, ttrafo, ylabel = "time / s",  group = ctg,
          bar_width = 0.67, legend = :none, ylims=(0,maxtime),
          lw = 0, framestyle = :box, size=(800,600), title = L"\textrm{NFFT}")

  p2 = groupedbar(name, tadjoint, ylabel = "time / s",  group = ctg,
          bar_width = 0.67, legend = :none, ylims=(0,maxtime),
          lw = 0, framestyle = :box, size=(800,600), title = L"\textrm{NFFT}^H")

  p3 = groupedbar(name, tpre, ylabel = "time / s",  group = ctg,
          bar_width = 0.67, ylims=(0,maxtime),
          lw = 0, framestyle = :box, size=(800,600), title = L"\textrm{Precompute}")

  p = plot(p1, p2, p3, layout=(3,1), size=(800,600), dpi=200)

  savefig(p, "../docs/src/assets/performance_mt_$(D)_$(N)_$(M).svg")
  return p
end

### run the code ###

if haskey(ENV, "NFFT_PERF")
  df = nfft_performance_comparison(4, 2.0)

  if isfile("performance_mt.csv")
    data, header = readdlm("performance_mt.csv", ',', header=true);
    df_ = DataFrame(data, vec(header))
    append!(df, df_)
  end

  writedlm("performance_mt.csv", Iterators.flatten(([names(df)], eachrow(df))), ',')

else
  rm("performance_mt.csv", force=true)
  ENV["NFFT_PERF"] = 1
  for t in threads
    cmd = `julia -t $t performance.jl`
    @info cmd
    run(cmd)

  end
  data, header = readdlm("performance_mt.csv", ',', header=true);
  df = DataFrame(data, vec(header))
  delete!(ENV, "NFFT_PERF")

  plot_performance(df, N=1024, M=1024*1024)
  #plot_performance(df, pre="LINEAR", N=1024, M=1024*1024*4)
  #plot_performance(df, pre="FULL")
  #plot_performance_serial(df)
end


text:
const TRY_BUT_ALLOW_FAILURES_URL_LIST = String[
    ]


text:
using DataFrames

const MODEL_MERGE = Dict{String, String}(
	"aae_vamp" 	=> "aae",
	"wae_vamp" 	=> "wae",
	"vae+ocsvm" => "vae_ocsvm"
)

const MODEL_TYPE = Dict{String, String}(
	"MAF" 				=> "flows",
	"RealNVP"			=> "flows",
	"sptn"				=> "flows",
	"MO_GAAL" 			=> "gans",
	"fmgan" 			=> "gans",
	"gan" 				=> "gans",
	"fAnoGAN" 			=> "gans",
	"fAnoGAN-GP" 		=> "gans",
	"abod" 				=> "classical",
	"hbos" 				=> "classical",
	"if" 				=> "classical",
	"knn" 				=> "classical",
	"loda" 				=> "classical",
	"lof" 				=> "classical",
	"ocsvm" 			=> "classical",
	"ocsvm_rbf"			=> "classical",
	"ocsvm_nu"			=> "classical",
	"pidforest" 		=> "classical",
	"GANomaly" 			=> "autoencoders",
	"aae" 				=> "autoencoders",
	"aae_vamp" 			=> "autoencoders",
	"aae_full" 			=> "autoencoders",
	"adVAE" 			=> "autoencoders",
	"vae" 				=> "autoencoders",
	"wae" 				=> "autoencoders",
	"wae_vamp" 			=> "autoencoders",
	"wae_full" 			=> "autoencoders",
	"vae_simple"		=> "autoencoders",
	"vae_full"			=> "autoencoders",
	"Conv-GANomaly" 	=> "autoencoders",
	"Conv-SkipGANomaly" => "autoencoders",
	"aae_ocsvm"			=> "two-stage",
	"vae_ocsvm" 		=> "two-stage",
	"vae_knn"	 		=> "two-stage",
	"DeepSVDD" 			=> "two-stage"
)

const MODEL_ALIAS = Dict{String, String}(
	"MAF" 				=> "maf",
	"RealNVP"			=> "rnvp",
	"sptn"				=> "sptn",
	"MO_GAAL" 			=> "mgal",
	"fmgan" 			=> "fmgn",
	"gan" 				=> "gan",
	"fAnoGAN" 			=> "fano",
	"fAnoGAN-GP" 		=> "fngp",
	"abod" 				=> "abod",
	"hbos" 				=> "hbos",
	"if" 				=> "if",
	"knn" 				=> "knn",
	"loda" 				=> "loda",
	"lof" 				=> "lof",
	"ocsvm" 			=> "osvm",
	"ocsvm_rbf"			=> "orbf",
	"ocsvm_nu"			=> "osnu",
	"pidforest" 		=> "pidf",
	"GANomaly" 			=> "gano",
	"aae" 				=> "aae",
	"aae_vamp" 			=> "aaev",
	"aae_full"			=> "aaef",
	"aae_ocsvm"			=> "aaeo",
	"adVAE" 			=> "avae",
	"vae" 				=> "vae",
	"wae" 				=> "wae",
	"wae_vamp" 			=> "waev",
	"wae_full"			=> "waef",
	"Conv-GANomaly" 	=> "gano",
	"Conv-SkipGANomaly" => "skip",
	"vae_simple" 		=> "vaes",
	"vae_full"			=> "vaef",
	"vae_ocsvm" 		=> "vaeo",
	"vae_knn"			=> "vaek",
	"DeepSVDD"			=> "dsvd"
)

const DATASET_ALIAS = Dict{String, String}(
	"abalone"					=> "aba",
	"annthyroid"				=> "ann",
	"arrhythmia"				=> "arr",
	"blood-transfusion"			=> "blt",
	"breast-cancer-wisconsin"	=> "bcw",
	"breast-tissue"				=> "bts",
	"cardiotocography"			=> "crd",
	"ecoli"						=> "eco",
	"glass"						=> "gls",
	"haberman"					=> "hab",
	"har"						=> "har",
	"htru2"						=> "htr",
	"ionosphere"				=> "ion",
	"iris"						=> "irs",
	"isolet"					=> "iso",
	"kdd99_small"				=> "kdd",
	"letter-recognition"		=> "ltr",
	"libras"					=> "lbr",
	"magic-telescope"			=> "mgc",
	"mammography"				=> "mam",
	"miniboone"					=> "mnb",
	"multiple-features"			=> "mlt",
	"page-blocks"				=> "pgb",
	"parkinsons"				=> "prk",
	"pendigits"					=> "pen",
	"pima-indians"				=> "pim",
	"seismic"					=> "sei",
	"sonar"						=> "snr",
	"spambase"					=> "spm",
	"spect-heart"				=> "sph",
	"statlog-satimage"			=> "sat",
	"statlog-segment"			=> "seg",
	"statlog-shuttle"			=> "sht",
	"statlog-vehicle"			=> "vhc",
	"synthetic-control-chart"	=> "scc",
	"wall-following-robot" 		=> "wrb",
	"waveform-1"				=> "wf1",
	"waveform-2"				=> "wf2",
	"wine"						=> "wne",
	"yeast"						=> "yst",
	"MNIST" 					=> "mnist",
	"FashionMNIST" 				=> "fmnist",
	"CIFAR10" 					=> "cifar10",
	"SVHN2" 					=> "svhn2"
)

const AC_CONVERSION = Dict(
	"mnist" 			=> ["5", "0", "4", "1", "9", "2", "3", "6", "7", "8"],
	"svhn2"				=> ["1", "9", "2", "3", "5", "8", "7", "4", "6", "0"],
	"fmnist" 			=> ["Ankle boot", "T-Shirt", "Dress", "Pullover", "Sneaker", "Sandal", "Trouser", "Shirt", "Coat", "Bag"],
	"cifar10"			=> ["frog", "truck", "deer", "automobile", "bird", "horse", "ship", "cat", "dog", "airplane"]
)

"""
    apply_aliases!(df; col="modelname", d=MODEL_ALIAS)

Renames entries of column `col` according to dictionary `d`.
"""
function apply_aliases!(df::D; col="modelname", d=MODEL_ALIAS) where {D <: AbstractDataFrame}
    !(col in names(df)) && error("$col is not in the DataFrame columns.")
    df[:,col] = map(x -> get(d, x, x), df[:, col])
    df
end

"""
	convert_anomaly_class(ac, dataset="MNIST")

Returns string representation of anomaly class id `ac`.
"""
function convert_anomaly_class(ac, dataset="mnist")
	getindex(AC_CONVERSION[dataset], ac)
end