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

using Markdown
using InteractiveUtils

# ╔═╡ 6f9a94bd-dd28-4110-a7ca-0ed84e9c7c3f
begin
	import Pkg
    # activate the shared project environment
    Pkg.activate(Base.current_project())
	using Interpolations
	using Plots
	using Dierckx
	using GasChromatographySimulator
	using GasChromatographyTools
	using Plots
	using PlutoUI
	TableOfContents()
end

# ╔═╡ 0b51a792-d9b8-4c29-928f-57aaf41f1c20
plotly()

# ╔═╡ 93ba6afc-4e9a-11ec-08d9-81c0a9bc502e
md"""
# Test of the RT lock functions
The test is necessary, because of failing (infinit loop in case of opt_ipt="linear", NaN abort in case of opt_itp="spline").
"""

# ╔═╡ 51ec0223-199a-4a25-8768-cd0b7e9f864d
begin
	function stretched_program(n::Float64, par::GasChromatographySimulator.Parameters)
		# stretch the temperature program in 'par' by a factor 'n'
		if isa(par, Array)==true
			error("Select an element of the array of GC-system parameters.")
		else
			new_tsteps = n.*par.prog.time_steps
			new_T_itp = GasChromatographySimulator.temperature_interpolation(new_tsteps, par.prog.temp_steps, par.prog.gf, par.col.L)
			new_pin_itp = GasChromatographySimulator.pressure_interpolation(new_tsteps, par.prog.pin_steps)
			new_pout_itp = GasChromatographySimulator.pressure_interpolation(new_tsteps, par.prog.pout_steps)
			new_prog = GasChromatographySimulator.Program(new_tsteps, par.prog.temp_steps, par.prog.pin_steps, par.prog.pout_steps, par.prog.gf, par.prog.a_gf, new_T_itp, new_pin_itp, new_pout_itp)
			new_par = GasChromatographySimulator.Parameters(par.col, new_prog, par.sub, par.opt)
			return new_par
		end
	end

	function initial_n(n::Float64, tR_lock::Float64, ii::Int, par::GasChromatographySimulator.Parameters)
	    par_n = stretched_program(n, par)
	    sol_n = GasChromatographySimulator.solve_system_multithreads(par_n)
	    tR_n = sol_n[ii].u[end][1]
	    if n>1
	        while tR_n-tR_lock<0 && n<130.0
	            n = n*2.0
	            par_n = stretched_program(n, par)
	            sol_n = GasChromatographySimulator.solve_system_multithreads(par_n)
	            tR_n = sol_n[ii].u[end][1]
	        end
	    elseif n<1
	        while tR_n-tR_lock>0 && n>0.01
	            n = n*0.5
	            par_n = stretched_program(n, par)
	            sol_n = GasChromatographySimulator.solve_system_multithreads(par_n)
	            tR_n = sol_n[ii].u[end][1]
	        end
	    end
	    if n>130.0
	        error("The choosen retention time for locking is to big.")
	    elseif n<0.01
	        error("The choosen retention time for locking is to small.")
	    else
	        return n
	    end
	end

	function RT_locking(par::GasChromatographySimulator.Parameters, tR_lock::Float64, tR_tol::Float64, solute_RT::String; opt_itp="linear")
		# estimate the factor 'n' for the temperature program to achieve the retention time 'tR_lock' for 'solute_RT' with the GC-system defined by 'par'
		if isa(par, Array)==true
			error("Select an element of the array of GC-systems.")
		else
			# find 'solute_RT' in the substances of 'par'
			name = Array{String}(undef, length(par.sub))
			for i=1:length(par.sub)
				name[i] = par.sub[i].name
			end
			ii = findfirst(name.==solute_RT)
			# calculate the retention time for the original (un-stretched) program
			sol₀ = GasChromatographySimulator.solving_odesystem_r(par.col, par.prog, par.sub[ii], par.opt)
			tR₀ = sol₀.u[end][1]
			# start value for the factor 'n'
			if tR₀-tR_lock<0
				n₁ = initial_n(2.0, tR_lock, ii, par)
			else
				n₁ = initial_n(0.5, tR_lock, ii, par)
			end
			# using a recursive function to estimate 'n'
			n = recur_RT_locking(n₁, [1.0], [tR₀], par, tR_lock, tR_tol, ii; opt_itp=opt_itp)
		end
		return n
	end

	function recur_RT_locking(n::Float64, n_vec::Array{Float64,1}, tR_vec::Array{Float64,1}, par::GasChromatographySimulator.Parameters, tR_lock::Float64, tR_tol::Float64, ii::Int64; opt_itp="linear")
		# recursive function to find the factor 'n' for the temperature program to achieve the retention time 'tR_lock' for solute index 'ii' with the GC-system defined by 'par'
		# calculate the retention time with the input guess 'n'
		par₁ = stretched_program(n, par)
		if par.opt.odesys==true
			sol₁ = GasChromatographySimulator.solving_odesystem_r(par₁.col, par₁.prog, par₁.sub[ii], par₁.opt)
			tR₁ = sol₁.u[end][1]
		else
			sol₁ = GasChromatographySimulator.solving_migration(par₁.col, par₁.prog, par₁.sub[ii], par₁.opt)
			tR₁ = sol₁.u[end]
		end
		if abs(tR₁-tR_lock)<tR_tol
			# if retention time is less than 'tR_tol' from 'tR_lock' we found the factor 'n'
			return n
		else
			# estimate a new factor 'new_n' by linear interpolation of the factors 'n_vec' + 'n' over the corresponding retention times
			new_n_vec = sort([n_vec; n])
			new_tR_vec = sort([tR_vec; tR₁])
			if opt_itp=="spline"
				# Dierckx.jl
				if length(new_tR_vec)<4
					k = length(new_tR_vec)-1
				else
					k = 3
				end
				itp = Spline1D(new_tR_vec, new_n_vec, k=k)
			else # opt_itp=="linear"
				# Interpolations.jl
				itp = LinearInterpolation(sort([tR_vec; tR₁]), sort([n_vec; n]))
			end
			new_n = itp(tR_lock)
			#println("new_n=$(new_n), tR₁=$(tR₁)")
			# use the new factor 'new_n' and call the recursive function again
			return recur_RT_locking(new_n, new_n_vec, new_tR_vec, par, tR_lock, tR_tol, ii; opt_itp=opt_itp)
		end
	end

	md"""
	## Copies of the RT lock functions
	In the package GasChromatographyTools.jl the functions where adapted to solve this problem.
	"""
end

# ╔═╡ 5e642f09-9a8f-4cca-b61f-b27c8433a2e5
begin
	opt = GasChromatographySimulator.Options(OwrenZen5(), 1e-6, 1e-3, "inlet", true)
	L = 4.0
	d = 0.1e-3
	df = 0.1e-6
	sp = "SLB5ms" # ["Rxi17SilMS" -> ok, "SLB5ms" -> ERR, "SPB50" -> ok, "Wax" -> ok, "DB5ms" -> ok, "Rxi5MS" -> ok, "genericLB", "genericJL"]
	gas = "He"
	col = GasChromatographySimulator.constructor_System(L, d, df, sp, gas)

	db_path = "/Users/janleppert/Documents/GitHub/Publication_GCsim/data/Databases/"
	db_file = "Database_append.csv"
	first_alkane = "C8"
	last_alkane = "C15"
	sub = GasChromatographySimulator.load_solute_database(db_path, db_file, sp, gas, [first_alkane, last_alkane], zeros(2), zeros(2))

	Tst = 273.15
	Tref = 150.0 + Tst
	pn = 101300.0
	pin = 300000.0 + pn
	pout = 0.0
	dimless_rate = 0.4
	Theat = 1000.0
	Tshift = 40.0
	tMref = GasChromatographySimulator.holdup_time(Tref, pin, pout, L, d, gas)
	rate = dimless_rate*30/tMref
	theat = Theat/rate
	Tstart = sub[1].Tchar-Tst-Tshift

	ΔT = 30.0
	α = -3.0
	prog0 = GasChromatographySimulator.constructor_Program([0.0, theat],[Tstart, Tstart+Theat], [pin, pin],[pout, pout],[ΔT, ΔT], [0.0, 0.0], [L, L], [α, α], opt.Tcontrol, L)

	par0 = GasChromatographySimulator.Parameters(col, prog0, sub, opt)

	tR_lock = 12*tMref
	tR_tol = 1e-3
	md"""
	## Settings
	The following settings produce the problem:
	"""
end

# ╔═╡ 3de48bb6-8eb7-4a33-9a98-d5fe3a19f6c6
md"""
## The ERROR
"""

# ╔═╡ 19d1abb6-32b6-414f-bc95-55635cbaa73a
n = RT_locking(par0, tR_lock, 1e-3, last_alkane; opt_itp="spline")

# ╔═╡ 2cad8185-4fc5-4009-98df-07a2be2133c6
md"""
## Reproducing the error by stepwise calculations
"""

# ╔═╡ 30207849-cafe-4ec6-af8d-ee7b2e2e6de0
begin
	# initial simulation
	sol₀ = GasChromatographySimulator.solving_odesystem_r(par0.col, par0.prog, par0.sub[2], par0.opt)
	tR₀ = sol₀.u[end][1]
end

# ╔═╡ 84e7e869-0fbf-4590-a8de-28855856661f
tR₀-tR_lock

# ╔═╡ c5ace7ce-cfa3-4c15-bcc1-8b66ad1c16e9
begin
	# first stretch of the program
	if tR₀-tR_lock<0
	    n₁ = initial_n(2.0, tR_lock, 2, par0)
	else
	    n₁ = initial_n(0.5, tR_lock, 2, par0)
	end
	# recur function
	par₁ = stretched_program(n₁, par0)
	sol₁ = GasChromatographySimulator.solving_odesystem_r(par₁.col, par₁.prog, par₁.sub[2], par₁.opt)
	tR₁ = sol₁.u[end][1]
	n_vec₁ = sort([1.0; n₁])
	tR_vec₁ = sort([tR₀; tR₁])
	# first interpolation
	itp₁ = Spline1D(tR_vec₁, n_vec₁, k=1)
end

# ╔═╡ c4d53222-03b1-4ce4-a49a-690945347432
tR₁-tR_lock

# ╔═╡ 6acfbf9b-f8ac-4483-8c93-17920d0d9f0e
begin
	# second stretch
	n₂ = itp₁(tR_lock)
	par₂ = stretched_program(n₂, par0)
	sol₂ = GasChromatographySimulator.solving_odesystem_r(par₂.col, par₂.prog, par₂.sub[2], par₂.opt)
	tR₂ = sol₂.u[end][1]
	n_vec₂ = sort([n_vec₁; n₂])
	tR_vec₂ = sort([tR_vec₁; tR₂])
	itp₂ = Spline1D(tR_vec₂, n_vec₂, k=2)

	p1 = plot([tR_lock, tR_lock], [0.95, 0.98], label="tR_lock")
	scatter!(p1, [tR₂, tR₂], [n₂, n₂], label="2")
	plot!(p1, 51.0:0.001:52.0, itp₂.(51.0:0.001:52.0), xlims=(51.0,52.0), ylims=(0.95, 0.98), label="itp₂")
	p1
end

# ╔═╡ 3b31f580-efe2-4a7a-89dc-228b38f2a71e
tR₂-tR_lock

# ╔═╡ b6c5126a-3d41-4809-ab9b-d554262e0668
begin
	n₃ = itp₂(tR_lock)
	par₃ = stretched_program(n₃, par0)
	sol₃ = GasChromatographySimulator.solving_odesystem_r(par₃.col, par₃.prog, par₃.sub[2], par₃.opt)
	tR₃ = sol₃.u[end][1]
	n_vec₃ = sort([n_vec₂; n₃])
	tR_vec₃ = sort([tR_vec₂; tR₃])
	k = length(tR_vec₃)-1
	itp₃ = Spline1D(tR_vec₃, n_vec₃, k=k)

	scatter!(p1, [tR₃, tR₃], [n₃, n₃], label="3")
	plot!(p1, 51.0:0.001:52.0, itp₃.(51.0:0.001:52.0), label="itp₃")
end

# ╔═╡ 6c5443bc-3dab-4317-99c9-1bb32b184fbc
tR₃-tR_lock

# ╔═╡ 3b1bb187-66f9-40bd-a0a5-a3c4e1a23819
begin
	n₄ = itp₃(tR_lock)
	par₄ = stretched_program(n₄, par0)
	sol₄ = GasChromatographySimulator.solving_odesystem_r(par₄.col, par₄.prog, par₄.sub[2], par₄.opt)
	tR₄ = sol₄.u[end][1]
	n_vec₄ = sort([n_vec₃; n₄])
	tR_vec₄ = sort([tR_vec₃; tR₄])
	itp₄ = Spline1D(tR_vec₄, n_vec₄, k=3)

	scatter!(p1, [tR₄, tR₄], [n₄, n₄], label="4")
	plot!(p1, 51.0:0.001:52.0, itp₄.(51.0:0.001:52.0), label="itp₄")
end

# ╔═╡ b81b2905-6756-445b-b1be-c54298df8b3f
tR₄-tR_lock

# ╔═╡ eb2d9f6a-8384-427a-a036-b4f0019d8251
md"""
The proposed settings seam to lead to jumping estimations around the searched retention time.
"""

# ╔═╡ 734bae50-12fa-4717-8446-addb963b8673
begin
	n₅ = itp₄(tR_lock)
	par₅ = GasChromatographyTools.stretched_program(n₅, par0)
	sol₅ = GasChromatographySimulator.solving_odesystem_r(par₅.col, par₅.prog, par₅.sub[2], par₅.opt)
	tR₅ = sol₅.u[end][1]
	n_vec₅ = sort([n_vec₄; n₅])
	tR_vec₅ = sort([tR_vec₄; tR₅])
	itp₅ = Spline1D(tR_vec₅, n_vec₅, k=3)

	scatter!(p1, [tR₅, tR₅], [n₅, n₅], label="5")
	plot!(p1, 51.0:0.001:52.0, itp₅.(51.0:0.001:52.0), xlims=(51.45,51.55), ylims=(0.969, 0.972), label="itp₅")
end

# ╔═╡ 8de716b9-ca56-42f5-aebf-39ad467f4613
tR₅-tR_lock

# ╔═╡ dadbefbf-a107-4c89-a4ef-0eb756517c1e
begin
	n₆ = itp₅(tR_lock)
	par₆ = GasChromatographyTools.stretched_program(n₆, par0)
	sol₆ = GasChromatographySimulator.solving_odesystem_r(par₆.col, par₆.prog, par₆.sub[2], par₆.opt)
	tR₆ = sol₆.u[end][1]
	n_vec₆ = sort([n_vec₅; n₆])
	tR_vec₆ = sort([tR_vec₅; tR₆])
	tR_vec₆.-tR_lock
	itp₆ = Spline1D(tR_vec₆, n_vec₆, k=3)

	scatter!(p1, [tR₆, tR₆], [n₆, n₆], label="6")
	plot!(p1, 51.0:0.001:52.0, itp₆.(51.0:0.001:52.0), xlims=(51.45,51.55), ylims=(0.969, 0.972), label="itp₆")
end

# ╔═╡ 313d1431-ec39-4d1e-91fd-e025efc1f5c3
tR₆-tR_lock

# ╔═╡ fab59b1e-ba32-49d3-b0f1-564db037400c
begin
	n₇ = itp₆(tR_lock)
	par₇ = GasChromatographyTools.stretched_program(n₇, par0)
	sol₇ = GasChromatographySimulator.solving_odesystem_r(par₇.col, par₇.prog, par₇.sub[2], par₇.opt)
	tR₇ = sol₇.u[end][1]
	n_vec₇ = sort([n_vec₆; n₇])
	tR_vec₇ = sort([tR_vec₆; tR₇])
	tR_vec₇.-tR_lock
	itp₇ = Spline1D(tR_vec₇, n_vec₇, k=3)

	scatter!(p1, [tR₇, tR₇], [n₇, n₇], label="7")
	plot!(p1, 51.0:0.001:52.0, itp₇.(51.0:0.001:52.0), xlims=(51.47,51.49), ylims=(0.9699, 0.9701), label="itp₇")
end

# ╔═╡ 25725568-ae9a-4ab0-8f98-87fec12c867a
begin
	n₈ = itp₇(tR_lock)
	par₈ = GasChromatographyTools.stretched_program(n₈, par0)
	sol₈ = GasChromatographySimulator.solving_odesystem_r(par₈.col, par₈.prog, par₈.sub[2], par₈.opt)
	tR₈ = sol₈.u[end][1]
	n_vec₈ = sort([n_vec₇; n₈])
	tR_vec₈ = sort([tR_vec₇; tR₈])
	tR_vec₈.-tR_lock
	itp₈ = Spline1D(tR_vec₈, n_vec₈, k=3)

	scatter!(p1, [tR₈, tR₈], [n₈, n₈], label="8")
	plot!(p1, 51.0:0.001:52.0, itp₈.(51.0:0.001:52.0), xlims=(51.47,51.49), ylims=(0.9699, 0.9701), label="itp₈")
end

# ╔═╡ 404c1d12-5217-485e-b3a6-6f024d29b544
md"""
The estimates continue to jump around tR_lock. The last estimate is further away than the one before.

Also, there seems to be some form of discontiuity of the simulation result around tR_lock, which produces the problem in the first place.
"""

# ╔═╡ 57634f00-45d5-4cf9-94f5-76319d4e5436
begin
	n₉ = itp₈(tR_lock)
	par₉ = GasChromatographyTools.stretched_program(n₉, par0)
	sol₉ = GasChromatographySimulator.solving_odesystem_r(par₉.col, par₉.prog, par₉.sub[2], par₉.opt)
	tR₉ = sol₉.u[end][1]
	n_vec₉ = sort([n_vec₈; n₉])
	tR_vec₉ = sort([tR_vec₈; tR₉])
	tR_vec₉.-tR_lock
	itp₉ = Spline1D(tR_vec₉, n_vec₉, k=3)

	scatter!(p1, [tR₉, tR₉], [n₉, n₉], label="9")
	plot!(p1, 51.0:0.001:52.0, itp₉.(51.0:0.001:52.0), xlims=(51.47,51.49), ylims=(0.9699, 0.9701), label="itp₉")
end

# ╔═╡ ae803a26-c436-497a-b245-a0afec92e46f
md"""
Make simulations for a range of stretch factors ``n``.
"""

# ╔═╡ 1983258e-e84b-4f39-9ce8-0e20e78a0893
begin
	nn = [0.969970, 0.969971, 0.969972, 0.969973, 0.969974, 0.969975, 0.969976, 0.969977]
	ttR = Array{Float64}(undef, length(nn))
	for i=1:length(nn)
	    pars = stretched_program(nn[i], par0)
	    sols = GasChromatographySimulator.solving_odesystem_r(pars.col, pars.prog, pars.sub[2], pars.opt)
	    ttR[i] = sols.u[end][1]
	end
	p2 = plot([tR_lock, tR_lock], [0.95, 0.98], label="tR_lock")
	plot!(p2, ttR, nn, line=(2,:solid), markers=:square, xlims=(51.479, 51.484), ylims=(0.969969,0.96998))
end

# ╔═╡ 89cbd5c0-0df6-43d3-bd5c-d3d55d669f33
begin
	nnn = 0.969974.+collect(0.0000001:0.0000001:0.0000009)
	tttR = Array{Float64}(undef, length(nnn))
	for i=1:length(nnn)
	    pars = stretched_program(nnn[i], par0)
	    sols = GasChromatographySimulator.solving_odesystem_r(pars.col, pars.prog, pars.sub[2], pars.opt)
	    tttR[i] = sols.u[end][1]
	end
	plot!(p2, tttR, nnn, ylims=(0.9699739,0.9699751), markers=:circle)
end

# ╔═╡ f541b210-3ef8-4003-b358-065e5c5949ad
par8 = stretched_program(nnn[8], par0)

# ╔═╡ 524b256b-0c1e-48bc-ac0b-2ee7fc1eab2b
par9 = stretched_program(nnn[9], par0)

# ╔═╡ 6d15c2a4-6d1e-4b89-b9e8-ecff004c4730
sol8 = GasChromatographySimulator.solving_odesystem_r(par8.col, par8.prog, par8.sub[2], par8.opt)

# ╔═╡ 665f302d-204b-47ac-90df-c5979350707c
sol9 = GasChromatographySimulator.solving_odesystem_r(par9.col, par9.prog, par9.sub[2], par9.opt)

# ╔═╡ 56095d71-6169-44b3-89d1-7ea7f1b6ddfb
sol8.destats

# ╔═╡ a51cd1dc-50bb-4444-a0f7-c0e4229b1257
sol9.destats

# ╔═╡ 225434d5-c753-4476-8f68-f5761a454852
sol8.retcode

# ╔═╡ 9a6f436e-e50e-4697-a1d1-3d2d43fc62fc
sol9.retcode

# ╔═╡ 8a6eccf0-8740-4d69-b96b-1248557d4c4d
begin
	nnnn = sort!(rand(0.9699745:0.000000001:0.9699755, 100))
	ttttR = Array{Float64}(undef, length(nnnn))
	for i=1:length(nnnn)
	    pars = stretched_program(nnnn[i], par0)
	    sols = GasChromatographySimulator.solving_odesystem_r(pars.col, pars.prog, pars.sub[2], pars.opt)
	    ttttR[i] = sols.u[end][1]
	end
	plot!(p2, ttttR, nnnn, ylims=(0.9699739,0.9699756), markers=:circle)
	md"""
	In the range around tR_lock the solution seems to be unstable. A small change in the temperature program leads to a bigger change in the retention time. There seem to be two branches of n(tR), which does not meet.

	$(embed_display(p2))
	"""
end

# ╔═╡ 08a0972a-55ff-42c9-838b-1a649afe9a46
md"""
## Decreased relative tolerance

The problem is not originated in the RT_lock algorithm but in the simulation itself. For certain settings a small change (like small strectch in the program) can result in bigger changes in retention time.
"""

# ╔═╡ 2d7da55e-2711-44a4-b8b9-bda6321a4c48
begin
	# decrease the relative tolerance
	opt_1 = GasChromatographySimulator.Options(OwrenZen5(), 1e-6, 1e-4, "inlet", true)
	par_1 = GasChromatographySimulator.Parameters(col, prog0, sub, opt_1)
	# repeat the simulation from above
	nnnn_1 = sort!(rand(0.9699745:0.000000001:0.9699755, 100))
	ttttR_1 = Array{Float64}(undef, length(nnnn_1))
	for i=1:length(nnnn_1)
	    pars = stretched_program(nnnn_1[i], par_1)
	    sols = GasChromatographySimulator.solving_odesystem_r(pars.col, pars.prog, pars.sub[2], pars.opt)
	    ttttR_1[i] = sols.u[end][1]
	end
	plot!(p2, ttttR_1, nnnn_1, ylims=(0.9699739,0.9699756), markers=:v, label="reltol=1e-4")
end

# ╔═╡ 89fa3139-1d10-4dd3-90b1-d9c0f65745c6
begin
	n_1 = RT_locking(par_1, tR_lock, 1e-3, last_alkane; opt_itp="spline")
	par_1_s = stretched_program(n_1, par_1)
	sol_1_s = GasChromatographySimulator.solving_odesystem_r(par_1_s.col, par_1_s.prog, par_1_s.sub[2], par_1_s.opt)
	tR_1_s = sol_1_s.u[end][1]
	scatter!(p2, [tR_1_s, tR_1_s], [n_1, n_1], ylims=(0.9699739, n_1*1.000001))
end

# ╔═╡ e9accaac-587f-498d-ab0c-a8064f573678
begin
	# simulation arround n_1
	n_range = (n_1-0.0001):0.000001:(n_1+0.0001)
	tR_range = Array{Float64}(undef, length(n_range))
	for i=1:length(n_range)
	    pars = stretched_program(n_range[i], par_1)
	    sols = GasChromatographySimulator.solving_odesystem_r(pars.col, pars.prog, pars.sub[2], pars.opt)
	    tR_range[i] = sols.u[end][1]
	end
	plot!(p2, tR_range, n_range, ylims=(n_range[1], n_range[end]), xlims=(tR_range[1], tR_range[end]))
end

# ╔═╡ 4cd795b1-efeb-430b-a025-77d0ee9d2b1b
md"""
## Modified Functions with relative difference
"""

# ╔═╡ 0bf4903c-074d-41b8-a4e8-5834ce8659d5
function recur_RT_locking_rel(n::Float64, n_vec::Array{Float64,1}, tR_vec::Array{Float64,1}, par::GasChromatographySimulator.Parameters, tR_lock::Float64, tR_tol::Float64, ii::Int64; opt_itp="linear")
	# recursive function to find the factor 'n' for the temperature program to achieve the retention time 'tR_lock' for solute index 'ii' with the GC-system defined by 'par'
	# calculate the retention time with the input guess 'n'
	par₁ = stretched_program(n, par)
	if par.opt.odesys==true
		sol₁ = GasChromatographySimulator.solving_odesystem_r(par₁.col, par₁.prog, par₁.sub[ii], par₁.opt)
		tR₁ = sol₁.u[end][1]
	else
		sol₁ = GasChromatographySimulator.solving_migration(par₁.col, par₁.prog, par₁.sub[ii], par₁.opt)
		tR₁ = sol₁.u[end]
	end
	if abs(tR₁-tR_lock)/tR_lock<tR_tol
		# if retention time is less than 'tR_tol' from 'tR_lock' we found the factor 'n'
		return n
	else
		# estimate a new factor 'new_n' by linear interpolation of the factors 'n_vec' + 'n' over the corresponding retention times
		new_n_vec = sort([n_vec; n])
		new_tR_vec = sort([tR_vec; tR₁])
		if opt_itp=="spline"
			# Dierckx.jl
			if length(new_tR_vec)<4
				k = length(new_tR_vec)-1
			else
				k = 3
			end
			itp = Spline1D(new_tR_vec, new_n_vec, k=k)
		else # opt_itp=="linear"
			# Interpolations.jl
			itp = LinearInterpolation(sort([tR_vec; tR₁]), sort([n_vec; n]))
		end
		new_n = itp(tR_lock)
		#println("new_n=$(new_n), tR₁=$(tR₁)")
		# use the new factor 'new_n' and call the recursive function again
		return recur_RT_locking_rel(new_n, new_n_vec, new_tR_vec, par, tR_lock, tR_tol, ii; opt_itp=opt_itp)
	end
end

# ╔═╡ be7d3bf5-d682-40e7-8161-4243354373f8
function RT_locking_rel(par::GasChromatographySimulator.Parameters, tR_lock::Float64, tR_tol::Float64, solute_RT::String; opt_itp="linear")
	# estimate the factor 'n' for the temperature program to achieve the retention time 'tR_lock' for 'solute_RT' with the GC-system defined by 'par'
	if isa(par, Array)==true
		error("Select an element of the array of GC-systems.")
	elseif tR_tol<par.opt.reltol
		error("The relative tolerance for retention time locking tR_tol is to small. Use tR_tol > par.opt.reltol ($(par.opt.reltol)).")
	else
		# find 'solute_RT' in the substances of 'par'
		name = Array{String}(undef, length(par.sub))
		for i=1:length(par.sub)
			name[i] = par.sub[i].name
		end
		ii = findfirst(name.==solute_RT)
		# calculate the retention time for the original (un-stretched) program
		sol₀ = GasChromatographySimulator.solving_odesystem_r(par.col, par.prog, par.sub[ii], par.opt)
		tR₀ = sol₀.u[end][1]
		# start value for the factor 'n'
		if tR₀-tR_lock<0
			n₁ = initial_n(2.0, tR_lock, ii, par)
		else
			n₁ = initial_n(0.5, tR_lock, ii, par)
		end
		# using a recursive function to estimate 'n'
		n = recur_RT_locking_rel(n₁, [1.0], [tR₀], par, tR_lock, tR_tol, ii; opt_itp=opt_itp)
	end
	return n
end

# ╔═╡ 8a7a24e3-3896-4012-8e50-473a70ec3a44
n_rel = RT_locking_rel(par0, tR_lock, 1e-3, last_alkane; opt_itp="spline")

# ╔═╡ 4cbd8054-d1c7-4132-ae61-ef09ce6ba9e7
par0_rel = stretched_program(n_rel, par0)

# ╔═╡ b842960f-ef61-4a70-8a17-5405ec3e5ed3
sol_rel = GasChromatographySimulator.solving_odesystem_r(par0_rel.col, par0_rel.prog, par0_rel.sub[2], par0_rel.opt)

# ╔═╡ 24269e34-fc1d-4ee2-87b6-afb6a3081596
tR_rel = sol_rel.u[end][1]

# ╔═╡ 0f7e36cb-c4a7-43b5-928f-d364b82d218b
tR_rel-tR_lock

# ╔═╡ 3bdf621e-969b-42b8-85f8-5ecd261646d3
(tR_rel-tR_lock)/tR_lock

# ╔═╡ 42a1fccd-e782-42fb-91de-18c8b33d507d
md"""
# End
"""

# ╔═╡ Cell order:
# ╠═6f9a94bd-dd28-4110-a7ca-0ed84e9c7c3f
# ╠═0b51a792-d9b8-4c29-928f-57aaf41f1c20
# ╠═93ba6afc-4e9a-11ec-08d9-81c0a9bc502e
# ╠═51ec0223-199a-4a25-8768-cd0b7e9f864d
# ╠═5e642f09-9a8f-4cca-b61f-b27c8433a2e5
# ╠═3de48bb6-8eb7-4a33-9a98-d5fe3a19f6c6
# ╠═19d1abb6-32b6-414f-bc95-55635cbaa73a
# ╠═2cad8185-4fc5-4009-98df-07a2be2133c6
# ╠═30207849-cafe-4ec6-af8d-ee7b2e2e6de0
# ╠═84e7e869-0fbf-4590-a8de-28855856661f
# ╠═c5ace7ce-cfa3-4c15-bcc1-8b66ad1c16e9
# ╠═c4d53222-03b1-4ce4-a49a-690945347432
# ╠═6acfbf9b-f8ac-4483-8c93-17920d0d9f0e
# ╠═3b31f580-efe2-4a7a-89dc-228b38f2a71e
# ╠═b6c5126a-3d41-4809-ab9b-d554262e0668
# ╠═6c5443bc-3dab-4317-99c9-1bb32b184fbc
# ╠═3b1bb187-66f9-40bd-a0a5-a3c4e1a23819
# ╠═b81b2905-6756-445b-b1be-c54298df8b3f
# ╟─eb2d9f6a-8384-427a-a036-b4f0019d8251
# ╠═734bae50-12fa-4717-8446-addb963b8673
# ╠═8de716b9-ca56-42f5-aebf-39ad467f4613
# ╠═dadbefbf-a107-4c89-a4ef-0eb756517c1e
# ╠═313d1431-ec39-4d1e-91fd-e025efc1f5c3
# ╠═fab59b1e-ba32-49d3-b0f1-564db037400c
# ╠═25725568-ae9a-4ab0-8f98-87fec12c867a
# ╠═404c1d12-5217-485e-b3a6-6f024d29b544
# ╠═57634f00-45d5-4cf9-94f5-76319d4e5436
# ╠═ae803a26-c436-497a-b245-a0afec92e46f
# ╠═1983258e-e84b-4f39-9ce8-0e20e78a0893
# ╠═89cbd5c0-0df6-43d3-bd5c-d3d55d669f33
# ╠═f541b210-3ef8-4003-b358-065e5c5949ad
# ╠═524b256b-0c1e-48bc-ac0b-2ee7fc1eab2b
# ╠═6d15c2a4-6d1e-4b89-b9e8-ecff004c4730
# ╠═665f302d-204b-47ac-90df-c5979350707c
# ╠═56095d71-6169-44b3-89d1-7ea7f1b6ddfb
# ╠═a51cd1dc-50bb-4444-a0f7-c0e4229b1257
# ╠═225434d5-c753-4476-8f68-f5761a454852
# ╠═9a6f436e-e50e-4697-a1d1-3d2d43fc62fc
# ╠═8a6eccf0-8740-4d69-b96b-1248557d4c4d
# ╟─08a0972a-55ff-42c9-838b-1a649afe9a46
# ╠═2d7da55e-2711-44a4-b8b9-bda6321a4c48
# ╠═89fa3139-1d10-4dd3-90b1-d9c0f65745c6
# ╠═e9accaac-587f-498d-ab0c-a8064f573678
# ╠═4cd795b1-efeb-430b-a025-77d0ee9d2b1b
# ╠═be7d3bf5-d682-40e7-8161-4243354373f8
# ╠═0bf4903c-074d-41b8-a4e8-5834ce8659d5
# ╠═8a7a24e3-3896-4012-8e50-473a70ec3a44
# ╠═4cbd8054-d1c7-4132-ae61-ef09ce6ba9e7
# ╠═b842960f-ef61-4a70-8a17-5405ec3e5ed3
# ╠═24269e34-fc1d-4ee2-87b6-afb6a3081596
# ╠═0f7e36cb-c4a7-43b5-928f-d364b82d218b
# ╠═3bdf621e-969b-42b8-85f8-5ecd261646d3
# ╠═42a1fccd-e782-42fb-91de-18c8b33d507d

ext: jl


text:
using POMDPs, Crux, Flux, POMDPGym, BSON
import POMDPPolicies:FunctionPolicy
import Distributions:Uniform
using Random
using Distributions
using Plots

## Pendulum
mdp = PendulumPOMDP()
as = [actions(mdp)...]
amin = [-2f0]
amax = [2f0]
rand_policy = FunctionPolicy((s) -> Float32.(rand.(Uniform.(amin, amax))))
S = state_space(mdp, σ=[6.3f0, 8f0])

# get expert trajectories
expert_trajectories = BSON.load("/home/anthonycorso/.julia/dev/Crux/examples/il/expert_data/pendulum.bson")[:data]
expert_perf = sum(expert_trajectories[:r]) / length(episodes(expert_trajectories))
expert_trajectories[:r] .= 1

# Define the networks we will use
QSA() = ContinuousNetwork(Chain(Dense(3, 64, relu), Dense(64, 64, relu), Dense(64, 1)))
QSA_SN(output=1) = ContinuousNetwork(Chain(DenseSN(3, 64, relu), DenseSN(64, 64, relu), DenseSN(64, 2), Dense(2,output)))
V() = ContinuousNetwork(Chain(Dense(2, 64, relu), Dense(64, 64, relu), Dense(64, 1)))
A() = ContinuousNetwork(Chain(Dense(2, 64, relu, init=Flux.orthogonal), Dense(64, 64, relu, init=Flux.orthogonal), Dense(64, 1, tanh, init=Flux.orthogonal), x -> 2f0 * x), 1)
G() = GaussianPolicy(A(), zeros(Float32, 1))

D_SN(output=1) = ContinuousNetwork(Chain(DenseSN(3, 12, relu), DenseSN(12, 2), Dense(2,output)))

function SAC_A()
    base = Chain(x -> x ./ [6.3f0, 8f0], Dense(2, 64, relu), Dense(64, 64, relu))
    mu = ContinuousNetwork(Chain(base..., Dense(64, 1)))
    logΣ = ContinuousNetwork(Chain(base..., Dense(64, 1)))
    SquashedGaussianPolicy(mu, logΣ)
end

## On-Policy GAIL - This currently doesn't work for some reason
𝒮_gail_on = OnPolicyGAIL(D=QSA_SN(), gan_loss=GAN_BCELoss(), 𝒟_demo=expert_trajectories, solver=REINFORCE, π=ActorCritic(G(), V()), S=S, N=1000000, ΔN=1024)
solve(𝒮_gail_on, mdp)

## Off-Policy GAIL
𝒮_gail = OffPolicyGAIL(D=D_SN(2),
                       𝒟_demo=expert_trajectories,
                       solver=TD3,
                       π=ActorCritic(A(), DoubleNetwork(QSA(), QSA())),
                       S=S,
                       ΔN=50,
                       N=30000,
                       buffer_size=Int(1e4),
                       c_opt=(batch_size=512, optimizer=ADAM(1e-3)),
                       a_opt=(optimizer=ADAM(1e-3),),
                       d_opt=(epochs=1, batch_size=256),
                       π_explore=GaussianNoiseExplorationPolicy(0.2f0, a_min=[-2.0], a_max=[2.0])
                       )
solve(𝒮_gail, mdp)

## Behavioral Cloning
𝒮_bc = BC(π=A(), 𝒟_demo=expert_trajectories, S=S, opt=(epochs=100,), log=(period=100,))
solve(𝒮_bc, mdp)

## Advil
𝒮_advil = AdVIL(π=ActorCritic(A(),QSA()), 𝒟_demo=expert_trajectories, S=S, a_opt=(epochs=1000, optimizer=ADAM(8f-4), batch_size=1024), c_opt=(optimizer=ADAM(8e-4),), max_steps=100, log=(period=100,))
solve(𝒮_advil, mdp)

## SQIL
𝒮_sqil = SQIL(π=ActorCritic(A(), DoubleNetwork(QSA(), QSA())),
              S=S,
              𝒟_demo=expert_trajectories,
              max_steps=100,
              N=30000,
              buffer_size=Int(1e4),
              c_opt=(batch_size=128, optimizer=ADAM(1e-3)),
              a_opt=(batch_size=128, optimizer=ADAM(1e-3)),
              solver=TD3,
              π_explore=GaussianNoiseExplorationPolicy(0.2f0, a_min=[-2.0], a_max=[2.0]))
solve(𝒮_sqil, mdp)

## Adril
Crux.set_crux_warnings(false)
𝒮_adril = AdRIL(π=ActorCritic(SAC_A(), DoubleNetwork(QSA(), QSA())),
              S=S,
              𝒟_demo=expert_trajectories,
              max_steps=100,
              N=30000,
              buffer_size=Int(1e4),
              c_opt=(batch_size=128, optimizer=ADAM(1e-3)),
              a_opt=(batch_size=128, optimizer=ADAM(1e-3)),
              π_explore=GaussianNoiseExplorationPolicy(0.2f0, a_min=[-2.0], a_max=[2.0]))
solve(𝒮_adril, mdp)

## ASAF
𝒮_ASAF = ASAF(π=G(),
              S=S,
              ΔN=2000,
              𝒟_demo=expert_trajectories,
              N=50000,
              max_steps=100,
              a_opt=(batch_size=256, optimizer=Flux.Optimise.Optimiser(Flux.ClipValue(1f0), ADAM(1e-3)), epochs=10))
solve(𝒮_ASAF, mdp)

p = plot_learning([𝒮_gail_on, 𝒮_gail, 𝒮_bc, 𝒮_advil, 𝒮_sqil, 𝒮_adril,𝒮_ASAF], title="Pendulum Swingup Imitation Learning Curves", labels=["On Policy GAIL", "Off-Policy GAIL", "BC", "AdVIL", "SQIL", "AdRIL", "ASAF"], legend=:right)
plot!(p, [1,50000], [expert_perf, expert_perf], color=:black, label="expert")

savefig("examples/il/pendulum_benchmark.pdf")


text:
module IterableTables

using Requires, IteratorInterfaceExtensions, TableTraits, TableTraitsUtils
using DataValues

# include("integrations/gadfly.jl")
# include("integrations/juliadb.jl")
include("integrations/generators.jl")

function __init__()
    @require DataFrames="a93c6f00-e57d-5684-b7b6-d8193f3e46c0" if !isdefined(DataFrames, :Tables)
        include("integrations/dataframes.jl")
    end
    @require StatsModels="3eaba693-59b7-5ba5-a881-562e759f1c8d" if !isdefined(StatsModels, :Tables)
        include("integrations/statsmodels.jl")
    end
    @require TimeSeries="9e3dc215-6440-5c97-bce1-76c03772f85e" if !isdefined(TimeSeries, :Tables)
        include("integrations/timeseries.jl")
    end
    @require Temporal="a110ec8f-48c8-5d59-8f7e-f91bc4cc0c3d" include("integrations/temporal.jl")
end

end # module


text:
 @testset "Nodes initialization algorithm n=3 entry = 3 output = 2" begin
    M = 1 # Constante para la función rampa
    # Bien definido para tamaño n = 2 y salida de dimensión 1
    f_regression(x,y,z)=[x*y-z,x]
    data_set_size =  6
    entry_dimension = 3
    output_dimension = 2
    # Número de neuronas
    n = data_set_size # Debe de ser mayor que 1 para que no de error
    X_train= rand(Float32, data_set_size, entry_dimension)
    Y_train::Matrix =  mapreduce(permutedims, vcat, map(x->f_regression(x...), eachrow(X_train)))

    h = nn_from_data(X_train, Y_train, n, M)

    # veamos que el tamaño de la salida es la adecuada
    @test size(h.W1) == (n,entry_dimension+1)
    @test size(h.W2) == (output_dimension,n)

    # Si ha sido bien construida:
    # Evaluar la red neuronal en los datos con los que se construyó
    # debería de resultar el valor de Y_train respectivo
    evaluar(x)=forward_propagation(h,
     RampFunction,x)

    for i in 1:n
        @test evaluar(X_train[i,:]) ≈ Y_train[i,:]
    end

end


text:
#############################################
"""
    eta_uniformgrid = interpsurface(amrgrid::Vector{VisClaw.AMRGrid}, topo::VisClaw.Topo)
    eta_uniformgrid = interpsurface(amr::VisClaw.AMR, topo::VisClaw.Topo)

Convert AMR tile data into uniform grid data using the SciPy scattered interpolation.\n
Interpolation of inundation height (land) is not supported. \n
"""
function interpsurface(amrgrid::Vector{VisClaw.AMRGrid}, topo::VisClaw.Topo)
    ## import
    scipyinterpolate = PyCall.pyimport("scipy.interpolate")

    ## make grid data from topo
    xmesh_all = repeat(topo.x', inner=(topo.nrows,1))
    ymesh_all = repeat(topo.y , inner=(1,topo.ncols))
    land = topo.elevation .>= 0.0

    ## get lims of the designated region
    xmin, xmax = extrema(topo.x)
    ymin, ymax = extrema(topo.y)

    ## scattered x y data
    ntile = length(amrgrid)
    x_all = empty([0.0])
    y_all = empty([0.0])
    z_all = empty([0.0])
    for i_tile = 1:ntile
        tile = amrgrid[i_tile]

        x1, x2, y1, y2 = VisClaw.getlims(tile)
        x2 < xmin && (continue)
        xmax < x1 && (continue)
        y2 < ymin && (continue)
        ymax < y1 && (continue)

        xmesh, ymesh = VisClaw.meshtile(tile)
        ind = isnan.(tile.eta)
        push!(x_all, xmesh[.!ind]...)
        push!(y_all, ymesh[.!ind]...)
        push!(z_all, tile.eta[.!ind]...)
    end

    ## nodata
    if length(z_all) == 0
        v_all = zeros(topo.nrows,topo.ncols)
    else
        v_all = scipyinterpolate.griddata([x_all y_all], z_all, (xmesh_all,ymesh_all), method="cubic")
    end

    v_all[land] .= NaN

    return v_all

end
#############################################

#############################################
function interpsurface(amrall::VisClaw.AMR, topo::VisClaw.Topo; timestep=1:amrall.nstep)
    eta_uniformgrid = [interpsurface(amrall.amr[k], topo) for k=timestep]
    eta_uniformgrid = cat(eta_uniformgrid...; dims=3)
    return eta_uniformgrid
end


text:
## header to provide surface and context
using Cairo
c = CairoRGBSurface(256,256);
cr = CairoContext(c);

save(cr);
set_source_rgb(cr,0.8,0.8,0.8);    # light gray
rectangle(cr,0.0,0.0,256.0,256.0); # background
fill(cr);
restore(cr);

save(cr);
## original example, following here
image = read_from_png("data/mulberry.png");
w = image.width;
h = image.height;

pattern = CairoPattern(image);
pattern_set_extend(pattern, Cairo.EXTEND_REPEAT);

translate(cr, 128.0, 128.0);
rotate(cr, pi / 4);
scale(cr, 1 / sqrt(2), 1 / sqrt(2));
translate(cr, -128.0, -128.0);

m = CairoMatrix(w/256.0 * 5.0,0,0,h/256.0 * 5.0,0,0);
#matrix_init_scale (&matrix, w/256.0 * 5.0, h/256.0 * 5.0);

set_matrix(pattern, m);
set_source(cr, pattern);

rectangle(cr, 0, 0, 256.0, 256.0);
fill(cr);

## mark picture with current date
restore(cr);
move_to(cr,0.0,12.0);
set_source_rgb(cr, 0,0,0);
show_text(cr,Libc.strftime(time()));
write_to_png(c,"sample_imagepattern.png");


text:
@defcomp emissions begin
    CCA     = Variable(index=[time])    # Cumulative indiustrial emissions
    E       = Variable(index=[time])    # Total CO2 emissions (GtC per year)
    EIND    = Variable(index=[time])    # Industrial emissions (GtC per year)

    etree   = Parameter(index=[time])   # Emissions from deforestation
    MIU     = Parameter(index=[time])   # Emission control rate GHGs
    sigma   = Parameter(index=[time])   # CO2-equivalent-emissions output ratio
    YGROSS  = Parameter(index=[time])   # Gross world product GROSS of abatement and damages (trillions 2005 USD per year)

    function run_timestep(p, v, d, t)
        # Define function for EIND
        v.EIND[t] = p.sigma[t] * p.YGROSS[t] * (1 - p.MIU[t])

        # Define function for E
        v.E[t] = v.EIND[t] + p.etree[t]

        # Define function for CCA
        if is_first(t)
            v.CCA[t] = v.EIND[t] * 10
        else
            v.CCA[t] = v.CCA[t - 1] + v.EIND[t] * 10
        end
    end
end


text:
#Qmodel keeps track of initial Q estimate (by logistic reg) and
#sequence of g estimates and fluctuations

module Qmodels

using ..Common, ..LReg
using StatsFuns

import ..LReg: linpred, predict
import StatsBase.nobs

export Qmodel, fluctuate!, defluctuate!, predict, linpred, nobs, weightedcovar, lastfluc, risk,
computefluc, valfluc

"""
A `Fluctuation` holds fluctuation covariates and weights. The weight times the fluctuation covariate is the so called "clever covariate" in the targeted
learning literature.
"""
type Fluctuation{T<:AbstractFloat}
    hA1::Vector{T}
    hA0::Vector{T}
    wts::Vector{T}
    epsilon::LR{T}
    weighted::Bool
    function Fluctuation(hA1, hA0, wts, epsilon, weighted)
        length(hA1) == length(hA0) == length(wts) || error(ArgumentError("lengths of hA1, hA0 and wts do not match"))
        new(hA1, hA0, wts, epsilon, weighted)
    end
end

Fluctuation{T<:AbstractFloat}(hA1::Vector{T}, hA0::Vector{T}, wts::Vector{T},
                              epsilon::LR{T}, weighted::Bool) =
    Fluctuation{T}(hA1, hA0, wts, epsilon, weighted)

nobs(fluc::Fluctuation) = length(fluc.hA1)

"""Computes the weight times flucutation covariate, sometimes called the "clever covariate"

** Arguments **

* `fluc` - fluctuation
* `a` - treatment vector
"""
function weightedcovar{T<:AbstractFloat}(fluc::Fluctuation{T}, a::Vector{T})
    ifelse(a.==1, fluc.hA1, fluc.hA0) .* fluc.wts
end

"""
Computes the fluctuated linear predictor given an offset

** Arguments **

* `fluc` - fluctuation
* `a` - treatment vector

** Keyword Arguments **

* `offset` - offset term, defaults to zeros, but the default should generally not be used
"""
function linpred{T<:AbstractFloat}(fluc::Fluctuation{T}, a::Vector{T}; offset::Vector{T}=zeros(T, 0))
    hAA = ifelse(a.==1, fluc.hA1, fluc.hA0)
    linpred(fluc.epsilon, hAA, offset=offset)
end

"""
The `Qmodel` type represents an estimate of \(\bar{Q}(a,w) = E(Y\mid A=a, W=w)\).

Initial estimate and fluctuations.
"""
type Qmodel{T<:AbstractFloat}
    logitQnA1::Vector{T} # \logit\bar{Q}_{n}(A=1, W)
    logitQnA0::Vector{T} # \logit\bar{Q}_{n}(A=0, W)
    flucseq::Vector{Fluctuation{T}} #flucseq[i] = ith fluctuation LR
    function Qmodel(logitQnA1, logitQnA0)
        length(logitQnA1) == length(logitQnA0) || error(ArgumentError("lengths of logitQnA1 and logitQnA0 do not match"))
        new(logitQnA1, logitQnA0, Fluctuation{T}[])
    end
end

Qmodel{T<:AbstractFloat}(logitQnA1::Vector{T}, logitQnA0::Vector{T}) = Qmodel{T}(logitQnA1, logitQnA0)

nobs(q::Qmodel) = length(q.logitQnA1)

"Returns the last fluctuation of `q`"
function lastfluc(q::Qmodel)
    length(q.flucseq) > 0 || error("q has not been fluctuated")
    q.flucseq[end]
end

# Qmodel{T<:AbstractFloat}(logitQnA1::Vector{T}, logitQnA0::Vector{T}, param::Parameter{T})=
#     Qmodel{T}(logitQnA1, logitQnA0, param, Vector{T}[], Vector{LR{Float64}}[])
"""Computes the linear predictor of an estimated `Qmodel` given `a`

** Arguments **

* `q` - `Qmodel`
* `a` - vector of treatments

"""
function linpred{T<:AbstractFloat}(q::Qmodel{T}, a::Vector{T})
    r = ifelse(a .==1, q.logitQnA1, q.logitQnA0)
    for fluc in q.flucseq
        r = linpred(fluc, a, offset=r)
    end
    r
end

"""Computes the prediction of an estimated `Qmodel` given `a`

** Arguments **

* `q` - `Qmodel`
* `a` - vector of treatments

"""
predict{T<:AbstractFloat}(q::Qmodel{T}, a::Vector{T}) = logistic(linpred(q,a))

function compute_h_wts(param::Parameter, gn1, A; weighted::Bool=false)
    n = length(A)
    n == length(gn1) || error()
    if weighted
        gna = ifelse(A.==1, gn1, (1 .- gn1))
        wts = 1 ./ gna
        hA1 = fluccovar(param, ones(n))
        hA0 = fluccovar(param, zeros(n))
    else
        wts = ones(n)
        hA1 = fluccovar(param, ones(n)) ./ gn1
        hA0 = fluccovar(param, zeros(n)) ./ (1 .- gn1)
    end
    (hA1, hA0, wts)
end

"""
Computes a new fluctuation for a (possible already fluctuated) initial `Qmodel`

** Arguments **

* `q` - `Qmodel`
* `param` - target parameter
* `gn1` - estimated propensity score
* `A` - observed treatments
* `Y` - observed outcomes

** Keyword Arguments **

* `weighted` - A boolean, `false` by default
If `false`, the "standard" fluctuation is performed, there the fluctuation covariate is divided
by \(g_n(A_i\mid W_i)\), and weights are set to 1.
If `true`, the fluctuation covariate is used directly, and weights are set to the reciprical
of \(g_n(A_i\mid W_i)\) where \(g_n(A_i \mid W_i)\) is computed as `gn1[i]` if `A[i]` is 1, or `1-gn1[i]` otherwise.
"""
function computefluc(q::Qmodel, param::Parameter, gn1, A, Y; weighted::Bool=false)
    #computes fluctuation for a given q and g
    hA1, hA0, wts = compute_h_wts(param, gn1, A, weighted=weighted)
    hAA = ifelse(A.==1, hA1, hA0)
    offset = linpred(q, A)
    epsilon = lreg(hAA, Y, offset=offset, wts=wts)
    Fluctuation(hA1, hA0, wts, epsilon, weighted)
end

function valfluc(fluc::Fluctuation, param::Parameter, gn1, A; weighted::Bool=false)
    hA1, hA0, wts = compute_h_wts(param, gn1, A, weighted=fluc.weighted)
    epsilon = fluc.epsilon
    Fluctuation(hA1, hA0, wts, epsilon, weighted)
end

function fluctuate!{T<:AbstractFloat}(q::Qmodel{T}, fluc::Fluctuation{T})
    nobs(q) == nobs(fluc) || error(ArgumentError("nobs of q and fluc do not match"))
    push!(q.flucseq, fluc)
    q
end

"""
Computes and adds a new fluctuation to a `Qmodel`

** Arguments **

* `q` - `Qmodel`
* `param` - target parameter
* `gn1` - estimated propensity score
* `A` - observed treatments
* `Y` - observed outcomes

** Keyword Arguments **

* `weighted` - A boolean, `false` by default.
See documentation for `computefluc` for details.
"""
function fluctuate!{T<:AbstractFloat}(q::Qmodel{T}, param::Parameter{T}, gn1::Vector{T},
                                      A::Vector{T}, Y::Vector{T};
                                      weighted::Bool=false)
    fluctuate!(q, computefluc(q, param, gn1, A, Y, weighted=weighted))
end

function defluctuate!(q::Qmodel)
    pop!(q.flucseq)
end

end


text:
using TopOpt.TopOptProblems: getE
using TopOpt.TrussTopOptProblems: truss_reinit!

@params mutable struct TrussElementKσ{T} <: AbstractFunction{T}
    problem::TrussProblem
    Kσes::AbstractVector{<:AbstractMatrix{T}}
    EALγ_s::AbstractVector{<:AbstractVector{T}}
    δmat_s::AbstractVector{<:AbstractMatrix{T}}
    L_s::AbstractVector{T}
    global_dofs::AbstractVector{<:Integer}
end

function Base.show(::IO, ::MIME{Symbol("text/plain")}, ::TrussElementKσ)
    return println("TopOpt element stress stiffness matrix (Kσ_e) construction function")
end

function TrussElementKσ(
    problem::TrussProblem{xdim,T}, solver::AbstractFEASolver
) where {xdim,T}
    Es = getE(problem)
    As = getA(problem)
    dh = problem.ch.dh

    # usually ndof_pc = xdim * n_basefuncs
    cellvalues = solver.elementinfo.cellvalues
    ndof_pc = ndofs_per_cell(dh)
    n_basefuncs = getnbasefunctions(cellvalues)
    @assert ndof_pc == xdim * n_basefuncs "$ndof_pc, $n_basefuncs"
    @assert n_basefuncs == 2

    global_dofs = zeros(Int, ndof_pc)
    δmat = zeros(T, ndof_pc, ndof_pc)
    Kσes = Matrix{T}[]
    EALγ_s = Vector{T}[]
    δmat_s = Matrix{T}[]
    L_s = T[]

    for (cellidx, cell) in enumerate(CellIterator(dh))
        truss_reinit!(cellvalues, cell, As[cellidx])
        celldofs!(global_dofs, dh, cellidx)
        E = Es[cellidx]
        A = As[cellidx]
        L = norm(cell.coords[1] - cell.coords[2])
        R = compute_local_axes(cell.coords[1], cell.coords[2])
        # local axial projection operator (local axis transformation)
        γ = vcat(-R[:, 1], R[:, 1])
        push!(EALγ_s, (E * A / L) * γ)

        fill!(δmat, 0.0)
        for i in 2:size(R, 2)
            δ = vcat(-R[:, i], R[:, i])
            # @assert δ' * γ ≈ 0
            δmat .+= δ * δ'
        end
        push!(δmat_s, copy(δmat))
        push!(L_s, T(L))

        push!(Kσes, zeros(T, ndof_pc, ndof_pc))
    end
    return TrussElementKσ(problem, Kσes, EALγ_s, δmat_s, L_s, global_dofs)
end

"""
    (eksig::TrussElementKσ)(u_e::AbstractVector, x_e::Number, ci::Integer)

Compute the stress (geometric) stiffness matrix for **truss element* index `ci` (axial bar element, no bending or torsion),
with nodal deformation `u_e` and density `x`.

See Matrix formulation defined in eq (3) and (4) in https://people.duke.edu/~hpgavin/cee421/truss-finite-def.pdf.
This matrix formulation is equivalent to the one used in

    M. Kočvara, “On the modelling and solving of the truss design problem with global stability constraints,” Struct Multidisc Optim, vol. 23, no. 3, pp. 189–203, Apr. 2002, doi: 10/br35mf.

Note:
Bar axial force is computed using first-order approx.
A better approx would be: EA/L * (u3-u1 + 1/(2*L0)*(u4-u2)^2) = EA/L * (γ'*u + 1/2*(δ'*u)^2)
see: https://people.duke.edu/~hpgavin/cee421/truss-finite-def.pdf
"""
function (eksig::TrussElementKσ)(u_e::AbstractVector, x_e::Number, ci::Integer)
    @unpack EALγ_s, δmat_s, L_s = eksig
    # x_e scales the cross section
    q_cell = x_e * EALγ_s[ci]' * u_e
    return q_cell / L_s[ci] * δmat_s[ci]
end

function (eksig::TrussElementKσ{T})(u::AbstractVector{T}, x::AbstractVector{T}) where {T}
    @unpack problem, Kσes, global_dofs = eksig
    dh = problem.ch.dh
    @assert getncells(dh.grid) == length(x)
    @assert ndofs(dh) == length(u)
    for ci in 1:length(x)
        celldofs!(global_dofs, dh, ci)
        Kσes[ci] = eksig(u[global_dofs], x[ci], ci)
    end
    return copy(Kσes)
end

function ChainRulesCore.rrule(eksig::TrussElementKσ{T}, u, x) where {T}
    @unpack problem, Kσes, global_dofs = eksig
    dh = problem.ch.dh
    Kσes = eksig(u, x)
    function pullback_fn(Δ)
        Δu = zeros(T, size(u))
        Δx = zeros(T, size(x))
        for ci in 1:length(x)
            celldofs!(global_dofs, dh, ci)
            function vec_eksig_fn(ux_vec)
                u_e = ux_vec[1:(end - 1)]
                x_e = ux_vec[end]
                return vec(eksig(u_e, x_e, ci))
            end
            jac_cell = ForwardDiff.jacobian(vec_eksig_fn, [u[global_dofs]; x[ci]])
            jtv = jac_cell' * vec(Δ[ci])
            Δu[global_dofs] += jtv[1:(end - 1)]
            Δx[ci] = jtv[end]
        end
        return Tangent{typeof(eksig)}(;
            problem=NoTangent(),
            Kσes=Δ,
            EALγ_s=NoTangent(),
            δmat_s=NoTangent(),
            L_s=NoTangent(),
            global_dofs=NoTangent(),
        ),
        Δu,
        Δx
    end
    return Kσes, pullback_fn
end

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

# TODO ElementKσ for volumetic cells


text:
# Translating
# https://proplot.readthedocs.io/en/latest/subplots.html#Automatic-figure-sizing
# to Julia

# Import proplot using PyCall
using PyCall
proplot = pyimport("proplot")

# Changing the reference subplot in the presence of unequal width/height ratios
suptitle = "Effect of reference subplot on figure size"
for ref in (1, 2)
    f, axs = proplot.subplots(
        ref=ref, nrows=3, ncols=3, wratios=(3, 2, 2),
        axwidth=1.1,
    )
    axs[ref].format(
        suptitle=suptitle,
        title="reference axes", titleweight="bold",
        titleloc="uc", titlecolor="red9"
    )
    f.savefig("figs/Automatic-figure-sizing-3a-$ref.svg", transparent=false)
end

# Changing the reference subplot in a complex grid
for (i,ref) in enumerate((3, 2))
    f, axs = proplot.subplots(
        [[1, 1, 2], [3, 4, 4]],
        hratios=(1, 1.5), wratios=(3, 2, 2),
        ref=ref, axwidth=1.1, span=false
    )
    axs[ref].format(
        suptitle=suptitle,
        title="reference axes", titleweight="bold",
        titleloc="uc", titlecolor="red9"
    )
    f.savefig("figs/Automatic-figure-sizing-3b-$i.svg", transparent=false)
end

max_stars_repo_name: natgeo-wong/JuliaProPlotExamples

max_issues_repo_name: natgeo-wong/JuliaProPlotExamples

max_forks_repo_name: natgeo-wong/JuliaProPlotExamples


text:
using VortexDistributions

# simple test field
Nx = 400; Ny = 400
Lx = 200; Ly = Lx
x = LinRange(-Lx/2,Lx/2, Nx+1)[1:end-1]; y = LinRange(-Ly/2,Ly/2, Ny+1)[1:end-1];
psi0 = one.(x*y') |> complex
psi = Torus(copy(psi0),x,y)
xp = 1.133
yp = 2.787
vp = PointVortex(xp,yp,1)
sp = ScalarVortex(vp)
vortex!(psi,sp)

ψ = psi.ψ
phase = angle.(ψ)

diffx,diffy = phase_jumps(phase, 1),phase_jumps(phase, 2)
a = zero(phase)
b = zero(phase)
phase_jumps!(a, phase)
@test a == diffx
phase_jumps!(b, phase, 2)
@test b == diffy


text:
# Basic scitypes

ST.scitype(::Integer,        ::DefaultConvention) = Count
ST.scitype(::AbstractFloat,  ::DefaultConvention) = Continuous
ST.scitype(::AbstractString, ::DefaultConvention) = Textual
ST.scitype(::TimeType,       ::DefaultConvention) = ScientificTimeType
ST.scitype(::Time    ,       ::DefaultConvention) = ScientificTime
ST.scitype(::Date,           ::DefaultConvention) = ScientificDate
ST.scitype(::DateTime,       ::DefaultConvention) = ScientificDateTime

# Image scitypes

ST.scitype(img::Arr{<:Gray,2}, ::DefaultConvention) = GrayImage{size(img)...}
ST.scitype(img::Arr{<:AbstractRGB,2}, ::DefaultConvention) =
ColorImage{size(img)...}

# Persistence diagrams

# ST.scitype(::PersistenceDiagram, ::DefaultConvention) = PersistenceDiagram

# CategoricalArray scitype

function ST.scitype(c::Cat, ::DefaultConvention)
    nc = length(levels(c.pool))
    return ifelse(c.pool.ordered, OrderedFactor{nc}, Multiclass{nc})
end

const CatArrOrSub{T, N} =
    Union{CategoricalArray{T, N}, SubArray{T, N, <:CategoricalArray}}

function ST.scitype(A::CatArrOrSub{T,N}, ::DefaultConvention) where {T,N}
    nlevels = length(levels(A))
    S = ifelse(isordered(A), OrderedFactor{nlevels}, Multiclass{nlevels})
    T >: Missing && (S = Union{S,Missing})
    return AbstractArray{S,N}
end

# Table scitype

function ST.scitype(X, ::DefaultConvention, ::Val{:table}; kw...)
    Xcol = Tables.columns(X)
    col_names = Tables.columnnames(Xcol)
    types = map(col_names) do name
        scitype(Tables.getcolumn(Xcol, name); kw...)
    end
    return Table{Union{types...}}
end

# Distributions

const Dist = Distributions

scalar_scitype(::Type) = Unknown
scalar_scitype(::Type{<:Dist.Discrete}) = Count
scalar_scitype(::Type{<:Dist.Continuous}) = Continuous
space_scitype(variate_form::Type{<:Dist.ArrayLikeVariate{0}},
              value_support) =
                  scalar_scitype(value_support)
space_scitype(variate_form::Type{<:Dist.ArrayLikeVariate{N}},
              value_support) where N =
                  AbstractArray{scalar_scitype(value_support),N}

ST.scitype(::Distributions.Sampleable{F,S}) where {F,S} =
    Sampleable{space_scitype(F,S)}

ST.scitype(::Distributions.Distribution{F,S}) where {F,S} =
    Density{space_scitype(F,S)}

# Scitype for fast array broadcasting

ST.Scitype(::Type{<:Integer},            ::DefaultConvention) = Count
ST.Scitype(::Type{<:AbstractFloat},      ::DefaultConvention) = Continuous
ST.Scitype(::Type{<:AbstractString},     ::DefaultConvention) = Textual
ST.Scitype(::Type{<:TimeType},           ::DefaultConvention) = ScientificTimeType
ST.Scitype(::Type{<:Date},               ::DefaultConvention) = ScientificDate
ST.Scitype(::Type{<:Time},               ::DefaultConvention) = ScientificTime
ST.Scitype(::Type{<:DateTime},           ::DefaultConvention) = ScientificDateTime
# ST.Scitype(::Type{<:PersistenceDiagram}, ::DefaultConvention) = PersistenceDiagram


text:
type DualAveragingMCTune <: MCTunerState
  step::Real # Stepsize of MCMC iteration (for ex leapfrog in HMC or drift stepsize in MALA)
  λ::Real
  μ::Real
  nleaps::Integer
  εbar::Real
  hbar::Real
  hweight::Real
  εweight::Real
  accepted::Integer # Number of accepted MCMC samples during current tuning period
  proposed::Integer # Number of proposed MCMC samples during current tuning period
  totproposed::Integer # Total number of proposed MCMC samples during burnin
  rate::Real # Observed acceptance rate over current tuning period

  function DualAveragingMCTune(
    step::Real,
    λ::Real,
    μ::Real,
    nleaps::Integer,
    εbar::Real,
    hbar::Real,
    hweight::Real,
    εweight::Real,
    accepted::Integer,
    proposed::Integer,
    totproposed::Integer,
    rate::Real
  )
    @assert (step > 0 || isnan(step)) "Stepsize of MCMC iteration should be positive or not a number (NaN)"
    @assert (εbar > 0 || isnan(εbar)) "εbar should be positive or not a number (NaN)"
    @assert accepted >= 0 "Number of accepted MCMC samples should be non-negative"
    @assert proposed >= 0 "Number of proposed MCMC samples should be non-negative"
    @assert totproposed >= 0 "Total number of proposed MCMC samples should be non-negative"
    if !isnan(rate)
      @assert 0 <= rate <= 1 "Observed acceptance rate should be in [0, 1]"
    end
    new(step, λ, μ, nleaps, εbar, hbar, hweight, εweight, accepted, proposed, totproposed, rate)
  end
end

DualAveragingMCTune(;
  step::Real=1.,
  λ::Real=2.,
  εbar::Real=1.,
  hbar::Real=0.,
  accepted::Integer=0,
  proposed::Integer=0,
  totproposed::Integer=0
) =
  DualAveragingMCTune(step, λ, NaN, 0, εbar, hbar, NaN, NaN, accepted, proposed, totproposed, NaN)

immutable DualAveragingMCTuner <: MCTuner
  targetrate::Real # Target acceptance rate
  nadapt::Integer
  ε0bar::Real
  h0bar::Real
  γ::Real
  t0::Int
  κ::Real
  period::Integer # Tuning period over which acceptance rate is reported in verbose mode
  verbose::Bool # If verbose=false then the tuner is silent, else it is verbose

  function DualAveragingMCTuner(
    targetrate::Real,
    nadapt::Integer,
    ε0bar::Real,
    h0bar::Real,
    γ::Real,
    t0::Int,
    κ::Real,
    period::Integer,
    verbose::Bool
  )
    @assert 0 < targetrate < 1 "Target acceptance rate should be between 0 and 1"
    @assert nadapt > 0 "Number of adaptation steps should be positive"
    @assert ε0bar > 0 "ε0bar should be positive"
    @assert period > 0 "Period over which acceptance rate is reported in verbose mode should be positive"
    @assert t0 > 0 "t0 should be positive"
    new(targetrate, nadapt, ε0bar, h0bar, γ, t0, κ, period, verbose)
  end
end

DualAveragingMCTuner(
  targetrate::Real,
  nadapt::Integer;
  ε0bar::Real=1.,
  h0bar::Real=0.,
  γ::Real=0.05,
  t0::Integer=10,
  κ::Real=0.75,
  period::Integer=100,
  verbose::Bool=false
) = DualAveragingMCTuner(targetrate, nadapt, ε0bar, h0bar, γ, t0, κ, period, verbose)

function tune!(tune::DualAveragingMCTune, tuner::DualAveragingMCTuner, count::Integer, a::Real)
  tune.hweight = 1/(count+tuner.t0)
  tune.hbar = (1-tune.hweight)*tune.hbar+tune.hweight*(tuner.targetrate-a)
  tune.step = exp(tune.μ-sqrt(count)*tune.hbar/tuner.γ)
  tune.εweight = count^(-tuner.κ)
  tune.εbar = exp((1-tune.εweight)*log(tune.εbar)+tune.εweight*log(tune.step))
end

Base.show(io::IO, tuner::DualAveragingMCTuner) =
  print(
    io,
    string(
      "DualAveragingMCTuner: target rate = ",
      tuner.targetrate,
      ", nadapt = ",
      tuner.nadapt,
      ", ε0bar = ",
      tuner.ε0bar,
      ", h0bar = ",
      tuner.h0bar,
      ", γ = ",
      tuner.γ,
      ", t0 = ",
      tuner.t0,
      ", κ = ",
      tuner.κ,
      ", period = ",
      tuner.period,
      ", verbose = ",
      tuner.verbose
    )
  )

Base.writemime(io::IO, ::MIME"text/plain", tuner::DualAveragingMCTuner) = show(io, tuner)


text:
module ValueHistories

using DataStructures
using RecipesBase
using Compat

export

    ValueHistory,
      UnivalueHistory,
        History,
        QHistory,
      MultivalueHistory,
        MVHistory,
        increment!,
    @trace

include("abstract.jl")
include("history.jl")
include("qhistory.jl")
include("mvhistory.jl")
include("recipes.jl")

end # module


text:
module Data

using YAML: load_file
using Images
using Images.FileIO
using LMDB
using LMDB: create, open, close, put!, start

struct MainDataset
    config::Dict
    num_frames::Int
    txn_map::Dict
    dbi_map::Dict
    idx_map::Dict
    yaw_map::Dict
    file_map::Dict

    function MainDataset(data_dir::AbstractString, config_path::AbstractString)
        # read the config file
        num_frames = 0
        config = load_file(config_path)
        txn_map = Dict()
        dbi_map = Dict()
        idx_map = Dict()
        yaw_map = Dict()
        file_map = Dict()

        for full_path in readdir(data_dir, join=true)
            if !isdir(full_path)
                continue
            end
            env = create()
            open(env, full_path)
            txn = start(env)

            dbi = open(txn)

            n = parse(Int64, LMDB.get(txn, dbi, "len", String)) - config["num_plan"] # last num_plan labels are missing from the dataset
            if n < config["num_plan"]+1
                print(full_path, " is too small. consider deleting")
                close(env, dbi)
            else
                offset = num_frames
                for i in 0:n-config["num_plan"]-1
                    num_frames+=1
                    for j in 0:length(config["camera_yaws"])-1
                        txn_map[(offset+i)*length(config["camera_yaws"])+j] = txn
                        dbi_map[(offset+i)*length(config["camera_yaws"])+j] = dbi
                        idx_map[(offset+i)*length(config["camera_yaws"])+j] = i
                        yaw_map[(offset+i)*length(config["camera_yaws"])+j] = j
                        file_map[(offset+i)*length(config["camera_yaws"])+j] = full_path
                    end
                end
            end
        end
        print(data_dir,": ", num_frames, " frames ", "x ", length(config["camera_yaws"]))
        new(config, num_frames, txn_map, dbi_map, idx_map, yaw_map, file_map)
    end
end

function access(tag::AbstractString, txn::Transaction, dbi::DBI, index::Int, T::Int, dtype::Type; preprocess = x->x)
    try
        reduce(hcat, preprocess.([LMDB.get(txn, dbi, "$(tag)_$(lpad(t, 5, '0'))", dtype) for t in index:index+(T-1)]))
    catch
        @show tag
    end
end

function get_item(d::MainDataset, state::Int)
    T = d.config["num_plan"]
    idx = state-1
    if !d.config["multi_cam"]
        idx *= length(d.config["camera_yaws"])
    end

    lmdb_txn = d.txn_map[idx]
    lmdb_dbi = d.dbi_map[idx]
    index = d.idx_map[idx]
    cam_index = d.yaw_map[idx]

    locs = access("loc", lmdb_txn, lmdb_dbi, index, T+1, Vector{Float32})
    rots = access("rot", lmdb_txn, lmdb_dbi, index, T, Vector{Float32})
    spds = access("spd", lmdb_txn, lmdb_dbi, index, T, Vector{Float32})

    decode = (x -> FileIO.load(IOBuffer(x)))
    lbls = [access("lbl_$(lpad(l, 2, '0'))", lmdb_txn, lmdb_dbi, index+1, T, Vector{UInt8}; preprocess=decode) |> channelview for l in 0:11]
    wide_rgb = permuteddimsview(access("wide_$(cam_index)", lmdb_txn, lmdb_dbi, index, 1, Vector{UInt8}; preprocess=decode) |> channelview, (2,3,1))
    wide_sem = access("wide_sem_$(cam_index)", lmdb_txn, lmdb_dbi, index, 1, Vector{UInt8}; preprocess=decode) |> channelview
    narr_rgb = access("narr_$(cam_index)", lmdb_txn, lmdb_dbi, index, 1, Vector{UInt8}; preprocess=decode) |> channelview

    cmd = access("cmd", lmdb_txn, lmdb_dbi, index, 1, Float32)
    wide_rgb, wide_sem, narr_rgb, lbls, locs, rots, spds, trunc.(Int, cmd)
end

Base.length(d::MainDataset) = (if (d.config["multi_cam"]) d.num_frames*length(d.config["camera_yaws"]) else d.num_frames end)
Base.iterate(d::MainDataset, state=1) = state > length(d) ? nothing : (get_item(d, state), state+1)

end


text:
export time_modeling

# Setup time-domain linear or nonlinear foward and adjoint modeling and interface to OPESCI/devito
function time_modeling(model_full::Model, srcGeometry, srcData, recGeometry, recData, dm, srcnum::Int64, op::Char, mode::Int64, options)
    pm = load_pymodel()

    # Load full geometry for out-of-core geometry containers
    typeof(recGeometry) == GeometryOOC && (recGeometry = Geometry(recGeometry))
    typeof(srcGeometry) == GeometryOOC && (srcGeometry = Geometry(srcGeometry))
    length(model_full.n) == 3 ? dims = [3,2,1] : dims = [2,1]   # model dimensions for Python are (z,y,x) and (z,x)

    # limit model to area with sources/receivers
    if options.limit_m == true
        model = deepcopy(model_full)
        if op=='J' && mode==1
            model,dm = limit_model_to_receiver_area(srcGeometry,recGeometry,model,options.buffer_size;pert=dm)
        else
            model = limit_model_to_receiver_area(srcGeometry,recGeometry,model,options.buffer_size)
        end
    else
        model = model_full
    end

    # Set up Python model structure
    if op=='J' && mode == 1
        modelPy = pm."Model"(origin=(0.,0.,0.), spacing=model.d, shape=model.n, vp=process_physical_parameter(sqrt.(1f0./model.m), dims), nbpml=model.nb,
            rho=process_physical_parameter(model.rho, dims), dm=process_physical_parameter(reshape(dm,model.n), dims), space_order=options.space_order)
    else
        modelPy = pm."Model"(origin=(0.,0.,0.), spacing=model.d, shape=model.n, vp=process_physical_parameter(sqrt.(1f0./model.m), dims), nbpml=model.nb,
            rho=process_physical_parameter(model.rho, dims), space_order=options.space_order)
    end

    # Load shot record if stored on disk
    if recData != nothing
        if typeof(recData[1]) == SegyIO.SeisCon
            recDataCell = Array{Any}(undef, 1); recDataCell[1] = convert(Array{Float32,2},recData[1][1].data); recData = recDataCell
        elseif typeof(recData[1]) == String
            recData = load(recData[1])."d".data
        end
    end

    # Remove receivers outside the modeling domain (otherwise leads to segmentation faults)
    if mode==1 && recGeometry != nothing
        recGeometry = remove_out_of_bounds_receivers(recGeometry, model)
    elseif mode==-1 && recGeometry != nothing
        recGeometry, recData = remove_out_of_bounds_receivers(recGeometry, recData, model)
    end

    # Devito interface
    argout = devito_interface(modelPy, model.o, srcGeometry, srcData, recGeometry, recData, dm, options)

    # Extend gradient back to original model size
    if op=='J' && mode==-1 && options.limit_m==true
        argout = vec(extend_gradient(model_full, model, reshape(argout, model.n)))
    end

    return argout
end

# Function instance without options
time_modeling(model::Model, srcGeometry, srcData, recGeometry, recData, perturbation, srcnum::Int64, op::Char, mode::Int64) =
    time_modeling(model, srcGeometry, srcData, recGeometry, recData, perturbation, srcnum, op, mode, Options())

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

# d_obs = Pr*F*Ps'*q
function devito_interface(modelPy::PyCall.PyObject, origin, srcGeometry::Geometry, srcData::Array, recGeometry::Geometry, recData::Nothing, dm::Nothing, options::Options)
    ac = load_acoustic_codegen()

    # Interpolate input data to computational grid
    dtComp = modelPy.critical_dt
    qIn = time_resample(srcData[1],srcGeometry,dtComp)[1]
    ntComp = size(qIn,1)
    ntRec = Int(trunc(recGeometry.t[1]/dtComp + 1))

    # Set up coordinates with devito dimensions
    src_coords = setup_grid(srcGeometry, modelPy.shape, origin)
    rec_coords = setup_grid(recGeometry, modelPy.shape, origin)

    # Devito call
    dOut = get(pycall(ac."forward_modeling", PyObject, modelPy, PyReverseDims(copy(transpose(src_coords))), PyReverseDims(copy(transpose(qIn))), PyReverseDims(copy(transpose(rec_coords))),
                  space_order=options.space_order, nb=modelPy.nbpml, free_surface=options.free_surface), 0)
    ntRec > ntComp && (dOut = [dOut zeros(size(dOut,1), ntRec - ntComp)])
    dOut = time_resample(dOut,dtComp,recGeometry)

    # Output shot record as judiVector
    if options.save_data_to_disk
        container = write_shot_record(srcGeometry,srcData,recGeometry,dOut,options)
        return judiVector(container)
    else
        return judiVector(recGeometry,dOut)
    end
end

# q_ad = Ps*F'*Pr'*d_obs
function devito_interface(modelPy::PyCall.PyObject, origin, srcGeometry::Geometry, srcData::Nothing, recGeometry::Geometry, recData::Array, dm::Nothing, options::Options)
    ac = load_acoustic_codegen()

    # Interpolate input data to computational grid
    dtComp = modelPy.critical_dt
    dIn = time_resample(recData[1],recGeometry,dtComp)[1]
    ntComp = size(dIn,1)
    ntSrc = Int(trunc(srcGeometry.t[1]/dtComp + 1))

    # Set up coordinates with devito dimensions
    src_coords = setup_grid(srcGeometry, modelPy.shape, origin)
    rec_coords = setup_grid(recGeometry, modelPy.shape, origin)

    # Devito call
    qOut = pycall(ac."adjoint_modeling", Array{Float32,2}, modelPy, PyReverseDims(copy(transpose(src_coords))), PyReverseDims(copy(transpose(rec_coords))), PyReverseDims(copy(transpose(dIn))),
                  space_order=options.space_order, nb=modelPy.nbpml, free_surface=options.free_surface)
    ntSrc > ntComp && (qOut = [qOut zeros(size(qOut), ntSrc - ntComp)])
    qOut = time_resample(qOut,dtComp,srcGeometry)

    # Output adjoint data as judiVector
    return judiVector(srcGeometry,qOut)
end

# u = F*Ps'*q
function devito_interface(modelPy::PyCall.PyObject, origin, srcGeometry::Geometry, srcData::Array, recGeometry::Nothing, recData::Nothing, dm::Nothing, options::Options)
    ac = load_acoustic_codegen()

    # Interpolate input data to computational grid
    dtComp = modelPy.critical_dt
    qIn = time_resample(srcData[1],srcGeometry,dtComp)[1]
    ntComp = size(qIn,1)

    # Set up coordinates with devito dimensions
    src_coords = setup_grid(srcGeometry, modelPy.shape, origin)

    # Devito call
    u = pycall(ac."forward_modeling", PyObject, modelPy, PyReverseDims(copy(transpose(src_coords))), PyReverseDims(copy(transpose(qIn))), nothing,
               space_order=options.space_order, nb=modelPy.nbpml, free_surface=options.free_surface)

    # Output forward wavefield as judiWavefield
    options.save_wavefield_to_disk && (u = dump_wavefield(u))
    return judiWavefield(Info(prod(modelPy.shape), 1, ntComp), dtComp, u)
end

# v = F'*Pr'*d_obs
function devito_interface(modelPy::PyCall.PyObject, origin, srcGeometry::Nothing, srcData::Nothing, recGeometry::Geometry, recData::Array, dm::Nothing, options::Options)
    ac = load_acoustic_codegen()

    # Interpolate input data to computational grid
    dtComp = modelPy.critical_dt
    dIn = time_resample(recData[1],recGeometry,dtComp)[1]
    ntComp = size(dIn,1)

    # Set up coordinates with devito dimensions
    rec_coords = setup_grid(recGeometry, modelPy.shape, origin)

    # Devito call
    v = pycall(ac."adjoint_modeling", PyObject, modelPy, nothing, PyReverseDims(copy(transpose(rec_coords))), PyReverseDims(copy(transpose(dIn))),
               space_order=options.space_order, nb=modelPy.nbpml, free_surface=options.free_surface)

    # Output adjoint wavefield as judiWavefield
    options.save_wavefield_to_disk && (v = dump_wavefield(v))
    return judiWavefield(Info(prod(modelPy.shape), 1, ntComp), dtComp, v)
end

# d_obs = Pr*F*u
function devito_interface(modelPy::PyCall.PyObject, origin, srcGeometry::Nothing, srcData::Array, recGeometry::Geometry, recData::Nothing, dm::Nothing, options::Options)
    ac = load_acoustic_codegen()

    # Interpolate input data to computational grid
    dtComp = modelPy.critical_dt
    ntRec = Int(trunc(recGeometry.t[1]/dtComp + 1))

    # Set up coordinates with devito dimensions
    rec_coords = setup_grid(recGeometry, modelPy.shape, origin)

    # Devito call
    dOut = get(pycall(ac."forward_modeling", PyObject, modelPy, nothing, srcData[1], PyReverseDims(copy(transpose(rec_coords))),
                  space_order=options.space_order, nb=modelPy.nbpml, free_surface=options.free_surface), 0)
    #ntRec > ntComp && (dOut = [dOut zeros(size(dOut,1), ntRec - ntComp)])
    dOut = time_resample(dOut,dtComp,recGeometry)

    # Output shot record as judiVector
    if options.save_data_to_disk
        throw("Writing shot record to SEG-Y file not supported for modeling with wavefield as right-hand-side.")
    else
        return judiVector(recGeometry,dOut)
    end
end

# q_ad = Ps*F'*v
function devito_interface(modelPy::PyCall.PyObject, origin, srcGeometry::Geometry, srcData::Nothing, recGeometry::Nothing, recData::Array, dm::Nothing, options::Options)
    ac = load_acoustic_codegen()

    # Interpolate input data to computational grid
    dtComp = modelPy.critical_dt
    ntSrc = Int(trunc(srcGeometry.t[1]/dtComp + 1))

    # Set up coordinates with devito dimensions
    src_coords = setup_grid(srcGeometry, modelPy.shape, origin)

    # Devito call
    qOut = pycall(ac."adjoint_modeling", Array{Float32,2}, modelPy, PyReverseDims(copy(transpose(src_coords))), nothing, recData[1],
                  space_order=options.space_order, nb=modelPy.nbpml, free_surface=options.free_surface)
    #ntSrc > ntComp && (qOut = [qOut zeros(size(qOut), ntSrc - ntComp)])
    qOut = time_resample(qOut,dtComp,srcGeometry)

    # Output adjoint data as judiVector
    return judiVector(srcGeometry,qOut)
end

# u_out = F*u_in
function devito_interface(modelPy::PyCall.PyObject, origin, srcGeometry::Nothing, srcData::Array, recGeometry::Nothing, recData::Nothing, dm::Nothing, options::Options)
    ac = load_acoustic_codegen()

    # Interpolate input data to computational grid
    dtComp = modelPy.critical_dt
    ntComp = srcData[1].shape[1]

    # Devito call
    u = pycall(ac."forward_modeling", PyObject, modelPy, nothing, srcData[1], nothing,
                  space_order=options.space_order, nb=modelPy.nbpml, free_surface=options.free_surface)

    # Output forward wavefield as judiWavefield
    options.save_wavefield_to_disk && (u = dump_wavefield(u))
    return judiWavefield(Info(prod(modelPy.shape), 1, ntComp), dtComp, u)
end

# v_out = F'*v_in
function devito_interface(modelPy::PyCall.PyObject, origin, srcGeometry::Nothing, srcData::Nothing, recGeometry::Nothing, recData::Array, dm::Nothing, options::Options)
    ac = load_acoustic_codegen()

    # Interpolate input data to computational grid
    dtComp = modelPy.critical_dt
    ntComp = recData[1].shape[1]

    # Devito call
    v = pycall(ac."adjoint_modeling", PyObject, modelPy, nothing, nothing, recData[1],
                  space_order=options.space_order, nb=modelPy.nbpml, free_surface=options.free_surface)

    # Output adjoint wavefield as judiWavefield
    options.save_wavefield_to_disk && (v = dump_wavefield(v))
    return judiWavefield(Info(prod(modelPy.shape), 1, ntComp), dtComp, v)
end

# d_lin = J*dm
function devito_interface(modelPy::PyCall.PyObject, origin, srcGeometry::Geometry, srcData::Array, recGeometry::Geometry, recData::Nothing, dm::Array, options::Options)
    ac = load_acoustic_codegen()

    # Interpolate input data to computational grid
    dtComp = modelPy.critical_dt
    tmaxRec = recGeometry.t[1]
    qIn = time_resample(srcData[1],srcGeometry,dtComp)[1]
    ntComp = size(qIn,1)
    ntRec = Int(trunc(tmaxRec/dtComp + 1))

    # Set up coordinates with devito dimensions
    #origin = get_origin(modelPy)
    src_coords = setup_grid(srcGeometry, modelPy.shape, origin)
    rec_coords = setup_grid(recGeometry, modelPy.shape, origin)

    # Devito call
    dOut = pycall(ac."forward_born", Array{Float32,2}, modelPy, PyReverseDims(copy(transpose(src_coords))), PyReverseDims(copy(transpose(qIn))), PyReverseDims(copy(transpose(rec_coords))),
                  space_order=options.space_order, nb=modelPy.nbpml, isic=options.isic)
    ntRec > ntComp && (dOut = [dOut zeros(size(dOut,1), ntRec - ntComp)])
    dOut = time_resample(dOut,dtComp,recGeometry)

    # Output linearized shot records as judiVector
    if options.save_data_to_disk
        container = write_shot_record(srcGeometry,srcData,recGeometry,dOut,options)
        return judiVector(container)
    else
        return judiVector(recGeometry,dOut)
    end
end

# dm = J'*d_lin
function devito_interface(modelPy::PyCall.PyObject, origin, srcGeometry::Geometry, srcData::Array, recGeometry::Geometry, recData::Array, dm::Nothing, options::Options)
    ac = load_acoustic_codegen()

    # Interpolate input data to computational grid
    dtComp = modelPy.critical_dt
    qIn = time_resample(srcData[1],srcGeometry,dtComp)[1]
    dIn = time_resample(recData[1],recGeometry,dtComp)[1]

    # Set up coordinates with devito dimensions
    src_coords = setup_grid(srcGeometry, modelPy.shape, origin)
    rec_coords = setup_grid(recGeometry, modelPy.shape, origin)

    if options.optimal_checkpointing == true
        op_F = pycall(ac."forward_modeling", PyObject, modelPy, PyReverseDims(copy(transpose(src_coords))), PyReverseDims(copy(transpose(qIn))), PyReverseDims(copy(transpose(rec_coords))),
                      op_return=true, space_order=options.space_order, nb=modelPy.nbpml)
        grad = pycall(ac."adjoint_born", Array{Float32, length(modelPy.shape)}, modelPy, PyReverseDims(copy(transpose(rec_coords))), PyReverseDims(copy(transpose(dIn))), op_forward=op_F, space_order=options.space_order,
            nb=modelPy.nbpml, is_residual=true, isic=options.isic, n_checkpoints=options.num_checkpoints, maxmem=options.checkpoints_maxmem)
    elseif ~isempty(options.frequencies)    # gradient in frequency domain
        typeof(options.frequencies) == Array{Any,1} && (options.frequencies = options.frequencies[1])
        d_pred, uf_real, uf_imag = pycall(ac."forward_freq_modeling", PyObject, modelPy, PyReverseDims(copy(transpose(src_coords))), PyReverseDims(copy(transpose(qIn))), PyReverseDims(copy(transpose(rec_coords))),
                                          options.frequencies, space_order=options.space_order, nb=modelPy.nbpml, factor=options.dft_subsampling_factor)
        grad = pycall(ac."adjoint_freq_born", Array{Float32, length(modelPy.shape)}, modelPy, PyReverseDims(copy(transpose(rec_coords))), PyReverseDims(copy(transpose(dIn))),
                      options.frequencies, uf_real, uf_imag, space_order=options.space_order, nb=modelPy.nbpml, isic=options.isic, factor=options.dft_subsampling_factor)
    else
        u0 = get(pycall(ac."forward_modeling", PyObject, modelPy, PyReverseDims(copy(transpose(src_coords))), PyReverseDims(copy(transpose(qIn))), PyReverseDims(copy(transpose(rec_coords))),
                    space_order=options.space_order, nb=modelPy.nbpml, save=true, tsub_factor=options.subsampling_factor), 1)
        grad = pycall(ac."adjoint_born", Array{Float32, length(modelPy.shape)}, modelPy, PyReverseDims(copy(transpose(rec_coords))), PyReverseDims(copy(transpose(dIn))), u=u0,
                      space_order=options.space_order, tsub_factor=options.subsampling_factor, nb=modelPy.nbpml, isic=options.isic)
    end

    # Remove PML and return gradient as Array
    grad = remove_padding(grad,modelPy.nbpml, true_adjoint=options.sum_padding)
    return vec(grad)
end


text:
# functions that might need to be loaded multiple times (dont add docs here)

function plotFactorMeasurements(dfg::AbstractDFG,
                                fctsym::Symbol,
                                fct::ExtendedPose2Pose2Types;
                                hdl=[],
                                dist::Vector{Float64}=[0.0;]  )
  #
  asMeas, asPred = solveFactorMeasurements(dfg, fctsym)

  plotFactorValues(asMeas, asPred, Pose2Pose2, fctsym=fctsym,hdl=hdl,dist=dist)
end

function plotFactor(dfg::AbstractDFG,
                    fctsym::Symbol,
                    fct::ExtendedPose2Pose2Types;
                    hdl=[],
                    dist::Vector{Float64}=Float64[0.0;]  )

  # variables
  fct = getFactor(dfg, fctsym)
  vars = fct._variableOrderSymbols

  pv1 = plotPose(dfg, vars[1], hdl=hdl)
  pv2 = plotPose(dfg, vars[2], hdl=hdl)

  pv12 = plotFactorMeasurements(dfg, fctsym, hdl=hdl, dist=dist)

  vstack(pv1, pv2, pv12)
end

function reportFactors(dfg::AbstractDFG,
                       T::PlotTypesPose2,
                       fcts::Vector{Symbol}=ls(dfg, T);
                       prefix::AbstractString="",
                       filepath=joinpath(getSolverParams(dfg).logpath, "$prefix"*getTimeEasy()*"_$T.pdf"),
                       show::Bool=true,
                       pdfWidth=20cm,
                       pdfHeight=30cm)
  #
  ss = split(filepath, '/')
  path = joinpath("/", joinpath(ss[1:(end-1)]...), "tmp")
  mkpath(path)
  alldists= Vector{Float64}()

  files = String[]
  ndist = Float64[0.0;]
  for fc in fcts
    if isMultihypo(getFactor(dfg, fc))
      # skip this factor
      continue
    end
    file = joinpath(path,"$fc.pdf")
    ndist[1] = 0.0
    plotFactor(dfg, fc, dist=ndist) |> PDF(file, pdfWidth, pdfHeight)
    push!(files, file)
    push!(alldists, ndist[1])
  end
  fileord = sortperm(alldists, rev=true)
  files = files[fileord]
  push!(files, filepath)

  2 < length(files) ? run(`pdfunite $files`) : nothing
  !show ? nothing : (@async run(`evince $filepath`))
  return filepath
end

#


text:
#--------------------------------------------------------------------
# Spacetime Discretization methods in Julia
# Soham 02-2018
#--------------------------------------------------------------------
"TODO: Fix convergence for n > 243
       Test convergence for analytic potential
       Test self-convergence with arbitrary potentials"

function testpconvforprojection(fn::Function, Nx::Int, Ny::Int)::Float64
    rPatch = Patch([1,1], projectonPatchbyRestriction(fn, Nx, Ny, 1, [1,1]))
    sPatch  = Float64[fn(x,y) for x in chebgrid(4Nx), y in chebgrid(4Ny)]
    return L2norm(sPatch, interpolatePatch(rPatch, 4Nx, 4Ny).value, chebweights(4Nx), chebweights(4Ny))
end

function testhconvforprojection(fn::Function, Nx::Int, Ny::Int, M::Int)::Float64
    L2err  = 0.0
    for m in 1:M, n in 1:M
        rPatch = Patch([1,1], projectonPatchbyRestriction(fn, Nx, Ny, M, [m,n]))
        sPatch = Float64[fn(x,y) for x in chebgrid(4Nx, M, m), y in chebgrid(4Ny, M, n)]
        L2err += (L2norm(sPatch, interpolatePatch(rPatch, 4Nx, 4Ny).value, chebweights(4Nx)/M, chebweights(4Ny)/M))^2
    end
    return sqrt(L2err)
end

# test p-refinement
showconv(x->sin(pi*x), y->sin(pi*y),
        (x,y)-> 0,
        (x,y)-> sin(pi*x) + sin(pi*y),
        20,  # maxmodes
        1,   # maxlevels
        1)   # h-factor

# test h-refinement
showconv(x->sin(pi*x), y->sin(pi*y),
        (x,y)-> 0,
        (x,y)-> sin(pi*x) + sin(pi*y),
        3,  # maxmodes
        4,  # maxlevels
        3)  # h-factor


text:
#=
The file contains methods which can override file operations to ensure that
the files are closed properly. This needs to be included in the test case with
DEBUG=true. default is false.

This is kept external to the regular library to ensure normal library operations
are not affected.
=#

const DEBUG=false

@static if DEBUG

import PDFIO.Common: get_tempfilepath,
                     util_open, util_close, utilPrintOpenFiles

import Base: close

IODebug=[0, Vector{Tuple{AbstractString, IO}}()]

function get_tempfilepath()
    global IODebug
    IODebug[1] += 1
    path = joinpath(get_tempdir(), string(IODebug[1]))
    return (path, util_open(path, "w"))
end

function util_open(filename, mode)
    global IODebug
    io=open(filename, mode)
    @printf "Opening file: %s\n" filename
    push!(IODebug[2], (filename, io))
    return io
end

function util_close(handle::IOStream)
    global IODebug
    idx=1
    for file in IODebug[2]
        if (handle === file[2])
            @printf "Closing file: %s, %d\n" file[1] idx
            close(handle)
            deleteat!(IODebug[2],idx)
            return
        end
        idx+=1
    end
    error("IO handle not found")
end

util_close(handle::IOBuffer) = close(handle)

function utilPrintOpenFiles()
    global IODebug
    println("The following files are opened currently:")
    println("-----------------------------------------")
    for f in IODebug[2]
        println(f)
    end
    println("-----------------------------------------")
    return IODebug[2]
end

end

function files_equal(f1, f2)
    io1 = util_open(f1, "r"); io2 = util_open(f2, "r")
    buf1 = read(io1); buf2 = read(io2)
    util_close(io1); util_close(io2)
    return buf1 == buf2
end

function extract_text(io, doc)
    npage = pdDocGetPageCount(doc)
    for i=1:npage
        page = pdDocGetPage(doc, i)
        if pdPageIsEmpty(page) == false
            pdPageGetContentObjects(page)
            pdPageExtractText(io, page)
            println(io)
        end
    end
end

function testfiles(filename)
    name, ext = splitext(filename)
    return (name*".res", joinpath(@__DIR__, "files", name*".txt"))
end

function local_testfiles(filename)
    name, ext = splitext(filename)
    return (name*".res", joinpath(@__DIR__, "files", name*".txt"),
            joinpath(@__DIR__, "files", filename))
end


text:
import HorizonSideRobots
using HorizonSideRobots
a = [Nord, Ost, Sud, West]
nx = [0, 1, 0, -1]
ny = [1, 0, -1, 0]
x = 25
y = 25
spawn = 25
used = [Int[0 for i in 1:50] for i in 1:50]
max_x = 25
min_x = 25
max_y = 25
min_y = 25

function clear()
    for i in 1:50
        for j in 1:50
            used[i][j] = 0
        end
    end
end

function inverse(side)
    return HorizonSide((Int(side) + 2) % 4)
end

function try_move(r, side)
    if (!isborder(r, side))
        move!(r, side)
    end
end

function moves(side, cnt)
    for i in 1:cnt
        try_move(r, side)
        if (!isborder(r, Nord))
            return false
        end
    end
    return true
end

function main(r)
    cnt = 1
    flag = true
    side = West
    while (flag)
        print(1)
        flag = moves(side, cnt)
        side = inverse(side)
        cnt *= 2
    end
end

r = Robot(animate = false, "mirea-progs/Task8/temp.sit")
main(r)
show(r)

#=
ДАНО: Робот - в произвольной клетке ограниченного прямоугольного поля (без внутренних перегородок)
РЕЗУЛЬТАТ: Робот - в исходном положении, в клетке с роботом стоит маркер, и все остальные клетки поля промаркированы в шахматном порядке
=#

max_stars_repo_name: Sergej-Kiryakov/mirea-progs

max_issues_repo_name: Sergej-Kiryakov/mirea-progs

max_forks_repo_name: Sergej-Kiryakov/mirea-progs


text:
struct PerfectMagneticConductor <: AbstractBoundaryCondition end
const PMC = PerfectMagneticConductor()


text:
"""
Contains functions that take in time-invariant network structure parameters,
such as the coupling matrix J and gain parameter g, and run simulations of
network activity. Includes functions that numerically calculate network
properties (such as Lyapunov spectrum) in parallel with the simulation.
"""

using LinearAlgebra

#-------------------------------------------------------------------------------
"""
	sim_activity(J::Matrix{T}, g::T, h0::Vector{T}; dt=0.1, t_bound=1000.0,
		s=0.0, savestep=1, minimum_norm=0) where T <: AbstractFloat

Integrates a tanh rate network with coupling matrix `J` and strength scaling `g`
from initial condition `h0`.

When `minimum_norm > 0` the integration will be cut off when the network state
falls too close to the origin to prevent floating point underflow.
"""
function sim_activity(J::Matrix{T}, g::T, h0::Vector{T}; dt=0.1, t_bound=1000.0,
		s=0.0, savestep=1, minimum_norm=0) where T <: AbstractFloat

	N = size(J)[1]
	# Check network size and type consistency
	@assert N == size(J)[2] == length(h0)
	@assert T == typeof(dt) == typeof(t_bound) == typeof(s)

	maxstep = round(Int, t_bound/dt) # Discretized simulation length

	# Allocate memory for output
	out = zeros((N, Int(maxstep/savestep)))

	h = copy(h0) # Set initial condition
	r = Vector{T}(undef, N)

	J = g .* J + s * I

	step = 0
	while step < maxstep
		step += 1
		map!(tanh, r, h)
		BLAS.gemv!('N', dt, J, r, 1 - dt, h)
		if step % savestep == 0 # Save every specified time interval
			savept = step ÷ savestep
			out[:,savept] .= h
		end
		if norm(h) < minimum_norm
			D = -I + J .* ((sech.(h)).^2)'
			if all(real(eigvals(D)) .< 0)
				break
			end
		end
	end

	lastindex = step ÷ savestep
	return out[:,1:lastindex]
end

"""
	sim_lyapunov(J::Matrix{T}, g::T, h0::Vector{T}, dh::Vector{T},
		Q0::Matrix{T}; s=0.0, dt=0.1, t_bound=1000.0, t_ONS=1.0, wait=0.0,
		return_max_LE=false, printstuff=false) where T <: AbstractFloat

Integrates a tanh network with a given weight matrix and gain parameter from a
specified initial condition, and returns a Lyapunov spectrum calculated along
the trajectory.

Optionally returns the maximum Lyapunov exponent calculated using a direct
numerical simulation of a perturbed trajectory.
"""
function sim_lyapunov(J::Matrix{T}, g::T, h0::Vector{T}, dh::Vector{T},
		Q0::Matrix{T}; s=0.0, dt=0.1, t_bound=1000.0, t_ONS=1.0, wait=0.0,
		return_max_LE=false, printstuff=false) where T <: AbstractFloat

	N = size(J)[1]
	# Check dimensional consistency
	@assert size(J) == (length(h0), length(dh)) == size(Q0) == (N,N)
	@assert T == typeof(s) == typeof(dt) == typeof(t_bound) == typeof(t_ONS) == typeof(wait)

	# Allocate memory for output
	spectrum = zeros(T, N)

	maxstep = Int(t_bound / dt) # When to stop the simulation
	startstep = Int(wait / dt) # When to start tracking Lypunov exponents
	tONSstep = Int(t_ONS / dt) # Interval between re-orthonormalization
	try
		Int(t_bound/t_ONS)
		Int(wait/t_ONS)
	catch e
		if isa(e, InexactError)
			msg = "t_ONS must divide t_bound and wait"
			throw(ArgumentError(msg))
		end
	end

	# Allocate vector for activity and firing rate of current state
	h = copy(h0)
	r = Vector{T}(undef, N)

	# Allocate memory for Jacobian calculation and orthonormal basis
	sech2(x::T)::T = (sech(x))^2
	sech2h = Matrix{T}(undef, 1, N)
	D = Matrix{T}(undef, N, N)
	Q = copy(Q0)
	dQ = Matrix{T}(undef, N, N)

	if return_max_LE
		max_LE = 0.0
		p_factor = T(1e-6)  # Order of the perturbation handled by the sim
		dh = copy(dh) ./ norm(dh)
		h2 = h .+ p_factor .* dh
		r2 = Vector{T}(undef, N)
	end

	# The effective coupling matrix J
	J = g .* J + s * I

	step = 0
	while step < maxstep
		step += 1

		#Evolve the trajectory and the orthonormal basis in the tangent space
		map!(tanh, r, h)
		BLAS.gemv!('N', dt, J, r, 1 - dt, h)
		map!(sech2, sech2h, h)
		copyto!(D, J)
		D .*= sech2h
		D[diagind(D)] .-= 1
		BLAS.gemm!('N', 'N', dt, D, Q, T(0), dQ)
		Q .+= dQ

		#Evolve the perturbed trajectory and calculate how much it has grown
		if return_max_LE
			map!(tanh, r2, h2)
			BLAS.gemv!('N', dt, J, r2, 1 - dt, h2)
			dh .= (h2 .- h) ./ p_factor
			scale = norm(dh)

			# If the perturbation scale has grown or shrunk too much, then
			# reset it to 1 and track the cumulative sum of log scales
			if scale > 100 || scale < 0.01 || step == startstep
				h2 .= h .+ dh ./ scale .* p_factor
				if step > startstep
					max_LE += log(scale)
				end
			end
		end

		# At each re-orthonormalization step, the log diagonal entries of R are
		# saved in a cumulative sum
		if step % tONSstep == 0
			printstuff && println("Time step: ", step)
			copyto!(dQ, Q)
			F = qr!(dQ)
			diags = view(F.factors, diagind(F.factors))
			if step > startstep  # forcing positive R-diagonal elements
				spectrum .+= log.(abs.(diags))
			end
			fill!(Q, 0)
			map!(sign, view(Q, diagind(Q)), diags)
			lmul!(F.Q, Q)
		end
	end

	# Update the maximum Lyapunov exponent after the simulation has ended so
	# that the full simulation time is used for the calculation
	if return_max_LE
		max_LE += log(norm((h2 .- h) ./ p_factor))
		return (spectrum ./ (t_bound - wait), max_LE / (t_bound - wait))
	else
		return spectrum ./ (t_bound - wait)
	end
end

"""
	sim_direct_max_LE(J::Matrix{Float64}, g::Float64, h0::Vector{Float64},
		dh::Vector{Float64}; dt=0.1, t_bound=1000.0, wait=0.0)

Calculates the maximum Lyapunov exponent of a tanh rate network given by
coupling matrix J, dependent on a given initial condition and initial
perturbation, using a direct method that tracks the growth of the perturbed
trajectory relative to the unperturbed one.
"""
function sim_direct_max_LE(J::Matrix{T}, g::T, h0::Vector{T}, dh::Vector{T};
		s=0.0, dt=0.1, t_bound=1000.0, wait=0.0, printstuff=false) where T <: AbstractFloat

	N = size(J)[1]
	# Check dimensional consistency
	@assert size(J) == (length(h0), length(dh)) == (N,N)
	p_factor = 10^-10  # Order of the perturbation handled by the sim

	# Allocate memory for output
	max_LE = 0.0

	maxstep = Int(t_bound/dt) # When to stop the simulation
	startstep = Int(wait/dt) # When to start tracking Lypunov exponents

	# Set the initial condition and perturbed IC
	h = copy(h0)
	dh = copy(dh) ./ norm(dh)
	h2 = h .+ p_factor .* dh

	r = Vector{Float64}(undef, N)
	r2 = Vector{Float64}(undef, N)

	step = 0
	while step < maxstep
		step += 1

		# Evolve the trajectory and its perturbation
		map!(tanh, r, h)
		map!(tanh, r2, h2)
		BLAS.gemv!('N', dt*g, J, r, 1.0-dt, h)
		BLAS.gemv!('N', dt*g, J, r2, 1.0-dt, h2)
		copyto!(dh, h2)
		dh .-= h
		dh ./= p_factor
		scale = norm(dh)

		# If the perturbation scale has grown or shrunk too much, then
		# reset it to 1 and track the cumulative sum of log scales
		if scale > 10 || scale < 0.1 || step == startstep
			printstuff && println("Time step: ", step)
			copyto!(h2, h)
			dh ./= (scale / p_factor)
			h2 .+= dh
			if step > startstep
				max_LE += log(scale)
			end
		end
	end

	# Update the maximum Lyapunov exponent after the simulation has ended so
	# that the full simulation time is used for the calculation
	max_LE += log(norm((h2 .- h) ./ p_factor))
	return max_LE / (t_bound - wait)
end

"""
    sim_lyapunov(J::Matrix{Float64}, g::Float64, A::Matrix{Float64},
        onb::Matrix{Float64}; dt=0.1, t_bound=1000.0, t_ONS=1.0, wait=0.0,
		printstuff=false)

Calculates the Lyapunov spectrum of a tanh rate network with coupling matrix
`J` and gain parameter `g` by integrating an orthonormal system according to
the Jacobians in the tangent space along a trajectory given by activity
matrix `A`.
"""
function sim_lyapunov(J::Matrix{Float64}, g::Float64, A::Matrix{Float64},
		onb::Matrix{Float64}; dt=0.1, t_bound=1000.0, t_ONS=1.0, wait=0.0,
		printstuff=false)

	N = size(J)[1]
	# Check dimensional consistency
	@assert(size(J)==(size(A,1),N)==size(onb)==(N,N))

	# Ranges that index the diagonal entries of useful matrix sizes
	diagonals = diagind(J)

	# Allocate memory for output
	spectrum = zeros(N)

	maxstep = Int(t_bound/dt) # When to stop the simulation
	startstep = Int(wait/dt) # When to start tracking Lypunov exponents
	tONSstep = Int(t_ONS/dt) # Interval between re-orthonormalization
	@assert maxstep==size(A,2) "Simulation length must match time parameters"
	try
		Int(t_bound/t_ONS)
		Int(wait/t_ONS)
	catch e
		if isa(e, InexactError)
			msg = "t_ONS must divide t_bound and wait"
			throw(ArgumentError(msg))
		end
	end

	# Allocate memory for Jacobian calculation and orthonormal basis
	sech2(x::Float64)::Float64 = (sech(x))^2
	sech2h = Matrix{Float64}(undef, 1, N)
	D = Matrix{Float64}(undef, N, N)
	Q = copy(onb)
	dQ = Matrix{Float64}(undef, N, N)

	J = g .* J

	step = 0
	while step < maxstep
		step += 1

		#Evolve the trajectory and the orthonormal basis in the tangent space
		map!(sech2, sech2h, view(A, :, step))
		copyto!(D, J)
		D .*= sech2h
		D[diagonals] .-= 1.0
		BLAS.gemm!('N', 'N', dt, D, Q, 0.0, dQ)
		Q .+= dQ

		# At each re-orthonormalization step, the log diagonal entries of R are
		# saved in a cumulative sum
		if step % tONSstep == 0
			printstuff && println("Time step: ", step)
			copyto!(dQ, Q)
			F = qr!(dQ)
			diags = view(F.factors, diagonals)
			if step > startstep  # forcing positive R-diagonal elements
				spectrum .+= log.(abs.(diags))
			end
			fill!(Q, 0.0)
			map!(sign, view(Q, diagonals), diags)
			lmul!(F.Q, Q)
		end
	end

	# Update the maximum Lyapunov exponent after the simulation has ended so
	# that the full simulation time is used for the calculation
	return spectrum ./ (t_bound - wait)
end

"""
	sim_lyapunov_halfspec(J::Matrix{Float64}, g::Float64,
		A::Matrix{Float64}, onb::Matrix{Float64}; dt=0.1, t_bound=1000.0,
		t_ONS=1.0, wait=0.0, printstuff=false)

Calculates the Lyapunov spectrum of a tanh rate network with coupling matrix
`J` and gain parameter `g` by integrating an orthonormal system according to
the Jacobians in the tangent space along a trajectory given by activity
matrix `A`. Saves time by only calculating half of the Lyapunov spectrum
"""
function sim_lyapunov_halfspec(J::Matrix{Float64}, g::Float64,
		A::Matrix{Float64}, onb::Matrix{Float64}; dt=0.1, t_bound=1000.0,
		t_ONS=1.0, wait=0.0, printstuff=false)

	N = size(J)[1]
	@assert N%2==0 "Calculating half the spectrum only works for even N"
	Nh = Int(N/2)
	# Check dimensional consistency
	@assert(size(J)==(size(A,1),N)==size(onb)==(N,N))

	# Ranges that index the diagonal entries of useful matrix sizes
	diagonals = diagind(J)
	diags_1halfrng = range(1, step=N+1, length=Nh)
	diags_2halfrng = range(1 + Nh * (N + 1), step=N+1, length=Nh)

	# Allocate memory for output
	spectrum = zeros(Nh)

	maxstep = Int(t_bound/dt) # When to stop the simulation
	startstep = Int(wait/dt) # When to start tracking Lypunov exponents
	tONSstep = Int(t_ONS/dt) # Interval between re-orthonormalization
	@assert maxstep==size(A,2) "Simulation length must match time parameters"
	try
		Int(t_bound/t_ONS)
		Int(wait/t_ONS)
	catch e
		if isa(e, InexactError)
			msg = "t_ONS must divide t_bound and wait"
			throw(ArgumentError(msg))
		end
	end

	# Allocate memory for Jacobian calculation and orthonormal basis
	sech2(x::Float64)::Float64 = (sech(x))^2
	sech2h = Matrix{Float64}(undef, 1, N)
	D = Matrix{Float64}(undef, N, N)
	Q = copy(onb)
	dQ = Matrix{Float64}(undef, N, N)
	Qh = Matrix{Float64}(undef, N, Nh)

	J = g .* J

	step = 0
	while step < maxstep
		step += 1

		#Evolve the orthonormal basis in the tangent space
		map!(sech2, sech2h, view(A, :, step))
		copyto!(D, J)
		D .*= sech2h
		D[diagonals] .-= 1.0
		BLAS.gemm!('N', 'N', dt, D, Q, 0.0, dQ)
		Q .+= dQ

		# At each re-orthonormalization step, the log diagonal entries of R are
		# saved in a cumulative sum
		if step % tONSstep == 0
			printstuff && println("Time step: ", step)
			copyto!(Qh, view(Q, :, 1:Nh))
			F = qr!(Qh)
			diags = view(F.factors, diags_1halfrng)
			if step > startstep  # forcing positive R-diagonal elements
				spectrum .+= log.(abs.(diags))
			end
			fill!(Q, 0.0)
			map!(sign, view(Q, diags_1halfrng), diags)
			Q[diags_2halfrng] .= 1.0
			lmul!(F.Q, Q)
		end
	end

	# Update the maximum Lyapunov exponent after the simulation has ended so
	# that the full simulation time is used for the calculation
	return spectrum ./ (t_bound - wait)
end


text:
"""
LinSolvers (C) 2020, Martin Gurtner
Solvers of linear systems of equations tailored to the needs of
DistCtrl4DistMan.jl package
"""
module LinSolvers

using LinearAlgebra

export kernel_solve, kernel_solve!, LinSystemData

struct LinSystemData{T<:Real, I<:Int}
    n::I;
    m::I;
    A::Matrix{T};
    b::Vector{T};
    R::Matrix{T};
    x_sol::Vector{T};

    LinSystemData{T}(n::I, m::I) where {I<:Int, T<:Real} = new{T, I}(
            n, m,
            zeros(T, m+n,m),                   # A
            zeros(T, m),                       # b
            zeros(T, m,m),    # R
            zeros(T, n)
    );

    LinSystemData(n::I, m::I, x_sol::Vector{T}) where {I<:Int, T<:Real} = new{T, I}(
        n, m,
        zeros(T, m+n,m),  # A
        zeros(T, m),      # b
        zeros(T, m,m),    # R
        x_sol
    );
end # struct

"""
    update!(lsd::LinSystemData, Dᵀ::Matrix, λ<:Real, b::Vector)

Update 'lsd' structure storing data for solving matrix equation (DᵀD + 1/λI)x=Dᵀb.
"""
function update!(lsd::LinSystemData{T}, Dᵀ::Matrix{T}, λ::T, b::Vector{T}) where T<:Real
    @assert (size(Dᵀ, 1) == lsd.n) && (size(Dᵀ, 2) == lsd.m) "LinSystemData does not have compatible dimensions with matrix Dᵀ."
    @assert length(b) == lsd.m "LinSystemData does not have compatible dimensions with vector b."

    @inbounds for i in 1:lsd.n, j in 1:lsd.m
        lsd.A[i,j] = Dᵀ[i,j];
    end

    diag_val = 1/sqrt(λ);
    @inbounds for i in 1:lsd.m, j in 1:lsd.m
        lsd.A[i+lsd.n,j] = i==j ? diag_val : 0;
    end

    copyto!(lsd.b, b);

    nothing
end

"""
    kernel_solve(Dᵀ::Matrix, λ<:Real, b::Vector)

Solve a matrix equation in form (DᵀD + 1/λ)x=Dᵀb.

# Examples
```julia-repl
julia> n, m = 10, 3;
julia> for i in 1:10
           x_sol = kernel_solve(rand(n,m), 1.0, rand(m));
           print("Solution is \$x_sol\\n");
       end
Solution is ...
```
"""
function kernel_solve(Dt::Matrix{T}, λ::T, b::Vector{T}) where T<:Real
    n, m = size(Dt);
    A = zeros(m+n, m);
    b_tmp = zeros(m);
    x_sol = zeros(n);
    R = zeros(m, m);

    A[1:n,1:m] = Dt;
    @simd for i in (n+1):(n+m)
        A[i,i-n] = 1/sqrt(λ);
    end

    qr_r!(A, R);

    copyto!(b_tmp, b);
    forward_substitution!(R, b_tmp)
    backward_substitution!(R, b_tmp)

    mul!(x_sol, Dt, b_tmp);
    return x_sol;
end

"""
    kernel_solve!(lsd::LinSystemData, Dᵀ::Matrix, λ<:Real, b::Vector)

Solve a matrix equation in form (DᵀD + 1/λ)x=Dᵀb using the pre-allocated data structure 'lsd'.

# Example
```julia-repl
julia> n, m = 10, 3;
julia> lsd = LinSystemData{Float64}(n, m);
julia> for i in 1:10
           kernel_solve!(lsd, rand(n,m), 1.0, rand(m));
           print("Solution is \$(lsd.x_sol)\\n");
       end
Solution is ...
```
"""
function kernel_solve!(lsd::LinSystemData{T}, Dt::Matrix{T}, λ::T, b::Vector{T}) where T<:Real
    update!(lsd, Dt, λ, b);

    qr_r!(lsd.A, lsd.R);

    forward_substitution!(lsd.R, lsd.b);
    backward_substitution!(lsd.R, lsd.b);

    mul!(lsd.x_sol, Dt, lsd.b);

    nothing
end

"""
    kernel_solve_prealloc!(Dt::Matrix{T}, λ::T, Fi::Vector{T}, x_sol::Vector{T},
            A::Matrix{T},
            b_tmp::Vector{T},
            R::Matrix{T}) where T<:Real

Solve a matrix equation in form (DᵀD + 1/λ)x=Dᵀb with preallocated vectors and matrices provided.

# Examples
```julia-repl
julia> n, m = 10, 3;
julia> for i in 1:10
           x_sol = kernel_solve(rand(n,m), 1.0, rand(m));
           print("Solution is \$x_sol\\n");
       end
Solution is ...
```
"""
function kernel_solve_prealloc!(Dt::Matrix{T}, λ::T, b::Vector{T}, x_sol::Vector{T},
        A::Matrix{T},
        b_tmp::Vector{T},
        R::Matrix{T}) where T<:Real
    n, m = size(Dt);

    @inbounds for i in 1:n, j in 1:m
        A[i,j] = Dt[i,j];
    end

    diag_val = 1/sqrt(λ);
    @inbounds for i in 1:m, j in 1:m
        A[i+n,j] = i==j ? diag_val : 0;
    end

    qr_r!(A, R);

    copyto!(b_tmp, b);
    forward_substitution!(R, b_tmp)
    backward_substitution!(R, b_tmp)

    mul!(x_sol, Dt, b_tmp);

    nothing
end

# Forward substitution
function forward_substitution!(R::Matrix, b::Vector)
    n = size(R,1);
    @assert n==length(b) "R and b must have compatible dimensions!"

    @fastmath @inbounds for i in 1:n
        @simd for j in 1:(i-1)
            b[i] -= b[j]*R[j,i];
        end
        b[i] /= R[i,i];
    end

    nothing
end

# Backward substitution
function backward_substitution!(R::Matrix, b::Vector)
    n = size(R,1);
    @assert n==length(b) "R and b must have compatible dimensions!"

    @fastmath @inbounds for i in n:-1:1
        @simd for j in (i+1):n
            b[i] -= b[j]*R[i,j];
        end
        b[i] /= R[i,i];
    end

    nothing
end

function qr_r!(A::Matrix, R::Matrix)
    m, n = size(A);
    @assert n==size(R,1) "Dimensions of A and R must be compatible!"

    @fastmath @inbounds for i in 1:n
        R[i,i] = sqrt(colTimesCol(A, i, i));
        for j in range(i+1, stop=n)
            R[i,j] = colTimesCol(A, i , j)/R[i,i];
            α = R[i,j]/R[i,i];
            for k in 1:m
                A[k,j] -= A[k,i]*α;
            end
        end
    end

    nothing
end

function colTimesCol(A::Matrix, i::Int, j::Int)
    res = 0.0;
    n = size(A,1);
    @inbounds @simd for k in 1:n
        res += A[k,i]*A[k,j];
    end

    return res
end

end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    LSIF(σ=2.0, b=100, λ=0.001)

Least squares importance fitting.

## Parameters

* `σ` - Bandwidth of Gaussian kernel (default to `2.0`)
* `b` - Maximum number of radial basis functions (default to `10`)
* `λ` - Regularization parameter (default to `0.001`)

## References

* Kanamori et al. 2009. A Least-squares Approach to Direct
Importance Estimation

## Author

* Júlio Hoffimann (julio.hoffimann@gmail.com)
"""
@with_kw struct LSIF{T} <: DensityRatioEstimator
  σ::T=2.0
  b::Int=10
  λ::T=0.001
end

default_optlib(dre::Type{<:LSIF}) = OptimLib

available_optlib(dre::Type{<:LSIF}) = [OptimLib, JuMPLib]

function _densratio(x_nu, x_de, dre::LSIF,
                    optlib::Type{<:OptimizationLibrary})
  K_de, H, h, x_ba = _lsif_consts(x_nu, x_de, dre)
  α = _lsif_coeffs(H, h, dre, optlib)
  K_de*α
end

function _densratiofunc(x_nu, x_de, dre::LSIF,
                        optlib::Type{<:OptimizationLibrary})
  K_de, H, h, x_ba = _lsif_consts(x_nu, x_de, dre)
  α = _lsif_coeffs(H, h, dre, optlib)
  function r(x)
    K = gaussian_gramian([x], x_ba, σ=dre.σ)
    dot(K, α)
  end
end

function _lsif_consts(x_nu, x_de, dre)
  x_ba = sample(x_nu, min(length(x_nu), dre.b), replace=false)
  K_nu = gaussian_gramian(x_nu, x_ba, σ=dre.σ)
  K_de = gaussian_gramian(x_de, x_ba, σ=dre.σ)

  b = length(x_ba)
  Φ = Matrix{eltype(K_de)}(undef, b, b)
  for l′ in 1:b
    φ′ = view(K_de, :, l′)
    for l in 1:l′
      φ = view(K_de, :, l)
      Φ[l,l′] = mean(φ .* φ′)
    end
  end
  H = Symmetric(Φ)
  h = vec(mean(K_nu, dims=1))

  K_de, H, h, x_ba
end

"""
    _lsif_coeffs(x_nu, x_de, basis, dre, optlib)

Return the coefficients of LSIF basis expansion.
"""
_lsif_coeffs(H, h, dre::LSIF, optlib::Type{<:OptimizationLibrary}) =
  @error "not implemented"