text:
#========================================================================

	Steady state solution for an Aiyagari model with a diffusion process

	Translated julia code based on code from Ben Moll's Website:
		http://www.princeton.edu/%7Emoll/HACTproject.htm

		For Julia 1.0.0
========================================================================#
using LinearAlgebra, SparseArrays, Plots

# Model parameters
γ = 2 # CRRA utility
α = 1/3 # parameter for cobb-douglas production function
δ = 0.05 # captal depreciation
ρ = 0.05 # discount rate
A_prod =0.1 # term in production function

# Parameters for productivity and the poisson process q
z_1 = 1
z_2 = 2*z_1
z = [z_1 z_2]
λ_1 = 1/3 # prob of transitioning from state 1 to 2
λ_2 = 1/3  # state 2 to 1
λ = [λ_1 λ_2]
z_ave = (z_1*λ_2 + z_2*λ_1)/(λ_1 + λ_2) # Average value of z

T =200 #timing setup for shock
N = 400
dt = T/N
time = (0:N-1)*dt
max_price_it = 300
con_crit=10^(-5)
relax =0.1

# Construct the TFP sequence
corr =0.8
ν=1-corr
A_prod_t = zeros(N,1)
A_prod_t[1] = .97*A_prod

for n =1:N-1
	A_prod_t[n+1] = dt*ν*(A_prod-A_prod_t[n]) + A_prod_t[n]
end

plot(time, A_prod_t, xlims=(0,40))

# Now to set up the other grid spaces
H = 1000 # Number of elements in grid space for a
a_min = -0.8
a_max = 20.0

# Grid spaces for variables
a = LinRange(a_min, a_max, H)
a = convert(Array,a)
da = (a_max-a_min)/(H-1)

aa = [a a]
zz = ones(H,1)*z #check later

# Paramters for the simulation
maxit = 100
Ir = 40
crit = 10^(-6)
critS = 10^(-5)
Δ = 1000 #step size for HJB

# Spaces for the finite difference terms
Vaf, Vab, Vzf, Vzb, Vzz, c= [zeros(H,2) for i in 1:6]
V_r = zeros(H,2,H)

# Spaces for other terms
r_r =[]
r_min_r =[]
r_max_r =[]
KD=[]
KS=[]
S=[]
adot=[]
V_n=[]
g_r=[]
gg_s=[]

# Construct the matrix to summarize the evolution of z over time
# this matrix will not change over the finite differences and is based on the poisson process

Aswitch = [-sparse(I,H,H)*λ[1] sparse(I,H,H)*λ[2] ; sparse(I,H,H)*λ[1] -sparse(I,H,H)*λ[2]]

# Prices
r =0.04 #initial r
w = 0.05

r0=0.03
global r_min=0.01
global r_max=.99*ρ

dist = []
v0=zeros(H,2)

v0[:,1] = (w*z[1].+r.*a).^(1-γ)/(1-γ)/ρ
v0[:,2] = (w*z[2].+r.*a).^(1-γ)/(1-γ)/ρ

# The main loop, this makes sure markets clear

for ir in 1:Ir
	global r =r
	global r_min = r_min
	global r_max = r_max

	println("Main Loop iteration $(ir)")

	push!(r_r, r)
	push!(r_min_r, r_min)
	push!(r_max_r, r_max )

	push!(KD, (α*A_prod/(r+δ))^(1/(1-α))*z_ave)
	global w = (1-α)*A_prod*KD[ir].^α * z_ave^(-α)

	if ir >1
		global v0 = V_r[:,:,ir-1]
	end

	global v=v0

# Inner loop, value function iteration

	for n in 1:maxit
		global V=v
		push!(V_n, V)
		# Forward difference
		Vaf[1:H-1,:] = (V[2:H,:]-V[1:H-1,:])/da
		Vaf[H,:] = (w*z .+ r.*a_max).^(-γ)

		# Backward difference
		Vab[2:H,:] = (V[2:H,:]-V[1:H-1,:])/da
		Vab[1,:] = (w*z .+ r.*a_min).^(-γ)

		# caculate consumption and savings
		# First the forward difference case
		cf = Vaf.^(-1/γ)
		sf = w*zz + r.*aa -cf

		# consumption and savings for the backward difference
		cb = Vab.^(-1/γ)
		sb = w*zz + r.*aa -cb

		#Consumption and V' at the steady state
		c0 = w*zz + r.*aa
		Va0 = c0.^(-γ)

		# Setup to the Upwind scheme
		If = sf .> 0
		Ib = sb .< 0
		I0 = (1 .- If .- Ib)

		#Now implement the upwind scheme
		Va_Upwind = Vaf.*If + Vab.*Ib + Va0.*I0

		global c= Va_Upwind.^(-1/γ)
		u = c.^(1-γ)/(1-γ)

		X = -min.(sb,0)/da
		Y = -max.(sf,0)/da + min.(sb,0)/da
		Z = max.(sf,0)/da

		A1 = spdiagm(0=>Y[:,1], -1 => X[2:H,1], 1=> Z[1:H-1,1])
		A2 = spdiagm(0=>Y[:,2], -1 => X[2:H,2], 1=> Z[1:H-1,2])
		AA = [A1 spzeros(H,H); spzeros(H,H) A2]

		global A= AA + Aswitch

		B = (1/Δ + ρ)*sparse(I,2*H,2*H)-A

		u_stacked = [u[:,1]; u[:,2]]
		V_stacked = [V[:,1]; V[:,2]]

		b = u_stacked + V_stacked/Δ

		V_stacked = B\b

		global V = [V_stacked[1:H] V_stacked[H+1:2*H]]

		V_change = V-v

		global v = V

		push!(dist,findmax(abs.(V_change))[1])

		if n>1 && dist[n] < crit
			println("Value Function Converge Iteration = $(n)")
			break
		end
	end
	# Now solve the Fokker-Planck equation
	AT = A'
	b = zeros(2*H,1)

	# Dirty fix to avoid matrix singularities
	i_fix=1
	b[i_fix] =.1
	row = [zeros(1,i_fix-1) 1 zeros(1,H*2-i_fix)]
	AT[i_fix,:] = row

	# Now solve the linear system
	gg = AT\b
	g_sum = gg'*ones(2*H,1).*da
	gg = gg./g_sum

	global g = [gg[1:H] gg[H+1:2*H]]

	push!(gg_s,gg)
	push!(g_r, g)
	push!(adot, w*zz+r.*aa-c)
	V_r[:,:,ir] = V

	push!(KS, g[:,1]'*a*da + g[:,2]'*a*da)
	push!(S, KS[ir]-KD[ir])

	# Update aggregate capital
	if S[ir]>critS
		println("Excess Supply")
		global r_max =r
		global r = 0.5*(r+r_min)
	elseif S[ir] < -critS
		println("Excess Demand")
		global r_min =r
		global r = 0.5*(r+r_max)
	elseif abs(S[ir]) < critS
		println("Equilibrium Found!!!")
		break
	end
end

# Save some objects
v_st = v
gg_st = gg_s[end]
K_st = KS[end]
w_st =w
r_st =r
C_st = gg_s[end]'*reshape(c,H*2,1).*da
A_prod_st = A_prod

# Transition Dynamics
gg0 = gg_st

# Set up locations for preallocting values
gg_br =[]
c_t=[]
A_t=[]
C_t=[]
K_t, K_out, r_t, w_t = [zeros(N,1) for i in 1:4]
dist_it=[]

K_t = K_st*ones(N,1)

for it in 1:max_price_it
	println("Price Iteration")
	println(it)

	global w_t = (1-α)*A_prod_t.*K_t.^α .*z_ave^(-α)
	global r_t = (α*A_prod_t.*K_t.^(α-1) .*z_ave^(1-α)) .- δ

	global V=v_st

	for n in N:-1:1
		# Forward difference
		Vaf[1:H-1,:] = (V[2:H,:]-V[1:H-1,:])/da
		Vaf[H,:] = (w_t[n]*z .+ r_t[n].*a_max).^(-γ)

		# Backward difference
		Vab[2:H,:] = (V[2:H,:]-V[1:H-1,:])/da
		Vab[1,:] = (w_t[n]*z .+ r_t[n].*a_min).^(-γ)

		# caculate consumption and savings
		# First the forward difference case
		cf = Vaf.^(-1/γ)
		sf = w_t[n]*zz + r_t[n].*aa -cf

		# consumption and savings for the backward difference
		cb = Vab.^(-1/γ)
		sb = w_t[n]*zz + r_t[n].*aa -cb

		#Consumption and V' at the steady state
		c0 = w_t[n]*zz + r_t[n].*aa
		Va0 = c0.^(-γ)

		# Setup to the Upwind scheme
		If = sf .> 0
		Ib = sb .< 0
		I0 = (1 .- If .- Ib)

		#Now implement the upwind scheme
		Va_Upwind = Vaf.*If + Vab.*Ib + Va0.*I0

		global c= Va_Upwind.^(-1/γ)
		u = c.^(1-γ)/(1-γ)
		push!(c_t,c)

		X = -min.(sb,0)/da
		Y = -max.(sf,0)/da + min.(sb,0)/da
		Z = max.(sf,0)/da

		A1 = spdiagm(0=>Y[:,1], -1 => X[2:H,1], 1=> Z[1:H-1,1])
		A2 = spdiagm(0=>Y[:,2], -1 => X[2:H,2], 1=> Z[1:H-1,2])
		AA = [A1 spzeros(H,H); spzeros(H,H) A2]

		global A= AA + Aswitch

		push!(A_t,A)

		B = (1/dt + ρ)*sparse(I,2*H,2*H)-A

		u_stacked = [u[:,1]; u[:,2]]
		V_stacked = [V[:,1]; V[:,2]]

		b = u_stacked + V_stacked/dt

		V_stacked = B\b

		global V = [V_stacked[1:H] V_stacked[H+1:2*H]]
	end
	# Now update the distribution
	push!(gg_br,gg0)
	for n = 1:N
		AT = A_t[end-(n-1)]'
		push!(gg_br,(sparse(I,2*H,2*H)-AT*dt)\gg_br[n])
		global K_out[n] = (gg_br[n][1:H])'*a*da + gg_br[n][H+1:2*H]'*a*da
		push!(C_t,gg_br[n]'*reshape(c_t[n],2*H,1).*da)
	end

	push!(dist_it, findmax(abs.(K_out-K_t))[1])

	global K_t = relax.*K_out + (1-relax).*K_t

	if dist_it[it]<con_crit
		println("Equilibrium Found")
		break
	end

	global gg_br=[]
end

# Inequality measurements

Wealth_neg_t,Wealth_Gini_t,CapInc_Gini_t,Inc_Gini_t,top_inc_t,top_wealth_t,neg_frac_t = [zeros(1,400) for i in 1:7]

for n in 1:N
	Wealth_neg_t[n] = gg_br[n][1:H]'*min.(a,0).*da + gg_br[n][H+1:2*H]'*min.(a,0).*da
	g_a_count = gg_br[n][1:H] + gg_br[n][H+1:2*H]

	# Discrete Gini for checking
	g_a = g_a_count*da
	S_a = cumsum(g_a.*a,dims=1)/sum(g_a.*a)
	trapez_a = .5 .* (S_a[1]*g_a[1]+sum((S_a[2:H]+S_a[1:H-1]).*g_a[2:H]))
	Wealth_Gini_t[n] = 1 .-2*trapez_a

	# Gini of Capital income
	yk = r_t[n].*a
	g_yk=g_a

	S_yk = cumsum(g_yk.*yk,dims=1)/sum(g_yk.*yk)
	trapez_yk = .5 .* (S_yk[1]*g_yk[1]+sum((S_yk[2:H]+S_yk[1:H-1]).*g_yk[2:H]))
	CapInc_Gini_t[n] = 1 -2*trapez_yk

	# Gini of total income
	y = w_t[n]*zz + r_t[n]*aa
	Ny = 2*H
	yy =reshape(y,Ny,1)
	index = sortperm(yy[:])
	yy=yy[index]
	g_y = gg_br[n][index].*da

	S_y = cumsum(g_y.*yy,dims=1)/sum(g_y.*yy)
	trapez_y = .5 .* (S_y[1]*g_y[1]+sum((S_y[1:H-1]+S_y[1:H-1]).*g_y[2:H]))
	Inc_Gini_t[n] = 1 -2*trapez_y

	G_y =cumsum(g_y,dims=2)
	G_a = cumsum(g_a,dims=2)

	#Top 10% income share
	p1=0.1
	obj, index = findmin(abs.((1 .-G_y) .-p1))
	top_inc_t[n]=1 .-S_y[index]

	#Top 10% Wealth share
	obj, index = findmin(abs.((1 .-G_a) .-p1))
	top_wealth_t[n]=1 .-S_a[index]

	neg_frac_t[n]=findmax(G_a[a.<=0])[1]

end

p1 = plot(time, A_prod_t,xlims=(0,100), title="Aggregate Prod.",legend=false)
plot!(time,A_prod_st*ones(N,1),line=:dash,color=:black)

p2 = plot(time, K_t,xlims=(0,100),ylims=(0.293,0.3),title="Aggregate Capital St.",legend=false )
plot!(p2,time,K_st.*ones(N,1),line=:dash,color=:black)

p3 = plot(time,w_t, xlims=(0,100),title="Wage",legend=false)
plot!(p3,time,w_st.*ones(N,1),line=:dash,color=:black)

p4=plot(time,r_t,xlims=(0,100),title="Interest Rate",legend=false)
plot!(p4,time,r_st.*ones(N,1),line=:dash,color=:black)

p5 =plot(time,Wealth_Gini_t',xlims=(0,100),title="Wealth Gini",legend=false)
plot!(p5,time,Wealth_Gini_t[N].*ones(N,1),line=:dash,color=:black)

p6 =plot(time,Inc_Gini_t',xlims=(0,100),title="Income Gini",legend=false)
plot!(p6,time,Inc_Gini_t[N].*ones(N,1),line=:dash,color=:black)

plot(p1,p2,p3,p4,p5,p6,layout=(2,3),legend=false)
png("Inequality_Measures")


text:
"""
$(SIGNATURES)

Helper function for `convert_inter_html` that processes an extracted block given a latex context
`lxc` and returns the processed html that needs to be plugged in the final html.
"""
function convert_block(β::AbstractBlock, lxcontext::LxContext)::AS
    # case for special characters / html entities
    β isa HTML_SPCH     && return ifelse(isempty(β.r), β.ss, β.r)
    # Return relevant interpolated string based on case
    βn = β.name
    βn ∈ MD_HEADER         && return convert_header(β, lxcontext.lxdefs)
    βn == :CODE_INLINE     && return html_code_inline(content(β) |> htmlesc)
    βn == :CODE_BLOCK_LANG && return convert_code_block(β.ss)
    βn == :CODE_BLOCK_IND  && return convert_indented_code_block(β.ss)
    βn == :CODE_BLOCK      && return html_code(strip(content(β)), "{{fill lang}}")
    βn == :ESCAPE          && return chop(β.ss, head=3, tail=3)
    βn == :FOOTNOTE_REF    && return convert_footnote_ref(β)
    βn == :FOOTNOTE_DEF    && return convert_footnote_def(β, lxcontext.lxdefs)
    βn == :LINK_DEF        && return ""

    # Math block --> needs to call further processing to resolve possible latex
    βn ∈ MATH_BLOCKS_NAMES && return convert_math_block(β, lxcontext.lxdefs)

    # Div block --> need to process the block as a sub-element
    if βn == :DIV
        raw_ct = ignore_starting_line_spaces(content(β))
        ct, _ = convert_md(raw_ct * EOS, lxcontext.lxdefs;
                           isrecursive=true, has_mddefs=false)
        name = chop(otok(β).ss, head=2, tail=0)
        return html_div(name, ct)
    end

    # default case, ignore block (should not happen)
    return ""
end
convert_block(β::LxCom, λ::LxContext) = resolve_lxcom(β, λ.lxdefs)

"""
MATH_BLOCKS_PARENS

Dictionary to keep track of how math blocks are fenced in standard LaTeX and how these fences need
to be adapted for compatibility with KaTeX. Each tuple contains the number of characters to chop
off the front and the back of the maths expression (fences) as well as the KaTeX-compatible
replacement.
For instance, `\$ ... \$` will become `\\( ... \\)` chopping off 1 character at the front and the
back (`\$` sign).
"""
const MATH_BLOCKS_PARENS = LittleDict{Symbol, Tuple{Int,Int,String,String}}(
    :MATH_A     => ( 1,  1, "\\(", "\\)"),
    :MATH_B     => ( 2,  2, "\\[", "\\]"),
    :MATH_C     => ( 2,  2, "\\[", "\\]"),
    :MATH_ALIGN => (13, 11, "\\[\\begin{aligned}", "\\end{aligned}\\]"),
    :MATH_EQA   => (16, 14, "\\[\\begin{array}{c}", "\\end{array}\\]"),
    :MATH_I     => ( 4,  4, "", "")
    )

"""
$(SIGNATURES)

Helper function for the math block case of `convert_block` taking the inside of a math block,
resolving any latex command in it and returning the correct syntax that KaTeX can render.
"""
function convert_math_block(β::OCBlock, lxdefs::Vector{LxDef})::String
    # try to find the block out of `MATH_BLOCKS_PARENS`, if not found, error
    pm = get(MATH_BLOCKS_PARENS, β.name) do
        throw(MathBlockError("Unrecognised math block name."))
    end

    # convert the inside, decorate with KaTex and return, also act if
    # if the math block is a "display" one (with a number)
    inner = chop(β.ss, head=pm[1], tail=pm[2])
    htmls = IOBuffer()
    if β.name ∉ [:MATH_A, :MATH_I]
        # NOTE: in the future if allow equation tags, then will need an `if`
        # here and only increment if there's no tag. For now just use numbers.

        # increment equation counter
        PAGE_EQREFS[PAGE_EQREFS_COUNTER] += 1

        # check if there's a label, if there is, add that to the dictionary
        matched = match(r"\\label{(.*?)}", inner)

        if !isnothing(matched)
            name   = refstring(matched.captures[1])
            write(htmls, "")
            inner  = replace(inner, r"\\label{.*?}" => "")
            # store the label name and associated number
            PAGE_EQREFS[name] = PAGE_EQREFS[PAGE_EQREFS_COUNTER]
        end
    end
    # assemble the katex decorators, the resolved content etc
    write(htmls, pm[3], convert_md_math(inner * EOS, lxdefs, from(β)), pm[4])
    return String(take!(htmls))
end

"""
$(SIGNATURES)

Helper function for the case of a header block (H1, ..., H6).
"""
function convert_header(β::OCBlock, lxdefs::Vector{LxDef})::String
    hk       = lowercase(string(β.name)) # h1, h2, ...
    title, _ = convert_md(content(β) * EOS, lxdefs;
                          isrecursive=true, has_mddefs=false)
    rstitle  = refstring(title)
    # check if the header has appeared before and if so suggest
    # an altered refstring; if that altered refstring also exist
    # (pathological case, see #241), then extend it with a letter
    # if that also exists (which would be really crazy) add a random
    # string
    keys_headers = keys(PAGE_HEADERS)
    if rstitle in keys_headers
        title, n, lvl = PAGE_HEADERS[rstitle]
        # then update the number of occurence
        PAGE_HEADERS[rstitle] = (title, n+1, lvl)
        # update the refstring, note the double _ (see refstring)
        rstitle *= "__$(n+1)"
        PAGE_HEADERS[rstitle] = (title, 0, parse(Int, hk[2]))
    else
        PAGE_HEADERS[rstitle] = (title, 1, parse(Int, hk[2]))
    end
    # return the title
    return html_hk(hk, html_ahref_key(rstitle, title); id=rstitle)
end

"""Convenience function to increment the eval' code block counter"""
increment_code_head() = (LOCAL_PAGE_VARS["jd_code_head"].first[] += 1)

"""Convenience function to mark `jd_code_eval` as true (subsequent blocks will be evaled)"""
toggle_jd_code_eval() = (LOCAL_PAGE_VARS["jd_code_eval"].first[] = true)

"""
$SIGNATURES

Helper function to eval a code block, write it where appropriate, and finally return
a resolved block that can be displayed in the html.
"""
function eval_and_resolve_code(code::AS, rpath::AS;
                               eval::Bool=true, nopush::Bool=false)::String
    # Here we have a julia code block that was provided with a script path
    # It will consequently be
    #  1. written to script file unless it's already there
    #  2. eval with redirect (unless file+output already there)
    #  3. inserted after cleaning out lines (see resolve_lx_input)

    # start by adding it to code scope
    nopush || push!(LOCAL_PAGE_VARS["jd_code_scope"].first, rpath, code)

    # form the path
    path = resolve_assets_rpath(rpath; canonical=true, code=true)

    # lazy names are allowed without extensions, add one if that's the case
    endswith(path, ".jl") || (path *= ".jl")

    # output directory etc
    out_path, fname = splitdir(path)
    out_path = mkpath(joinpath(out_path, "output"))
    out_name = splitext(fname)[1] * ".out"
    res_name = splitext(fname)[1] * ".res"
    res_path = joinpath(out_path, res_name)
    out_path = joinpath(out_path, out_name)

    # if we're in the no-eval case, check that the relevant files are
    # there otherwise do re-eval with re-write
    if !eval && isfile(path) && isfile(out_path) && isfile(res_path)
        # just return the resolved code block
        return resolve_lx_input_hlcode(rpath, "julia")
    end

    write(path, MESSAGE_FILE_GEN_JMD * code)
    SILENT_MODE[] || print(rpad("\r→ evaluating code [...] ($(CUR_PATH[]), $rpath)", 79) * "\r")
    # - execute the code while redirecting stdout to file
    Logging.disable_logging(Logging.LogLevel(3_000))
    res = nothing
    open(out_path, "w") do outf        # for stdout
        redirect_stdout(outf) do
            res = try
                include(path)
            catch e
                print("There was an error running the code:\n$(e.error)")
            end
        end
    end
    open(res_path, "w") do outf
        redirect_stdout(outf) do
            show(stdout, "text/plain", res)
        end
    end
    Logging.disable_logging(Logging.Debug)
    SILENT_MODE[] || print(rpad("\r→ evaluating code [✓]", 79) * "\r")

    # resolve the code block (highlighting) and return it
    return resolve_lx_input_hlcode(rpath, "julia")
end

"""
$(SIGNATURES)

Helper function for the code block case of `convert_block`.
"""
function convert_code_block(ss::SubString)::String
    fencer = ifelse(startswith(ss, "`````"), "`````", "```")
    reg    = Regex("$fencer([a-z-]*)(\\:[a-zA-Z\\\\\\/-_\\.]+)?\\s*\\n?((?:.|\\n)*)$fencer")
    m      = match(reg, ss)
    lang   = m.captures[1]
    rpath  = m.captures[2]
    code   = strip(m.captures[3])

    # if no rpath is given, the code is not eval'
    isnothing(rpath) && return html_code(code, lang)

    # if the code is not in julia, it's not eval'ed
    if lang!="julia"
        @warn "Eval of non-julia code blocks is not yet supported."
        return html_code(code, lang)
    end

    # path currently has an indicative `:` we don't care about
    rpath = rpath[2:end]

    # extract handles of relevant local page variables
    reeval = LOCAL_PAGE_VARS["reeval"].first         # full page re-eval
    eval   = LOCAL_PAGE_VARS["jd_code_eval"].first[] # eval toggle from given point
    freeze = LOCAL_PAGE_VARS["freezecode"].first
    scope  = LOCAL_PAGE_VARS["jd_code_scope"].first
    head   = increment_code_head()

    # In the case of forced re-eval, we don't care about the
    # code scope just force-reeval everything sequentially
    if FORCE_REEVAL[] || reeval || eval
        length(scope.codes) ≥ head && purgefrom!(scope, head)
        return eval_and_resolve_code(code, rpath)
    end

    # Now we have Julia code with a path; check briefly if there
    # are explicit flags indicating we should not eval; if that's
    # the case then there will be a check to see if the relevant
    # files exist, if they don't exist the code *will* be eval'ed
    # (see `eval_and_resolve_code`)
    if FULL_PASS[] || freeze
        length(scope.codes) ≥ head && purgefrom!(scope, head)
        return eval_and_resolve_code(code, rpath, eval=false)
    end

    # Here we're either in
    # A. full pass but with forced-reeval
    # B. local pass with non-frozen code

    # check if the page we're looking at is in scope
    if CUR_PATH[] != CUR_PATH_WITH_EVAL[]
        # we're necessarily at the first code block of the page.
        # need to re-instantiate a code scope; note that if we
        # are here then necessarily a def_LOCAL_PAGE_VARS was
        # called, so LOCAL_PAGE_VARS["jd_code_head"] points to 1
        reset!(scope)
        # keep track that the page is now in scope
        CUR_PATH_WITH_EVAL[] = CUR_PATH[]
        # flag rest of page as to be eval-ed (might be stale)
        toggle_jd_code_eval()
        # eval and resolve code
        return eval_and_resolve_code(code, rpath)
    end

    # we're in scope, compare the code block with the
    # current code scope and act appropriately
    ncodes = length(scope.codes)

    # there is only one case where we might not add and eval
    # --> if c ≤ length(code_dict)  -- code block may be among seen ones
    # --> code == code_dict[rpath]  -- the content matches exactly
    if (head ≤ ncodes)
        if (scope.rpaths[head] == rpath && code == scope.codes[head])
            # resolve with no eval (the function will check if the files are
            # present and if they're not, there will be an evaluation)
            return eval_and_resolve_code(code, rpath; eval=false, nopush=true)
        else
            # purge subsequent code blocks as stale
            purgefrom!(scope, head)
            # flag rest of page as to be eval-ed (stale)
            toggle_jd_code_eval()
        end
    end
    return eval_and_resolve_code(code, rpath)
end

"""
$(SIGNATURES)

Helper function for the indented code block case of `convert_block`.
"""
function convert_indented_code_block(ss::SubString)::String
    # 1. decrease indentation of all lines (either frontal \n\t or \n⎵⎵⎵⎵)
    code = replace(ss, r"\n(?:\t| {4})" => "\n")
    # 2. return; lang is a LOCAL_PAGE_VARS that is julia by default and can be set
    sc = strip(code)
    isempty(sc) && return ""
    return html_code(sc, "{{fill lang}}")
end

"""
$(SIGNATURES)

Helper function to convert a `[^1]` into a html sup object with appropriate ref and backref.
"""
function convert_footnote_ref(β::Token)::String
    # β.ss is [^id]; extract id
    id = string(match(r"\[\^(.*?)\]", β.ss).captures[1])
    # add it to the list of refs unless it's been seen before
    pos = 0
    for (i, pri) in enumerate(PAGE_FNREFS)
        if pri == id
            pos = i
            break
        end
    end
    if pos == 0
        push!(PAGE_FNREFS, id)
        pos = length(PAGE_FNREFS)
    end
    return html_sup("fnref:$id", html_ahref(url_curpage() * "#fndef:$id", "[$pos]"; class="fnref"))
end

"""
$(SIGNATURES)

Helper function to convert a `[^1]: ...` into a html table for the def.
"""
function convert_footnote_def(β::OCBlock, lxdefs::Vector{LxDef})::String
    # otok(β) is [^id]:
    id = match(r"\[\^(.*?)\]:", otok(β).ss).captures[1]
    pos = 0
    for (i, pri) in enumerate(PAGE_FNREFS)
        if pri == id
            pos = i
            break
        end
    end
    if pos == 0
        # this was never referenced before, so probably best not to show it
        return ""
    end
    # need to process the content which could contain stuff
    ct, _ = convert_md(content(β) * EOS, lxdefs;
                       isrecursive=true, has_mddefs=false)
    """

            $(html_ahref(url_curpage() * "#fnref:$id", "[$pos]"))
            $(ct)

    """
end

ext: jl


text:
struct ReverseDiffBackend <: AbstractDiffBackend end

function Manifolds._gradient(f, p, ::ReverseDiffBackend)
    return ReverseDiff.gradient(f, p)
end

function Manifolds._gradient!(f, X, p, ::ReverseDiffBackend)
    return ReverseDiff.gradient!(X, f, p)
end

if default_differential_backend() === NoneDiffBackend()
    set_default_differential_backend!(ReverseDiffBackend())
end


text:
struct SOL

    inner_boundary :: Vector{Complex{Float64}}
    outer_boundary :: Vector{Complex{Float64}}
    psi :: Function
    iota :: Function
    bvec :: Function

end

function SOL(Bpol, Bz, Nc, r0, Nr, Nt, Nz, ds, divs::Matrix{Float64})

    # compute inner and outer boundaries
    outer_chain = Complex{Float64}[]
    inner_chain = Vector{Complex{Float64}}(undef, Nc)
    thetas = LinRange(0, 2*pi, Nc+1)[1:Nc]
    for i=1:Nc
        if !any([divs[1,j] < thetas[i] < divs[2,j] for j=1:size(divs)[2]])
            outer_chain = [outer_chain; [0.95*exp(im*thetas[i])]]
        end
        inner_chain[i] = r0*exp(im*thetas[i])
    end

    # compute boundary value for psi
    psib = Vector{Float64}(undef, 1001)
    ts = LinRange(0, 2*pi, 1002)[1:1001]
    for i=1:1001
        for j=1:size(divs)[2]
            if divs[1,j] < ts[i] < divs[2,j]
                u = abs(2*(ts[i] - (divs[1,j] + divs[2,j])/2) / (divs[2,j] - divs[1,j]))
                psib[i] = divs[3,j]*(1 - u)
            else
                psib[i] = 0.0
            end
        end
    end

    # find fourier coefficients for flux function
    C = fft(psib)
    c0 = real(C[1]) / 1001
    a = 2*real.(C[2:501]) / 1001
    b = -2*imag.(C[2:501]) / 1001

    # define flux function
    psi(z) = begin
        r = abs(z)
        t = angle(z)

        if r < 0.9
            kmax = minimum([500, abs(Int64(ceil(8 / log10(r/0.9))))])
            c0 + Bpol*log(r/0.9) + sum([(r/0.9)^k * (a[k]*cos(k*t) + b[k]*sin(k*t)) for k=1:kmax])
        else
            c0 + Bpol*log(r/0.9) + sum([(r/0.9)^(-k) * (a[k]*cos(k*t) + b[k]*sin(k*t)) for k=1:500])
        end
    end

    # define poloidal and total magnetic field
    gradPsi(r) = ForwardDiff.gradient(r -> psi(r[1] + im*r[2]), r)
    Bp(r) = begin dpsi = gradPsi(r); [-dpsi[2], dpsi[1]] end
    B(r) = [Bp(r)..., Bz]
    bvec(r) = begin Bvec = B(r); Bvec / norm(Bvec) end

    # define 1/q
    iota(z) = 2*pi*norm(Bp([real(z), imag(z)])) / Bz

    return SOL(inner_chain, outer_chain, psi, iota, bvec)
end

function fluxmap(z::ComplexF64, dPhi::Float64, ds::Float64, sol::SOL)

    points, dS = trace_fieldline([real(z), imag(z), 0.0], sol.bvec, r -> sol.psi(r[1] + im*r[2]), ds) do r
        abs(r[3]) > dPhi
    end
    z_final = Complex(points[1:2,end]...)

    return z_final, dS
end


text:
using DataFrames, Plots, NumericalMethodsforEngineers
pyplot()

old = pwd()
ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  f(x::Float64, y::Vector{Float64}) = (x + y[1])^2
  steps = 25
  h = 0.02

  x = 0.0
  y = [1.0]

  # y(3.0) = 4.2165

  r = Array{Float64,2}[]
  push!(r, euler(f, x, y, steps, h))
  push!(r, modified_euler(f, x, y, steps, h))
  push!(r, mid_point_euler(f, x, y, steps, h))
  push!(r, runga_kutta_4(f, x, y, steps, h))
  rr = hcat(r[1], r[2], r[3], r[4])
  println()
  rr[21:26,:] |> display
  println()

  xs = Float64[(i-1)*h for i in 1:steps+1]
  titles = ["Euler", "Modified_Euler", "Mid_Point_Euler", "Runga_Kutta_4"]

  p1 = plot(xs, r[1][:, 2])
  plot!(p1, xs, r[2][:, 2])
  plot!(p1, xs, r[3][:, 2])
  plot!(p1, xs, r[4][:, 2])
  show(p1)
  savefig("Ex7.2.png")
  gui()

end


text:
using ADI
using BenchmarkTools
using CSV
using DataFrames
using HCIDatasets: BetaPictoris, HR8799
using HCIToolbox
using InvertedIndices
using PyCall

results = []

# parameters
nbranch = 3
fwhm = 8

# julia benchmarks

alg = PCA(20)
psf = construct(Kernels.Normal(5), (21, 21))

for dataset in (BetaPictoris, HR8799)
    @info "Benchmarking - $dataset Contrast Curve"
    cube = dataset[:cube]
    angles = dataset[:pa]

    time_elapsed = @belapsed contrast_curve($alg, $cube, $angles, $psf; fwhm=fwhm, nbranch=nbranch)

    @info "Julia" time=time_elapsed
    push!(results, (framework="ADI.jl", alg="pca_20", N=length(cube), time=time_elapsed))
end

# python benchmarks

vip = pyimport("vip_hci")

for dataset in (BetaPictoris, HR8799)
    @info "Benchmarking - $dataset Contrast Curve"
    cube = dataset[:cube]
    angles = dataset[:pa]

    time_elapsed = @belapsed vip.metrics.contrast_curve($cube, $angles,  $psf, fwhm=fwhm, pxscale=1e-3, starphot=Metrics.estimate_starphot($cube, fwhm), algo=vip.pca.pca, nbranch=nbranch, verbose=false, plot=false, ncomp=20)

    @info "Python" time=time_elapsed
    push!(results, (framework="VIP", alg="pca_20", N=length(cube), time=time_elapsed))
end

# load results file and update
path = joinpath(@__DIR__, "contrast_benchmarks.csv")
df = isfile(path) ? CSV.File(path) |> DataFrame : DataFrame()
out = append!(DataFrame(results), df)
unique!(out, Not(:time)) |> CSV.write(path)


text:
module Gurobi

const _DEPS_FILE = joinpath(dirname(@__FILE__),"..","deps","deps.jl")
if isfile(_DEPS_FILE)
    include(_DEPS_FILE)
else
    error("""
        Gurobi not properly installed. Please run Pkg.build(\"Gurobi\"). For
        more information go to https://github.com/JuliaOpt/Gurobi.jl
    """)
end

import SparseArrays, LinearAlgebra

### imports

import Base.show, Base.copy, Base.read

### exports
export

# grb_env
free_env,

# grb_params
getparam, setparam!, setparams!,

# grb_model
set_sense!, update_model!, reset_model!, get_tune_result!,
read_model, write_model, tune_model, presolve_model, fixed_model,
copy, read,

# grb_attrs
model_name, model_sense, model_type,
num_vars, num_constrs, num_sos, num_qconstrs,
num_cnzs, num_qnzs, num_qcnzs,
is_qp, is_qcp, is_mip,

lowerbounds, upperbounds, objcoeffs, set_objcoeffs!,

# grb_vars
GRB_CONTINUOUS, GRB_BINARY, GRB_INTEGER,
add_var!, add_vars!, add_cvar!, add_cvars!,
add_bvar!, add_bvars!, add_ivar!, add_ivars!,
del_vars!,

# grb_constrs
add_constr!, add_constrs!, add_constrs_t!,
add_rangeconstr!, add_rangeconstrs!, add_rangeconstrs_t!,
get_constrmatrix, add_sos!, del_constrs!, chg_coeffs!,

# grb_quad
add_qpterms!, add_qconstr!,

# higher level
gurobi_model,

# grb_solve
optimize, computeIIS, get_solution,
get_status, OptimInfo, get_optiminfo, get_objval

### include source files

include("grb_common.jl")
include("grb_env.jl")

include("grb_model.jl")
include("grb_params.jl")
include("grb_vars.jl")
include("grb_attrs.jl")
include("grb_constrs.jl")
include("grb_quad.jl")
include("grb_highlevel.jl")

include("grb_solve.jl")
include("grb_callbacks.jl")

include("MPB_wrapper.jl")
include("MOI_wrapper.jl")
include("MOI_callbacks.jl")
include("MOI_multi_objective.jl")

end


text:
# # Эксперимент №1

# ## Набор параметров
α       = 0.004;        # Параметр регуляризации
w       = 0.0005;       # Эмпирический параметр регуляризации
S       = 35000;        # Количество итераций
β       = 0.003;        # Шаг минимизации
#
x_tp    = 0.25;         # Стартовое местоположение фронта
T_end   = 0.2;          # Регулируем конечное местоположение фронта
ε       = 10^(-1.5);    # Крутизна фронта
Nx      = 500;          # Число интервалов по ``X``
Mt      = 1000;         # Число интервалов по ``T``

# -----------------------------------------------------------------------------
using NonLinearReactionAdvectionDiffusionWithFrontData
using NonLinearReactionAdvectionDiffusionWithFrontData: heterogeneity_map;
using Serialization;
using Plots; pyplot();
# -----------------------------------------------------------------------------

# ## Решение на точных данных

# ### Решение прямой задачи для генерирования экспериментальной информации
# -----------------------------------------------------------------------------
a, b, t₀, T, N, M, ε, Xₙ, Tₘ, qₙ, ulₘ, urₘ, u₀ = dparams(x_tp = x_tp,
                                                         ε = ε,
                                                         Nx = Nx,
                                                         Mt = Mt,
                                                         T_end = T_end);
u, XX, TP = solve(u₀, Xₙ, N, Tₘ, M, ε, ulₘ, urₘ, qₙ);
ϕl, ϕr, ϕ, f1_data, f2_data = generate_obs_data(u, Xₙ, N, Tₘ, M, qₙ, ulₘ, urₘ);
directP = draft(u, Xₙ, N, Tₘ, M, title = "Эскиз прямого решения")
savefig(directP, "direct1.png");
# -----------------------------------------------------------------------------

# ### Найдем начальное приближение
# -----------------------------------------------------------------------------
q_guess = initial_guess(f1_data, Xₙ, N, Tₘ, M, ulₘ, urₘ, α);
plot(Xₙ, qₙ, label="Истинное");
plot!(Xₙ, q_guess, label="Найденное")
# -----------------------------------------------------------------------------

# -----------------------------------------------------------------------------
# В следующей строке, мы записываем все параметры в latex формате,
# не забываем экранировать все спецсимволы.
# Дальше, мы отобразим все эти параметры на графике.
using Printf;
params = "\$w = $(w), \\varepsilon = $(ε), "*
"\\beta = $(β), f_1 \\in $(@sprintf("[%.2f, %.2f]", extrema(f1_data)...)), "*
"\\alpha = $(α)\$"
nothing #hide
# -----------------------------------------------------------------------------

# ### Старт с найденного приближения на точных данных
# -----------------------------------------------------------------------------
q₀ = q_guess;
@time qs, Js, Qs = minimize(q₀, u₀, ulₘ, urₘ, Xₙ, N, Tₘ, M, ε, f1_data, f2_data,
                            S = S, β = β, w = w, showProgress = true)
serialize("qs.jld", qs);
serialize("Js.jld", Js);
serialize("Qs.jld", Qs);

# ### Эскиз процесса минимизации
# -----------------------------------------------------------------------------
a, b, c = minimization_draft(qₙ, Qs, Xₙ, N, Js, zoom = true,
                              annotate_string = params, zoom_chunk = 17//18)
plot(a, size = (800, 800))
plot(b)
plot(c)
withguessP = plot(a, b, c, layout = (3,1), size = (800, 2400));
savefig(withguessP, "withguess1.png")
nothing; #hide
# -----------------------------------------------------------------------------

# ## Слабо зашумленные данные
# -----------------------------------------------------------------------------
using Random
using NonLinearReactionAdvectionDiffusionWithFrontData: front_velocity;
rng = MersenneTwister(1234)

# Создадим мелкий шум.
f1_data_noised = f1_data + randn(rng, length(f1_data)) .* 0.00008
f2_data_noised = f2_data + randn(rng, length(f2_data)) .* 0.00008
plot(Tₘ, f1_data_noised, label="Зашумленные \$f_1\$")
plot!(Tₘ, f2_data_noised, label="Зашумленные \$f_2\$")

# Убедимся, что численное дифференцирование зашумленных функций --- это плохо.
v_f1 = front_velocity(f1_data_noised, Tₘ, M);
plot(Tₘ, v_f1[5:end-5], label="Скорость зашумленного \$f_1\$")

# Регуляризация во время поиска начального приближения справляется
q_guess = initial_guess(f1_data_noised, Xₙ, N, Tₘ, M, ulₘ, urₘ, α*3);
plot(Xₙ, qₙ, label="Истинное");
plot!(Xₙ, q_guess, label="Найденное")

# ### Старт с найденного приближения на зашумленных данных
# -----------------------------------------------------------------------------
q₀ = q_guess;
@time qs, Js, Qs = minimize(q₀, u₀, ulₘ, urₘ, Xₙ, N, Tₘ, M, ε, f1_data, f2_data,
                            S = S, β = β, w = w, showProgress = true)
serialize("qs_noised.jld", qs);
serialize("Js_noised.jld", Js);
serialize("Qs_noised.jld", Qs);

# ### Эскиз процесса минимизации
a, b, c = minimization_draft(qₙ, Qs, Xₙ, N, Js, zoom = true,
                              annotate_string = params, zoom_chunk = 17//18)
plot(a)
plot(b)
plot(c)
noisedP = plot(a, b, c);
savefig(noisedP, "noised.png");
nothing; #hide
# -----------------------------------------------------------------------------

# ### Старт с найденного приближения на точных данных
# -----------------------------------------------------------------------------
q₀ = q_guess;
@time qs, Js, Qs = minimize(q₀, u₀, ulₘ, urₘ, Xₙ, N, Tₘ, M, ε, f1_data, f2_data,
                            S = S, β = β, w = w, showProgress = true)
serialize("qs_noised.jld", qs);
serialize("Js_noised.jld", Js);
serialize("Qs_noised.jld", Qs);

# ### Эскиз процесса минимизации
a, b, c = minimization_draft(qₙ, Qs, Xₙ, N, Js, zoom = true,
                              annotate_string = params, zoom_chunk = 17//18)
plot(a)
plot(b)
plot(c)
noisedP = plot(a, b, c);
savefig(noisedP, "noised.png");
nothing; #hide
# -----------------------------------------------------------------------------


text:
"""
    expected_shortfall(returns, confidence, method; multiplier=1.0)

Computes the Expected Shortfall (ES), also known as Conditional Value-at-Risk (CVaR), Average Value-at-Risk (AVaR) or Expected Tail Loss (ETL). The ES is the expected return on the asset in the worst `α%` of cases, therefore quantifies the tail-risk of an asset. It is calculated by averaging all of the returns in the distribution that are worse than the VaR of the portfolio at a given significance level `α`. For instance, for a 5% significance level, the expected shortfall is calculated by taking the average of returns in the worst 5% of cases.

Expected Shortfall puts emphasis on the tail of the loss distribution, whereas Value-at-risk neglects this aspect.

# Arguments
- `returns`:        Vector of asset returns.
- `α`:              Significance level, e.g. use `0.05` for 95% confidence, or `0.01` for 99% confidence.
- `method`:         Distribution estimation method: `:historical`, `:gaussian` or `:cornish_fisher`.
- `multiplier`:     Optional scalar multiplier, i.e. use `12` to annualize monthly returns, and use `252` to annualize daily returns.

# Methods
- `:historical`:        Historical based on empirical distribution of returns.
- `:gaussian`:          Gaussian distribution based on parametric fit (mean, variance).
- `:cornish_fisher`:    Cornish-Fisher based on Gaussian parametric distribution fit adjusted for third and fourth moments (skewness, kurtosis). Cornish-Fisher expansion aims to approximate the quantile of a true distribution by using higher moments (skewness and kurtosis) of that distribution to adjust for its non-normality. See https://thema.u-cergy.fr/IMG/pdf/2017-21.pdf for details.

# Sources
- Amédée-Manesme, Charles-Olivier and Barthélémy, Fabrice and Maillard, Didier (2017). Computation of the Corrected Cornish–Fisher Expansion using the Response Surface Methodology: Application to VaR and CVaR. THEMA Working Paper n°2017-21, Université de Cergy-Pontoise, France.
"""
function expected_shortfall(returns, α, method::Symbol; multiplier=1.0)
    if method == :historical
        # average return below significance level (quantile)
        sorted = sort(returns)
        idx = floor(Int64, length(sorted) * α)
        return mean(sorted[1:idx]) * sqrt(multiplier)
    elseif method == :gaussian
        # derivation: http://blog.smaga.ch/expected-shortfall-closed-form-for-normal-distribution/
        q = quantile(Normal(), α)
        μ = mean(returns)
        σ = std(returns; corrected=false)
        return (μ - σ*pdf(Normal(), q)/α) * sqrt(multiplier)
    elseif method == :cornish_fisher
        # third/fourth moment adjusted Gaussian distribution fit
        # https://papers.ssrn.com/sol3/papers.cfm?abstract_id=1024151
        q = quantile(Normal(), α)
        S = skewness(returns)
        K = kurtosis(returns; method=:excess)
        g = q + 1/6*(q^2-1)S + 1/24*(q^3-3q)*K - 1/36*(2q^3-5q)*S^2
        ϕ = pdf(Normal(), g)
        EG2 = -1/α*ϕ * (1 + 1/6*(g^3)*S + 1/72*(g^6 - 9g^4 + 9g^2 + 3)*S^2 + 1/24*(g^4 - 2g^2 - 1)*K)
        μ = mean(returns)
        σ = std(returns; corrected=false)
        return (μ + σ*EG2) * sqrt(multiplier)
    end

    throw(ArgumentError("Passed method parameter '$(method)' is invalid, must be one of :historical, :gaussian, :cornish_fisher."))
end


text:
#####
##### `EDF.Signal` label handling
#####

# This function:
# - ensures the given label is whitespace-stripped, lowercase, and parens-free
# - strips trailing generic EDF references (e.g. "ref", "ref2", etc.)
# - replaces all references with the appropriate name as specified by `canonical_names`
# - replaces `+` with `_plus_` and `/` with `_over_`
# - returns the initial reference name (w/o prefix sign, if present) and the entire label;
#   the initial reference name should match the canonical channel name,
#   otherwise the channel extraction will be rejected.
function _normalize_references(original_label, canonical_names)
    label = replace(_safe_lowercase(original_label), r"\s"=>"")
    label = replace(replace(label, '('=>""), ')'=>"")
    label = replace(label, r"\*$"=>"")
    label = replace(label, '-'=>'…')
    label = replace(label, '+'=>"…+…")
    label = replace(label, '/'=>"…/…")
    m = match(r"^\[(.*)\]$", label)
    if m !== nothing
        label = only(m.captures)
    end
    parts = split(label, '…'; keepempty=false)
    final = findlast(part -> replace(part, r"\d" => "") != "ref", parts)
    parts = parts[1:something(final, 0)]
    isempty(parts) && return ("", "")
    for n in canonical_names
        if n isa Pair
            primary, alternatives = n
            primary = string(primary)
            for alternative in (string(a) for a in alternatives)
                for i in 1:length(parts)
                    if parts[i] == alternative
                        parts[i] = primary
                    end
                end
            end
        end
    end
    recombined = '-'^startswith(original_label, '-') * join(parts, '-')
    recombined = replace(recombined, "-+-"=>"_plus_")
    recombined = replace(recombined, "-/-"=>"_over_")
    return first(parts), recombined
end

_safe_lowercase(c::Char) = isvalid(c) ? lowercase(c) : c

# malformed UTF-8 chars are a choking hazard
_safe_lowercase(s::AbstractString) = map(_safe_lowercase, s)

function match_edf_label(label, signal_names, channel_name, canonical_names)
    label = _safe_lowercase(label)
    for signal_name in signal_names
        # match exact STANDARD (or custom) signal types at beginning of label, ignoring case
        # possibly bracketed by or prepended with `[`, `]`, `,` or whitespace
        # everything after is included in the spec a.k.a. label
        m = match(Regex("[\\s\\[,\\]]*$(signal_name)[\\s,\\]]*\\s+(?<spec>.+)", "i"), label)
        if !isnothing(m)
            label = m[:spec]
        end
        # if signal type does not match, use the entire label
    end
    label = replace(label, r"\s*-\s*" => "-")
    initial, normalized_label = _normalize_references(label, canonical_names)
    initial == channel_name && return normalized_label
    return nothing
end

#####
##### encodings
#####

struct MismatchedSampleRateError <: Exception
    sample_rates
end

function Base.showerror(io::IO, err::MismatchedSampleRateError)
    print(io, """
              found mismatched sample rate between channel encodings: $(err.sample_rates)

              OndaEDF does not currently automatically resolve mismatched sample rates;
              please preprocess your data before attempting `import_edf` so that channels
              of the same signal share a common sample rate.
              """)
end

# I wasn't super confident that the `sample_offset_in_unit` calculation I derived
# had the correctness/symmetry I was hoping it had, so basic algebra time (written
# so that you can try each step out in the REPL with random values):
#
#   res = (pmax - pmin) / (dmax - dmin)
#   pmax - (sample_resolution_in_unit * dmax) ≈ pmin - (sample_resolution_in_unit * dmin)
#   pmax - ((pmax * dmax - pmin * dmax) / (dmax - dmin)) ≈ pmin - ((pmax * dmin - pmin * dmin) / (dmax - dmin))
#   pmax - ((pmax * dmax - pmin * dmax) / (dmax - dmin)) + ((pmax * dmin - pmin * dmin) / (dmax - dmin)) ≈ pmin
#   pmax + ((pmin * dmax - pmax * dmax) / (dmax - dmin)) + ((pmax * dmin - pmin * dmin) / (dmax - dmin)) ≈ pmin
#   pmax + (pmin * dmax - pmax * dmax + pmax * dmin - pmin * dmin) / (dmax - dmin) ≈ pmin
#   pmax + (pmin * dmax + pmax * (-dmax) + pmax * dmin + pmin * (-dmin)) / (dmax - dmin) ≈ pmin
#   pmax + (pmax*(dmin - dmax) + pmin*(dmax - dmin)) / (dmax - dmin) ≈ pmin
#   pmax + pmin + (pmax*(dmin - dmax)/(dmax - dmin)) ≈ pmin
#   pmax + pmin + (-pmax) ≈ pmin
#   pmin ≈ pmin
function edf_signal_encoding(edf_signal_header::EDF.SignalHeader,
                             edf_seconds_per_record)
    dmin, dmax = edf_signal_header.digital_minimum, edf_signal_header.digital_maximum
    pmin, pmax = edf_signal_header.physical_minimum, edf_signal_header.physical_maximum
    sample_resolution_in_unit = (pmax - pmin) / (dmax - dmin)
    sample_offset_in_unit = pmin - (sample_resolution_in_unit * dmin)
    sample_rate = edf_signal_header.samples_per_record / edf_seconds_per_record
    return (sample_resolution_in_unit=Float64(sample_resolution_in_unit),
            sample_offset_in_unit=Float64(sample_offset_in_unit),
            sample_rate=Float64(sample_rate))
end

# TODO: implement a more disciplined widening protocol for `sample_type`
function promote_encodings(encodings; pick_offset=(_ -> 0.0), pick_resolution=minimum)
    sample_type = Int16

    sample_rates = [e.sample_rate for e in encodings]
    if all(==(first(sample_rates)), sample_rates)
        sample_rate = first(sample_rates)
    else
        throw(MismatchedSampleRateError(sample_rates))
    end

    offsets = [e.sample_offset_in_unit for e in encodings]
    if all(==(first(offsets)), offsets)
        sample_offset_in_unit = first(offsets)
    else
        sample_type = Int32
        sample_offset_in_unit = pick_offset(offsets)
    end

    resolutions = [e.sample_resolution_in_unit for e in encodings]
    if all(==(first(resolutions)), resolutions)
        sample_resolution_in_unit = first(resolutions)
    else
        sample_type = Int32
        sample_resolution_in_unit = pick_resolution(resolutions)
    end

    return (sample_type=sample_type,
            sample_offset_in_unit=sample_offset_in_unit,
            sample_resolution_in_unit=sample_resolution_in_unit,
            sample_rate=sample_rate)
end

#####
##### `EDF.Signal`s -> `Onda.Samples`
#####

function extract_channels(edf_signals, channel_matchers)
    extracted_channels = Pair{String,EDF.Signal}[]
    for channel_matcher in channel_matchers
        for edf_signal in edf_signals
            edf_signal isa EDF.Signal || continue
            any(x -> last(x) === edf_signal, extracted_channels) && continue
            channel_name = channel_matcher(edf_signal)
            channel_name === nothing && continue
            push!(extracted_channels, channel_name => edf_signal)
        end
    end
    return extracted_channels
end

"""
    edf_signals_to_samplesinfo(edf, edf_signals, kind, channel_names, samples_per_record; unit_alternatives=STANDARD_UNITS)

Generate a single `Onda.SamplesInfo` for the given collection of `EDF.Signal`s
corresponding to the channels of a single Onda signal.  Sample units are
converted to Onda units and checked for consistency, and a promoted encoding
(resolution, offset, and sample type/rate) is generated.

No conversion of the actual signals is performed at this step.
"""
function edf_signals_to_samplesinfo(edf::EDF.File, edf_signals::Vector{<:EDF.Signal}, kind, channel_names; unit_alternatives=STANDARD_UNITS)
    onda_units = map(s -> edf_to_onda_unit(s.header.physical_dimension, unit_alternatives), edf_signals)
    onda_sample_unit = first(onda_units)

    edf_encodings = map(s -> edf_signal_encoding(s.header, edf.header.seconds_per_record), edf_signals)
    onda_encoding = promote_encodings(edf_encodings)

    info = SamplesInfo(; kind=kind, channels=channel_names,
                       sample_unit=string(onda_sample_unit),
                       sample_resolution_in_unit=onda_encoding.sample_resolution_in_unit,
                       sample_offset_in_unit=onda_encoding.sample_offset_in_unit,
                       sample_type=onda_encoding.sample_type,
                       sample_rate=onda_encoding.sample_rate)
    return info
end

struct SamplesInfoError <: Exception
    msg::String
    cause::Exception
end

function Base.showerror(io::IO, e::SamplesInfoError)
    print(io, "SamplesInfoError: ", e.msg, " caused by: ")
    Base.showerror(io, e.cause)
end

function groupby(f, list)
    d = Dict()
    for v in list
        push!(get!(d, f(v), Vector{eltype(list)}()), v)
    end
    return d
end

"""
    extract_channels_by_label(edf::EDF.File, signal_names, channel_names;
                              unit_alternatives=STANDARD_UNITS, preprocess_labels=(l,t) -> l)

For one or more signal names and one or more channel names,
return a list of `[(infos, errors)...]` where `errors` is a vector of errors that occurred,
and `infos` is a vector of `(si::Onda.SamplesInfo, edf_signals::Vector{EDF.Signal})`,
where `edf_signals` align with `si.channel`. This list can have more than one pair
if channels with the same signal kind/type have different sample rates or
physical units.

`(label, transducer_type)` pairs will be transformed into labels by `preprocess_labels`
(default preprocessor returns the original label). An alternate to `STANDARD_UNITS`
for mapping different spellings of units to their canonical unit name can be passed in
as `unit_alternatives`.

`errors` contains `SamplesInfoError`s thrown if channels corresponding to a signal
were extracted but an error occured while interpreting physical units,
while promoting sample encodings, or otherwise constructing a `SamplesInfo`.

`signal_names` should be an iterable of `String`s naming the signal types to
extract (e.g., `["ecg", "ekg"]`; `["eeg"]`).

`channel_names` should be an iterable of channel specifications, each of which
can be either a `String` giving the generated channel name, or a `Pair` mapping
a canonical name to a list of alternatives that it should be substituted for
(e.g., `"canonical_name" => ["alt1", "alt2", ...]`).

`unit_alternatives` lists standardized unit names and alternatives that map to them.
See `OndaEDF.STANDARD_UNITS` for defaults.

`preprocess_labels(label::String, transducer_type::String)` is applied to raw edf signal header labels
beforehand; defaults to returning `label`.

See `OndaEDF.STANDARD_LABELS` for the labels (`signal_names => channel_names`
`Pair`s) that are used to extract EDF signals by default.

"""
function extract_channels_by_label(edf::EDF.File, signal_names, channel_names;
                                   unit_alternatives=STANDARD_UNITS, preprocess_labels=(l,t) -> l)
    matcher = x -> begin
        # yo I heard you like closures
        # x is either a channel name (string), or a channel_name => alternatives Pair.
        this_channel_name = x isa Pair ? first(x) : x
        return s -> match_edf_label(preprocess_labels(s.header.label, s.header.transducer_type),
                                    signal_names,
                                    this_channel_name,
                                    channel_names)
    end
    edf_channels = extract_channels(edf.signals, (matcher(x) for x in channel_names))

    # place channels with different physical units or sample rates in separate Onda signals
    grouped = groupby(edf_channels) do p
        channel_name, edf_signal = p
        return ((edf_signal.header.samples_per_record / edf.header.seconds_per_record), edf_signal.header.physical_dimension)
    end

    results = map(values(grouped)) do pairs
        # pairs is a vector of `standard_onda_name::String => EDF.Signal` pairs
        edf_channel_names, edf_channels = zip(pairs...)
        try
            edf_channels = collect(edf_channels)
            info = edf_signals_to_samplesinfo(edf, edf_channels, first(signal_names), collect(edf_channel_names))
            return info, edf_channels
        catch e
            # do not throw, but return any errors
            units = [s.header.label => s.header.physical_dimension for s in edf_channels]
            msg ="""Skipping signal: error while processing units and encodings
                    for $(first(signal_names)) signal with units $units"""
            return SamplesInfoError(msg, e)
        end
    end

    return ([info_channels for info_channels in results if !isa(info_channels, Exception)],
            [e for e in results if isa(e, Exception)])
end

#####
##### `import_edf!`
#####

function onda_samples_from_edf_signals(target::Onda.SamplesInfo, edf_signals,
                                       edf_seconds_per_record)
    sample_count = length(first(edf_signals).samples)
    if !all(length(s.samples) == sample_count for s in edf_signals)
        error("mismatched sample counts between `EDF.Signal`s: ", [length(s.samples) for s in edf_signals])
    end
    sample_data = Matrix{sample_type(target)}(undef, length(target.channels), sample_count)
    for (i, edf_signal) in enumerate(edf_signals)
        edf_encoding = edf_signal_encoding(edf_signal.header, edf_seconds_per_record)
        if target.sample_rate != edf_encoding.sample_rate
            throw(MismatchedSampleRateError((target.sample_rate, edf_encoding.sample_rate)))
        end
        if (target.sample_resolution_in_unit != edf_encoding.sample_resolution_in_unit ||
            target.sample_offset_in_unit != edf_encoding.sample_offset_in_unit ||
            sample_type(target) != eltype(edf_signal.samples))
            decoded_samples = Onda.decode(edf_encoding.sample_resolution_in_unit,
                                          edf_encoding.sample_offset_in_unit,
                                          edf_signal.samples)
            encoded_samples = Onda.encode(sample_type(target), target.sample_resolution_in_unit,
                                          target.sample_offset_in_unit, decoded_samples,
                                          missing)
        else
            encoded_samples = edf_signal.samples
        end
        copyto!(view(sample_data, i, :), encoded_samples)
    end
    return Samples(sample_data, target, true)
end

"""
    store_edf_as_onda(path, edf::EDF.File, uuid::UUID=uuid4();
                      custom_extractors=STANDARD_EXTRACTORS, import_annotations::Bool=true,
                      postprocess_samples=identity,
                      signals_prefix="edf", annotations_prefix=signals_prefix)

Convert an EDF.File to `Onda.Samples` and `Onda.Annotation`s, store the samples
in `\$path/samples/`, and write the Onda signals and annotations tables to
`\$path/\$(signals_prefix).onda.signals.arrow` and
`\$path/\$(annotations_prefix).onda.annotations.arrow`.  The default prefix is
"edf", and if a prefix is provided for signals but not annotations both will use
the signals prefix.  The prefixes cannot reference (sub)directories.

Returns `uuid => (signals, annotations)`.

Samples are extracted with [`edf_to_onda_samples`](@ref), and EDF+ annotations are
extracted with [`edf_to_onda_annotations`](@ref) if `import_annotations==true`
(the default).

Collections of `EDF.Signal`s are mapped as channels to `Onda.Signal`s via simple
"extractor" callbacks of the form:

    edf::EDF.File -> (samples_info::Onda.SamplesInfo,
                      edf_signals::Vector{EDF.Signal})

`store_edf_as_onda` automatically uses a variety of default extractors derived
from the EDF standard texts; see `src/standards.jl` and
[`extract_channels_by_label`](@ref) for details. The caller can provide
alternative extractors via the `custom_extractors` keyword argument, and the
[`edf_signals_to_samplesinfo`](@ref) utility can be used to extract a common
`Onda.SamplesInfo` from a collection of EDF.Signals.

`EDF.Signal` labels that are converted into Onda channel names undergo the
following transformations:

- the label's prepended signal type is matched against known types, if present
- the remainder of the label is whitespace-stripped, parens-stripped, and lowercased
- trailing generic EDF references (e.g. "ref", "ref2", etc.) are dropped
- any instance of `+` is replaced with `_plus_` and `/` with `_over_`
- all component names are converted to their "canonical names" when possible
  (e.g. for an EOG matched channel, "eogl", "loc", "lefteye", etc. are converted
  to "left").

`EDF.Signal`s which get extracted into more than one `Onda.Samples` are removed
and an `AmbiguousChannelError` displayed as a warning.

See the OndaEDF README for additional details regarding EDF formatting expectations.
"""
function store_edf_as_onda(path, edf::EDF.File, uuid::UUID=uuid4();
                           custom_extractors=STANDARD_EXTRACTORS, import_annotations::Bool=true,
                           postprocess_samples=identity,
                           signals_prefix="edf", annotations_prefix=signals_prefix)
    EDF.read!(edf)
    file_format = "lpcm.zst"

    signals = Any[]
    edf_samples, diagnostics = edf_to_onda_samples(edf; custom_extractors=custom_extractors)
    for e in diagnostics.errors
        @warn sprint(showerror, e)
    end
    edf_samples = postprocess_samples(edf_samples)
    for samples in edf_samples
        file_path = joinpath(path, "samples", string(uuid, "_", samples.info.kind, ".", file_format))
        signal = rowmerge(store(file_path, file_format, samples, uuid, Second(0)); file_path=string(file_path))
        push!(signals, signal)
    end

    signals_path = joinpath(path, "$(validate_arrow_prefix(signals_prefix)).onda.signals.arrow")
    write_signals(signals_path, signals)
    if import_annotations
        annotations = edf_to_onda_annotations(edf, uuid)
        if !isempty(annotations)
            annotations_path = joinpath(path, "$(validate_arrow_prefix(annotations_prefix)).onda.annotations.arrow")
            write_annotations(annotations_path, annotations)
        else
            @warn "No annotations found in $path"
        end
    else
        annotations = Annotation[]
    end
    return uuid => (signals, annotations)
end

function validate_arrow_prefix(prefix)
    prefix == basename(prefix) || throw(ArgumentError("prefix \"$prefix\" is invalid: cannot contain directory separator"))
    pm = match(r"(.*).onda.(signals|annotations).arrow", prefix)
    if pm !== nothing
        @warn "Extracting prefix \"$(pm.captures[1])\" from provided prefix \"$prefix\""
        prefix = pm.captures[1]
    end
    return prefix
end

# the same EDF.Signal was extracted into multiple `Onda.SampleInfos`.
struct AmbiguousChannelError <: Exception
    summary
end

Base.showerror(io::IO, e::AmbiguousChannelError) = print(io, "AmbiguousChannelError: the same `EDF.Signal` was extracted into multiple `Onda.SampleInfo`s\n", e.summary)

# return a list of tuples describing sample info channels that are ambiguous
# because they were extracted from the same edf_signal
# modifies `infos` and associated `edf_signals::Vector{EDF.Signal}` in place to remove ambiguous signals
function _ambiguous_channels!(extracted, edf_signal)
    ambiguous = []
    for (info, edf_signals) in extracted
        i = findfirst(==(edf_signal), edf_signals)
        isnothing(i) && continue
        channel_name = info.channels[i]
        summary = (info.kind, channel_name, info.sample_unit)
        push!(ambiguous, (info, edf_signals) => summary)
    end
    length(ambiguous) < 2 && return nothing
    for ((info, edf_signals), (_, channel_name, _)) in ambiguous
        filter!(!=(channel_name), info.channels)
        filter!(!=(edf_signal), edf_signals)
    end
    return map(last, ambiguous)
end

"""
    edf_to_onda_samples(edf::EDF.File; custom_extractors=())

Read signals from an `EDF.File` into a vector of `Onda.Samples`,
which are returned along with a NamedTuple with diagnostic information
(the same info returned by [`edf_header_to_onda_samples_info`](@ref)).

Collections of `EDF.Signal`s are mapped as channels to `Onda.Signal`s via simple
"extractor" callbacks of the form:

    edf::EDF.File -> (samples_info::Onda.SamplesInfo,
                      edf_signals::Vector{EDF.Signal})

`edf_to_onda_samples` automatically uses a variety of default extractors derived from
the EDF standard texts; see `src/standards.jl` for details. The caller can also
provide additional extractors via the `custom_extractors` keyword argument.

`EDF.Signal` labels that are converted into Onda channel names undergo the
following transformations:

- the label is whitespace-stripped, parens-stripped, and lowercased
- trailing generic EDF references (e.g. "ref", "ref2", etc.) are dropped
- any instance of `+` is replaced with `_plus_` and `/` with `_over_`
- all component names are converted to their "canonical names" when possible
  (e.g. "m1" in an EEG-matched channel name will be converted to "a1").

See the OndaEDF README for additional details regarding EDF formatting expectations.
"""
function edf_to_onda_samples(edf::EDF.File; custom_extractors=STANDARD_EXTRACTORS)
    EDF.read!(edf)
    info_map, nt = edf_header_to_onda_samples_info(edf; custom_extractors=custom_extractors)
    edf_samples = [onda_samples_from_edf_signals(info,
                                                 edf_signals,
                                                 edf.header.seconds_per_record)
                   for (info, edf_signals) in info_map if !isempty(info.channels)]
    return edf_samples, nt
end

"""
    edf_header_to_onda_samples_info(edf::EDF.File; custom_extractors=STANDARD_EXTRACTORS)

Read edf header, return a mapping from `Onda.SamplesInfo`s to the vector
of `EDF.Signals` it was extracted from, along with a `NamedTuple` containing
diagnostic information in fields:
- `header_map`: a vector of corresponding `Onda.SamplesInfo`, `Vector{EDF.SignalHeader}` pairs.
- `unextracted_edf_headers`: a vector of EDF signal headers that could not be extracted.
- `errors`: vector of exceptions thrown while attempting to extract header info.

The `NamedTuple` can be pretty-printed and used to compare outputs with `SamplesInfo`s
previously extracted from the same data, for testing purposes.

`EDF.read!` does not get called, this function will work
with only the first few bytes--30k should be enough--of the edf file,
which is convenient for developping custom extractors for a dataset without
reading and converting all the samples data.
"""
function edf_header_to_onda_samples_info(edf::EDF.File; custom_extractors=STANDARD_EXTRACTORS)
    info_map, errors = try
        matched = []
        errors = Exception[]
        for extractor in custom_extractors
            extracteds, errs = extractor(edf)
            append!(errors, errs)
            for extracted in extracteds
                push!(matched, extracted)
            end
        end
        # each edf_signal should get extracted into at most one SamplesInfo, otherwise it is ambiguous
        matched_signals = Set(Iterators.flatten(map(last, matched)))
        ambiguous_edf_signals = []
        for edf_signal in matched_signals
            ambiguous_channels = _ambiguous_channels!(matched, edf_signal)
            if !isnothing(ambiguous_channels)
                edf_signal_summary = (edf_signal.header.label,
                                      edf_signal.header.transducer_type,
                                      edf_signal.header.physical_dimension)
                push!(errors, AmbiguousChannelError(edf_signal_summary => ambiguous_channels))
            end
        end
        matched = [(info, edf_signals) for (info, edf_signals) in matched if !isempty(info.channels)]
        sort!(errors; by=string)
        matched, errors
    catch e
        [], [e]
    end
    matched_edf_headers = mapreduce(last, union, info_map; init=[])
    unextracted = [s.header for s in edf.signals if isa(s, EDF.Signal) && s ∉ matched_edf_headers]
    header_map = [info => [s.header for s in edf_signals if isa(s, EDF.Signal)] for (info, edf_signals) in info_map]
    return info_map, (header_map=header_map,
                      unextracted_edf_headers=unextracted,
                      errors=errors)
end

"""
    edf_to_onda_annotations(edf::EDF.File, uuid::UUID)

Extract EDF+ annotations from an `EDF.File` for recording with ID `uuid` and
return them as a vector of `Onda.Annotation`s.  Each returned annotation has
a  `value` field that contains the string value of the corresponding EDF+
annotation.

If no EDF+ annotations are found in `edf`, then an empty `Vector{Annotation}` is
returned.
"""
function edf_to_onda_annotations(edf::EDF.File, uuid::UUID)
    EDF.read!(edf)
    annotations = Annotation[]
    for annotation_signal in edf.signals
        annotation_signal isa EDF.AnnotationsSignal || continue
        for record in annotation_signal.records
            for tal in record
                start_nanosecond = Nanosecond(round(Int, 1e9 * tal.onset_in_seconds))
                if tal.duration_in_seconds === nothing
                    stop_nanosecond = start_nanosecond
                else
                    stop_nanosecond = start_nanosecond + Nanosecond(round(Int, 1e9 * tal.duration_in_seconds))
                end
                for annotation_string in tal.annotations
                    isempty(annotation_string) && continue
                    annotation = Annotation(; recording=uuid, id=uuid4(),
                                            span=TimeSpan(start_nanosecond, stop_nanosecond),
                                            value=annotation_string)
                    push!(annotations, annotation)
                end
            end
        end
    end
    return annotations
end


text:
@testset "easy test" begin
    apex = SVector(15.0,-7, 0.1)
    axis = SVector(0,0,-1.)
    opang = deg2rad(20)
    h = 10.
    sizet = (10,10)

    cps, cns = RANSACBenchmark.samplecone(apex, axis, opang, h, sizet)

    rrr = [25, 19, 97]
    hicone = RANSAC.fit3pointcone(cps[rrr], cns[rrr])
    rrr2 = [25, 97, 19]
    hicone2 = RANSAC.fit3pointcone(cps[rrr2], cns[rrr2])

    @test hicone.apex-apex |> norm < 1e-13
    @test hicone2.apex-apex |> norm < 1e-13

    @test hicone.axis+[0,0,1] |> norm < 1e-13
    @test hicone2.axis+[0,0,1] |> norm < 1e-13

    @test isapprox(hicone.opang, opang)
    @test isapprox(hicone2.opang, opang)

    p = RANSACParameters{Float64}()
    hcone = RANSAC.fitcone(cps[rrr], cns[rrr], p)
    hcone2 = RANSAC.fitcone(cps[rrr2], cns[rrr2], p)

    @test hcone isa FittedCone
    @test hcone2 isa FittedCone
    @test hcone.outwards == true
    @test hcone2.outwards == true
end


text:
using EmpiricalBayes
using Test

lfsr_target =  PosteriorTarget(LFSRNumerator(2.0))
mean_target = PosteriorTarget(PosteriorMeanNumerator(2.0))

EmpiricalBayes.pretty_label(lfsr_target)
EmpiricalBayes.pretty_label(mean_target)


text:
@with_kw struct Data_il
  protein::String
  MMP::Float64

  cation::String
  anion::String

  MM::Float64
  c::Float64

  cat_atomtypes::String = "$(cation)_atomtypes_VSIL.itp"
  cat_par::String = "$(cation)_VSIL.itp"

  an_atomtypes::String = "$(anion)_atomtypes_VSIL.itp"
  an_par::String = "$(anion)_VSIL.itp"
end

@with_kw struct Data_2il
  protein::String
  MMP::Float64

  cation::String
  anion1::String
  anion2::String

  MM1::Float64
  MM2::Float64
  c::Float64

  cat_atomtypes::String = "$(cation)_atomtypes_VSIL.itp"
  cat_par::String = "$(cation)_VSIL.itp"

  an1_atomtypes::String = "$(anion1)_atomtypes_VSIL.itp"
  an1_par::String = "$(anion1)_VSIL.itp"

  an2_atomtypes::String = "$(anion2)_atomtypes_VSIL.itp"
  an2_par::String = "$(anion2)_VSIL.itp"

end

@with_kw struct Data_elec

  cation::String
  anion::String

  MM::Float64
  c::Float64

  cat_atomtypes::String = "$(cation)_atomtypes_VSIL.itp"
  cat_par::String = "$(cation)_VSIL.itp"

  an_atomtypes::String = "$(anion)_atomtypes_VSIL.itp"
  an_par::String = "$(anion)_VSIL.itp"

end

@with_kw struct Data_twoelec

  cation::String
  anion1::String
  anion2::String

  MM1::Float64
  MM2::Float64

  c::Float64

  cat_atomtypes::String = "$(cation)_atomtypes_VSIL.itp"
  cat_par::String = "$(cation)_VSIL.itp"

  an1_atomtypes::String = "$(anion1)_atomtypes_VSIL.itp"
  an1_par::String = "$(anion1)_VSIL.itp"

  an2_atomtypes::String = "$(anion2)_atomtypes_VSIL.itp"
  an2_par::String = "$(anion2)_VSIL.itp"

end


text:
module Acorn

import REPL
using Printf

include("terminal.jl")

# Configuration submodule
include("EditorConfig.jl")
using .EditorConfig

include("row.jl")
include("cmds/Command.jl")
include("editor.jl")

# Load commands
include("cmds/open.jl")
include("cmds/save.jl")
include("cmds/quit.jl")
include("cmds/find.jl")
include("cmds/help.jl")
include("cmds/bind.jl")
include("cmds/set.jl")
include("cmds/echo.jl")

function acorn(filename::String; rel::Bool=true)
    ed = Editor()

    editorOpen(ed, filename)

    setStatusMessage(ed, "HELP: ctrl-p: command mode | ctrl-q: quit | ctrl-s: save")

    REPL.Terminals.raw!(ed.term, true)

    try
        while !ed.quit
            refreshScreen(ed)
            processKeypress(ed)
        end
    catch ex
        editorQuit(ed, force=true)
        rethrow(ex) # Don't reset stacktrace
    end

    REPL.Terminals.raw!(ed.term, false)

    return nothing
end

function acorn()
    # If a file is given, open it
    if length(ARGS) > 0
        filename = ARGS[1]
        acorn(filename, rel=false)
    else
        println("No filename detected.")
    end
end

export
    acorn

end # module


text:
export Header,
SegyToSeis,
SeisToSegy,
SeisCopy,
SeisHeaderInfo,
SeisRead,
SeisReadHeaders,
SeisRemove,
SeisWrite,
SeisWriteHeaders,
ExtractSegyHeader,
ReadSegyHeader,
SegyHeaderInfo
include("Header.jl")
include("SegyStruct.jl")
include("SegyToSeis.jl")
include("SeisToSegy.jl")
include("SeisCopy.jl")
include("SeisHeaderInfo.jl")
include("SeisRead.jl")
include("SeisReadHeaders.jl")
include("SeisRemove.jl")
include("SeisWrite.jl")
include("SeisWriteHeaders.jl")
include("ExtractSegyHeader.jl")
include("ReadSegyHeader.jl")
include("SegyHeaderInfo.jl")


text:
function AbstractTimeScales()
    return AbstractTimeScales(())
end

function create_besselian_epoch(obj::AbstractTimeScales, arg0::jdouble)
    return jcall(obj, "createBesselianEpoch", AbsoluteDate, (jdouble,), arg0)
end

function create_julian_epoch(obj::AbstractTimeScales, arg0::jdouble)
    return jcall(obj, "createJulianEpoch", AbsoluteDate, (jdouble,), arg0)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_bdt(obj::TimeScales)
    return jcall(obj, "getBDT", BDTScale, ())
end

function get_beidou_epoch(obj::AbstractTimeScales)
    return jcall(obj, "getBeidouEpoch", AbsoluteDate, ())
end

function get_ccsds_epoch(obj::AbstractTimeScales)
    return jcall(obj, "getCcsdsEpoch", AbsoluteDate, ())
end

function get_fifties_epoch(obj::AbstractTimeScales)
    return jcall(obj, "getFiftiesEpoch", AbsoluteDate, ())
end

function get_future_infinity(obj::AbstractTimeScales)
    return jcall(obj, "getFutureInfinity", AbsoluteDate, ())
end

function get_galileo_epoch(obj::AbstractTimeScales)
    return jcall(obj, "getGalileoEpoch", AbsoluteDate, ())
end

function get_glonass(obj::TimeScales)
    return jcall(obj, "getGLONASS", GLONASSScale, ())
end

function get_glonass_epoch(obj::AbstractTimeScales)
    return jcall(obj, "getGlonassEpoch", AbsoluteDate, ())
end

function get_gmst(obj::AbstractTimeScales, arg0::IERSConventions, arg1::jboolean)
    return jcall(obj, "getGMST", GMSTScale, (IERSConventions, jboolean), arg0, arg1)
end

function get_gps(obj::TimeScales)
    return jcall(obj, "getGPS", GPSScale, ())
end

function get_gps_epoch(obj::AbstractTimeScales)
    return jcall(obj, "getGpsEpoch", AbsoluteDate, ())
end

function get_gst(obj::TimeScales)
    return jcall(obj, "getGST", GalileoScale, ())
end

function get_irnss(obj::TimeScales)
    return jcall(obj, "getIRNSS", IRNSSScale, ())
end

function get_irnss_epoch(obj::AbstractTimeScales)
    return jcall(obj, "getIrnssEpoch", AbsoluteDate, ())
end

function get_j2000_epoch(obj::AbstractTimeScales)
    return jcall(obj, "getJ2000Epoch", AbsoluteDate, ())
end

function get_java_epoch(obj::AbstractTimeScales)
    return jcall(obj, "getJavaEpoch", AbsoluteDate, ())
end

function get_julian_epoch(obj::AbstractTimeScales)
    return jcall(obj, "getJulianEpoch", AbsoluteDate, ())
end

function get_modified_julian_epoch(obj::AbstractTimeScales)
    return jcall(obj, "getModifiedJulianEpoch", AbsoluteDate, ())
end

function get_past_infinity(obj::AbstractTimeScales)
    return jcall(obj, "getPastInfinity", AbsoluteDate, ())
end

function get_qzss(obj::TimeScales)
    return jcall(obj, "getQZSS", QZSSScale, ())
end

function get_qzss_epoch(obj::AbstractTimeScales)
    return jcall(obj, "getQzssEpoch", AbsoluteDate, ())
end

function get_tai(obj::TimeScales)
    return jcall(obj, "getTAI", TAIScale, ())
end

function get_tcb(obj::TimeScales)
    return jcall(obj, "getTCB", TCBScale, ())
end

function get_tcg(obj::TimeScales)
    return jcall(obj, "getTCG", TCGScale, ())
end

function get_tdb(obj::TimeScales)
    return jcall(obj, "getTDB", TDBScale, ())
end

function get_tt(obj::TimeScales)
    return jcall(obj, "getTT", TTScale, ())
end

function get_ut1(obj::AbstractTimeScales, arg0::IERSConventions, arg1::jboolean)
    return jcall(obj, "getUT1", UT1Scale, (IERSConventions, jboolean), arg0, arg1)
end

function get_utc(obj::TimeScales)
    return jcall(obj, "getUTC", UTCScale, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end


text:
using ITensors, ITensorsGPU

function my_polar( T::ITensor, inds...; kwargs... )
    U, S, V, spec, u, v = svd(T, inds; kwargs...)
    replaceind!(U, u, v)
    return U*V
end

function initQs( A::fPEPS, col::Int, next_col::Int; kwargs...)
    Ny, Nx = size(A)
    maxdim::Int = get(kwargs, :maxdim, 1)
    Q         = MPO(deepcopy(A[:, col]), 0, Ny+1)
    prev_col  = next_col > col ? col - 1 : col + 1
    A_r_inds  = [commonind(A[row, col], A[row, next_col]) for row in 1:Ny]
    QR_inds   = [Index(ITensors.dim(A_r_inds[row]), "Site,QR,c$col,r$row") for row in 1:Ny]
    A_up_inds = [commonind(A[row, col], A[row+1, col]) for row in 1:Ny-1]
    Q_up_inds = [Index(ITensors.dim(A_up_inds[row]), "Link,u,Qup$row") for row in 1:Ny-1]
    next_col_inds = [commonind(A[row, col], A[row, next_col]) for row in 1:Ny]
    prev_col_inds = 0 < prev_col < Nx ? [commonind(A[row, col], A[row, prev_col]) for row in 1:Ny] : Vector{Index}(undef, Ny)
    for row in 1:Ny
        row < Ny && replaceind!(Q[row], A_up_inds[row], Q_up_inds[row])
        row > 1  && replaceind!(Q[row], A_up_inds[row-1], Q_up_inds[row-1])
        replaceind!(Q[row], next_col_inds[row], QR_inds[row])
    end
    return Q, QR_inds, next_col_inds
end

function initQs( A::fPEPS, Q::MPO, col::Int, next_col::Int; kwargs...)
    Ny, Nx = size(A)
    maxdim::Int = get(kwargs, :maxdim, 1)
    Q         = MPO(data(Q), 0, Ny+1)
    prev_col  = next_col > col ? col - 1 : col + 1
    A_r_inds  = [commonind(A[row, col], A[row, next_col]) for row in 1:Ny]
    QR_inds   = [Index(ITensors.dim(A_r_inds[row]), "Site,QR,c$col,r$row") for row in 1:Ny]
    A_up_inds = [commonind(A[row, col], A[row+1, col]) for row in 1:Ny-1]
    Q_up_inds = [Index(ITensors.dim(A_up_inds[row]), "Link,u,Qup$row") for row in 1:Ny-1]
    next_col_inds = [commonind(A[row, col], A[row, next_col]) for row in 1:Ny]
    prev_col_inds = 0 < prev_col < Nx ? [commonind(A[row, col], A[row, prev_col]) for row in 1:Ny] : Vector{Index}(undef, Ny)
    for row in 1:Ny
        #row < Ny && replaceind!(Q[row], A_up_inds[row], Q_up_inds[row])
        #row > 1  && replaceind!(Q[row], A_up_inds[row-1], Q_up_inds[row-1])
        #replaceind!(Q[row], next_col_inds[row], QR_inds[row])
        QR_inds[row] = firstind(Q[row], "QR,Site")
        if ((row == 1 || row == Ny) && ndims(Q[row]) == 4) || ndims(Q[row]) == 5
            prev_ind = firstind(Q[row], "Link,r")
            prev_ind_A = firstind(A[row, col], tags(prev_ind))
            replaceind!(Q[row], prev_ind, prev_ind_A)
        end
    end
    return Q, QR_inds, next_col_inds
end

function initQsHorizontal( A::fPEPS, col::Int, next_col::Int; kwargs...)
    maxdim::Int = get(kwargs, :maxdim, 1)
    Ny        = length(A)
    Q         = MPO(deepcopy(A), 0, Ny+1)
    prev_col  = next_col > col ? col - 1 : col + 1
    A_r_inds  = [commonind(A[row, col], A[row, next_col]) for row in 1:Ny]
    QR_inds   = [Index(ITensors.dim(A_r_inds[row]), "Site,QR,c$col,r$row") for row in 1:Ny]
    A_up_inds = [commonind(A[row, col], A[row+1, col]) for row in 1:Ny-1]
    Q_up_inds = [Index(ITensors.dim(A_up_inds[row]), "Link,u,Qup$row") for row in 1:Ny-1]
    next_col_inds = [commonind(A[row, col], A[row, next_col]) for row in 1:Ny]
    prev_col_inds = 0 < prev_col < Nx ? [commonind(A[row, col], A[row, prev_col]) for row in 1:Ny] : Vector{Index}(undef, Ny)
    for row in 1:Ny
        row < Ny && replaceind!(Q[row], A_up_inds[row], Q_up_inds[row])
        row > 1  && replaceind!(Q[row], A_up_inds[row-1], Q_up_inds[row-1])
        replaceind!(Q[row], next_col_inds[row], QR_inds[row])
    end
    return Q, QR_inds, next_col_inds
end

function buildEnvs(A::fPEPS, Q::MPO, next_col_inds, QR_inds, col)
    Ny, Ny = size(A)
    thisTerm  = Vector{ITensor}(undef, Ny)
    thisfTerm = Vector{ITensor}(undef, Ny)
    for row in 1:Ny
        Ap = dag(deepcopy(A[row, col]))'
        Ap = setprime(Ap, 0, next_col_inds[row]')
        Qp = dag(deepcopy(Q[row]))'
        Qp = setprime(Qp, 0, QR_inds[row]')
        thisTerm[row]  = A[row, col] * Ap * Qp
        thisfTerm[row] = thisTerm[row] * Q[row]
    end
    Envs = thisTerm
    fF   = cumprod(thisfTerm)
    rF   = reverse(cumprod(reverse(thisfTerm)))
    for row in 1:Ny
        if row > 1
            Envs[row] *= fF[row - 1]
        end
        if row < Ny
            Envs[row] *= rF[row + 1]
        end
    end
    return Envs
end

function buildEnvsHorizontal(A::Vector{ITensor}, Q::MPO, next_col_inds, QR_inds, col, ncol)
    Ny = length(A)
    thisTerm  = Vector{ITensor}(undef, Ny)
    thisfTerm = Vector{ITensor}(undef, Ny)
    for row in 1:Ny
        Ap = dag(deepcopy(A[row]))'
        Ap = setprime(Ap, 0, next_col_inds[row]')
        Ap = setprime(Ap, 0, firstind(A[row], "Site,c$ncol,r$row")')
        Qp = dag(deepcopy(Q[row]))'
        Qp = setprime(Qp, 0, QR_inds[row]')
        thisTerm[row]  = A[row] * Ap * Qp
        thisfTerm[row] = thisTerm[row] * Q[row]
    end
    Envs = thisTerm
    fF   = cumprod(thisfTerm)
    rF   = reverse(cumprod(reverse(thisfTerm)))
    for row in 1:Ny
        if row > 1
            Envs[row] *= fF[row - 1]
        end
        if row < Ny
            Envs[row] *= rF[row + 1]
        end
    end
    return Envs
end

function buildEnvsHorizontalVariant(A::Vector{ITensor}, Q::MPO, next_col_inds, QR_inds, col, ncol)
    Ny    = length(A)
    is_cu = is_gpu(A)
    Envs  = Vector{ITensor}(undef, Ny)
    Env   = is_cu ? cuITensor(1.0) : ITensor(1.0)
    for row in reverse(1:Ny)
        Ap = dag(deepcopy(A[row]))'
        Ap = setprime(Ap, 0, next_col_inds[row]')
        Ap = setprime(Ap, 0, firstind(A[row], "Site,c$ncol,r$row")')
        Qp = dag(deepcopy(Q[row]))'
        Qp = setprime(Qp, 0, QR_inds[row]')
        Env *= A[row]
        Env *= Ap
        Env *= Qp
        if row > 1
            Env *= Q[row]
        end
        Envs[row] = Env
    end
    return Envs
end

function compute_overlap(Ampo, Q, R, is_cu)
    Ny = length(Q)
    aqr_overlap = is_cu ? cuITensor(1.0) : ITensor(1.0)
    a_norm      = is_cu ? cuITensor(1.0) : ITensor(1.0)
    cis         = [commonind(Ampo[row], Ampo[row+1]) for row in 1:Ny-1]
    Ampo_       = deepcopy(Ampo)
    for row in 1:Ny
        if row < Ny
            Ampo_[row] = prime(Ampo_[row], cis[row])
        end
        if row > 1
            Ampo_[row] = prime(Ampo_[row], cis[row-1])
        end
        aqr_overlap *= Ampo[row] * Q[row] * R[row]
        a_norm      *= dag(Ampo_[row]) * Ampo[row]
    end
    ratio = abs(scalar(collect(aqr_overlap)))/abs(scalar(collect(a_norm)))
    return ratio
end

function orthogonalize_Q!(Envs, Q, A, QR_inds, next_col_inds, dummy_nexts, col, two_sided; kwargs...)
    Ny    = length(Q)
    Ampo  = MPO(Ny)
    cmb_l = Vector{ITensor}(undef, Ny)
    for row in reverse(1:Ny)
        if row < Ny
            Q_ = my_polar(Envs[row], QR_inds[row], commonind(Q[row], Q[row+1]); kwargs...)
            Q[row] = deepcopy(noprime(Q_))
        else
            Q_ = my_polar(Envs[row], QR_inds[row]; kwargs...)
            Q[row] = deepcopy(noprime(Q_))
        end
        AQinds     = IndexSet(firstind(A[row, col], "Site"))
        if two_sided
            AQinds = IndexSet(AQinds..., commonind(inds(A[row, col], "Link"), inds(Q[row], "Link"))) # prev_col_ind
        end
        cmb_l[row] = combiner(AQinds, tags="Site,AQ,r$row")
        Ampo[row]  = A[row, col] * cmb_l[row]
        Ampo[row]  = replaceind!(Ampo[row], next_col_inds[row], dummy_nexts[row])
        Q[row]    *= cmb_l[row]
    end
    return Q, Ampo, cmb_l
end

function orthogonalize_QHorizontal!(Envs, Q, A, Af, QR_inds, next_col_inds, dummy_nexts, col, ncol, two_sided; kwargs...)
    Ny, Ny = size(Af)
    Ampo  = MPO(Ny)
    cmb_l = Vector{ITensor}(undef, Ny)
    for row in 1:Ny
        if row < Ny
            Q_ = my_polar(Envs[row], QR_inds[row], commonind(Q[row], Q[row+1]); kwargs...)
            Q[row] = deepcopy(noprime(Q_))
        else
            Q_ = my_polar(Envs[row], QR_inds[row]; kwargs...)
            Q[row] = deepcopy(noprime(Q_))
        end
        AQinds     = IndexSet(firstind(A[row], "Site,c$col"))
        if two_sided
            AQinds = IndexSet(AQinds..., commonind(inds(A[row], "Link"), inds(Q[row], "Link"))) # prev_col_ind
        end
        cmb_l[row] = combiner(AQinds, tags="Site,AQ,r$row")
        Ampo[row]  = A[row] * cmb_l[row]
        Ampo[row] *= dummy_nexts[row]
        Q[row]    *= cmb_l[row]
    end
    return Q, Ampo, cmb_l
end

function orthogonalize_QHorizontalRebuild!(Envs, Q, A, Af, QR_inds, next_col_inds, dummy_nexts, col, ncol, two_sided; kwargs...)
    Ny, Ny = size(Af)
    is_cu  = is_gpu(Af)
    Ampo   = MPO(Ny)
    cmb_l  = Vector{ITensor}(undef, Ny)
    for row in 1:Ny
        if row < Ny
            Q_ = my_polar(Envs[row], QR_inds[row], commonind(Q[row], Q[row+1]); kwargs...)
            Q[row] = deepcopy(noprime(Q_))
        else
            Q_ = my_polar(Envs[row], QR_inds[row]; kwargs...)
            Q[row] = deepcopy(noprime(Q_))
        end
        Ap = dag(deepcopy(A[row]))'
        Ap = setprime(Ap, 0, next_col_inds[row]')
        Ap = setprime(Ap, 0, firstind(A[row], "Site,c$ncol,r$row")')
        Qp = dag(deepcopy(Q[row]))'
        Qp = setprime(Qp, 0, QR_inds[row]')
        thisTerm  = is_cu ? cuITensor(1.0) : ITensor(1.0)
        if row > 1
            thisTerm *= Envs[row-1]
        end
        thisTerm = A[row] * Ap * Qp * Q[row]
        Envs[row] = row > 1 ? thisTerm*Envs[row-1] : thisTerm
        if row < Ny
            Ap = dag(deepcopy(A[row+1]))'
            Ap = setprime(Ap, 0, next_col_inds[row+1]')
            Ap = setprime(Ap, 0, firstind(A[row+1], "Site,c$ncol,r$(row+1)")')
            Qp = dag(deepcopy(Q[row+1]))'
            Qp = setprime(Qp, 0, QR_inds[row+1]')
            thisTerm  = Envs[row]
            thisTerm *= A[row+1] * Ap * Qp
            Envs[row+1] = thisTerm
            if row + 1 < Ny
                Envs[row+1] *= Envs[row+2]
            end
        end
        AQinds     = IndexSet(firstind(A[row], "Site,c$col"))
        if two_sided
            AQinds = IndexSet(AQinds..., commonind(inds(A[row], "Link"), inds(Q[row], "Link"))) # prev_col_ind
        end
        cmb_l[row] = combiner(AQinds, tags="Site,AQ,r$row")
        Ampo[row]  = A[row] * cmb_l[row]
        Ampo[row] *= dummy_nexts[row]
        Q[row]    *= cmb_l[row]
    end
    return Q, Ampo, cmb_l
end

function gaugeQR(A::fPEPS, col::Int, side::Symbol; kwargs...)
    overlap_cutoff::Real = get(kwargs, :overlap_cutoff, 1e-4)
    max_gauge_iter::Int  = get(kwargs, :max_gauge_iter, 100)
    Ny, Nx   = size(A)
    prev_col_inds = Vector{Index}(undef, Ny)
    is_cu    = is_gpu(A)
    next_col = side == :left ? col - 1 : col + 1
    prev_col = side == :left ? col + 1 : col - 1
    Q, QR_inds, next_col_inds = initQs(A, col, next_col; kwargs...)
    left_edge     = col == 1
    right_edge    = col == Nx
    ratio_history = Vector{Float64}()
    ratio         = 0.0
    best_overlap  = 0.0
    best_Q        = MPO(Ny)
    best_R        = MPO(Ny)
    iter          = 1
    dummy_nexts   = [Index(ITensors.dim(next_col_inds[row]), "DM,Site,r$row") for row in 1:Ny]
    iter          = 0
    while best_overlap < overlap_cutoff
        @timeit "build envs" begin
            Envs = buildEnvs(A, Q, next_col_inds, QR_inds, col)
        end
        @timeit "polar decomp" begin
            two_sided      = (side == :left && !right_edge) || (side == :right && !left_edge)
            Q, Ampo, cmb_l = orthogonalize_Q!(Envs, Q, A, QR_inds, next_col_inds, dummy_nexts, col, two_sided; kwargs...)
        end
        @timeit "Q*A -> R" begin
            R = multMPO(dag(Q), Ampo; kwargs...)
        end
        @timeit "compute overlap" begin
            ratio = compute_overlap(Ampo, Q, R, is_cu)
        end
        for row in 1:Ny
            Q[row]      *= cmb_l[row]
        end
        push!(ratio_history, ratio)
        if ratio > best_overlap || iter == 0
            best_Q = deepcopy(Q)
            best_R = deepcopy(R)
            best_overlap = ratio
        end
        ratio > overlap_cutoff && break
        iter += 1
        iter > max_gauge_iter && break
        if (iter > 10 && best_overlap < 0.5) || (iter > 20 && mod(iter, 20) == 0)
            Q_, QR_inds_, next_col_inds_ = initQs(A, col, next_col; kwargs...)
            for row in 1:Ny
                if row < Ny
                    old_u = commonind(Q[row], Q[row+1])
                    new_u = commonind(Q_[row], Q_[row+1])
                    replaceind!(Q_[row], new_u, old_u)
                    replaceind!(Q_[row+1], new_u, old_u)
                end
                replaceind!(Q_[row], QR_inds_[row], QR_inds[row])
            end
            for row in 1:Ny
                Q[row]  = Q_[row]
                salt    = is_cu ? cuITensor(randomITensor(inds(Q[row])))/100.0 : randomITensor(inds(Q[row]))/100.0
                salt_d  = ratio < 0.5 ? norm(salt) : 10.0*norm(salt)
                Q[row] += salt/salt_d
                Q[row] /= sqrt(norm(Q[row]))
            end
        end
    end
    println( "\tbest overlap: ", best_overlap)#, "\nratio_history: $ratio_history")
    return best_Q, best_R, next_col_inds, QR_inds, dummy_nexts
end

function gaugeQRHorizontal(A::fPEPS, Aj::Vector{ITensor}, col::Int, ncol::Int, side::Symbol; kwargs...)
    overlap_cutoff::Real = get(kwargs, :overlap_cutoff, 1e-4)
    max_gauge_iter::Int  = get(kwargs, :max_gauge_iter, 100)
    Ny, Nx = size(A)
    is_cu  = is_gpu(A)
    prev_col_inds = Vector{Index}(undef, Ny)
    next_col = side == :left ? col - 1 : col + 1
    prev_col = side == :left ? col + 1 : col - 1
    #Q, QR_inds, next_col_inds = initQs(A, col, ncol; kwargs...)
    Q = MPO()
    QR_inds = Vector{Index}(undef, Ny)
    next_col_inds = Vector{Index}(undef, Ny)
    #=if haskey(q_dict, col=>side)
        Q = q_dict[col=>side]
        Q, QR_inds, next_col_inds = initQs(A, Q, col, ncol; kwargs...)
    else=#
    old_Q, old_R, old_ncis, old_QRinds, old_dns = gaugeQR(A, col, side; max_gauge_iter=10, overlap_cutoff=0.999)
    Q = old_Q
    QR_inds = old_QRinds
    next_col_inds = old_ncis
    #end
    left_edge     = col == 1
    right_edge    = col == Nx
    ratio_history = Vector{Float64}()
    ratio         = 0.0
    best_overlap  = 0.0
    best_Q        = MPO(Ny)
    best_R        = MPO(Ny)
    iter          = 1
    dummy_nexts   = Vector{ITensor}(undef, Ny) #[Index(ITensors.dim(next_col_inds[row]), "DM,Site,r$row") for row in 1:Ny]
    for row in 1:Ny
        cis = IndexSet(firstind(Aj[row], "Site,c$ncol"))
        if ncol > col && ncol <= Nx - 1
            anci = commonindex(A[row, ncol], A[row, ncol+1])
            next_col_inds[row] = anci
            cis = IndexSet(cis..., anci)
        elseif ncol < col && ncol >= 2
            anci = commonindex(A[row, ncol], A[row, ncol-1])
            next_col_inds[row] = anci
            cis = IndexSet(cis..., anci)
        end
        dummy_nexts[row] = combiner(cis, tags="DM,Site,r$row")
    end
    iter          = 0
    downward      = 0
    below         = 0
    while best_overlap < overlap_cutoff
        @timeit "build envs" begin
            #Envs = buildEnvsHorizontal(Aj, Q, next_col_inds, QR_inds, col, ncol)
            Envs = buildEnvsHorizontalVariant(Aj, Q, next_col_inds, QR_inds, col, ncol)
        end
        @timeit "polar decomp" begin
            two_sided      = (side == :left && !right_edge) || (side == :right && !left_edge)
            Q, Ampo, cmb_l = orthogonalize_QHorizontalRebuild!(Envs, Q, Aj, A, QR_inds, next_col_inds, dummy_nexts, col, ncol, two_sided; kwargs...)
        end
        @timeit "Q*A -> R" begin
            R = multMPO(dag(Q), Ampo; kwargs...)
        end
        @timeit "compute overlap" begin
            ratio = compute_overlap(Ampo, Q, R, is_cu)
        end
        for row in 1:Ny
            Q[row]      *= cmb_l[row]
        end
        push!(ratio_history, ratio)
        if iter > 1
            if ratio < ratio_history[end-1]
                downward += 1
            else
                downward = 0
            end
            if ratio < 0.5
                below += 1
            else
                below = 0
            end
        end
        if ratio > best_overlap || iter == 0
            best_Q = deepcopy(Q)
            best_R = deepcopy(R)
            best_overlap = ratio
            #q_dict[col=>side] = Q
        end
        ratio > overlap_cutoff && break
        iter += 1
        iter > max_gauge_iter && break
        #if (iter > 10 && best_overlap < 0.5) || (iter > 20 && mod(iter, 40) == 0)
        if downward > 20 || below > 20
            Q_, QR_inds_, next_col_inds_ = initQs(A, col, next_col; kwargs...)
            for row in 1:Ny
                if row < Ny
                    old_u = commonind(Q[row], Q[row+1])
                    new_u = commonind(Q_[row], Q_[row+1])
                    replaceind!(Q_[row], new_u, old_u)
                    replaceind!(Q_[row+1], new_u, old_u)
                end
                replaceind!(Q_[row], QR_inds_[row], QR_inds[row])
            end
            for row in 1:Ny
                Q[row]  = Q_[row]
                salt    = is_cu ? cuITensor(randomITensor(inds(Q[row])))/100.0 : randomITensor(inds(Q[row]))/100.0
                salt_d  = ratio < 0.5 ? norm(salt) : 10.0*norm(salt)
                Q[row] += salt/salt_d
                Q[row] /= sqrt(norm(Q[row]))
            end
            push!(ratio_history, -1000.0)
            downward = 0
            below = 0
        end
    end
    println( "best overlap: ", best_overlap)
    println( "ratio_history:")
    display(ratio_history)
    println()
    println()
    #=for row in 1:Ny
        A_r_inds  = [commonind(A[row, col], A[row, next_col]) for row in 1:Ny]
        QR_inds  .= [Index(ITensors.dim(A_r_inds[row]), "Site,QR,c$col,r$row") for row in 1:Ny]
        replaceind!(best_Q[row], next_col_inds[row], QR_inds[row])
    end
    q_dict[col=>side] = best_Q=#
    return best_Q, best_R, next_col_inds, QR_inds, dummy_nexts
end

function gaugeColumn( A::fPEPS, col::Int, side::Symbol; kwargs...)
    Ny, Nx = size(A)

    prev_col_inds = Vector{Index}(undef, Ny)
    next_col_inds = Vector{Index}(undef, Ny)

    next_col   = side == :left ? col - 1 : col + 1
    prev_col   = side == :left ? col + 1 : col - 1
    left_edge  = col == 1
    right_edge = col == Nx
    is_cu      = is_gpu(A)

    @timeit "gauge QR" begin
        Q, R, next_col_inds, QR_inds, dummy_next_inds = gaugeQR(A, col, side; kwargs...)
    end
    @timeit "update next col" begin
        cmb_r = Vector{ITensor}(undef, Ny)
        cmb_u = Vector{ITensor}(undef, Ny - 1)
        if (side == :left && col > 1 ) || (side == :right && col < Nx)
            next_col_As = MPO(deepcopy(A[:, next_col]), 0, Ny+1)
            nn_col = side == :left ? next_col - 1 : next_col + 1
            cmb_inds_ = [IndexSet(firstind(A[row, next_col], "Site")) for row in 1:Ny]
            for row in 1:Ny
                cmb_inds_r = cmb_inds_[row]
                if 0 < nn_col < Nx + 1
                    cmb_inds_r = IndexSet(cmb_inds_[row]..., commonind(A[row, next_col], A[row, nn_col]))
                end
                cmb_r[row] = combiner(cmb_inds_r, tags="Site,CMB,c$col,r$row")
                next_col_As[row] *= cmb_r[row]
                if (side == :left && !left_edge) || (side == :right && !right_edge)
                    next_col_As[row] = replaceind!(next_col_As[row], next_col_inds[row], dummy_next_inds[row])
                end
            end
        end
        maxdim::Int  = get(kwargs, :maxdim, 1)
        result       = multMPO(R, next_col_As; kwargs...)
        a_norm      = is_cu ? cuITensor(1.0) : ITensor(1.0)
        q_norm      = is_cu ? cuITensor(1.0) : ITensor(1.0)
        r_norm      = is_cu ? cuITensor(1.0) : ITensor(1.0)
        for row in 1:Ny
            a_norm      *= dag(result[row]) * result[row]
            q_norm      *= dag(Q[row]) * Q[row]
            r_norm      *= dag(R[row]) * R[row]
        end
        true_QR_inds = [Index(ITensors.dim(QR_inds[row]), "Link,r,r$row" * (side == :left ? ",c$(col-1)" : ",c$col")) for row in 1:Ny]
        for row in 1:Ny
            A[row, col] = Q[row]
        end
        cUs = [commonind(A[row, col], A[row+1, col]) for row in 1:Ny-1]
        true_U_inds = [Index(ITensors.dim(cUs[row]), "Link,u,r$row,c$col") for row in 1:Ny-1]
        A[:, col] = [replaceind!(A[row, col], QR_inds[row], true_QR_inds[row]) for row in 1:Ny]
        A[:, col] = vcat([replaceind!(A[row, col], cUs[row], true_U_inds[row]) for row in 1:Ny-1], A[Ny, col])
        A[:, col] = vcat(A[1, col], [replaceind!(A[row, col], cUs[row-1], true_U_inds[row-1]) for row in 2:Ny])

        for row in 1:Ny
            A[row, next_col] = result[row]
        end
        A[:, next_col] = [replaceind!(A[row, next_col], QR_inds[row], true_QR_inds[row]) for row in 1:Ny]
        cUs            = [commonind(A[row, next_col], A[row+1, next_col]) for row in 1:Ny-1]
        true_nU_inds   = [Index(ITensors.dim(cUs[row]), "Link,u,r$row,c" * string(next_col)) for row in 1:Ny-1]
        A[:, next_col] = vcat([replaceind!(A[row, next_col], cUs[row], true_nU_inds[row]) for row in 1:Ny-1], A[Ny, next_col])
        A[:, next_col] = vcat(A[1, next_col], [replaceind!(A[row, next_col], cUs[row-1], true_nU_inds[row-1]) for row in 2:Ny])
        A[:, next_col] = [A[row, next_col] * cmb_r[row] for row in 1:Ny]
    end
    return A
end

function gaugeColumnHorizontal( A::fPEPS, Aj::Vector{ITensor}, col::Int, ncol::Int, side::Symbol; kwargs...)
    Ny, Nx = size(A)

    prev_col_inds = Vector{Index}(undef, Ny)
    next_col_inds = Vector{Index}(undef, Ny)

    next_col   = side == :left ? col - 1 : col + 1
    prev_col   = side == :left ? col + 1 : col - 1
    left_edge  = col == 1
    right_edge = col == Nx
    is_cu      = is_gpu(A)

    @timeit "gauge QR" begin
        Q, R, next_col_inds, QR_inds, dummy_next_inds = gaugeQRHorizontal(A, Aj, col, ncol, side; kwargs...)
    end
    @timeit "update next col" begin
        cmb_r = dummy_next_inds
        cmb_u = Vector{ITensor}(undef, Ny - 1)
        maxdim::Int  = get(kwargs, :maxdim, 1)
        true_QR_inds = [Index(ITensors.dim(QR_inds[row]), "Link,r,r$row" * (side == :left ? ",c$(col-1)" : ",c$col")) for row in 1:Ny]
        for row in 1:Ny
            A[row, col] = Q[row]
        end
        cUs = [commonind(A[row, col], A[row+1, col]) for row in 1:Ny-1]
        true_U_inds = [Index(ITensors.dim(cUs[row]), "Link,u,r$row,c$col") for row in 1:Ny-1]
        A[:, col] = [replaceind!(A[row, col], QR_inds[row], true_QR_inds[row]) for row in 1:Ny]
        A[:, col] = vcat([replaceind!(A[row, col], cUs[row], true_U_inds[row]) for row in 1:Ny-1], A[Ny, col])
        A[:, col] = vcat(A[1, col], [replaceind!(A[row, col], cUs[row-1], true_U_inds[row-1]) for row in 2:Ny])

        for row in 1:Ny
            A[row, next_col] = R[row]
        end
        A[:, next_col] = [replaceind!(A[row, next_col], QR_inds[row], true_QR_inds[row]) for row in 1:Ny]
        cUs            = [commonind(A[row, next_col], A[row+1, next_col]) for row in 1:Ny-1]
        true_nU_inds   = [Index(ITensors.dim(cUs[row]), "Link,u,r$row,c" * string(next_col)) for row in 1:Ny-1]
        A[:, next_col] = vcat([replaceind!(A[row, next_col], cUs[row], true_nU_inds[row]) for row in 1:Ny-1], A[Ny, next_col])
        A[:, next_col] = vcat(A[1, next_col], [replaceind!(A[row, next_col], cUs[row-1], true_nU_inds[row-1]) for row in 2:Ny])
        A[:, next_col] = [A[row, next_col] * cmb_r[row] for row in 1:Ny]
    end
    return A
end

function gaugeColumnForInsert( A::fPEPS, col::Int, side::Symbol; kwargs...)
    Ny, Nx = size(A)

    prev_col_inds = Vector{Index}(undef, Ny)
    next_col_inds = Vector{Index}(undef, Ny)

    next_col   = side == :left ? col - 1 : col + 1
    prev_col   = side == :left ? col + 1 : col - 1
    left_edge  = col == 1
    right_edge = col == Nx
    is_cu      = is_gpu(A)

    @timeit "gauge QR" begin
        Q, R, next_col_inds, QR_inds, dummy_next_inds = gaugeQR(A, col, side; kwargs...)
    end
    return A, Q, R, dummy_next_inds
end


text:
__precompile__(true)

module Single
export single

"""
    single(x [,default])

Returns the one and single element of collection `x`, and throws an error if the collection
has zero or multiple elements. If `default` is provided, `default` is returned in the case of
an empty iterators instead of throwing an error. For iterators with more then one element,
an error is still thrown.
"""
function single end

# https://github.com/JuliaLang/julia/pull/25078
Base.@propagate_inbounds function single(iter)
    i0 = start(iter)
    @boundscheck if done(iter, i0)
        throw(ArgumentError("Collection is empty, must contain exactly one element"))
    end
    (ret, i1) = next(iter, i0)
    @boundscheck if !done(iter, i1)
        throw(ArgumentError("Collection has multiple elements, must contain exactly one element"))
    end
    ret
end

Base.@propagate_inbounds function single(iter, default)
    i0 = start(iter)
    @boundscheck if done(iter, i0)
        return default
    end
    (ret, i1) = next(iter, i0)
    @boundscheck if !done(iter, i1)
        throw(ArgumentError("Collection has multiple elements, must contain exactly one element"))
    end
    return ret
end

# Collections of known size
single(::Tuple{}, default) = default
single(x::Tuple{Any}) = x[1]
single(x::Tuple{Any}, default) = x[1]

single(a::AbstractArray{<:Any, 0}) = @inbounds return a[]

end # module


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule MiniFB_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("MiniFB")
JLLWrappers.@generate_main_file("MiniFB", UUID("b93afbca-07d8-5443-9979-bb4221978155"))
end  # module MiniFB_jll