text:
"""
    cleanNames!(mydata::Vector{String})

Replaces problematic characters in column names, avoids duplicate names, and
prefixes an '_' if the name starts with a number.

# Arguments:
- `mydata`: vector of names (gets modified)
"""
function cleanNames!(mydata::Vector{String})
    # replace problematic characters,
    # put "_" in front of colname in case it starts with a number
    # avoid duplicate names (add suffixes _2, _3, ...)
    usedNames = Set{String}()
    for j in eachindex(mydata)
        mydata[j] = replace(mydata[j], "-" => "_")
        if isnumeric(first(mydata[j]))
            mydata[j] = "_" * mydata[j]
        end
        # avoid duplicate names
        if mydata[j] in usedNames
            idx = 2
            while "$(mydata[j])_$idx" in usedNames
                idx += 1
            end
            mydata[j] *= "_$idx"
        end
        push!(usedNames, mydata[j])
    end
end

"""
    getMetaData(f)
Collect the meta data information in a more user friendly format.

# Arguments:
- `f`: input structure with `.params` and `.data` fields
"""
function getMetaData(meta::Dict{String,String})::DataFrame

    # declarations and initializations
    metaKeys = keys(meta)
    channel_properties = []
    defaultValue = ""

    # determine the number of channels
    pars = parse(Int, strip(join(meta["\$PAR"])))

    # determine the available channel properties
    for (key,) in meta
        if length(key) >= 2 && key[1:2] == "\$P"
            i = 3
            while i <= length(key) && isdigit(key[i])
                i += 1
            end
            if i <= length(key) && !in(key[i:end], channel_properties)
                push!(channel_properties, key[i:end])
            end
        end
    end

    # create a data frame for the results
    df = Matrix{String}(undef, pars, length(channel_properties))
    df .= defaultValue
    df = DataFrame(df)
    rename!(df, Symbol.(channel_properties))

    # collect the data from params
    for ch = 1:pars
        for p in channel_properties
            if "\$P$ch$p" in metaKeys
                df[ch, Symbol(p)] = meta["\$P$ch$p"]
            end
        end
    end

    return df
end

"""
    getMarkerNames(meta::DataFrame)::Tuple{Vector{String}, Vector{String}}

Extract suitable raw names (useful for selecting columns) and pretty readable
names (useful for humans) from FCS file metadata.

"""
function getMarkerNames(meta::DataFrame)::Tuple{Vector{String},Vector{String}}
    orig = Array{String}(meta[:, :N])
    nice = copy(orig)
    if hasproperty(meta, :S)
        for i = 1:size(meta, 1)
            if strip(meta[i, :S]) != ""
                nice[i] = meta[i, :S]
            end
        end
    end
    return (orig, nice)
end

"""
    compensate!(data::Matrix{Float64}, spillover::Matrix{Float64}, cols::Vector{Int})

Apply a compensation matrix in `spillover` (the individual columns of which
describe, in order, the spillover of `cols` in `data`) to the matrix `data`
in-place.
"""
function compensate!(data::Matrix{Float64}, spillover::Matrix{Float64}, cols::Vector{Int})
    data[:, cols] = data[:, cols] * inv(spillover)
end

"""
    parseSpillover(str::String)::Union{Tuple{Vector{String},Matrix{Float64}}, Nothing}

Parses the spillover matrix from the string from FCS parameter value.
"""
function parseSpillover(str::String)::Tuple{Vector{String},Matrix{Float64}}
    fields = split(str, ',')
    n = parse(Int, fields[1])
    if length(fields) != 1 + n + n * n
        @error "Spillover matrix of size $n expects $(1+n+n*n) fields, got $(length(fields)) instead."
        error("Invalid spillover matrix")
    end

    names = fields[2:(1+n)]
    spill = collect(transpose(reshape(parse.(Float64, fields[(2+n):length(fields)]), n, n)))

    names, spill
end

"""
    getSpillover(params::Dict{String, String})::Union{Tuple{Vector{String},Matrix{Float64}}, Nothing}

Get a spillover matrix from FCS `params`. Returns a pair with description of
columns to be applied, and with the actual spillover matrix. Returns `nothing`
in case spillover is not present.
"""
function getSpillover(
    params::Dict{String,String},
)::Union{Tuple{Vector{String},Matrix{Float64}},Nothing}
    spillNames = ["\$SPILL", "\$SPILLOVER", "SPILL", "SPILLOVER"]
    for i in spillNames
        if in(i, keys(params))
            return parseSpillover(params[i])
        end
    end
    return nothing
end


text:
import SymbolicRegression: L2DistLoss

maximum_residual = 1e-2
custom_cos(x) = cos(x)

default_params = (
    binary_operators=(/, +, *),
    unary_operators=(exp, custom_cos),
    constraints=nothing,
    loss=L2DistLoss(),
    ns=10,
    topn=10,
    parsimony=0.000100f0,
    alpha=0.100000f0,
    maxsize=20,
    maxdepth=nothing,
    fast_cycle=false,
    migration=true,
    hofMigration=true,
    fractionReplacedHof=0.1f0,
    shouldOptimizeConstants=true,
    hofFile=nothing,
    npopulations=nothing,
    perturbationFactor=1.000000f0,
    annealing=true,
    batching=false,
    batchSize=50,
    mutationWeights=[
        10.000000, 1.000000, 1.000000, 3.000000, 3.000000, 0.010000, 1.000000, 1.000000
    ],
    crossoverProbability=0.0f0,
    warmupMaxsizeBy=0.0f0,
    useFrequency=false,
    npop=1000,
    ncyclesperiteration=300,
    fractionReplaced=0.1f0,
    verbosity=convert(Int, 1e9),
    probNegate=0.01f0,
    seed=nothing,
    bin_constraints=nothing,
    una_constraints=nothing,
    progress=false,
    terminal_width=nothing,
    warmupMaxsize=nothing,
    optimizer_algorithm="NelderMead",
    optimizer_nrestarts=3,
    optimize_probability=0.1f0,
    optimizer_iterations=100,
    recorder=nothing,
    recorder_file="pysr_recorder.json",
    probPickFirst=1.0,
    earlyStopCondition=nothing,
    stateReturn=false,
    timeout_in_seconds=nothing,
    skip_mutation_failures=false,
)

ext: jl


text:
#-----------------------------------------------------------------------------# find_matches
function find_matches(left_row, right)
    findall(x -> cmp(x, left_row), right)
end

cmp(x, y) = x == y

function check_join_fields(left, right, on)
    on = on isa Col ? field_name(left, on.sel) : on
    on_keys = names(select(left, on))
    lkeys = setdiff(names(left), on_keys)
    rkeys = setdiff(names(right), on_keys)
    for k in lkeys
        k in rkeys && error("left and right both contain non-key field: $k")
    end
end

#-----------------------------------------------------------------------------# hash_join
function hash_join(left::CTable, right::CTable, on)
    left_on = select(left, on)
    right_on = select(right, on)

    out = merge(
        OrderedDict(k => similar(v, 0) for (k,v) in pairs(cols(left))),
        OrderedDict(k => similar(v, 0) for (k,v) in pairs(cols(right))),
    )
    left_hash_table = OrderedDict(i => hash(row) for (i,row) in enumerate(left_on))
    right_hash_table = OrderedDict(hash(row) => i for (i,row) in enumerate(right_on))

    for (i,h) in pairs(left_hash_table)
        j = get(right_hash_table, h, nothing)
        if !isnothing(j)
            for (k,v) in pairs(cols(left))
                push!(out[k], v[i])
            end
            for (k,v) in pairs(cols(right))
                if on isa Col
                    k != on.sel && push!(out[k], v[j])
                elseif !(k in names(left_on))
                    push!(out[k], v[j])
                end
            end
        end
    end
    return CTable(namedtuple(out))
end

#-----------------------------------------------------------------------------# inner_join
function inner_join(left::CTable, right::CTable, on)
    check_join_fields(left, right, on)
    left_out = CTable(map(x -> similar(x, 0), cols(left)))
    right_out = CTable(map(x -> similar(x, 0), cols(right)))

    left_sub = select(left, on)
    right_sub = select(right, on)
    for (i,left_sub_row) in enumerate(left_sub)
        idx = find_matches(left_sub_row, right_sub)
        for j in idx
            push!(left_out, left[i])
            push!(right_out, right[j])
        end
    end
    return CTable(merge(cols(left_out), cols(right_out)))
end

#-----------------------------------------------------------------------------# left_join
function left_join(left::CTable, right::CTable, on)
    check_join_fields(left, right, on)
    out = CTable(merge(
        map(x -> similar(x, 0), cols(left)),
        map(x -> similar(x, Union{Missing,eltype(x)}, 0), cols(right))
    ))
    keys2add = setdiff(names(right), names(left))
    left_sub = select(left, on)
    right_sub = select(right, on)
    for left_row in left
        idx = find_matches(left_row, right_sub)
        if isempty(idx)
            push!(out, merge(left_row, namedtuple(keys2add, fill(missing, length(keys2add)))))
        else
            for j in idx
                push!(out, merge(left_row, right[j]))
            end
        end
    end
    return out
end


text:
using Revise
using OptiMode
using GeometryPrimitives
using LinearAlgebra
using StaticArrays
using LoopVectorization
using ChainRules
using Zygote
using ForwardDiff
using UnicodePlots
using FFTW
using OhMyREPL
using Crayons.Box       # for color printing
using Rotations: RotY, MRP
using RuntimeGeneratedFunctions
RuntimeGeneratedFunctions.init(@__MODULE__)
LNx = rotate(MgO_LiNbO₃,Matrix(MRP(RotY(π/2))),name=:LiNbO₃)
AD_style = BOLD*BLUE_FG #NEGATIVE*BOLD*BLUE_FG      # defined in Crayons.Box
FD_style = BOLD*RED_FG
MAN_style = BOLD*GREEN_FG
AD_style_N = NEGATIVE*BOLD*BLUE_FG #NEGATIVE*BOLD*BLUE_FG      # defined in Crayons.Box
FD_style_N = NEGATIVE*BOLD*RED_FG
MAN_style_N = NEGATIVE*BOLD*GREEN_FG

Δx,Δy,Δz,Nx,Ny,Nz = 6.0, 4.0, 1.0, 128, 128, 1;
grid = Grid(Δx,Δy,Nx,Ny)
rwg(x) = ridge_wg_partial_etch(x[1],x[2],x[3],x[4],0.5,LNx,SiO₂,Δx,Δy) # partially etched ridge waveguide with dispersive materials, x[3] is partial etch fraction of top layer, x[3]*x[2] is etch depth, remaining top layer thickness = x[2]*(1-x[3]).
p = [
       1.7,                #   top ridge width         `w_top`         [μm]
       0.7,                #   ridge thickness         `t_core`        [μm]
       0.5,                #   ridge thickness         `t_core`        [μm]
       π / 14.0,           #   ridge sidewall angle    `θ`             [radian]
               ];
geom = rwg(p)
##
ms = ModeSolver(1.45, geom, grid, nev=4)
ωs = range(0.6,0.8,length=16) |> collect
using OptiMode: _solve_n_serial
##

function E_relpower_xyz(ms::ModeSolver{ND,T},ω²H) where {ND,T<:Real}
    E = 1im * ε⁻¹_dot( fft( kx_tc( reshape(ω²H[2],(2,size(ms.grid)...)),mn(ms),ms.M̂.mag), (2:1+ND) ), flat( ms.M̂.ε⁻¹ ))
    Es = reinterpret(reshape, SVector{3,Complex{T}},  E)
    Pₑ_xyz_rel = normalize([mapreduce((ee,epss)->(abs2(ee[a])*inv(epss)[a,a]),+,Es,ms.M̂.ε⁻¹) for a=1:3],1)
    return Pₑ_xyz_rel
end
# E_relpower_xyz(ms,(real(ms.ω²[eigind]),ms.H⃗[:,eigind]))
TE_filter = (ms,ω²H)->E_relpower_xyz(ms,ω²H)[1]>0.7
TM_filter = (ms,ω²H)->E_relpower_xyz(ms,ω²H)[2]>0.7
oddX_filter = αX->sum(abs2,𝓟x̄(ms.grid)*αX[2])>0.9
evenX_filter = αX->sum(abs2,𝓟x(ms.grid)*αX[2])>0.9
ns_TE,ngs_TE,gvds_TE,Evs_TE = _solve_n_serial(ms,ωs,geom;f_filter=TE_filter)
ns_TM,ngs_TM,gvds_TM,Evs_TM = _solve_n_serial(ms,ωs,geom;f_filter=TM_filter)
##
using GLMakie, AbstractPlotting
using AbstractPlotting: lines, lines!, scatterlines, scatterlines!, GeometryBasics, Point, PointBased
using Colors
import Colors: JULIA_LOGO_COLORS
logocolors = JULIA_LOGO_COLORS

##
# interactions

function indicator(ax::Axis,ob)
    register_interaction!(ax, :indicator) do event::MouseEvent, axis
    if event.type === MouseEventTypes.over
        ob[] = event.data
    end
    end
end
function indicator(grid::GridLayout,ob)
    foreach(Axis,grid;recursive=true) do ax
    indicator(ax,ob)
    end
end
function indicator(grid::GridLayout)
    ob = Observable(Point2f0(0.,0.))
    indicator(grid,ob)
    ob
end
function indicator(fig,args...; tellwidth=false, kwargs...)
    Label(
        fig,
        lift(ind->"x: $(ind[1])  y: $(ind[2])",indicator(fig.layout)),
        args...; tellwidth=tellwidth, kwargs...
    )
end

##

ωind = 5
Es = [ view(real(Evs_TE),i,:,:,ωind) for i=1:3 ]
Emagmax = [sqrt(maximum(abs2.(Evs_TE[:,:,:,ωind]))) for i=1:3]

fig = Figure()

# dispersion plots
disp = ( (ns_TE,ns_TM), (ngs_TE,ngs_TM), (gvds_TE,gvds_TM) )
# labels_disp = [la*lb for la in ["neff", "ng", "gvd"], lb in [" TE₀₀", " TM₀₀"]]
labels_disp = [" TE₀₀", " TM₀₀"]
# labels_disp = ["n", "ng" ] .* [",TE",",TM"]

n_disp = length(disp)
n_modes = 2
λs = inv.(ωs)
ax_disp = fig[1:n_disp,1] = [Axis(fig) for i=1:n_disp]
colors_disp = [logocolors[:red],logocolors[:blue]]
sls_disp 	= 	[scatterlines!(
					ax_disp[i],
					λs,
					disp[i][j],
					color=colors_disp[j],
					markercolor=colors_disp[j],
					strokecolor=colors_disp[j],
					markersize=2,
					linewidth=2,
					label=labels_disp[j],
				) for i=1:n_disp, j=1:n_modes ]
disp_models = [:n,:ng,:gvd]
lns_bulk = [plot_model!(ax_disp[i],[LNx,];model=disp_models[i],linewidth=3) for i=1:3]

fig[4,1] = Legend(fig,ax_disp[1],orientation=:horizontal)

hidexdecorations!.(ax_disp[1:end-1])
[axx.xlabel= "λ [μm]" for axx in ax_disp[end,:]]
[axx.ylabel= lbl for (axx,lbl) in zip(ax_disp,["neff","ng","gvd"])]
linkxaxes!(ax_disp...)

# spatial plots
xs = x(ms.grid)
ys = y(ms.grid)
ax_geom = fig[1,2] = Axis(fig)
plot_shapes(geom,ax_geom)

ax_n	= fig[1,3] = Axis(fig)
ax_E 	= fig[2,2:3] = [Axis(fig) for i=1:2]

# poly!(ax_geom,geom.shapes[1],color=:red,strokecolor=:black,strokewidth=1)
# poly!(ax_geom,geom.shapes[2],color=:blue,strokecolor=:black,strokewidth=1)

label_n = "nₓ"
cmap_n 	  = :viridis
n_spatial = sqrt.(getindex.(inv.(ms.M̂.ε⁻¹),1,1))
n_spatial_max = maximum(n_spatial)
heatmap_n = heatmap!(ax_n, xs, ys, n_spatial, colormap=cmap_n,label=label_n,colorrange=(1,n_spatial_max);overdraw=false)
# wf_n = wireframe!(ax_n, xs, ys, n_spatial, colormap=cmap_n,linewidth=0.1,color=:white)
cbar_n = Colorbar(fig[1,4], heatmap_n,  width=20 )
text!(ax_n,label_n,position=(-1.4,1.1),textsize=0.7,color=:white)

cmap_E = :diverging_bkr_55_10_c35_n256
label_base = ["x","y"]
labels_E = "E".*label_base
heatmaps_E = [heatmap!(ax_E[j], xs, ys, Es[j],colormap=cmap_E,label=labels_E[j],colorrange=(-Emagmax[j],Emagmax[j])) for j=1:2]
cbars_E = Colorbar(fig[2,4], heatmaps_E[2],  width=20 )
# wfs_E = [wireframe!(ax_E[j], xs, ys, Es[j], colormap=cmap_E,linewidth=0.02,color=:white) for j=1:2]
map( (axx,ll)->text!(axx,ll,position=(-1.4,1.1),textsize=0.7,color=:white), ax_E, labels_E )

ax_spatial = vcat(ax_geom, ax_n, ax_E)
hidexdecorations!.(ax_spatial[1:end-1,:])
hideydecorations!.(ax_spatial[:,2:end])
[axx.xlabel= "x [μm]" for axx in ax_spatial[end,:]]
[axx.ylabel= "y [μm]" for axx in ax_spatial[:,1]]
[ axx.aspect=DataAspect() for axx in ax_spatial ]
linkaxes!(ax_spatial...)
txt= fig[4,2] = indicator(fig)
# for axx in ax_all
# 	on(mouseposition) do mpos
#
# end

fig

##

using Distributed
pids = addprocs(8)
@show wp = CachingPool(workers()) #default_worker_pool()
@everywhere begin
	using LinearAlgebra, Statistics, FFTW, StaticArrays, HybridArrays, ChainRules, Zygote, ForwardDiff, GeometryPrimitives, OptiMode
	using Zygote: dropgrad, @ignore
	using Rotations: RotY, MRP
	LNx = rotate(MgO_LiNbO₃,Matrix(MRP(RotY(π/2))))
	Δx,Δy,Δz,Nx,Ny,Nz = 6.0, 4.0, 1.0, 128, 128, 1;
	# Δx,Δy,Δz,Nx,Ny,Nz = 6.0, 4.0, 1.0, 256, 256, 1;
	gr = Grid(Δx,Δy,Nx,Ny)
	p_pe = [
	       1.7,                #   top ridge width         `w_top`         [μm]
	       0.7,                #   ridge thickness         `t_core`        [μm]
	       0.5,                #   top layer etch fraction `etch_frac`     [1]
	       π / 14.0,           #   ridge sidewall angle    `θ`             [radian]
	               ];
	rwg_pe(x) = ridge_wg_partial_etch(x[1],x[2],x[3],x[4],0.5,LNx,SiO₂,Δx,Δy) # partially etched ridge waveguide with dispersive materials, x[3] is partial etch fraction of top layer, x[3]*x[2] is etch depth, remaining top layer thickness = x[2]*(1-x[3]).

	geom_pe = rwg_pe(p_pe)
	ms = ModeSolver(1.45, rwg_pe(p_pe), gr)
	p_pe_lower = [0.4, 0.3, 0., 0.]
	p_pe_upper = [2., 2., 1., π/4.]

	λs = collect(reverse(1.45:0.02:1.65))
	ωs = 1 ./ λs

	geom1 = rwg_pe(p_pe)
	# numngs1(oms) = sum(solve_n(oms,rwg_pe(p_pe),gr)[2])
	# numngs2(oms,x) = sum(solve_n(oms,rwg_pe(x),gr)[2])
	# numngs3(oms) = sum(pmap(x->solve_n(x,rwg_pe(p_pe),gr)[2]),oms))

	# numngs1([0.6,0.7])
	# n1F,ng1F = solve_n(ms,ωs,rwg_pe(p_pe)); n1S,ng1S = solve_n(ms,2*ωs,rwg_pe(p_pe))
end

##
@everywhere function foo1(x)
	Zygote.@ignore begin
	Δx,Δy,Δz,Nx,Ny,Nz = 6.0, 4.0, 1.0, 128, 128, 1;
	p_pe = [
	      1.7,                #   top ridge width         `w_top`         [μm]
	      0.7,                #   ridge thickness         `t_core`        [μm]
	      0.5,                #   top layer etch fraction `etch_frac`     [1]
	      π / 14.0,           #   ridge sidewall angle    `θ`             [radian]
	               ];
	ms = ModeSolver(1.45, rwg_pe(p_pe), gr)
	end
	return x^2
end

@everywhere function foo2(x)
	Δx,Δy,Δz,Nx,Ny,Nz = 6.0, 4.0, 1.0, 128, 128, 1;
	p_pe = [
	      1.7,                #   top ridge width         `w_top`         [μm]
	      0.7,                #   ridge thickness         `t_core`        [μm]
	      0.5,                #   top layer etch fraction `etch_frac`     [1]
	      π / 14.0,           #   ridge sidewall angle    `θ`             [radian]
	               ];
	ms = ModeSolver(1.45, rwg_pe(p_pe), gr)
	return x^2
end

@everywhere function foo3(x)
	sumeps = Zygote.@ignore begin
	Δx,Δy,Δz,Nx,Ny,Nz = 6.0, 4.0, 1.0, 128, 128, 1;
	p_pe = [
	      1.7,                #   top ridge width         `w_top`         [μm]
	      0.7,                #   ridge thickness         `t_core`        [μm]
	      0.5,                #   top layer etch fraction `etch_frac`     [1]
	      π / 14.0,           #   ridge sidewall angle    `θ`             [radian]
	               ];
	ms = ModeSolver(1.45, rwg_pe(p_pe), gr)
	sum(sum(ms.M̂.ε⁻¹))
	end
	return sumeps*x^2
end

function _solve_n_parallel1(ωs::Vector{T},geom::Vector{<:Shape},gr::Grid;nev=1,eigind=1,maxiter=3000,tol=1e-8,log=false,ω²_tol=tol,wp=nothing) where {ND,T<:Real}
	wp = CachingPool(workers())
	nω = length(ωs)
	ind0 = Int(ceil(nω/2))
	ω0 = ωs[ind0]
	ms = @ignore(ModeSolver(kguess(ω0,geom), geom, gr))
	nng0 = solve_n(ms,ω0,geom)
	ms_copies = [ deepcopy(ms) for om in 1:nω ]
	geoms = [deepcopy(geom) for om in 1:nω ]
	nng = pmap(wp,ms_copies,ωs,geoms) do m,om,s
		@ignore( replan_ffts!(m) );
		solve_n(m,om,s)
	end
	n = [res[1] for res in nng]; ng = [res[2] for res in nng]
	return n, ng
end

function _solve_n_parallel2(ωs::Vector{T},geom::Vector{<:Shape},gr::Grid;nev=1,eigind=1,maxiter=3000,tol=1e-8,log=false,ω²_tol=tol,wp=nothing) where {ND,T<:Real}
	wp = CachingPool(workers())
	nω = length(ωs)
	# ind0 = Int(ceil(nω/2))
	# ω0 = ωs[ind0]
	# ms = @ignore(ModeSolver(kguess(ω0,geom), geom, gr))
	# nng0 = solve_n(ms,ω0,geom)
	# ms_copies = [ deepcopy(ms) for om in 1:nω ]
	geoms = [deepcopy(geom) for om in 1:nω ]
	nng = pmap(wp,ωs,geoms) do om,s
		@ignore( replan_ffts!(m) );
		solve_n(om,s,gr)
	end
	n = [res[1] for res in nng]; ng = [res[2] for res in nng]
	return n, ng
end

n1,ng1 = _solve_n_parallel1(ωs,rwg_pe(p_pe),gr)
gr2 = Grid(6.0, 4.0, 256, 256)
n2,ng2 = _solve_n_parallel1(ωs,rwg_pe(p_pe),gr2)

n1,ng1 = _solve_n_parallel2(ωs,rwg_pe(p_pe),gr)

##
gradient(var_ng2,ωs,p)

##

solve_n(ms,ωs,shapes;n_procs=10)
var_ng2(ωs,p)
var_ng1(ωs,p)

gradient(var_ng2,ωs,p)

@everywhere begin
	function var_ng1(ms,ωs,p,procs)
		ng = solve_n(ms,ωs,rwg(p);n_procs=procs)[2]
		mean( abs2.( ng ) ) - abs2(mean(ng))
	end

	function var_ng_pb(ms,ωs,p)
		var_ng,var_ng_pb = Zygote.pullback(ωs,p) do ωs,p
			ng = solve_n(ms,ωs,rwg(p))[2]
			mean( abs2.( ng ) ) - abs2(mean(ng))
		end
		return (var_ng, var_ng_pb(1))
	end
end

var_ng1(ms,ωs,p,10)
gradient(p) do p
	var_ng1(ms,ωs,p,10)
end

gradient((om,x)->var_ng1(ms,om,x,10), ωs, p)

gradient((om,x)->var_ng1(ms,om,x,1), ωs, p)

ωs2 = collect(0.55:0.05:1.0)
(n2,ng2), nng_pb = pullback(ωs2,p) do om, p
	s = rwg(p)
	solve_n(ms,om,s;n_procs=10)
end

plot_nng(ωs2,solve_n(ms,ωs2,shapes;n_procs=10)...;c_ng=:green,m=".")

# n,ng = solve_n(ms,ωs,shapes;n_procs=10)
using Plots
function plot_nng(ωs,n,ng;c_n=:blue,c_ng=:red,ls_n=:solid,ls_ng=:dash,
		legend=:bottomleft,m=nothing,xlabel="λ (μm)",ylabel="n, ng")
	p_nng = plot(
		(1 ./ ωs),
		n;
		ls=ls_n,
		label="n",
		color=c_n,
		legend,
		xlabel,
		ylabel,
		m
		)
	plot!(p_nng,
		(1 ./ ωs),
		ng;
		ls=ls_ng,
		label="ng",
		color=c_ng,
		m
		)
	return p_nng
end

plot_nng(ωs,n,ng;c_ng=:green)

vng,∂vng = var_ng(ms,ωs,p)
ε⁻¹ = make_εₛ⁻¹(rwg(p),dropgrad(ms))
nω = length(ωs)
ms_copies = [ deepcopy(ms) for om in 1:nω ]
vmap(ms_copies,ωs) do m,om
	solve_n(m,om,ε⁻¹)
end

solve_ω²(ms,ks)
solve_k(ms,ωs)
solve_k(ms,ωs,rwg(p))
solve_n(ms,ωs,rwg(p))

@btime solve_ω²($ms,$ks) # 1.982 s (26246 allocations: 19.68 MiB)
@btime solve_k($ms,$ωs) # 4.873 s (57515 allocations: 279.31 MiB)
@btime solve_k($ms,$ωs,rwg($p)) # 5.074 s (3335372 allocations: 450.57 MiB)
@btime solve_n($ms,$ωs,rwg($p)) # 5.147 s (3337492 allocations: 514.33 MiB)

@everywhere using OptiMode, ArrayInterface, StaticArrays, HybridArrays, ChainRules, Zygote

@everywhere p = [
    # 1.45,               #   propagation constant    `kz`            [μm⁻¹]
    1.7,                #   top ridge width         `w_top`         [μm]
    0.7,                #   ridge thickness         `t_core`        [μm]
    π / 14.0,           #   ridge sidewall angle    `θ`             [radian]
    2.4,                #   core index              `n_core`        [1]
    1.4,                #   substrate index         `n_subs`        [1]
    0.5,                #   vacuum gap at boundaries `edge_gap`     [μm]
]

m.M̂.𝓕! = plan_fft!(randn(ComplexF64, (3,m.M̂.Nx,m.M̂.Ny,m.M̂.Nz)),(2:4),flags=FFTW.PATIENT);
m.M̂.𝓕⁻¹! = plan_bfft!(randn(ComplexF64, (3,m.M̂.Nx,m.M̂.Ny,m.M̂.Nz)),(2:4),flags=FFTW.PATIENT);
m.M̂.𝓕 = plan_fft(randn(ComplexF64, (3,m.M̂.Nx,m.M̂.Ny,m.M̂.Nz)),(2:4),flags=FFTW.PATIENT);
m.M̂.𝓕⁻¹ = plan_bfft(randn(ComplexF64, (3,m.M̂.Nx,m.M̂.Ny,m.M̂.Nz)),(2:4),flags=FFTW.PATIENT);

@everywhere rwg2(p) = ridge_wg(p[1],p[2],p[3],p[6],p[4],p[5],6.0,4.0)

pmap(ms_copies,ωs) do m,om
	solve_n(m,om,rwg2(p))
end

n,ng = solve_n(ms,ωs,rwg(p))
var(ng)
mean( abs2.( ng .- mean(ng) ) )
mean( abs2.( ng ) ) - abs2(mean(ng))
var_ng
# solve_ω²(ms,1.5,ε⁻¹)
# (ω²,H⃗), ω²H⃗_pb = Zygote.pullback(1.5,ε⁻¹) do x,y
# 	solve_ω²(ms,x,y)
# end

function calc_ng(p)
	# ε⁻¹ = make_εₛ⁻¹(ridge_wg(p[2],p[3],p[4],p[7],p[5],p[6],Δx,Δy),g)
	ε⁻¹ = HybridArray{Tuple{3,3,Dynamic(),Dynamic(),Dynamic()},Float64,5,5,Array{Float64,5}}(make_εₛ⁻¹(ridge_wg(p[2],p[3],p[4],p[7],p[5],p[6],Δx,Δy),ms))
	solve_nω(ms,p[1],ε⁻¹;eigind=1)[2]
end

function calc_ng(ms,p)
	ng,ng_pb = Zygote.pullback(p) do p
		solve_nω(ms,p[1],ridge_wg(p[2],p[3],p[4],p[7],p[5],p[6],ms.M̂.Δx,ms.M̂.Δy);eigind=1)[2]
	end
	return (ng, real(ng_pb(1)[1]))
end

function calc_ng(p)
	solve_n(ms,p[1],ridge_wg(p[2],p[3],p[4],p[7],p[5],p[6],ms.M̂.Δx,ms.M̂.Δy);eigind=1)[2]
end

function calc_ng_pb(ms,p)
	ng,ng_pb = Zygote.pullback(p) do p
		solve_n(ms,p[1],ridge_wg(p[2],p[3],p[4],p[7],p[5],p[6],ms.M̂.Δx,ms.M̂.Δy);eigind=1)[2]
	end
	return (ng, real(ng_pb(1)[1]))
end

p̄ω_FD3 = FiniteDifferences.jacobian(central_fdm(3,1),x->calc_ng2(x),pω)[1][1,:]
p̄_FD3 = FiniteDifferences.jacobian(central_fdm(3,1),x->calc_ng(x),p)[1][1,:]

## minimal working pmap + Zygote example from:
# https://discourse.julialang.org/t/passing-constructed-closures-to-child-processes/34723
# where issues with closures are discussed
using Distributed
addprocs(8)
@everywhere begin
  using OptiMode, ChainRules, Zygote
  function f_pmap_zygote_solve(ωs, p)
	  Δx,Δy,Δz,Nx,Ny,Nz = 6.0, 4.0, 1.0, 128, 128, 1;
  	# Δx,Δy,Δz,Nx,Ny,Nz = 6.0, 4.0, 1.0, 256, 256, 1;
  	gr = Grid(Δx,Δy,Nx,Ny)
  	rwg_pe(x) = ridge_wg_partial_etch(x[1],x[2],x[3],x[4],0.5,MgO_LiNbO₃,SiO₂,Δx,Δy)
	geom = rwg_pe(p)
	nng = pmap(ωs) do om
  	  # @ignore( replan_ffts!(m) );
  	  solve_n(om,geom,gr)
    end
    # return sum(sum(xs))
	# n = [res[1] for res in nng]
	ng = [res[2] for res in nng]
	# return n, ng
	return sum(ng)
  end

end
##
# wp = default_worker_pool()
# A = randn(8,8) #sprand(200, 200, 0.01) + 200*I
# b0s = [randn(8) for i=1:8]
# Zygote.gradient(f_pmap_zygote_solve, A, b0s)
ωs = collect(range(0.6,0.7,length=10))
Zygote.gradient(f_pmap_zygote_solve, ωs, p_pe)


text:
using FileIO, LASArrays
using Test

@testset "LasArrays.jl" begin
    # Write your tests here.

end


text:
# unit under test
include("../../src/image/mask.jl")

# test data

array1 = ones(3, 3)
array2 = ones(5, 5)

mask1 = BitArray(
[0 0 0;
 0 1 0;
 0 0 0])

mask1Frisket = BitArray(
 [1 1 1;
  1 0 1;
  1 1 1])

mask2 = BitArray(
[0 0 0 0 0;
 0 0 0 0 0;
 0 0 1 0 0;
 0 0 0 0 0;
 0 0 0 0 0])

mask2Frisket = BitArray(
[0 0 0 0 0;
 0 1 1 1 0;
 0 1 0 1 0;
 0 1 1 1 0;
 0 0 0 0 0])

 mask2Frisket2 = BitArray(
 [1 1 1 1 1;
  1 1 1 1 1;
  1 1 0 1 1;
  1 1 1 1 1;
  1 1 1 1 1])

result = frisketMaskAroundMask(array1, mask1, 1)
@assert result == mask1Frisket

result = frisketMaskAroundMask(array2, mask2, 1)
@assert result == mask2Frisket

result = frisketMaskAroundMask(array2, mask2, 2)
@assert result == mask2Frisket2

# frisketDepth out of bounds of array
result = frisketMaskAroundMask(array2, mask2, 3)
@assert result == mask2Frisket2

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
using CEnum
using SQLParser
using SQLParser: NodeTag

const IGNORE_LIST = [
    :T_Invalid,
    :T_TIDBitmap, #Cvoid
    :T_List,
    :T_SeqScan, #duplicate of Scan
    :T_DistinctExpr, #duplicate of OpExpr
    :T_NullIfExpr, #duplicate of OpExpr
    :T_Value, # special union struct
    :T_String,
    :T_Null,
    :T_Integer,
    :T_Float,
]

function generate_code_for_abstract_tree_interface()
    open(joinpath(@__DIR__, "abstract_trees_ext.jl"), "w") do io
        SQLParser.generate_code_for_abstract_tree_interface(io)
    end
end

function generate_code_for_abstract_tree_interface(io)

    Base.print(io,
    """
    # This is an auto-generated file.
    # You can regenerate this file using the following command:
    # `julia --project=. ./src/abstract_trees_codegen.jl`
    """)

    node_tag_enum_names = map(x -> x[1], CEnum.name_value_pairs(NodeTag))

    target_node_tag_enum_names = filter(x -> !(x in IGNORE_LIST), node_tag_enum_names)

    for enum_name in target_node_tag_enum_names
        enum_name_str = string(enum_name)
        @assert startswith(enum_name_str, "T_")
        tp_name_str = SubString(enum_name_str, 3)
        tp_name = Symbol(tp_name_str)
        try
            tp = eval(Meta.parse("SQLParser.$tp_name"))
            tp == Nothing && continue
            local printnode_body =
            """
            hasfield($tp_name_str, :type) && $(enum_name != :T_Expr) && @assert x.type == $enum_name "x.type is \$(x.type), but expected to be $enum_name"
                Base.print(io, "$tp_name_str")
                Base.print(io, "(")
            """
            local children_body =
            """
            hasfield($tp_name_str, :type) && $(enum_name != :T_Expr) && @assert x.type == $enum_name "x.type is \$(x.type), but expected to be $enum_name"
                children = []
            """
            is_first_simple_field = true
            for (i, (fld, fld_tp)) in enumerate(zip(fieldnames(tp), tp.types))
                fld == :type && continue
                fld_name = string(fld)
                if is_simple_type(fld_tp)
                    printnode_body *=
                    """
                        $(!is_first_simple_field) && Base.print(io, ", ")
                        Base.print(io, "$fld_name")
                        Base.print(io, "=")
                        Base.print(io, simple_type_value(x.$fld_name))
                    """
                    is_first_simple_field = false
                else
                    children_body *=
                    """
                        push!(children, FieldNameAndValuePair("$fld_name", x.$fld_name))
                    """
                end
            end

            q =
            """
            function convert_to_node_type(::Val{$enum_name})
                return $tp_name
            end
            function AbstractTrees.printnode(io::IO, x::$tp_name)
                $printnode_body
                Base.print(io, ")")
            end

            function AbstractTrees.children(x::$tp_name)
                $children_body
                return children
            end
            """
            Base.print(io, q)
        catch e
            #swallow the error
        end
    end
end


text:
# Validated from http://www.ams.jhu.edu/~castello/625.414/Handouts/FractionalProg.pdf

@testset "Simple example" begin
    lfp = LinearFractionalModel(with_optimizer(Clp.Optimizer))
    x1 = @variable(lfp, lower_bound=0, base_name="x1")
    x2 = @variable(lfp, lower_bound=0, upper_bound=6, base_name="x2")
    @constraint(lfp, -x1 + x2 <= 4)
    @constraint(lfp, 2x1 + x2 <= 14)

    LinearFractional.set_objective(lfp, MOI.MIN_SENSE,
              @expression(lfp, -2x1 + x2 + 2),
              @expression(lfp, x1 + 3x2 + 4))

    optimize!(lfp)
    @test termination_status(lfp) === MOI.OPTIMAL
    @test value(x1) ≈ 7.0
    @test value(x2) ≈ 0.0
end


text:
using AlphaGo

set_all_params(6)

cur_nn = NeuralNet(; tower_height = 1)
prev_nn = deepcopy(cur_nn)

NUM_GAMES = 32

MEM_SIZE = 50000
BATCH_SIZE = 32
EVAL_FREQ = 4

pos_buffer = Vector{AlphaGo.go.Position}()
π_buffer = Vector{Vector{Float32}}()
res_buffer = Vector{Int}()

function get_replay_batch(pos_buffer, π_buffer, res_buffer)
  idxs = rand(1:length(pos_buffer), BATCH_SIZE)
  pos_replay = pos_buffer[idxs]
  π_replay = hcat(π_buffer[idxs]...)
  res_replay = res_buffer[idxs]

  pos_replay, π_replay, res_replay
end

for i = 1:NUM_GAMES
  player = selfplay(cur_nn, 8)
  p, π, v = extract_data(player)

  pos_buffer = vcat(pos_buffer, p)
  π_buffer = vcat(π_buffer, π)
  res_buffer = vcat(res_buffer, v)

  if length(pos_buffer) > MEM_SIZE
    pos_buffer = pos_buffer[end-MEM_SIZE+1:end]
    π_buffer = π_buffer[end-MEM_SIZE+1:end]
    res_buffer = res_buffer[end-MEM_SIZE+1:end]
  end

  replay_pos, replay_π, replay_res = get_replay_batch(pos_buffer, π_buffer, res_buffer)
  loss = train!(cur_nn, (replay_pos, replay_π, replay_res))

  print("Episode $i over. Loss: $loss")
  if i % EVAL_FREQ == 0
    cur_is_winner = evaluate(cur_nn, prev_nn; num_games = 3, ro = 64)
    if cur_is_winner
      prev_nn = deepcopy(cur_nn)
    else
      cur_nn = deepcopy(prev_nn)
    end
    print("Evaluated")
  end
  println()
end


text:
const var = :value

type GetterSpec
    name
    argnames
    argtypes

    attr_idx
    index_idx

    result_specs

    # Simple validating constructor
    function GetterSpec(name, argnames, argtypes, attr_idx = 2, index_idx = 3;
        result_specs = Dict())
        namelen = length(argnames)
        typelen = length(argtypes)
        @assert namelen == typelen
        # Support for the following call signatures:
        # getter(attr, value)
        # getter(obj, attr, value)
        # getter(obj, attr, index, value)
        @assert namelen >= 2 && namelen <= 4

        if namelen <= 3
            index_idx = nothing
        end

        new(name, argnames, argtypes, attr_idx, index_idx, result_specs)
    end
end

type ResultSpec
    eltype
    count
    indexed

    function ResultSpec(eltype, count = 1;
        indexed = false
        )
        return new(eltype, count, indexed)
    end
end

type Getter
    signature
    setup_code
    c_call
    convert_code

    function Getter(getter_spec, key)
        signature = getter_signature(getter_spec, key)
        setup_code, convert_code = getter_code(getter_spec, key)
        c_call = getter_ccall(getter_spec, key)

        return new(signature, setup_code, c_call, convert_code)
    end
end

import Base.convert

function convert(::Type{Expr}, g :: Getter)
    return Expr(:function,
        g.signature,
        quote
            $(g.setup_code)
            err = $(g.c_call)

            test_status(err)

            $(g.convert_code)
        end)
end

function getter_name(discrim :: Symbol)
    replace(lowercase(string(discrim)), "hsa_", "")
end

function getter_code(getter_spec, key)
    result_spec = getter_spec.result_specs[key]

    getter_code(getter_spec, Val{result_spec.eltype}, result_spec.count)
end

# max size Array
function getter_code(getter_spec, eltype, count :: Int)
    t = eltype.parameters[1]

    if count == 1
        return :($var = Ref{$t}(Base.zero($t))),
            :($var = $var.x)
    else
        return :(v_arr = Array($t, $count);
                 $var = pointer(v_arr)),
            :($var = v_arr)
    end
end

function getter_code(getter_spec, eltype :: Type{Val{AbstractString}}, count :: Int)
    if count == 1
        # not a single character but a string constant
        return :($var = Ptr{UInt8}[0]),
            :($var = bytestring($var[1]))
    else
        # max size string
        return :($var = Array(UInt8, $count)),
            :($var = strip(ascii($var), '\0'))
    end
end

function getter_code(getter_spec, eltype :: Type{Val{AbstractString}}, count :: Symbol)
    # AbstractString with separate length getter
    len_sig = getter_signature(getter_spec, count)

    return :(len = $len_sig;
        $var = Array(UInt8, len)),
        :($var = strip(ascii($var), '\0'))
end

function getter_signature(getter_spec, key)
    name = getter_name(key)
    argnames = getter_spec.argnames
    result_spec = getter_spec.result_specs[key]

    signature = Expr(:call, symbol(name))

    # last one is presumed to be the output arg
    for i in 1:length(argnames) - 1
        arg_name = argnames[i]

       if (i != getter_spec.index_idx || result_spec.indexed) && # only emit index argument if it is necessary
           (i != getter_spec.attr_idx) # attr argument is superfluous

           push!(signature.args, arg_name)
       end
    end

    return signature
end

function getter_ccall(getter_spec, key)
    c_name = getter_spec.name
    argnames = getter_spec.argnames
    argtypes = getter_spec.argtypes
    result_spec = getter_spec.result_specs[key]

    ccall_args = [
        Expr(:tuple, Expr(:quote, c_name), :libhsa),
        :hsa_status_t,
        Expr(:tuple, argtypes...)
    ]

    for i in 1:length(argnames) - 1
        arg_name = argnames[i]

        if i == getter_spec.index_idx
            if result_spec.indexed
                push!(ccall_args, arg_name)
            else
                arg_type = argtypes[i]
                push!(ccall_args, :(Base.zero($arg_type)))
            end
        elseif i == getter_spec.attr_idx
            push!(ccall_args, key)
        else
            push!(ccall_args, arg_name)
        end
    end

    # output argument
    push!(ccall_args, var)

    return Expr(:ccall, ccall_args...)
end

function is_type_tuple(t)
    if !isa(t, Tuple)
        return false
    end

    for x in t
        if !isa(x, DataType)
            return false
        end
    end
    return true
end

function getter(c_name::Symbol, argnames, argtypes, map)
   for key in keys(map)
        jl_name = getter_name(key)
        return_type = map[key]

        setup_code, convert_code = getter_code_for(:value, return_type)

        getter_args = []

        push!(ccall_args, key, :value)

        getter_ccall = Expr(:ccall, ccall_args...)

        getter_code = Expr(:function,
                Expr(:call, symbol(jl_name), getter_args...),
                quote
                    local value

                    $setup_code

                    err = $getter_ccall

                    test_status(err)

                    $convert_code

                    return value
                end
            )

        # Debug output of generated code
        if false && endswith(jl_name, "max_dim")
            println(getter_code)
        end

        push!(obuf, getter_code)
    end
end

const getters = []

push!(getters,
    GetterSpec(
        :hsa_system_get_info,
        (:info, :data),
        (:hsa_system_info_t, Ptr{Void}),
        1, nothing,
        result_specs = Dict(
            :HSA_SYSTEM_INFO_VERSION_MAJOR => ResultSpec(Cushort),
            :HSA_SYSTEM_INFO_VERSION_MINOR => ResultSpec(Cushort),
            :HSA_SYSTEM_INFO_TIMESTAMP => ResultSpec(Culonglong),
            :HSA_SYSTEM_INFO_TIMESTAMP_FREQUENCY => ResultSpec(Cushort),
            :HSA_SYSTEM_INFO_SIGNAL_MAX_WAIT => ResultSpec(Culonglong),
            :HSA_SYSTEM_INFO_ENDIANNESS => ResultSpec(:hsa_endianness_t),
            :HSA_SYSTEM_INFO_MACHINE_MODEL => ResultSpec(:hsa_machine_model_t),
            :HSA_SYSTEM_INFO_EXTENSIONS => ResultSpec(UInt8, 128)
        )
    )
)

push!(getters,
    GetterSpec(
        :hsa_agent_get_info,
        (:agent, :info, :data),
        (:hsa_agent_t, :hsa_agent_info_t, Ptr{Void}),
        result_specs = Dict(
            :HSA_AGENT_INFO_NAME => ResultSpec(AbstractString,64),
            :HSA_AGENT_INFO_VENDOR_NAME => ResultSpec(AbstractString,64),
            :HSA_AGENT_INFO_FEATURE => ResultSpec(:hsa_agent_feature_t),
            :HSA_AGENT_INFO_MACHINE_MODEL => ResultSpec(:hsa_machine_model_t),
            :HSA_AGENT_INFO_PROFILE => ResultSpec(:hsa_profile_t),
            :HSA_AGENT_INFO_DEFAULT_FLOAT_ROUNDING_MODE => ResultSpec(:hsa_default_float_rounding_mode_t),
            :HSA_AGENT_INFO_BASE_PROFILE_DEFAULT_FLOAT_ROUNDING_MODES => ResultSpec(:hsa_default_float_rounding_mode_t),
            :HSA_AGENT_INFO_FAST_F16_OPERATION => ResultSpec(Bool),
            :HSA_AGENT_INFO_WAVEFRONT_SIZE => ResultSpec(UInt32),
            :HSA_AGENT_INFO_WORKGROUP_MAX_DIM => ResultSpec(Tuple{UInt16, UInt16, UInt16}),
            :HSA_AGENT_INFO_WORKGROUP_MAX_SIZE => ResultSpec(UInt32),
            :HSA_AGENT_INFO_GRID_MAX_DIM => ResultSpec(:hsa_dim3_t),
            :HSA_AGENT_INFO_GRID_MAX_SIZE => ResultSpec(UInt32),
            :HSA_AGENT_INFO_FBARRIER_MAX_SIZE => ResultSpec(UInt32),
            :HSA_AGENT_INFO_QUEUES_MAX => ResultSpec(UInt32),
            :HSA_AGENT_INFO_QUEUE_MIN_SIZE => ResultSpec(UInt32),
            :HSA_AGENT_INFO_QUEUE_MAX_SIZE => ResultSpec(UInt32),
            :HSA_AGENT_INFO_QUEUE_TYPE => ResultSpec(:hsa_queue_type_t),
            :HSA_AGENT_INFO_NODE => ResultSpec(UInt32),
            :HSA_AGENT_INFO_DEVICE => ResultSpec(:hsa_device_type_t),
            :HSA_AGENT_INFO_CACHE_SIZE => ResultSpec(Tuple{UInt32, UInt32, UInt32, UInt32}),
            :HSA_AGENT_INFO_ISA => ResultSpec(:hsa_isa_t),
            :HSA_AGENT_INFO_EXTENSIONS => ResultSpec(UInt8, 128),
            :HSA_AGENT_INFO_VERSION_MAJOR => ResultSpec(UInt16),
            :HSA_AGENT_INFO_VERSION_MINOR => ResultSpec(UInt16),
            :HSA_EXT_AGENT_INFO_IMAGE1D_MAX_DIM => ResultSpec(:hsa_dim3_t),
            :HSA_EXT_AGENT_INFO_IMAGE2D_MAX_DIM => ResultSpec(:hsa_dim3_t),
            :HSA_EXT_AGENT_INFO_IMAGE3D_MAX_DIM => ResultSpec(:hsa_dim3_t),
            :HSA_EXT_AGENT_INFO_IMAGE_ARRAY_MAX_SIZE => ResultSpec(UInt32),
            :HSA_EXT_AGENT_INFO_IMAGE_RD_MAX => ResultSpec(UInt32),
            :HSA_EXT_AGENT_INFO_IMAGE_RDWR_MAX => ResultSpec(UInt32),
            :HSA_EXT_AGENT_INFO_SAMPLER_MAX => ResultSpec(UInt32),
        )
    )
)

push!(getters,
    GetterSpec(
        :hsa_region_get_info,
        (:region, :attribute, :value),
        (:hsa_region_t, :hsa_region_info_t, Ptr{Void}),
        result_specs = Dict(
            :HSA_REGION_INFO_SEGMENT => ResultSpec(:hsa_region_segment_t),
            :HSA_REGION_INFO_GLOBAL_FLAGS => ResultSpec(:hsa_region_global_flag_t),
            :HSA_REGION_INFO_SIZE => ResultSpec(Csize_t),
            :HSA_REGION_INFO_ALLOC_MAX_SIZE => ResultSpec(Csize_t),
            :HSA_REGION_INFO_RUNTIME_ALLOC_ALLOWED => ResultSpec(Bool),
            :HSA_REGION_INFO_RUNTIME_ALLOC_GRANULE => ResultSpec(Csize_t),
            :HSA_REGION_INFO_RUNTIME_ALLOC_ALIGNMENT => ResultSpec(Csize_t),
        )
    )
)

push!(getters,
    GetterSpec(
        :hsa_isa_get_info,
        (:isa, :attribute, :index, :value),
        (:hsa_isa_t, :hsa_isa_info_t, UInt32, Ptr{Void}),
        result_specs = Dict(
            :HSA_ISA_INFO_NAME_LENGTH => ResultSpec(UInt32),
            :HSA_ISA_INFO_NAME => ResultSpec(AbstractString, :HSA_ISA_INFO_NAME_LENGTH)
        )
    )
)

push!(getters,
    GetterSpec(
        :hsa_executable_get_info,
        (:executable, :attribute, :value),
        (:hsa_executable_t, :hsa_executable_info_t, Ptr{Void}),
        result_specs = Dict(
            :HSA_EXECUTABLE_INFO_PROFILE => ResultSpec(:hsa_profile_t),
            :HSA_EXECUTABLE_INFO_STATE => ResultSpec(:hsa_executable_state_t)
        )
    )
)

push!(getters,
    GetterSpec(
        :hsa_executable_symbol_get_info,
        (:symbol, :attribute, :value),
        (:hsa_executable_symbol_t, :hsa_executable_symbol_info_t, Ptr{Void}),
        result_specs = Dict(
            :HSA_EXECUTABLE_SYMBOL_INFO_TYPE=> ResultSpec(:hsa_symbol_kind_t),
            :HSA_EXECUTABLE_SYMBOL_INFO_NAME_LENGTH => ResultSpec(UInt32),
            :HSA_EXECUTABLE_SYMBOL_INFO_NAME=> ResultSpec(AbstractString, :HSA_EXECUTABLE_SYMBOL_INFO_NAME_LENGTH),
            :HSA_EXECUTABLE_SYMBOL_INFO_MODULE_NAME_LENGTH => ResultSpec(UInt32),
            :HSA_EXECUTABLE_SYMBOL_INFO_MODULE_NAME=> ResultSpec(AbstractString, :HSA_EXECUTABLE_SYMBOL_INFO_MODULE_NAME_LENGTH),
            :HSA_EXECUTABLE_SYMBOL_INFO_AGENT => ResultSpec(:hsa_agent_t),
            :HSA_EXECUTABLE_SYMBOL_INFO_VARIABLE_ADDRESS=> ResultSpec(UInt64),
            :HSA_EXECUTABLE_SYMBOL_INFO_LINKAGE=> ResultSpec(:hsa_symbol_kind_linkage_t),
            :HSA_EXECUTABLE_SYMBOL_INFO_IS_DEFINITION => ResultSpec(Bool),
            :HSA_EXECUTABLE_SYMBOL_INFO_VARIABLE_ALLOCATION => ResultSpec(:hsa_variable_allocation_t),
            :HSA_EXECUTABLE_SYMBOL_INFO_VARIABLE_SEGMENT => ResultSpec(:hsa_variable_segment_t),
            :HSA_EXECUTABLE_SYMBOL_INFO_VARIABLE_ALIGNMENT => ResultSpec(UInt32),
            :HSA_EXECUTABLE_SYMBOL_INFO_VARIABLE_SIZE=> ResultSpec(UInt32),
            :HSA_EXECUTABLE_SYMBOL_INFO_VARIABLE_IS_CONST => ResultSpec(Bool),
            :HSA_EXECUTABLE_SYMBOL_INFO_KERNEL_OBJECT => ResultSpec(UInt64),
            :HSA_EXECUTABLE_SYMBOL_INFO_KERNEL_KERNARG_SEGMENT_SIZE => ResultSpec(UInt32),
            :HSA_EXECUTABLE_SYMBOL_INFO_KERNEL_KERNARG_SEGMENT_ALIGNMENT => ResultSpec(UInt32),
            :HSA_EXECUTABLE_SYMBOL_INFO_KERNEL_GROUP_SEGMENT_SIZE=> ResultSpec(UInt32),
            :HSA_EXECUTABLE_SYMBOL_INFO_KERNEL_PRIVATE_SEGMENT_SIZE=> ResultSpec(UInt32),
            :HSA_EXECUTABLE_SYMBOL_INFO_KERNEL_DYNAMIC_CALLSTACK=> ResultSpec(Bool),
            :HSA_EXECUTABLE_SYMBOL_INFO_INDIRECT_FUNCTION_OBJECT => ResultSpec(UInt32),
            :HSA_EXECUTABLE_SYMBOL_INFO_INDIRECT_FUNCTION_CALL_CONVENTION => ResultSpec(UInt32)
        )
    )
)

function gen_getters(obuf)
    for spec in getters
        for m_key in keys(spec.result_specs)
            method = Getter(spec, m_key)

            code = convert(Expr, method)

            push!(obuf, code)
        end
    end
end


text:
# # Constrained optimization

#=
In this example, we show how to differentiate through the solution of the following constrained optimization problem:
```math
\hat{y}(x) = \min_{y \in \mathcal{C}} f(x, y)
```
where ``\mathcal{C}`` is a closed convex set.
The optimal solution can be found as the fixed point of the projected gradient algorithm for any step size ``\eta``. This insight yields the following optimality conditions:
```math
F(x, \hat{y}(x)) = 0 \quad \text{with} \quad F(x,y) = \mathrm{proj}_{\mathcal{C}}(y - \eta \nabla_2 f(x, y)) - y
```
=#

using ImplicitDifferentiation
using Ipopt
using JuMP
using Zygote

using ChainRulesTestUtils  #src
using Test  #src

# ## Projecting onto the simplex

#=
We focus on minimizing ``f(x,y) = \lVert x - y \rVert_2^2``.
We also assume that ``\mathcal{C} = \Delta^n`` is the ``n``-dimensional probability simplex, because we know an exact procedure to compute the projection *and* its Jacobian.

Both of these procedures are outlined in [_From Softmax to Sparsemax: A Sparse Model of Attention and Multi-Label Classification_](https://arxiv.org/abs/1602.02068).
You can find an implementation in the function [`simplex_projection`](@ref).
Because this function involves a call to `sort`, standard AD backends cannot differentiate through it, which is why we also had to define a chain rule for it.
=#

# ## Implicit function wrapper

# We now wrap a black box optimizer inside an `ImplicitFunction` to compare its implicit differentiation with the explicit procedure given above.

function forward(x)
    n = length(x)
    optimizer = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)
    model = Model(optimizer)
    @variable(model, y[1:n] >= 0)
    @constraint(model, sum(y) == 1)
    @objective(model, Min, sum((y .- x) .^ 2))
    optimize!(model)
    return value.(y)
end;

conditions(x, y) = simplex_projection(y - 0.1 * 2(y - x)) - y;

implicit = ImplicitFunction(forward, conditions);

# ## Testing

# Let us study the behavior of our implicit function.

x = rand(5)

# We can see that the forward pass computes the projection correctly, at least up to numerical precision.

hcat(simplex_projection(x), implicit(x))

# And the same goes for the Jacobian.

cat(Zygote.jacobian(simplex_projection, x)[1], Zygote.jacobian(implicit, x)[1]; dims=3)

# The following tests are not included in the docs.  #src

@test implicit(x) ≈ simplex_projection(x) atol = 1e-5  #src

@testset verbose = true "ChainRules" begin  #src
    test_rrule(simplex_projection, x; check_inferred=false)  #src
    test_rrule(implicit, x; atol=1e-2, check_inferred=false)  #src
end  #src


text:
function ts_hbt(x::Neuro1DRealSignal)::Neuro1DAnalyticSignal
    r, c = size(x.signal)

    result::Array{ComplexF64} = zeros(Float64, r, c) .+ zeros(Float64, r, c)*im

    for i = 1:c
        result[:, i] = hilbert(x.signal[:, i])
    end

    Neuro1DAnalyticSignal(
        signal = result,
        sample_rate = x.sample_rate
    )
end


text:
function initialize!(integrator, cache::ABDF2ConstantCache)
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator, cache::ABDF2ConstantCache, repeat_step=false)
  @unpack t,f,p = integrator
  @unpack dtₙ₋₁,nlsolver = cache
  alg = unwrap_alg(integrator, true)
  dtₙ, uₙ, uₙ₋₁, uₙ₋₂ = integrator.dt, integrator.u, integrator.uprev, integrator.uprev2

  if integrator.iter == 1 && !integrator.u_modified
    cache.dtₙ₋₁ = dtₙ
    perform_step!(integrator, cache.eulercache, repeat_step)
    cache.fsalfirstprev = integrator.fsalfirst
    return
  end

  # precalculations
  fₙ₋₁ = integrator.fsalfirst
  ρ = dtₙ/dtₙ₋₁
  d = 2//3
  ddt = d*dtₙ
  dtmp = 1//3*ρ^2
  d1 = 1+dtmp
  d2 = -dtmp
  d3 = -(ρ-1)*1//3

  # calculate W
  update_W!(integrator, cache, ddt, repeat_step)

  zₙ₋₁ = dtₙ*fₙ₋₁
  # initial guess
  if alg.extrapolant == :linear
    z = dtₙ*fₙ₋₁
  else # :constant
    z = zero(uₙ)
  end
  nlsolver.z = z

  nlsolver.tmp = d1*uₙ₋₁ + d2*uₙ₋₂ + d3*zₙ₋₁
  z = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  uₙ = nlsolver.tmp + d*z
  integrator.fsallast = f(uₙ,p,t+dtₙ)
  integrator.destats.nf += 1

  if integrator.opts.adaptive
    tmp = integrator.fsallast - (1+dtₙ/dtₙ₋₁)*integrator.fsalfirst + (dtₙ/dtₙ₋₁)*cache.fsalfirstprev
    est = (dtₙ₋₁+dtₙ)/6 * tmp
    atmp = calculate_residuals(est, uₙ₋₁, uₙ, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

  ################################### Finalize

  if integrator.EEst < one(integrator.EEst)
    cache.fsalfirstprev = integrator.fsalfirst
    cache.dtₙ₋₁ = dtₙ
  end

  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = uₙ
  return
end

function initialize!(integrator, cache::ABDF2Cache)
  integrator.kshortsize = 2
  integrator.fsalfirst = cache.fsalfirst
  integrator.fsallast = cache.k
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # For the interpolation, needs k at the updated point
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator, cache::ABDF2Cache, repeat_step=false)
  @unpack t,dt,f,p = integrator
  @unpack z,k,b,J,W,tmp,atmp,dtₙ₋₁,zₙ₋₁,nlsolver = cache
  alg = unwrap_alg(integrator, true)
  uₙ,uₙ₋₁,uₙ₋₂,dtₙ = integrator.u,integrator.uprev,integrator.uprev2,integrator.dt

  if integrator.iter == 1 && !integrator.u_modified
    cache.dtₙ₋₁ = dtₙ
    perform_step!(integrator, cache.eulercache, repeat_step)
    cache.fsalfirstprev .= integrator.fsalfirst
    nlsolver.tmp = tmp
    return
  end

  # precalculations
  fₙ₋₁ = integrator.fsalfirst
  ρ = dtₙ/dtₙ₋₁
  d = 2//3
  ddt = d*dtₙ
  dtmp = 1//3*ρ^2
  d1 = 1+dtmp
  d2 = -dtmp
  d3 = -(ρ-1)*1//3

  update_W!(integrator, cache, ddt, repeat_step)

  @.. zₙ₋₁ = dtₙ*fₙ₋₁
  # initial guess
  if alg.extrapolant == :linear
    @.. z = dtₙ*fₙ₋₁
  else # :constant
    z .= zero(eltype(z))
  end

  @.. tmp = d1*uₙ₋₁ + d2*uₙ₋₂ + d3*zₙ₋₁
  z = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return

  @.. uₙ = tmp + d*z

  f(integrator.fsallast, uₙ, p, t+dtₙ)
  integrator.destats.nf += 1
  if integrator.opts.adaptive
    btilde0 = (dtₙ₋₁+dtₙ)*1//6
    btilde1 = 1+dtₙ/dtₙ₋₁
    btilde2 = dtₙ/dtₙ₋₁
    @.. tmp = btilde0*(integrator.fsallast - btilde1*integrator.fsalfirst + btilde2*cache.fsalfirstprev)
    calculate_residuals!(atmp, tmp, uₙ₋₁, uₙ, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

  ################################### Finalize

  if integrator.EEst < one(integrator.EEst)
    @.. cache.fsalfirstprev = integrator.fsalfirst
    cache.dtₙ₋₁ = dtₙ
  end
  return
end

# SBDF

function initialize!(integrator, cache::SBDFConstantCache)
  @unpack uprev, p, t = integrator
  @unpack f1, f2 = integrator.f
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  cache.du₁ = f1(uprev,p,t)
  cache.du₂ = f2(uprev,p,t)
  integrator.destats.nf += 1
  integrator.destats.nf2 += 1
  integrator.fsalfirst = cache.du₁ + cache.du₂

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

function perform_step!(integrator,cache::SBDFConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p,alg = integrator
  @unpack uprev2,uprev3,uprev4,du₁,du₂,k₁,k₂,k₃,nlsolver = cache
  @unpack f1, f2 = integrator.f
  cnt = cache.cnt = min(alg.order, integrator.iter+1)
  integrator.iter == 1 && !integrator.u_modified && ( cnt = cache.cnt = 1 )
  nlsolver.γ = γ = inv(γₖ[cnt])
  if cnt == 1
    tmp = uprev + dt*du₂
  elseif cnt == 2
    tmp = γ * (2*uprev - 1//2*uprev2 + dt*(2*du₂ - k₁))
  elseif cnt == 3
    tmp = γ * (3*uprev - 3//2*uprev2 + 1//3*uprev3 + dt*(3*(du₂ - k₁) + k₂))
  else
    tmp = γ * (4*uprev - 3*uprev2 + 4//3*uprev3 - 1//4*uprev4 + dt*(4*du₂ - 6*k₁ + 4*k₂ - k₃))
  end
  nlsolver.tmp = tmp

  # Implicit part
  # precalculations
  γdt = γ*dt
  update_W!(integrator, cache, γdt, repeat_step)

  # initial guess
  if alg.extrapolant == :linear
    z = dt*du₁
  else # :constant
    z = zero(u)
  end
  nlsolver.z = z

  z = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return
  u = nlsolver.tmp + γ*z

  cnt == 4 && ( cache.uprev4 = uprev3; cache.k₃ = k₂ )
  cnt >= 3 && ( cache.uprev3 = uprev2; cache.k₂ = k₁ )
              ( cache.uprev2 = uprev;  cache.k₁ = du₂ )
  cache.du₁ = f1(u, p, t+dt)
  cache.du₂ = f2(u, p, t+dt)
  integrator.destats.nf += 1
  integrator.destats.nf2 += 1
  integrator.fsallast = cache.du₁ + cache.du₂
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function initialize!(integrator, cache::SBDFCache)
  @unpack uprev, p, t = integrator
  @unpack f1, f2 = integrator.f
  integrator.kshortsize = 2
  integrator.fsalfirst = cache.fsalfirst
  integrator.fsallast = cache.k
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  f1(cache.du₁, uprev, p, t)
  f2(cache.du₂, uprev, p, t)
  integrator.destats.nf += 1
  integrator.destats.nf2 += 1
  @.. integrator.fsalfirst = cache.du₁ + cache.du₂
end

function perform_step!(integrator, cache::SBDFCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p,alg = integrator
  @unpack tmp,uprev2,uprev3,uprev4,k,k₁,k₂,k₃,du₁,du₂,z,nlsolver = cache
  @unpack f1, f2 = integrator.f
  cnt = cache.cnt = min(alg.order, integrator.iter+1)
  integrator.iter == 1 && !integrator.u_modified && ( cnt = cache.cnt = 1 )
  nlsolver.γ = γ = inv(γₖ[cnt])
  # Explicit part
  if cnt == 1
    @.. tmp = uprev + dt*du₂
  elseif cnt == 2
    @.. tmp = γ * (2*uprev - 1//2*uprev2 + dt*(2*du₂ - k₁))
  elseif cnt == 3
    @.. tmp = γ * (3*uprev - 3//2*uprev2 + 1//3*uprev3 + dt*(3*(du₂ - k₁) + k₂))
  else
    @.. tmp = γ * (4*uprev - 3*uprev2 + 4//3*uprev3 - 1//4*uprev4 + dt*(4*du₂ - 6*k₁ + 4*k₂ - k₃))
  end
  # Implicit part
  # precalculations
  γdt = γ*dt
  update_W!(integrator, cache, γdt, repeat_step)

  # initial guess
  if alg.extrapolant == :linear
    @.. z = dt*du₁
  else # :constant
    @.. z = zero(eltype(u))
  end

  z = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return
  @.. u = tmp + γ*z

  cnt == 4 && ( cache.uprev4 .= uprev3; cache.k₃ .= k₂ )
  cnt >= 3 && ( cache.uprev3 .= uprev2; cache.k₂ .= k₁ )
              ( cache.uprev2 .= uprev;  cache.k₁ .= du₂ )
  f1(du₁, u, p, t+dt)
  f2(du₂, u, p, t+dt)
  integrator.destats.nf += 1
  integrator.destats.nf2 += 1
  @.. k = du₁ + du₂
end

# QNDF1

function initialize!(integrator, cache::QNDF1ConstantCache)
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

function perform_step!(integrator,cache::QNDF1ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack uprev2,D,D2,R,U,dtₙ₋₁,nlsolver = cache
  alg = unwrap_alg(integrator, true)
  cnt = integrator.iter
  k = 1
  if cnt == 1
    κ = zero(alg.kappa)
  else
    κ = alg.kappa
    ρ = dt/dtₙ₋₁
    D[1] = uprev - uprev2   # backward diff
    if ρ != 1
      R!(k,ρ,cache)
      D[1] = D[1] * (R[1] * U[1])
    end
  end

  # precalculations
  γ₁ = 1//1
  γ = inv((1-κ)*γ₁)

  u₀ = uprev + D[1]
  ϕ = γ * (γ₁*D[1])
  nlsolver.tmp = u₀ - ϕ

  γdt = γ*dt
  update_W!(integrator, cache, γdt, repeat_step)

  # initial guess
  nlsolver.z = dt*integrator.fsalfirst
  nlsolver.γ = γ

  z = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return
  u = nlsolver.tmp + γ*z

  if integrator.opts.adaptive && integrator.success_iter > 0
    D2[1] = u - uprev
    D2[2] = D2[1] - D[1]
    utilde = (κ*γ₁ + inv(k+1)) * D2[2]
    atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
    if integrator.EEst > one(integrator.EEst)
      return
    end
  else
    integrator.EEst = one(integrator.EEst)
  end
  cache.dtₙ₋₁ = dt
  cache.uprev2 = uprev
  integrator.fsallast = f(u, p, t+dt)
  integrator.destats.nf += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function initialize!(integrator, cache::QNDF1Cache)
  integrator.kshortsize = 2
  integrator.fsalfirst = cache.fsalfirst
  integrator.fsallast = cache.k
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # For the interpolation, needs k at the updated point
  integrator.destats.nf += 1
end

function perform_step!(integrator,cache::QNDF1Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack uprev2,D,D2,R,U,dtₙ₋₁,tmp,z,W,utilde,atmp,nlsolver = cache
  alg = unwrap_alg(integrator, true)
  cnt = integrator.iter
  k = 1
  if cnt == 1
    κ = zero(alg.kappa)
  else
    κ = alg.kappa
    ρ = dt/dtₙ₋₁
    @.. D[1] = uprev - uprev2 # backward diff
    if ρ != 1
      R!(k,ρ,cache)
      @.. D[1] = D[1] * (R[1] * U[1])
    end
  end

  # precalculations
  γ₁ = 1//1
  nlsolver.γ = γ = inv((1-κ)*γ₁)
  @.. tmp = uprev + D[1] - γ * (γ₁*D[1])

  γdt = γ*dt
  update_W!(integrator, cache, γdt, repeat_step)

  # initial guess
  @.. z = dt*integrator.fsalfirst

  z = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return
  @.. u = tmp + γ*z

  if integrator.opts.adaptive && integrator.success_iter > 0
    @.. D2[1] = u - uprev
    @.. D2[2] = D2[1] - D[1]
    @.. utilde = (κ*γ₁ + inv(k+1)) * D2[2]
    calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
    if integrator.EEst > one(integrator.EEst)
      return
    end
  else
    integrator.EEst = one(integrator.EEst)
  end
  cache.dtₙ₋₁ = dt
  cache.uprev2 .= uprev
  f(integrator.fsallast, u, p, t+dt)
  integrator.destats.nf += 1
end

function initialize!(integrator, cache::QNDF2ConstantCache)
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

function perform_step!(integrator,cache::QNDF2ConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack uprev2,uprev3,dtₙ₋₁,dtₙ₋₂,D,D2,R,U,nlsolver = cache
  alg = unwrap_alg(integrator, true)
  cnt = integrator.iter
  k = 2
  if cnt == 1 || cnt == 2
    κ = zero(alg.kappa)
    γ₁ = 1//1
    γ₂ = 1//1
  elseif dtₙ₋₁ != dtₙ₋₂
    κ = alg.kappa
    γ₁ = 1//1
    γ₂ = 1//1 + 1//2
    ρ₁ = dt/dtₙ₋₁
    ρ₂ = dt/dtₙ₋₂
    D[1] = uprev - uprev2
    D[1] = D[1] * ρ₁
    D[2] = D[1] - ((uprev2 - uprev3) * ρ₂)
  else
    κ = alg.kappa
    γ₁ = 1//1
    γ₂ = 1//1 + 1//2
    ρ = dt/dtₙ₋₁
    # backward diff
    D[1] = uprev - uprev2
    D[2] = D[1] - (uprev2 - uprev3)
    if ρ != 1
      R!(k,ρ,cache)
      R .= R * U
      D[1] = D[1] * R[1,1] + D[2] * R[2,1]
      D[2] = D[1] * R[1,2] + D[2] * R[2,2]
    end
  end

  # precalculations
  nlsolver.γ = γ = inv((1-κ)*γ₂)
  u₀ = uprev + D[1] + D[2]
  ϕ = γ * (γ₁*D[1] + γ₂*D[2])
  nlsolver.tmp = u₀ - ϕ

  γdt = γ*dt
  update_W!(integrator, cache, γdt, repeat_step)

  # initial guess
  nlsolver.z = dt*integrator.fsalfirst

  z = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return
  u = nlsolver.tmp + γ*z

  if integrator.opts.adaptive
    if integrator.success_iter == 0
      integrator.EEst = one(integrator.EEst)
    elseif integrator.success_iter == 1
      utilde = (u - uprev) - ((uprev - uprev2) * dt/dtₙ₋₁)
      atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
      integrator.EEst = integrator.opts.internalnorm(atmp,t)
    else
      D2[1] = u - uprev
      D2[2] = D2[1] - D[1]
      D2[3] = D2[2] - D[2]
      utilde = (κ*γ₂ + inv(k+1)) * D2[3]
      atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
      integrator.EEst = integrator.opts.internalnorm(atmp,t)
    end
  end
  if integrator.EEst > one(integrator.EEst)
    return
  end

  cache.uprev3 = uprev2
  cache.uprev2 = uprev
  cache.dtₙ₋₂ = dtₙ₋₁
  cache.dtₙ₋₁ = dt
  integrator.fsallast = f(u, p, t+dt)
  integrator.destats.nf += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
  return
end

function initialize!(integrator, cache::QNDF2Cache)
  integrator.kshortsize = 2
  integrator.fsalfirst = cache.fsalfirst
  integrator.fsallast = cache.k
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # For the interpolation, needs k at the updated point
  integrator.destats.nf += 1
end

function perform_step!(integrator,cache::QNDF2Cache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack uprev2,uprev3,dtₙ₋₁,dtₙ₋₂,D,D2,R,U,tmp,utilde,atmp,W,nlsolver = cache
  alg = unwrap_alg(integrator, true)
  cnt = integrator.iter
  k = 2
  if cnt == 1 || cnt == 2
    κ = zero(alg.kappa)
    γ₁ = 1//1
    γ₂ = 1//1
  elseif dtₙ₋₁ != dtₙ₋₂
    κ = alg.kappa
    γ₁ = 1//1
    γ₂ = 1//1 + 1//2
    ρ₁ = dt/dtₙ₋₁
    ρ₂ = dt/dtₙ₋₂
    @.. D[1] = uprev - uprev2
    @.. D[1] = D[1] * ρ₁
    @.. D[2] = D[1] - ((uprev2 - uprev3) * ρ₂)
  else
    κ = alg.kappa
    γ₁ = 1//1
    γ₂ = 1//1 + 1//2
    ρ = dt/dtₙ₋₁
    # backward diff
    @.. D[1] = uprev - uprev2
    @.. D[2] = D[1] - (uprev2 - uprev3)
    if ρ != 1
      R!(k,ρ,cache)
      R .= R * U
      @.. D[1] = D[1] * R[1,1] + D[2] * R[2,1]
      @.. D[2] = D[1] * R[1,2] + D[2] * R[2,2]
    end
  end

  # precalculations
  nlsolver.γ = γ = inv((1-κ)*γ₂)
  @.. tmp = uprev + D[1] + D[2] - γ * (γ₁*D[1] + γ₂*D[2])

  γdt = γ*dt
  update_W!(integrator, cache, γdt, repeat_step)

  # initial guess
  @.. nlsolver.z = dt*integrator.fsalfirst

  z = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return
  @.. u = tmp + γ*z

  if integrator.opts.adaptive
    if integrator.success_iter == 0
      integrator.EEst = one(integrator.EEst)
    elseif integrator.success_iter == 1
      @.. utilde = (u - uprev) - ((uprev - uprev2) * dt/dtₙ₋₁)
      calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
      integrator.EEst = integrator.opts.internalnorm(atmp,t)
    else
      @.. D2[1] = u - uprev
      @.. D2[2] = D2[1] - D[1]
      @.. D2[3] = D2[2] - D[2]
      @.. utilde = (κ*γ₂ + inv(k+1)) * D2[3]
      calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
      integrator.EEst = integrator.opts.internalnorm(atmp,t)
    end
  end
  if integrator.EEst > one(integrator.EEst)
    return
  end

  cache.uprev3 .= uprev2
  cache.uprev2 .= uprev
  cache.dtₙ₋₂ = dtₙ₋₁
  cache.dtₙ₋₁ = dt
  f(integrator.fsallast, u, p, t+dt)
  integrator.destats.nf += 1
  return
end

function initialize!(integrator, cache::QNDFConstantCache)
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

function perform_step!(integrator,cache::QNDFConstantCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack udiff,dts,order,max_order,D,D2,R,U,nlsolver = cache
  k = order
  cnt = integrator.iter
  κ = integrator.alg.kappa[k]
  γ = inv((1-κ)*γₖ[k])
  flag = true
  for i in 2:k
    if dts[i] != dts[1]
      flag = false
      break
    end
  end
  if cnt > 2
    if flag
      ρ = dt/dts[1]
      # backward diff
      n = k+1
      if cnt == 3
        n = k
      end
      for i = 1:n
        D2[1,i] = udiff[i]
      end
      backward_diff!(cache,D,D2,k)
      if ρ != 1
        U!(k,U)
        R!(k,ρ,cache)
        R .= R * U
        reinterpolate_history!(cache,D,R,k)
      end
    else
      n = k+1
      if cnt == 3
        n = k
      end
      for i = 1:n
        D2[1,i] = udiff[i] * dt/dts[i]
      end
      backward_diff!(cache,D,D2,k)
    end
  else
    γ = 1//1
  end
  nlsolver.γ = γ
  # precalculations
  u₀ = uprev + sum(D)  # u₀ is predicted value
  ϕ = zero(γ)
  for i = 1:k
    ϕ += γₖ[i]*D[i]
  end
  ϕ *= γ
  nlsolver.tmp = u₀ - ϕ
  γdt = γ*dt
  update_W!(integrator, cache, γdt, repeat_step)
  # initial guess
  nlsolver.z = dt*integrator.fsalfirst

  z = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return
  u = nlsolver.tmp + γ*z

  if integrator.opts.adaptive
    if cnt == 1
      integrator.EEst = one(integrator.EEst)
    elseif cnt == 2
      utilde = (u - uprev) - (udiff[1] * dt/dts[1])
      atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
      integrator.EEst = integrator.opts.internalnorm(atmp,t)
    else
      δ = u - uprev
      for i = 1:k
        δ -= D[i]
      end
      utilde = (κ*γₖ[k] + inv(k+1)) * δ
      atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
      integrator.EEst = integrator.opts.internalnorm(atmp,t)
    end

    if cnt == 1
      cache.order = 1
    elseif cnt <= 3
      cache.order = 2
    else
      errm1 = 0
      if k > 1
        utildem1 = (κ*γₖ[k-1] + inv(k)) * D[k]
        atmpm1 = calculate_residuals(utildem1, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
        errm1 = integrator.opts.internalnorm(atmpm1, t)
      end
      backward_diff!(cache,D,D2,k+1,false)
      δ = u - uprev
      for i = 1:(k+1)
        δ -= D2[i,1]
      end
      utildep1 = (κ*γₖ[k+1] + inv(k+2)) * δ
      atmpp1 = calculate_residuals(utildep1, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
      errp1 = integrator.opts.internalnorm(atmpp1, t)
      pass = stepsize_and_order!(cache, integrator.EEst, errm1, errp1, dt, k)
      if pass == false
        cache.c = cache.c + 1
        fill!(D, zero(u)); fill!(D2, zero(u))
        fill!(R, zero(t)); fill!(U, zero(t))
        return
      end
      cache.c = 0
    end # cnt == 1
  end # integrator.opts.adaptive
  for i = 6:-1:2
    dts[i] = dts[i-1]
    udiff[i] = udiff[i-1]
  end
  dts[1] = dt
  udiff[1] = u - uprev
  fill!(D, zero(u)); fill!(D2, zero(u))
  fill!(R, zero(t)); fill!(U, zero(t))

  integrator.fsallast = f(u, p, t+dt)
  integrator.destats.nf += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function initialize!(integrator, cache::QNDFCache)
  integrator.kshortsize = 2
  integrator.fsalfirst = cache.fsalfirst
  integrator.fsallast = cache.k
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # For the interpolation, needs k at the updated point
  integrator.destats.nf += 1
end

function perform_step!(integrator,cache::QNDFCache,repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack udiff,dts,order,max_order,D,D2,R,U,tmp,utilde,atmp,W,nlsolver = cache
  cnt = integrator.iter
  k = order
  κ = integrator.alg.kappa[k]
  γ = inv((1-κ)*γₖ[k])
  flag = true
  for i in 2:k
    if dts[i] != dts[1]
      flag = false
      break
    end
  end
  if cnt > 2
    if flag
      ρ = dt/dts[1]
      # backward diff
      n = k+1
      if cnt == 3
        n = k
      end
      for i = 1:n
        D2[1,i] .= udiff[i]
      end
      backward_diff!(cache,D,D2,k)
      if ρ != 1
        U!(k,U)
        R!(k,ρ,cache)
        R .= R * U
        reinterpolate_history!(cache,D,R,k)
      end
    else
      n = k+1
      if cnt == 3
        n = k
      end
      for i = 1:n
        @.. D2[1,i] = udiff[i] * dt/dts[i]
      end
      backward_diff!(cache,D,D2,k)
    end
  else
    γ = one(γ)
  end
  nlsolver.γ = γ
  # precalculations
  ϕ = fill!(utilde, zero(eltype(u)))
  for i = 1:k
    @.. ϕ += γₖ[i]*D[i]
  end
  @.. ϕ *= γ
  tm = fill!(cache.b, zero(eltype(u)))
  for i = 1:k
    @.. tm += D[i]
  end
  @.. nlsolver.tmp = uprev + tm - ϕ

  γdt = γ*dt
  update_W!(integrator, cache, γdt, repeat_step)
  # initial guess
  @.. nlsolver.z = dt*integrator.fsalfirst

  z = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return
  @.. u = nlsolver.tmp + γ*z

  if integrator.opts.adaptive
    if cnt == 1
      integrator.EEst = one(integrator.EEst)
    elseif cnt == 2
      @.. utilde = (u - uprev) - (udiff[1] * dt/dts[1])
      calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
      integrator.EEst = integrator.opts.internalnorm(atmp,t)
    else
      @.. tmp = u - uprev
      for i = 1:k
        @.. tmp -= D[i]
      end
      @.. utilde = (κ*γₖ[k] + inv(k+1)) * tmp
      calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
      integrator.EEst = integrator.opts.internalnorm(atmp,t)
    end

    if cnt == 1
      cache.order = 1
    elseif cnt <= 3
      cache.order = 2
    else
      errm1 = 0
      if k > 1
        @.. utilde = (κ*γₖ[k-1] + inv(k)) * D[k]
        calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
        errm1 = integrator.opts.internalnorm(atmp,t)
      end
      backward_diff!(cache,D,D2,k+1,false)
      @.. tmp = u - uprev
      for i = 1:(k+1)
        @.. tmp -= D2[i,1]
      end
      @.. utilde = (κ*γₖ[k+1] + inv(k+2)) * tmp
      calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
      errp1 = integrator.opts.internalnorm(atmp,t)
      pass = stepsize_and_order!(cache, integrator.EEst, errm1, errp1, dt, k)
      if pass == false
        for i = 1:5
          fill!(D[i], zero(eltype(u)))
        end
        for i = 1:6
          for j = 1:6
            fill!(D2[i,j], zero(eltype(u)))
          end
        end
        fill!(R, zero(t)); fill!(U, zero(t))
        cache.c = cache.c + 1
        return
      end
      cache.c = 0
    end # cnt == 1
  end # integrator.opts.adaptive
  swap_tmp = udiff[6]
  for i = 6:-1:2
    dts[i] = dts[i-1]
    udiff[i] = udiff[i-1]
  end
  udiff[1] = swap_tmp
  dts[1] = dt
  @.. udiff[1] = u - uprev
  for i = 1:5
    fill!(D[i], zero(eltype(u)))
  end
  for i = 1:6, j = 1:6
    fill!(D2[i,j], zero(eltype(u)))
  end
  fill!(R, zero(t)); fill!(U, zero(t))

  f(integrator.fsallast, u, p, t+dt)
  integrator.destats.nf += 1
end

### MEBDF2
function initialize!(integrator, cache::MEBDF2ConstantCache)
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
  integrator.destats.nf += 1

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
end

@muladd function perform_step!(integrator, cache::MEBDF2ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  nlsolver = cache.nlsolver
  alg = unwrap_alg(integrator, true)
  update_W!(integrator, cache, dt, repeat_step)

  # initial guess
  if alg.extrapolant == :linear
    nlsolver.z = dt*integrator.fsalfirst
  else # :constant
    nlsolver.z = zero(u)
  end

### STEP 1
  nlsolver.tmp = uprev
  z = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return
  z₁ = nlsolver.tmp + z
### STEP 2
  nlsolver.tmp = z₁
  nlsolver.c = 2
  set_new_W!(nlsolver, false)
  z = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return
  z₂ = z₁ + z
### STEP 3
  tmp2 = 0.5uprev + z₁ - 0.5z₂
  nlsolver.tmp = tmp2
  nlsolver.c = 1
  set_new_W!(nlsolver, false)
  z = nlsolve!(integrator, cache)
  nlsolvefail(nlsolver) && return
  u = tmp2 + z

### finalize
  integrator.fsallast = f(u, p, t+dt)
  integrator.destats.nf += 1
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.u = u
end

function initialize!(integrator, cache::MEBDF2Cache)
  integrator.kshortsize = 2
  integrator.fsalfirst = cache.fsalfirst
  integrator.fsallast = cache.k
  resize!(integrator.k, integrator.kshortsize)
  integrator.k[1] = integrator.fsalfirst
  integrator.k[2] = integrator.fsallast
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # For the interpolation, needs k at the updated point
  integrator.destats.nf += 1
end

@muladd function perform_step!(integrator, cache::MEBDF2Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack z,z₁,z₂,tmp2,tmp,nlsolver = cache
  mass_matrix = integrator.f.mass_matrix
  alg = unwrap_alg(integrator, true)
  update_W!(integrator, cache, dt, repeat_step)

  # initial guess
  if alg.extrapolant == :linear
    @.. z = dt*integrator.fsalfirst
  else # :constant
    z .= zero(eltype(u))
  end

### STEP 1
 nlsolver.tmp = uprev
 z = nlsolve!(integrator, cache)
 nlsolvefail(nlsolver) && return
 @.. z₁ = uprev + z
### STEP 2
 nlsolver.tmp = z₁
 nlsolver.c = 2
 set_new_W!(nlsolver, false)
 z = nlsolve!(integrator, cache)
 nlsolvefail(nlsolver) && return
 @.. z₂ = z₁ + z
### STEP 3
 # z .= zero(eltype(u))
 @.. tmp2 = 0.5uprev + z₁ - 0.5z₂
 nlsolver.tmp = tmp2
 nlsolver.c = 1
 set_new_W!(nlsolver, false)
 z = nlsolve!(integrator, cache)
 nlsolvefail(nlsolver) && return
 @.. u = tmp2 + z

### finalize
 f(integrator.fsallast,u,p,t+dt)
 integrator.destats.nf += 1
end


text:
module utilities

export build_tree

#=
    build_tree()

    TODO: Make version for array and DataFrame
=#
"""
"""
function build_tree()

end

end # module

max_stars_repo_name: DopedCow/FPGrowth

max_issues_repo_name: DopedCow/FPGrowth

max_forks_repo_name: DopedCow/FPGrowth


text:
using Luxor

origin = Point(0, 0)
scale = 2
r = 90 * scale
c1 = Point(0, -r) + origin
c2 = Point(-r * 3^0.5 / 2, r/2) + origin
c3 = Point(0, r) + origin
c4 =  Point(r * 3^0.5 / 2, r/2) + origin
box_size = 15 * scale
circle_size = 30 * scale

box_node1 = Point(0,-box_size) + origin
box_node2 = Point(-box_size, box_size) + origin
box_node3 = Point(0, box_size) + origin
box_node4 = Point(box_size, box_size) + origin

c_node1 = c1 + Point(0, circle_size)
c_node2 = c2 + Point(circle_size*3^0.5/2, - circle_size / 2)
c_node3 = c3 + Point(0, -circle_size)
c_node4 = c4 + Point(-circle_size*3^0.5/2, - circle_size / 2)

arc_node1 = box_node2 + (c_node2 - box_node2) * 0.2
arc_node2 = box_node4 + (c_node4 - box_node4) * 0.2

arrow_width = 5 * scale
aha = pi/4

@png begin
    box(Point(-box_size, -box_size) + origin, Point(box_size, box_size) + origin, :fill)
    arc2r(origin, arc_node2, arc_node1, :stroke)
    setcolor(0.251, 0.388, 0.847)  # dark blue
    circle(c1, circle_size, :fill)
    arrow(c_node1, box_node1, linewidth=arrow_width, arrowheadangle =aha)
    setcolor(0.796, 0.235, 0.2)  # dark red
    circle(c2, circle_size, :fill)
    arrow(box_node2, c_node2, linewidth=arrow_width, arrowheadangle =aha)
    setcolor(0.22, 0.596, 0.149) # dark green
    circle(c3, circle_size, :fill)
    arrow(box_node3, c_node3, linewidth=arrow_width, arrowheadangle =aha)
    setcolor(0.584, 0.345, 0.698) # dark purple
    circle(c4, circle_size, :fill)
    arrow(box_node4, c_node4, linewidth=arrow_width, arrowheadangle =aha)
    setcolor(0., 0., 0.) # black

    end 640 640 "docs/src/assets/logo.png"


text:
# ---
# title: 581. Shortest Unsorted Continuous Subarray
# id: problem581
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/shortest-unsorted-continuous-subarray/description/>
# hidden: true
# ---
#
# Given an integer array `nums`, you need to find one **continuous subarray**
# that if you only sort this subarray in ascending order, then the whole array
# will be sorted in ascending order.
#
# Return  _the shortest such subarray and output its length_.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [2,6,4,8,10,9,15]
#     Output: 5
#     Explanation: You need to sort [6, 4, 8, 10, 9] in ascending order to make the whole array sorted in ascending order.
#
#
# **Example 2:**
#
#
#
#     Input: nums = [1,2,3,4]
#     Output: 0
#
#
# **Example 3:**
#
#
#
#     Input: nums = [1]
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 104`
#   * `-105 <= nums[i] <= 105`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using MethodAnalysis

# Analyze MethodInstance signatures and select those that seem at risk for being invalidated.
function atrisktype(@nospecialize(typ))
    # signatures like `convert(Vector, a)`, `foo(::Vararg{Synbol,N}) where N` do not seem to pose a problem
    isa(typ, TypeVar) && return false
    # isbits parameters are not a problem
    isa(typ, Type) || return false
    if isa(typ, UnionAll)
        typ = Base.unwrap_unionall(typ)
    end
    # Exclude signatures with Union{}
    typ === Union{} && return false
    isa(typ, Union) && return atrisktype(typ.a) | atrisktype(typ.b)
    # Type{T}: signatures like `convert(::Type{AbstractString}, ::String)` are not problematic, so mark Type as OK
    typ <: Type && return false
    if typ <: Tuple && length(typ.parameters) >= 1
        p1 = typ.parameters[1]
        # Constructor calls are not themselves a problem (any `convert`s they trigger might be, but those are covered)
        isa(p1, Type) && p1 <: Type && return false
        # convert(::Type{T}, ::S) where S<:T is not problematic
        if p1 === typeof(Base.convert) || p1 === typeof(Core.convert) || p1 === typeof(Core.Compiler.convert)
            p2, p3 = typ.parameters[2], typ.parameters[3]
            if isa(p2, Type)
                p2 = Base.unwrap_unionall(p2)
                if isa(p2, DataType) && length(p2.parameters) === 1
                    T = p2.parameters[1]
                    isa(p3, Type) && isa(T, Type) && p3 <: T && return false
                end
            end
        # `getindex`, `length`, etc are OK for various Tuple{T1,T2,...}
        elseif (p1 === typeof(Base.getindex) || p1 === typeof(Core.Compiler.getindex)) ||
               (p1 === typeof(Base.length)  || p1 === typeof(Core.Compiler.length)) ||
               (p1 === typeof(Base.isempty) || p1 === typeof(Core.Compiler.isempty)) ||
               (p1 === typeof(Base.iterate) || p1 === typeof(Core.iterate) || p1 === typeof(Core.Compiler.iterate))
            p2 = typ.parameters[2]
            if isa(p2, Type)
                p2 = Base.unwrap_unionall(p2)
                p2 <: Tuple && return false
            end
        # show(io::IO, x) is OK as long as typeof(x) is safe
        elseif p1 === typeof(Base.show)
            atrisktype(typ.parameters[2]) && return true
            length(typ.parameters) == 3 && atrisktype(typ.parameters[3]) && return true
            return false
        end
    end
    # Standard DataTypes
    isconcretetype(typ) && return false
    # ::Function args are excluded
    typ === Function && return false
    !isempty(typ.parameters) && (any(atrisktype, typ.parameters) || return false)
    return true
end

# A few tests
@assert  atrisktype(Tuple{typeof(==),Any,Any})
@assert  atrisktype(Tuple{typeof(==),Symbol,Any})
@assert  atrisktype(Tuple{typeof(==),Any,Symbol})
@assert !atrisktype(Tuple{typeof(==),Symbol,Symbol})
@assert !atrisktype(Tuple{typeof(convert),Type{Any},Any})
@assert !atrisktype(Tuple{typeof(convert),Type{AbstractString},AbstractString})
@assert !atrisktype(Tuple{typeof(convert),Type{AbstractString},String})
@assert  atrisktype(Tuple{typeof(convert),Type{String},AbstractString})
@assert !atrisktype(Tuple{typeof(map),Function,Vector{Any}})
@assert !atrisktype(Tuple{typeof(getindex),Dict{Union{String,Int},Any},Union{String,Int}})
@assert  atrisktype(Tuple{typeof(getindex),Dict{Union{String,Int},Any},Any})
@assert !atrisktype(Tuple{Type{BoundsError},Any,Any})
@assert  atrisktype(Tuple{typeof(sin),Any})
@assert !atrisktype(Tuple{typeof(length),Tuple{Any,Any}})

function collect_mis(m::Method)
    list = Core.MethodInstance[]
    visit(m) do item
        if isa(item, Core.MethodInstance)
            push!(list, item)
            return false
        end
        return true
    end
    return list
end

const mis = Dict{Method,Vector{Core.MethodInstance}}()
visit() do item
    if item isa Method
        m = item
        mis[m] = collect_mis(m)
        return false
    end
    return true
end

# Count # of backedges for MethodInstances with abstract types
const becounter = Dict{Core.MethodInstance,Int}()
visit() do item
    if item isa Core.MethodInstance
        if atrisktype(item.specTypes)
            becounter[item] = length(all_backedges(item))
        end
        return false
    end
    return true
end

prs = sort!(collect(becounter); by=last)

# Organize them by method instead

const mcounter = Dict{Method,Int}()
for (mi, c) in becounter
    oc = get(mcounter, mi.def, 0)
    mcounter[mi.def] = oc + c
end

mprs = sort!(collect(mcounter); by=last)

open("/tmp/methinstdata_$VERSION.log", "w") do io
    for (mi, c) in prs
        c == 0 && continue
        println(io, mi.specTypes=>c)
    end
end
open("/tmp/methdata_$VERSION.log", "w") do io
    for (m, c) in mprs
        c == 0 && continue
        println(io, m.sig=>c)
    end
end


text:
export shift, shift!

"""
    shift!(arr, shifts; kwargs...)

Shifts an array in-place. For real arrays it is based on `rfft`.
For complex arrays based on `fft`.
`shifts` can be non-integer, for integer shifts one should prefer
`circshift` or `ShiftedArrays.circshift` because a FFT-based methods
introduces numerical errors.

## kwargs...
* `fix_nyquist_frequency=false`: Fourier shifting of even-sized arrays is not revertible. However, if you did
    `shift(x, δ)` you can it revert by `shift(x, δ, fix_nyquist_frequency=true)`. This only works if `δ` is the same.
* `take_real=true`: For even-sized arrays we take by default the `real` part of the exponential phase at the Nyquist frequency.
    This satisfies the property of real valuedness and the aliasing of the Nyquist term.

## Memory Usage
Note that for complex arrays we can avoid any large memory allocations
because of `fft!`. For `rfft` there does not exist a usable implementation
yet, so for real arrays there might be a temporary larger memory usage.

# Examples
```jldoctest
julia> x = [1.0 2.0 3.0; 4.0 5.0 6.0]
2×3 Matrix{Float64}:
 1.0  2.0  3.0
 4.0  5.0  6.0

julia> shift!(x, (1, 2))
2×3 Matrix{Float64}:
 5.0  6.0  4.0
 2.0  3.0  1.0

julia> x = [0, 1.0, 0.0, 1.0]
4-element Vector{Float64}:
 0.0
 1.0
 0.0
 1.0

julia> shift!(x, 0.5)
4-element Vector{Float64}:
 0.49999999999999994
 0.5
 0.49999999999999994
 0.5
```
"""
function shift!(arr::AbstractArray{<:Complex, N}, shifts; soft_fraction=0, fix_nyquist_frequency=false, take_real=true) where {N}
    return shift_by_1D_FT!(arr, shifts; soft_fraction=soft_fraction, fix_nyquist_frequency=fix_nyquist_frequency, take_real=take_real)
end

function shift!(arr::AbstractArray{<:Real, N}, shifts; soft_fraction=0, fix_nyquist_frequency=false, take_real=true) where {N}
    return shift_by_1D_RFT!(arr, shifts; soft_fraction=soft_fraction, fix_nyquist_frequency=fix_nyquist_frequency, take_real=take_real)
end

"""
    shift(arr, shifts)

Returning a shifted array.
See [`shift!`](@ref shift!) for more details
"""
function shift(arr, shifts; soft_fraction=0, fix_nyquist_frequency=false, take_real=true)
    return shift!(copy(arr), shifts; soft_fraction=soft_fraction, fix_nyquist_frequency=fix_nyquist_frequency, take_real=take_real)
end

function soft_shift(freqs, shift, fraction=eltype(freqs)(0.1); corner=false)
    rounded_shift = round.(shift);
    if corner
        w = window_half_cos(size(freqs),border_in=2.0-2*fraction, border_out=2.0, offset=CtrCorner)
    else
        w = ifftshift_view(window_half_cos(size(freqs),border_in=1.0-fraction, border_out=1.0))
    end
    return cispi.(-freqs .* 2 .* (w .* shift + (1.0 .-w).* rounded_shift))
end

function shift_by_1D_FT!(arr::AbstractArray{<:Complex, N}, shifts; soft_fraction=0, take_real=false, fix_nyquist_frequency=false) where {N}
    for (d, shift) in pairs(shifts) # iterates of the dimension d using the corresponding shift
        if iszero(shift)
            continue
        end
        # better use reorient from NDTools here?
        freqs = reshape(fftfreq(size(arr, d)), ntuple(i -> 1, Val(d-1))..., size(arr,d))
        # allocates a 1D slice of exp values
        if iszero(soft_fraction)
            ϕ = cispi.(- freqs .* 2 .* shift)
        else
            ϕ = soft_shift(freqs, shift, soft_fraction)
        end
        # in even case, set one value to real
        if iseven(size(arr, d))
            s = size(arr, d) ÷ 2 + 1
            ϕ[s] = take_real ? real(ϕ[s]) : ϕ[s]
            invr = 1 / ϕ[s]
            invr = isinf(invr) ? 0 : invr
            ϕ[s] = fix_nyquist_frequency ? invr : ϕ[s]
        end
        # go to fourier space and apply ϕ
        fft!(arr, d)
        arr .*= ϕ
    end

    # do one single multi dimensional fft to go back to real space
    # but only over those dimensions where shift != 0
    dims = Int[]
    for (i, s) in pairs(shifts)
        if !iszero(s)
            push!(dims, i)
        end
    end
    ifft!(arr, dims)

    return arr
end

# the idea is the following:
# rfft(x, 1) -> exp shift -> fft(x, 2) -> exp shift ->  fft(x, 3) -> exp shift -> ifft(x, [2,3]) -> irfft(x, 1)
# So once we did a rft to shift something we can call the routine for complex arrays to shift
function shift_by_1D_RFT!(arr::AbstractArray{<:Real, N}, shifts; soft_fraction=0, fix_nyquist_frequency=false, take_real=true) where {T, N}
    for (d, shift) in pairs(shifts)
        if iszero(shift)
            continue
        end

        s = size(arr, d) ÷ 2 + 1
        freqs = reshape(fftfreq(size(arr, d))[1:s], ntuple(i -> 1, d-1)..., s)
        if iszero(soft_fraction)
            ϕ = cispi.(-freqs .* 2 .* shift)
        else
            ϕ = soft_shift(freqs, shift, soft_fraction, corner=true)
        end
        if iseven(size(arr, d))
            # take real and maybe fix nyquist frequency
            ϕ[s] = take_real ? real(ϕ[s]) : ϕ[s]
            invr = 1 / ϕ[s]
            invr = isinf(invr) ? 0 : invr
            ϕ[s] = fix_nyquist_frequency ? invr : ϕ[s]
        end
        p = plan_rfft(arr, d)

        arr_ft = p * arr
        arr_ft .*= ϕ
        # since we now did a single rfft dim, we can switch to the complex routine
        new_shifts = ntuple(i -> i ≤ d ? 0 : shifts[i], N)
         # workaround to mimic in-place rfft
        shift_by_1D_FT!(arr_ft, new_shifts; soft_fraction=soft_fraction, take_real=take_real, fix_nyquist_frequency=fix_nyquist_frequency)
        # go back to real space now and return because shift_by_1D_FT processed
        # the other dimensions already
        mul!(arr, inv(p), arr_ft)
        return arr # this breaks the for loop and finishes the algorithm
    end
    return arr
end

 #"""
 #    shift_by_ND(arr, shift_vector)
 #
 #Shifts an array by the shift_vector using sinc-interpolation via FFTs
 #`shift_vector`.
 #The vector by which to shift
 #```
 #"""
 #function shift_by_ND(arr::AbstractArray{T, N}, shift_vector) where {T, N}
 #    # for complex arrays we need a full FFT
 #    if T <: Complex
 #        arr_out = shift_by_FFT(arr, Tuple(shift_vector))
 #    else
 #        arr_out = shift_by_RFFT(arr, Tuple(shift_vector))
 #    end
 #    return arr_out
 #end
 #
 #
 #"""
 #    shift_by_RFFT(mat, new_size)
 #
 #Does shifting based on `rfft`. This function is called by `shift`.
 #"""
 #function shift_by_ND_RFT(mat, shift_vec) where {T}
 #    old_size=size(mat)
 #    irft(rft(mat).*exp.(
 #        dot.([shift_vec], idx(rft_size(mat),scale=ScaRFT,offset=CtrRFT)
 #        ).*(2im*pi)), size(mat,1))
 #end
 #
 #
 #function shift_by_ND_FT(mat, shift_vec) where {T}
 #    old_size=size(mat)
 #    ift(ft(mat).*exp.(
 #        dot.([shift_vec], idx(size(mat),scale=ScaFT,offset=CtrFT)
 #        ).*(2im*pi)))
 #end
 #


text:
using LPVSpectral
using Test, LinearAlgebra, Statistics, Random
using Plots, DSP
Random.seed!(0)
# write your own tests here
function generate_signal(f,w,N, modphase=false)
    x = sort(10rand(N)) # Sample points
    v = range(0, stop=1, length=N) # Scheduling variable

    # generate output signal
    # phase_matrix
    dependence_matrix = Float64[f[(i-1)%length(f)+1].(v) for v in v, i in eachindex(w)] # N x nw
    frequency_matrix = [cos(w*x -0.5modphase*(dependence_matrix[i,j])) for (i,x) in enumerate(x), (j,w) in enumerate(w)] # N x nw
    y = sum(dependence_matrix.*frequency_matrix,dims=2)[:] # Sum over all frequencies
    y += 0.1randn(size(y))
    y,v,x,frequency_matrix, dependence_matrix
end

@testset "LPV methods" begin

    N      = 500 # Number of training data points
    f      = [v->2v^2, v->2/(5v+1), v->3exp(-10*(v-0.5)^2),] # Functional dependences on the scheduling variable
    w      = 2pi*[2,10,20] # Frequency vector
    w_test = 2π*collect(2:2:25)
    Y,V,X,frequency_matrix, dependence_matrix = generate_signal(f,w,N,true)

    λ      = 0.02 # Regularization parmater
    normal = true # Use normalized basis functions
    Nv     = 50 # Number of basis functions

    se = ls_spectral_lpv(Y,X,V,w_test,Nv; λ = λ, normalize = normal) # Perform LPV spectral estimation
    windowpsd = ls_windowpsd_lpv(Y,X,V,w_test,Nv; λ = λ, normalize = normal)

    spectrum_lpv   = psd(se)

    si = sortperm(spectrum_lpv[:],rev=true)
    @test Set(si[1:3]) == Set([1,5,10])

    si = sortperm(windowpsd[:],rev=true)
    @test Set(si[1:3]) == Set([1,5,10])
end

@testset "detrend" begin
    tre = [1,2,3]
    @test detrend(tre) == [-1,0,1]
    @test tre == [1,2,3] # Shouldn't have changed
    detrend!(tre)
    @test tre == [-1,0,1] # Should have changed
end

## Test ComplexNormal
@testset "Complex Normal" begin
    n  = 10
    n2 = 5
    a  = randn(n)
    A  = randn(n,n)
    A  = A'A + I
    b  = randn(n2)
    B  = randn(n2,n2)
    B  = B'B + I
    X  = randn(n,n2)
    Y  = randn(n,n2)
    cn = ComplexNormal(X,Y)
    @test size(cn.m) == (n2,)
    @test size(cn.Γ) == (n2,n2)
    @test size(cn.C) == (n2,n2)
    # affine_transform(cn, B, b)
    pdf(cn,b)
    pdf(cn,im*b)

    @test isa(cn_Vxx(A,A), Cholesky)
    @test isa(cn_fVxx(A,A), Matrix)
    @test issymmetric(cn_fVxx(A,A))
    cn_V(A,0.1A)

    cn  = ComplexNormal(a,A)
    cn  = ComplexNormal(im*b,A)

    A   = randn(4,4);
    A   = A'A
    x   = randn(1000,3)
    y   = randn(1000,3)
    cn  = ComplexNormal(x,y)
    z   = rand(cn,1000000);
    cn2 = ComplexNormal(z)
    @test norm(Matrix(cn.Γ)-Matrix(cn2.Γ)) < 0.01
    @test norm(Matrix(cn.C)-Matrix(cn2.C)) < 0.01

end

@testset "ls methods" begin
    T = 1000
    t = 0:T-1
    f = LPVSpectral.default_freqs(t)
    @test f[1] == 0
    @test f[end] == 0.5-1/length(t)/2
    @test length(f) == T÷2

    # f2 = LPVSpectral.default_freqs(t,10)
    # @test f2[1] == 0
    # @test f2[end] == 0.5-1/length(t)/2
    # @test length(f2) == T÷2÷10

    @test LPVSpectral.check_freq(f) == 1
    @test_throws ArgumentError LPVSpectral.check_freq([1,0,2])
    A, z = LPVSpectral.get_fourier_regressor(t,f)
    @test size(A) == (T,2length(f)-1)

    Base.isapprox(t1::Tuple{Float64,Int64}, t2::Tuple{Float64,Int64}; atol) = all(t -> isapprox(t[1],t[2],atol=atol), zip(t1,t2))
    y = sin.(t)
    x,_ = ls_spectral(y,t)
    @test findmax(abs.(x)) ≈ (0.9999773730281, 160) atol=0.001

    W = ones(length(y))
    x,_ = ls_spectral(y,t,f,W)
    @test findmax(abs.(x)) ≈ (0.999977373027, 160) atol=0.001

    x,_ = tls_spectral(y,t)
    @test findmax(abs.(x)) ≈ (0.9999777508878254, 160) atol=0.001

    x,_ = ls_windowpsd(y,t; nw=20)
    @test findmax(abs.(x)) ≈ (1.000783557456378, 160) atol=0.001

    x,_ = ls_windowpsd(y,t)
    @test findmax(abs.(x)) ≈ (1.0011490769234443, 160) atol=0.001

    x,_ = ls_windowcsd(y,y,t)
    @test findmax(abs.(x)) ≈ (1.0011490769234443, 160) atol=0.001

    x,_ = ls_cohere(y,y,t)
    @test findmax(abs.(x))[1] ≈ 1 atol=0.001

    x,_ = ls_cohere(y,y .+ 10randn.(),t)
    @test mean(abs.(x)) ≈ 0.5 atol = 0.15

end

@testset "plots" begin

    y = randn(1000)
    plot(periodogram(y))
    plot(periodogram(filtfilt(ones(4), [4], y)))
    plot(welch_pgram(y))

end

include("test_lasso.jl")


text:
"""
What element(s) is this material suitable as a standard for?
"""
struct DBStandardFor
    database::SQLite.DB
    element::Element
    material::Material
end

function Base.write(db::SQLite.DB, ::Type{DBStandardFor}, elm::Element, mat::Union{String, Material, Integer})
    if mat isa String
        matkey = find(db, Material, mat)
        mat = read(db, Material, matkey)
    elseif mat isa Material
        matkey = write(db, mat)
    else
        matkey = mat
        mat = read(db, Material, matkey)
    end
    @assert haskey(mat, elm) "The material $mat does not contain the element $elm."
    stmt1 = SQLite.Stmt(db, "SELECT * FROM STANDARDFOR WHERE ELEMENT=? AND MATKEY=?;")
    q1 = DBInterface.execute(stmt1, ( z(elm), matkey, ) )
    if SQLite.done(q1)
        stmt1 = SQLite.Stmt(db, "INSERT INTO STANDARDFOR ( ELEMENT, MATKEY ) VALUES ( ?, ? );")
        DBInterface.execute(stmt1, ( z(elm), matkey, ) )
    end
    return lastrowid(q1)
end

function Base.findall(db::SQLite.DB, ::Type{DBStandardFor}, elm::Element)::Array{DBMaterial}
    stmt1 = SQLite.Stmt(db, "SELECT * FROM STANDARDFOR WHERE ELEMENT=?;")
    q1 = DBInterface.execute(stmt1, ( z(elm), ) )
    return [ read(db, r[:MATKEY]) for r in q1 ]
end


text:
function start!(r)
    move!(r,Nord)
    move!(r,West)
end
function main!(r::Robot)
    k = 0
    while !ismarker(r)
        start!(r)
        k = k + 2
        a = k
        while a > 0
            move!(r, Ost)
            if ismarker(r)
                break
            end
            a = a - 1
        end
        if ismarker(r)
            break
        end
        a = k
        while a > 0
            move!(r,Sud)
            if ismarker(r)
                break
            end
            a = a - 1
        end
        if ismarker(r)
            break
        end
        a = k
        while a > 0
            move!(r, West)
            if ismarker(r)
                break
            end
            a = a - 1
        end
        if ismarker(r)
            break
        end
        a = k
        while a > 0
            move!(r, Nord)
            if ismarker(r)
                break
            end
            a = a - 1
        end
        if ismarker(r)
            break
        end
    end
end


text:
using Coulter
using Test
using Dates
using Distributions
using KernelDensity
using StatsBase
using Random

@testset "Loading" begin
    data = loadZ2("testdata/b_0_1.=#Z2", "blank")

    @test data.timepoint == DateTime(2017, 12, 22, 17, 46, 58)
    @test all(data.binheights .== 0.0)

    run = loadZ2("testdata/lat0um_0_1.=#Z2")

    @test mode(run.data) == run.binvols[findmax(run.binheights)[2]]
    @test minimum(run.data) == run.binvols[findfirst(!iszero, run.binheights)]
    @test maximum(run.data) == run.binvols[findlast(!iszero, run.binheights)]

    run2 = loadZ2("testdata/wt_0min_0_1.=#Z2", "wt")

    # correct values from Coulter software

    # make sure the common volume is in position 76
    @test findmax(run2.binheights)[2] == 76
    @test all(run2.binheights[75:77] .== [281, 326, 289])

    # load in total volume mode where the y axis is now counts*volume
    # this prevents overweighting of small volume objects
    run3 = loadZ2("testdata/wt_0min_0_1.=#Z2", "wt"; yvariable=:volume)

    # the coulter rounds to 4 digits
    @test all(map(x->round(x, digits=-2), run3.binheights[75:77]) .== [113.7e3, 133.6e3, 120.0e3])

    # test loading whole folders
    runs = Coulter.load_folder("testdata/")

    @test "wt" in keys(runs)
    @test sum(runs["wt"][1].data) == sum(run3.data)
end

@testset "Analysis" begin
    @testset "Peak-finding" begin
        # flat line, no peaks
        ys = fill(0.1, 20)
        xs = collect(1.0:20)

        @test length(Coulter._find_peaks(xs, ys, minx=0.0)) == 0

        # add two peaks at 10.5 and 16
        ys[1] = 0.0
        ys[10] = 0.2
        ys[11] = 0.2
        ys[16] = 0.4
        ys[17] = 0.3
        ys[18] = 0.1
        ys[20] = 0.0

        @test all(Coulter._find_peaks(xs, ys, minx=0.0) .== [10.5, 16.0])

        # more realistic example
        dist = MixtureModel([Normal(8.2, 0.75), Normal(9.6, 0.5)], [0.6, 0.4])
        xs = range(7, stop=17, length=400)
        ys = pdf.(dist, xs)

        Random.seed!(1234)
        sim_data = volume.(rand(dist, 10000))
        Random.seed!()

        kd_est = kde(sim_data)
        peaks = Coulter._find_peaks(collect(kd_est.x), kd_est.density, minx=0.0)
        @test all((peaks .- [289.325, 446.655]) .< 0.01)

        @test Coulter.extract_peak(sim_data) ≈ 446.654739
    end
end

@testset "Misc" begin
    @test volume(0.0) == 0.0
    @test volume(10.0) ≈ 523.59878
    @test diameter(0.0) == 0.0
    @test diameter(523.59878) ≈ 10.0
end


text:
"""
    quick_two_sum(a, b)

Computes `s = fl(a+b)` and `e = err(a+b)`. Assumes `|a| ≥ |b|`.
"""
@inline function quick_two_sum(a::Float64, b::Float64)
    s = a + b
    e = b - (s - a)
    s, e
end

"""
    two_sum(a, b)

Computes `s = fl(a+b)` and `e = err(a+b)`.
"""
@inline function two_sum(a::Float64, b::Float64)
    s = a + b
    v = s - a
    e = (a - (s - v)) + (b - v)

    s, e
end

"""
    quick_two_diff(a, b)

Computes `s = fl(a-b)` and `e = err(a-b)`.  Assumes `|a| ≥ |b|`.
"""
@inline function quick_two_diff(a::Float64, b::Float64)
    s = a - b
    e = (a - s) - b
    s, e
end

"""
    two_diff(a, b)

Computes `s = fl(a-b)` and `e = err(a-b)`.
"""
@inline function two_diff(a::Float64, b::Float64)
    s = a - b
    v = s - a
    e = (a - (s - v)) - (b + v)

    s, e
end

"""
    two_prod(a, b)

Computes `s = fl(a*b)` and `e = err(a*b)`.
"""
@inline function two_prod(a::Float64, b::Float64)
    p = a * b
    e = fma(a, b, -p)
    p, e
end

"""
    two_square(a)

Computes `s = fl(a*a)` and `e = err(a*a)`. Faster than [`two_prod(a, a)`](@ref).
"""
@inline function two_square(a::Float64)
    p = a * a
    e = fma(a, a, -p)
    p, e
end