text:
module TestDataMatrixs
    using Test
    using DataArrays

    # Base.getindex(d::DataMatrix, i::SingleIndex, j::SingleIndex)
    dm = @data([1 2; missing 4])
    dm[1, 2]

    # Base.getindex(x::DataMatrix, i::SingleIndex, col_inds::AbstractDataVector{Bool})
    dm = @data([1 2; missing 4])
    dm[1, @data([true, false])]

    # Base.getindex(x::DataMatrix, i::SingleIndex, col_inds::AbstractDataVector)
    dm = @data([1 2; missing 4])
    dm[1, @data([1, 2])]

    # Base.getindex(x::DataMatrix, i::SingleIndex, col_inds::MultiIndex)
    dm = @data([1 2; missing 4])
    dm[1, [1, 2]]

    # Base.getindex(x::DataMatrix, row_inds::AbstractDataVector{Bool}, j::SingleIndex)
    dm = @data([1 2; missing 4])
    dm[@data([true, false]), 1]

    # Base.getindex(x::DataMatrix, row_inds::AbstractDataVector, j::SingleIndex)
    dm = @data([1 2; missing 4])
    dm[@data([1, 2]), 1]

    # Base.getindex(x::DataMatrix, row_inds::MultiIndex, j::SingleIndex)
    dm = @data([1 2; missing 4])
    dm[[1, 2], 1]

    # Base.getindex(x::DataMatrix, row_inds::AbstractDataVector{Bool}, col_inds::AbstractDataVector{Bool})
    dm = @data([1 2; missing 4])
    dm[@data([true, false]), @data([true, false])]

    # Base.getindex(x::DataMatrix, row_inds::AbstractDataVector{Bool}, col_inds::AbstractDataVector)
    dm = @data([1 2; missing 4])
    dm[@data([true, false]), @data([1, 2])]

    # Base.getindex(x::DataMatrix, row_inds::AbstractDataVector{Bool}, col_inds::MultiIndex)
    dm = @data([1 2; missing 4])
    dm[@data([true, false]), [1, 2]]

    # Base.getindex(x::DataMatrix, row_inds::AbstractDataVector, col_inds::AbstractDataVector{Bool})
    dm = @data([1 2; missing 4])
    dm[@data([1, 2]), @data([true, false])]

    # Base.getindex(x::DataMatrix, row_inds::AbstractDataVector, col_inds::AbstractDataVector)
    dm = @data([1 2; missing 4])
    dm[@data([1, 2]), @data([1, 2])]

    # Base.getindex(x::DataMatrix, row_inds::AbstractDataVector, col_inds::MultiIndex)
    dm = @data([1 2; missing 4])
    dm[@data([1, 2]), [1, 2]]

    # Base.getindex(x::DataMatrix, row_inds::MultiIndex, col_inds::AbstractDataVector{Bool})
    dm = @data([1 2; missing 4])
    dm[[1, 2], @data([true, false])]

    # Base.getindex(x::DataMatrix, row_inds::MultiIndex, col_inds::AbstractDataVector)
    dm = @data([1 2; missing 4])
    dm[[1, 2], @data([1, 2])]

    # Base.getindex(x::DataMatrix, row_inds::MultiIndex, col_inds::MultiIndex)
    dm = @data([1 2; missing 4])
    dm[[1, 2], [1, 2]]

    # Base.setindex!(dm::DataMatrix, val::Missing, i::SingleIndex, j::SingleIndex)
    dm = @data([1 2; missing 4])
    dm[1, 1] = missing

    # Base.setindex!(dm::DataMatrix, val::Any, i::SingleIndex, j::SingleIndex)
    dm = @data([1 2; missing 4])
    dm[1, 1] = 3

    # Base.setindex!(dm::DataMatrix, val::Missing, row_inds::MultiIndex, j::SingleIndex)
    dm = @data([1 2; missing 4])
    dm[[1, 2], 1] = 3

    # Base.setindex!{S, T}(dm::DataMatrix{S}, vals::Vector{T}, row_inds::MultiIndex, j::SingleIndex)
    dm = @data([1 2; missing 4])
    dm[[1, 2], 1] = [3, 4]

    # Base.setindex!(dm::DataMatrix, val::Any, row_inds::MultiIndex, j::SingleIndex)
    dm = @data([1 2; missing 4])
    dm[[1, 2], 1] = 3

    # Base.setindex!(dm::DataMatrix, val::Missing, i::SingleIndex, col_inds::MultiIndex)
    dm = @data([1 2; missing 4])
    dm[[1, 2], 1] = missing

    # Base.setindex!{S, T}(dm::DataMatrix{S}, vals::Vector{T}, i::SingleIndex, col_inds::MultiIndex)
    dm = @data([1 2; missing 4])
    dm[1, [1, 2]] = [3, 4]

    # Base.setindex!(dm::DataMatrix, val::Any, i::SingleIndex, col_inds::MultiIndex)
    dm = @data([1 2; missing 4])
    dm[1, [1, 2]] = 3

    # Base.setindex!(dm::DataMatrix, val::Missing, row_inds::MultiIndex, col_inds::MultiIndex)
    dm = @data([1 2; missing 4])
    dm[[1, 2], [1, 2]] = missing

    # Base.setindex!{S, T}(dm::DataMatrix{S}, vals::Vector{T}, row_inds::MultiIndex, col_inds::MultiIndex)
    dm = @data([1 2; missing 4])
    dm[[1, 2], [1, 2]] = [1, 2, 3, 4]

    # Base.setindex!(dm::DataMatrix, val::Any, row_inds::MultiIndex, col_inds::MultiIndex)
    dm = @data([1 2; missing 4])
    dm[[1, 2], [1, 2]] = 5

    # Base.diag{T}(dm::DataMatrix{T})
    dm = @data([1 2; missing 4])
    diag(dm)
end


text:
using DistributedReinforcementLearning
using Test

@testset "DistributedReinforcementLearning.jl" begin
    # Write your tests here.
end

ext: jl


text:
using SolverLogging
using Test

@testset "Basics" begin
    include("basic.jl")
    include("utils.jl")
end

@testset "Macros" begin
    include("macros.jl")
end


text:
__precompile__()

module AxisArrayVariables

using JuMP
using AxisArrays
using MacroTools

export @axis_variable

# Keyword arguments use :(=) on v0.6 and :kw on v0.5, and their escaping rules
# are slightly different:
@static if VERSION >= v"0.6-"
    make_kwarg(key, val) = Expr(:(=), esc(key), esc(val))
else
    make_kwarg(key, val) = Expr(:kw, key, esc(val))
end

macro axis_variable(m, varexpr, args...)
    jump_args = Expr[esc(e) for e in args]
    (var, axis_args, lb, ub) = @match varexpr begin
        (var_[axis_args__] |
        (var_[axis_args__] <= ub_) |
        (var_[axis_args__] >= lb_) |
        (lb_ <= var_[axis_args__] <= ub_)) => (var, axis_args, lb, ub)
    end
    if var === nothing || axis_args === nothing
        error("Unrecognized expression")
    end
    if lb !== nothing
        unshift!(jump_args, make_kwarg(:lowerbound, lb))
    end
    if ub !== nothing
        unshift!(jump_args, make_kwarg(:upperbound, ub))
    end
    axes = []
    domains = Expr[]
    for arg in axis_args
        if @capture(arg, name_ = domain_)
            push!(axes, Expr(:call, Expr(:curly, :Axis, QuoteNode(name)), esc(domain)))
            push!(domains, esc(domain))
        else
            push!(axes, esc(arg))
            push!(domains, esc(arg))
        end
    end
    range_exprs = [:(1:length($a)) for a in domains]
    quote
        ranges = $(Expr(:vect, range_exprs...))
        vars = let
            local $(esc(var))
            $(Expr(:macrocall, Symbol("@variable"), esc(m), Expr(:ref, esc(var), range_exprs...), jump_args...))
        end
        $(esc(var)) = $(Expr(:call, :AxisArray, :vars, axes...))
    end
end

end


text:
using Random
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities;

using OSQP

n = 20 # variable dimension
m = 15 # no of inequality constraints

x̂ = rand(n)
Q = rand(n, n)
Q = Q' * Q # ensure PSD
q = rand(n)
G = rand(m, n)
h = G * x̂ + rand(m);

model = MOI.instantiate(OSQP.Optimizer, with_bridge_type=Float64)
x = MOI.add_variables(model, n);

# define objective

quad_terms = MOI.ScalarQuadraticTerm{Float64}[]
for i in 1:n
    for j in i:n # indexes (i,j), (j,i) will be mirrored. specify only one kind
        push!(quad_terms, MOI.ScalarQuadraticTerm(Q[i,j],x[i],x[j]))
    end
end

objective_function = MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm.(q, x),quad_terms,0.0)
MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(), objective_function)
MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

# maintain constrain to index map - will be useful later
constraint_indices = [
    MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(G[i,:], x), 0.),MOI.LessThan(h[i]),
    ) for i in 1:m
]

MOI.optimize!(model)

@assert MOI.get(model, MOI.TerminationStatus()) in (MOI.LOCALLY_SOLVED, MOI.OPTIMAL)

x̄ = MOI.get(model, MOI.VariablePrimal(), x)

# objective value (predicted vs actual) sanity check
@test 0.5*x̄'*Q*x̄ + q'*x̄  <= 0.5*x̂'*Q*x̂ + q'*x̂

# NOTE: can't use Ipopt
# Ipopt.Optimizer doesn't supports accessing MOI.ObjectiveFunctionType

#
# Verifying KKT Conditions
#

# complimentary slackness  + dual feasibility
for i in 1:size(constraint_indices)[1]
    con_index = constraint_indices[i]
    μ_i = MOI.get(model, MOI.ConstraintDual(), con_index)

    # μ[i] * (G * x - h)[i] = 0
    @test abs(μ_i * (G[i,:]' * x̄ - h[i])) < 3e-2

    # μ[i] <= 0
    @test μ_i <= 1e-2
end

# checking stationarity
for j in 1:n
    G_mu_sum = 0

    for i in 1:size(constraint_indices)[1]
        con_index = constraint_indices[i]
        μ_i = MOI.get(model, MOI.ConstraintDual(), con_index)

        G_mu_sum += μ_i * G[i,j]
    end

    @test abs(G_mu_sum - (Q * x̄ + q)[j]) < 1e-2
end


text:
include("../vb-cb/plot_vb_results.jl")

using Distributed
addprocs(10)

@everywhere include("../vb-cb/plot_vb_results.jl")

# Path to simulation data
path_to_data = ARGS[1]
results_dir = ARGS[2]

@everywhere path_to_data = $(path_to_data)

# Read CB data
@everywhere y = Matrix.(BSON.load(path_to_data)[:y])

paths = [joinpath(root, file)
         for (root, dirs, files) in walkdir(results_dir)
         for file in files
         if file == "output.bson"]

@everywhere function makegraph(path)
  println("Processing: $(path)"); flush(stdout)
  makeplots(y=y, path_to_output=path, markernames=[], nlam=30)
  return
end

# Make graphs in parallel
status = pmap(makegraph, paths, on_error=identity)
println(status)

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch

max_forks_repo_name: luiarthur/CytofResearch


text:
struct QNVars{T, Ty}
    d::T # search direction
    s::T # change in x
    y::Ty # change in successive gradients
end
function QNVars(x, g)
    QNVars(copy(g), copy(x), copy(x))
end
function preallocate_qn_caches(mstyle, x0)
    if mstyle === InPlace()
        # Maintain gradient and state pairs in QNVars
        cache = QNVars(x0, x0)
        return cache
    else
        return nothing
    end
end

function solve(problem::OptimizationProblem, x0, approach::LineSearch, options::OptimizationOptions, cache=preallocate_qn_caches(mstyle(problem), x0))
    _solve(mstyle(problem), problem, (x0, nothing), approach, options, cache)
end
function solve(problem::OptimizationProblem, s0::Tuple, approach::LineSearch, options::OptimizationOptions, cache=preallocate_qn_caches(mstyle(problem), first(s0)))
    _solve(mstyle(problem), problem, s0, approach, options, cache)
end

function _solve(mstyle, problem::OptimizationProblem, s0::Tuple, approach::LineSearch, options::OptimizationOptions, cache)
    t0 = time()

    #==============
         Setup
    ==============#
    x0, B0 = s0
    T = eltype(x0)

    objvars = prepare_variables(problem, approach, x0, copy(x0), B0)
    P = initial_preconditioner(approach, x0)
    f0, ∇f0 = objvars.fz, norm(objvars.∇fz, Inf) # use user norm

    if any(initial_converged(approach, objvars, ∇f0, options))
        x, fx, ∇fx, z, fz, ∇fz, B, Pg = objvars
        return ConvergenceInfo(approach, (P=P, B=B, ρs=norm(x.-z), ρx=norm(x), minimizer=z, fx=fx, minimum=fz, ∇fz=∇fz, f0=f0, ∇f0=∇f0, iter=0, time=time()-t0), options)
    end
    qnvars = QNVars(copy(objvars.∇fz), copy(objvars.∇fz), copy(objvars.∇fz))

    #==============================
             First iteration
    ==============================#
    objvars, P, qnvars = iterate(mstyle, qnvars, objvars, P, approach, problem, options)
    iter = 1
    # Check for gradient convergence
    is_converged = converged(approach, objvars, ∇f0, options)
    print_trace(approach, options, iter, t0, objvars)
    while iter < options.maxiter && !any(is_converged)
        iter += 1
        #==============================
                     iterate
        ==============================#
        objvars, P, qnvars = iterate(mstyle, qnvars, objvars, P, approach, problem, options, false)
        #==============================
                check convergence
        ==============================#
        is_converged = converged(approach, objvars, ∇f0, options)
        print_trace(approach, options, iter, t0, objvars)
    end
    x, fx, ∇fx, z, fz, ∇fz, B, Pg = objvars
    return ConvergenceInfo(approach, (P=P, B=B, ρs=norm(x.-z), ρx=norm(x), minimizer=z, fx=fx, minimum=fz, ∇fz=∇fz, f0=f0, ∇f0=∇f0, iter=iter, time=time()-t0), options)
end
function print_trace(approach::LineSearch, options, iter, t0, objvars)
    if false
       println(@sprintf("iter: %d   time: %.4f   obj: %.4e   ||∇f||: %.4e    α: %.4e", iter, time()-t0, objvars.fz, norm(objvars.∇fz, Inf), objvars.α))
    end
end
function iterate(mstyle::InPlace, cache, objvars, P, approach::LineSearch, problem::OptimizationProblem, options::OptimizationOptions, is_first=nothing)
    # split up the approach into the hessian approximation scheme and line search
    x, fx, ∇fx, z, fz, ∇fz, B, Pg = objvars

    Tf = typeof(fx)
    scheme, linesearch = modelscheme(approach), algorithm(approach)
    y, d, s = cache.y, cache.d, cache.s

    # Move nexts into currs
    fx = fz
    copyto!(x, z)
    copyto!(∇fx, ∇fz)

    # Update preconditioner
    P = update_preconditioner(scheme, x, P)
    # Update current gradient and calculate the search direction
    d = find_direction!(d, B, P, ∇fx, scheme) # solve Bd = -∇fx
    # real is needed to convert complex dots to actually being real
    φ = _lineobjective(mstyle, problem, ∇fz, z, x, d, fx, real(dot(∇fx, d)))

    # Perform line search along d
    # Also returns final step vector and update the state
    α, f_α, ls_success = find_steplength(mstyle, linesearch, φ, Tf(1))

    @. s = α * d
    z = retract(problem, z, x, s)

    # Update approximation
    fz, ∇fz, B, s, y = update_obj!(problem, s, y, ∇fx, z, ∇fz, B, scheme, is_first)
    return (x=x, fx=fx, ∇fx=∇fx, z=z, fz=fz, ∇fz=∇fz, B=B, Pg=Pg, α=α), P, QNVars(d, s, y)
end
function print_trace(approach::LineSearch, options, iter, t0, objvars, Δ)
    if false
        println(@sprintf("iter: %d   time: %f   f: %.4e   ||∇f||: %.4e    Δ: %.4e", iter, time()-t0, objvars.fz, norm(objvars.∇fz, Inf), Δ))
    end
end

function iterate(mstyle::OutOfPlace, cache, objvars, P, approach::LineSearch, problem::OptimizationProblem, options::OptimizationOptions, is_first=nothing)
    # split up the approach into the hessian approximation scheme and line search
    x, fx, ∇fx, z, fz, ∇fz, B, Pg = objvars
    Tf = typeof(fx)
    scheme, linesearch = modelscheme(approach), algorithm(approach)
    # Move nexts into currs
    fx = fz
    x = copy(z)
    ∇fx = copy(∇fz)

    # Update preconditioner
    P = update_preconditioner(scheme, x, P)
    # Update current gradient and calculate the search direction
    d = find_direction(B, P, ∇fx, scheme) # solve Bd = -∇fx
    # real is needed to convert complex dots to actually being real
    φ = _lineobjective(mstyle, problem, ∇fz, z, x, d, fx, real(dot(∇fx, d)))

    # Perform line search along d
    α, f_α, ls_success = find_steplength(mstyle, linesearch, φ, Tf(1))

    # # Calculate final step vector and update the state
    s = @. α * d
    z = retract(problem, z, x, s)

    # Update approximation
    fz, ∇fz, B, s, y = update_obj(problem, s, ∇fx, z, ∇fz, B, scheme, is_first)

    return (x=x, fx=fx, ∇fx=∇fx, z=z, fz=fz, ∇fz=∇fz, B=B, Pg=Pg), P, QNVars(d, s, y)
end


text:
using Test, Colors

@testset "Algorithms" begin

    # issue #349
    msc_h_diff = 0
    for hsv_h in 0:0.1:360
        hsv = HSV(hsv_h, 1.0, 1.0) # most saturated
        lch = convert(LCHuv, hsv)
        msc = MSC(lch.h)
        msc_h_diff = max(msc_h_diff, colordiff(msc, lch))
    end
    @test msc_h_diff < 1

    @test MSC(123.45, 100) ≈ 0
    @test MSC(123.45, 0) ≈ 0

    msc_h_l_sat = 1
    for h = 0:0.1:359.9, l = 1:1:99
        c = MSC(h, l)
        hsv = convert(HSV, LCHuv(l, c, h))
        # When the color is saturated (except black/white), `s` or `v` is 1.
        msc_h_l_sat = min(msc_h_l_sat, max(hsv.s, hsv.v))
    end
    @test msc_h_l_sat > 1 - 1e-4

    # the linear interpolation introduces some approximation errors.(issue #349)
    @test MSC(0, 90, linear=true) > MSC(0, 90)
    @test MSC(280, 50, linear=true) < MSC(280, 50)

    @testset "find_maximum_chroma hsv_h=$hsv_h" for hsv_h in 0:60:300
        hsv = HSV(hsv_h, 1.0, 1.0) # corner
        lchab = convert(LCHab, hsv)
        lchuv = convert(LCHuv, hsv)
        lab = convert(Lab, hsv)
        luv = convert(Luv, hsv)
        @test Colors.find_maximum_chroma(lchab) ≈ lchab.c atol=0.01
        @test Colors.find_maximum_chroma(lchuv) ≈ lchuv.c atol=0.01
        @test Colors.find_maximum_chroma(lab) ≈ lchab.c atol=0.01
        @test Colors.find_maximum_chroma(luv) ≈ lchuv.c atol=0.01
    end

    # yellow in LCHab
    @test Colors.find_maximum_chroma(LCHab(94.2, 0, 100)) ≈ 93.749 atol=0.01
    @test Colors.find_maximum_chroma(LCHab(97.6, 0, 105)) ≈ 68.828 atol=0.01
end


text:
# ==========================================================================================
# COMPARISON OF FEATURE EXTRACTION BETWEEN 3 DIFFERENT LOCAL DISCRIMINANT BASIS (LDB)
# ALGORITHMS
#
# This is the source code used to compare the differences in feature extraction capabilities
# between 3 LDB algorithms:
#   1) LDBK:
#           LDB that relies on the Time-Frequency energy distribution of classes. When
#           computing the discriminant measures of each node in the decomposed binary tree,
#           only the top K largest coefficients from each subspace are used.
#
#   2) LDBKASH:
#           LDB that relies on the Empirical Probability Densities of each class. When
#           computing the discriminant measures of each node in the decomposed binary tree,
#           only the top K largest coefficients from each subspace are used.
#
#   3) LDBKEMD:
#           LDB that depends on the use of signatures and Earth Mover Distance (EMD). When
#           computing the discriminant measures of each node in the decomposed binary tree,
#           only the top K largest coefficients from each subspace are used.
#
# In this file, the default values used to compare the algorithms are as follows:
#   1) LDBK:
#           Energy map          : Time-Frequency
#           Discriminant measure: Asymmetric Relative Entropy (Kullback-Leibler divergence)
#           K                   : 10
#           Discriminant power  : Output from discriminant measure
#
#   2) LDBKASH:
#           Energy map          : Empirical Probability Density
#           Discriminant measure: Asymmetric Relative Entropy (Kullback-Leibler divergence)
#           K                   : 10
#           Discriminant power  : Output from discriminant measure
#
#   3) LDBKEMD:
#           Energy map          : Signatures
#           Discriminant measure: Earth Mover's Distance
#           K                   : 10
#           Discriminant power  : Output from discriminant measure
#
#   Miscellaneous:
#           Dataset             : Triangular waveform dataset
#           Wavelet filter      : 6-tap Coiflet
#           Decomposition levels: 5
#           Train size          : 33/class
#           Test size           : 333/class
#           Features extracted  : 10
#
# The areas of the code that was set to the default values can be switched out. Look for the
# comments labeled with
#           # TODO: ___________
# at the start.
#
# ------------------------------------------------------------------------------------------

using Wavelets,
      WaveletsExt,
      Statistics,
      Plots,
      MLJ,
      DataFrames

include("utils.jl")

## ========== Setup ==========
train_x, train_y = generateclassdata(ClassData(:tri, 33 , 33 , 33 ), false)
test_x , test_y  = generateclassdata(ClassData(:tri, 333, 333, 333), true)
# TODO: Change wavelet filter
wt = wavelet(WT.coif2)
# TODO: Change decomposition levels
L = 5
# TODO: Change the value of K
k = 10
# TODO: Change the number of features extracted
n_features = 10
# TODO： Change the type of discriminant power
dp = BasisDiscriminantMeasure()
# TODO: Change the discriminant measures
dm_ldbk    = AsymmetricRelativeEntropy()
dm_ldbkash = SymmetricRelativeEntropy()
dm_ldbkemd = EarthMoverDistance()
# Create LDB objects
ldbk    = LocalDiscriminantBasis(wt = wt,
                                 max_dec_level = L,
                                 dm = dm_ldbk,
                                 en = TimeFrequency(),
                                 dp = dp,
                                 top_k = k,
                                 n_features = n_features)
ldbkash = LocalDiscriminantBasis(wt = wt,
                                 max_dec_level = L,
                                 dm = dm_ldbkash,
                                 en = ProbabilityDensity(),
                                 dp = dp,
                                 top_k = k,
                                 n_features = n_features)
ldbkemd = LocalDiscriminantBasis(wt = wt,
                                 max_dec_level = L,
                                 dm = dm_ldbkemd,
                                 en = Signatures(),
                                 dp = dp,
                                 top_k = k,
                                 n_features = n_features)

## ========== Display sample dataset ==========
# TODO: Change indexing if train size is changed above
class1_mean = mean(train_x[:, 1:33], dims=2)
class2_mean = mean(train_x[:,34:66], dims=2)
class3_mean = mean(train_x[:,67:99], dims=2)
p01 = plot(class1_mean, linestyle=:solid, color=:black, label="Class 1")
plot!(p01, class2_mean, linestyle=:dash, color=:black, label="Class 2")
plot!(p01, class3_mean, linestyle=:dot, color=:black, label="Class 3")
plot!(p01, title="Mean training waveform")

## ========== Fit LDB ==========
# Fit and transform training data
ldbk_train_features    = WaveletsExt.fit_transform(ldbk   , train_x, train_y)
ldbkash_train_features = WaveletsExt.fit_transform(ldbkash, train_x, train_y)
ldbkemd_train_features = WaveletsExt.fit_transform(ldbkemd, train_x, train_y)
# Transform test data
ldbk_test_features    = WaveletsExt.transform(ldbk, test_x)
ldbkash_test_features = WaveletsExt.transform(ldbkash, test_x)
ldbkemd_test_features = WaveletsExt.transform(ldbkemd, test_x)

## ========== Display top N LDB vectors ==========
# TODO: Change the value of N (Default: 3)
N = 3
ldbk_vec = compute_ldb_vectors(ldbk, N)
ldbkash_vec = compute_ldb_vectors(ldbkash, N)
ldbkemd_vec = compute_ldb_vectors(ldbkemd, N)
p02 = plot(ldbk_vec, label="")
p03 = plot(ldbkash_vec, label="")
p04 = plot(ldbkemd_vec, label="")

## ========== Display best basis nodes from each LDB algorithm ==========
p05 = plot_tfbdry(ldbk.tree)
p06 = plot_tfbdry(ldbkash.tree)
p07 = plot_tfbdry(ldbkemd.tree)

## ===== Scatter plot of coefficients in the top 2 most discriminating LDB coordinates =====
p08 = plot_coefficients(ldbk_train_features) |> p -> plot!(p, title="LDBK")
p09 = plot_coefficients(ldbkash_train_features) |> p -> plot!(p, title="LDBKASH")
p10 = plot_coefficients(ldbkemd_train_features) |> p -> plot!(p, title="LDBKEMD")

## ========== Fit features in Linear Discriminant Analysis (LDA) classifier ==========
# Data wrangling to fit MLJ.jl syntax
original_train_features = DataFrame(train_x', :auto)
ldbk_train_features     = DataFrame(ldbk_train_features', :auto)
ldbkash_train_features  = DataFrame(ldbkash_train_features', :auto)
ldbkemd_train_features  = DataFrame(ldbkemd_train_features', :auto)
train_y = coerce(train_y, Multiclass)
original_test_features = DataFrame(test_x', :auto)
ldbk_test_features     = DataFrame(ldbk_test_features', :auto)
ldbkash_test_features  = DataFrame(ldbkash_test_features', :auto)
ldbkemd_test_features  = DataFrame(ldbkemd_test_features', :auto)
test_y = coerce(test_y, Multiclass)
# Model fitting
LDA = @load LDA pkg=MultivariateStats
original_classifier = machine(LDA(), original_train_features, train_y)
ldbk_classifier     = machine(LDA(), ldbk_train_features, train_y)
ldbkash_classifier  = machine(LDA(), ldbkash_train_features, train_y)
ldbkemd_classifier  = machine(LDA(), ldbkemd_train_features, train_y)
MLJ.fit!(original_classifier)
MLJ.fit!(ldbk_classifier)
MLJ.fit!(ldbkash_classifier)
MLJ.fit!(ldbkemd_classifier)
# Model predictions
original_train_ŷ = predict_mode(original_classifier, original_train_features)
ldbk_train_ŷ     = predict_mode(ldbk_classifier    , ldbk_train_features)
ldbkash_train_ŷ  = predict_mode(ldbkash_classifier , ldbkash_train_features)
ldbkemd_train_ŷ  = predict_mode(ldbkemd_classifier , ldbkemd_train_features)
original_test_ŷ = predict_mode(original_classifier, original_test_features)
ldbk_test_ŷ     = predict_mode(ldbk_classifier    , ldbk_test_features)
ldbkash_test_ŷ  = predict_mode(ldbkash_classifier , ldbkash_test_features)
ldbkemd_test_ŷ  = predict_mode(ldbkemd_classifier , ldbkemd_test_features)
# Evaluate Model
original_train_accuracy = Accuracy()(original_train_ŷ, train_y)
ldbk_train_accuracy     = Accuracy()(ldbk_train_ŷ    , train_y)
ldbkash_train_accuracy  = Accuracy()(ldbkash_train_ŷ , train_y)
ldbkemd_train_accuracy  = Accuracy()(ldbkemd_train_ŷ , train_y)
original_test_accuracy = Accuracy()(original_test_ŷ, test_y)
ldbk_test_accuracy     = Accuracy()(ldbk_test_ŷ    , test_y)
ldbkash_test_accuracy  = Accuracy()(ldbkash_test_ŷ , test_y)
ldbkemd_test_accuracy  = Accuracy()(ldbkemd_test_ŷ , test_y)


text:
@recipe function f(b::FullBinner{T}) where T<:Number
    framestyle --> :grid
    size --> (750, 500)
    plottype = get(plotattributes, :seriestype, :plot)

    if plottype == :histogram
        @series begin
            color --> :lightblue
            alpha --> 0.6
            label --> ""
            legend --> true
            b.x
        end
    elseif plottype == :binning
        bss, errors, error_means = BinningAnalysis.all_binning_errors(b)
        @series begin
            label --> "standard error"
            seriestype := :line
            markershape --> :circle
            markercolor --> :lightgreen
            markersize --> 2
            color --> :lightgreen
            legend --> true
            bss, errors
        end

        @series begin
            seriestype := :line
            xlabel --> "bin sizes"
            label --> "cumulative mean"
            ylabel --> "standard error of the mean"
            markershape --> :circle
            markercolor --> :green
            markersize --> 2
            color --> :green
            legend --> true
            bss, error_means
        end
    elseif plottype == :corrplot
        y = autocor(b.x)
        @series begin
            seriestype := :line
            markershape --> :circle
            markercolor --> :blue
            markersize --> 2
            color --> :blue
            y
        end

        @series begin
            seriestype := :line
            xlabel --> "lag"
            ylabel --> "autocorrelation"
            legend --> false
            color --> :grey
            linestyle --> :dash
            ylim --> (min(minimum(y),0), 1.1)
            fill(0, length(y))
        end
    elseif plottype == :tauplot
        bss, errors, error_means = BinningAnalysis.all_binning_errors(b)
        taus = tau.(Ref(b), errors)
        @series begin
            seriestype := :line
            xlabel --> "bin sizes"
            label --> "autocorrelation time"
            markershape --> :circle
            markercolor --> :lightblue
            markersize --> 2
            color --> :blue
            bss, taus
        end

        @series begin
            tau_means = [mean(@view taus[1:i]) for i in 1:length(taus)]
            seriestype := :line
            xlabel --> "bin sizes"
            label --> "cumulative mean"
            ylabel --> "autocorrelation time"
            markershape --> :circle
            markercolor --> :blue
            markersize --> 2
            color --> :green
            legend --> true
            bss, tau_means
        end
    else
        @series begin
            xlabel --> "time t"
            label --> "time series"
            markershape --> :circle
            markercolor --> :grey
            markersize --> 2
            color --> :grey
            legend --> true
            b.x
        end
    end

    μ = mean(b)
    if plottype == :histogram
        @series begin
            seriestype := :vline
            color --> :black
            label --> "mean"
            fill(μ, length(b))
        end
    elseif plottype in (:binning, :corrplot, :tauplot)
        nothing
    else
        @series begin
            seriestype := :hline
            color --> :purple
            label --> "mean"
            fill(μ, length(b))
        end
    end

    Δ = std_error(b)
    if !(plottype in (:binning, :corrplot, :tauplot))
        for i in (1, -1)
            @series begin
                seriestype := (plottype == :histogram ? :vline : :hline)
                color --> :violet
                label --> (i == 1 ? "std error" : "")
                fill(μ + i*Δ, length(b))
            end
        end
    end
end

@shorthands binning
@shorthands corrplot
@shorthands tauplot

# ---------------------- LogBinner -----------------------
@recipe function f(b::LogBinner{T}) where T<:Number
    framestyle --> :grid
    size --> (750, 500)
    plottype = get(plotattributes, :seriestype, :plot)

    if plottype == :tauplot
        taus = all_taus(b)
        bss = 2 .^ (2:length(taus)+1)
        idxs = findall(bs->bs<=length(b)/30, bss)
        taus = taus[idxs]
        bss = bss[idxs]
        @series begin
            seriestype := :line
            xlabel --> "bin sizes"
            label --> "autocorrelation time"
            markershape --> :circle
            markercolor --> :lightblue
            markersize --> 2
            color --> :blue
            bss, taus
        end

        @series begin
            tau_means = [mean(@view taus[1:i]) for i in 1:length(taus)]
            seriestype := :line
            xlabel --> "bin sizes"
            label --> "cumulative mean"
            ylabel --> "autocorrelation time"
            markershape --> :circle
            markercolor --> :blue
            markersize --> 2
            color --> :green
            legend --> true
            bss, tau_means
        end
    elseif plottype == :corrplot
        y = autocor(b.x_sum)
        @series begin
            seriestype := :line
            markershape --> :circle
            markercolor --> :blue
            markersize --> 2
            color --> :blue
            y
        end

        @series begin
            xlabel --> "lag"
            ylabel --> "autocorrelation (of log bin means!)"
            seriestype := :line
            color --> :grey
            linestyle --> :dash
            ylim --> (min(minimum(y),0), 1.1)
            legend --> false
            fill(0, length(y))
        end
    elseif plottype == :binning
        errors = all_std_errors(b)
        bss = 2 .^ (2:length(errors)+1)
        idxs = findall(bs->bs<=length(b)/30, bss)
        bss = bss[idxs]
        errors = errors[idxs]
        error_means = [mean(errors[1:i]) for i in 1:length(errors)] # TODO: optimize
        @series begin
            label --> "standard error"
            seriestype := :line
            markershape --> :circle
            markercolor --> :lightgreen
            markersize --> 2
            color --> :lightgreen
            legend --> true
            bss, errors
        end

        @series begin
            seriestype := :line
            xlabel --> "bin sizes"
            label --> "cumulative mean"
            ylabel --> "standard error of the mean"
            markershape --> :circle
            markercolor --> :green
            markersize --> 2
            color --> :green
            legend --> true
            bss, error_means
        end
    end
end


text:
#######################################################################
# An example of creating Excel Column charts with XlsxWriter
# Original Python Copyright 2013-2016, John McNamara, jmcnamara@cpan.org
# https://github.com/jmcnamara/XlsxWriter

using XlsxWriter

function test()

	wb = Workbook("chart_column.xlsx")
	ws = add_worksheet!(wb)
	bold = add_format!(wb, Dict("bold"=>1))

	# Add the worksheet data that the charts will refer
	headings = ["Number" "Batch 1" "Batch 2"]
	data = [
	 [2 3 4 5 6 7];
	 [10 40 50 20 10 50];
	 [30 60 70 50 40 30];
	]

	 write_row!(ws, "A1", headings, fmt=bold)
	 write_column!(ws, "A2", data[1,:])
	 write_column!(ws, "B2", data[2,:])
	 write_column!(ws, "C2", data[3,:])

	#######################################################################
	#
	# Create a new column
	#
	chart1 =  add_chart!(wb, Dict("type"=>"column"))

	# Configure the first
	 add_series!(chart1, Dict(
		"name"=>"=Sheet1!\$B\$1",
		"categories"=>"=Sheet1!\$A\$2:\$A\$7",
		"values"=>"=Sheet1!\$B\$2:\$B\$7",
	))

	# Configure a second   Note use of alternative syntax to define
	 add_series!(chart1, Dict(
		"name"=>["Sheet1", 0, 2],
		"categories"=>["Sheet1", 1, 0, 6, 0],
		"values"=>["Sheet1", 1, 2, 6, 2],
	))

	# Add a chart title and some axis
	 set_title!(chart1, Dict("name"=>"Results of sample analysis"))
	 set_x_axis!(chart1, Dict("name"=>"Test number"))
	 set_y_axis!(chart1, Dict("name"=>"Sample length (mm)"))

	# Set an Excel chart
	 set_style!(chart1, 11)

	# Insert the chart into the worksheet (with an offset).
	 insert_chart!(ws, "D2", chart1, Dict("x_offset"=>25, "y_offset"=>10))

	#######################################################################
	#
	# Create a stacked chart sub-
	#
	chart2 =  add_chart!(wb, Dict("type"=>"column", "subtype"=>"stacked"))

	# Configure the first
	 add_series!(chart2, Dict(
		"name"=>"=Sheet1!\$B\$1",
		"categories"=>"=Sheet1!\$A\$2:\$A\$7",
		"values"=>"=Sheet1!\$B\$2:\$B\$7",
	))

	# Configure second
	 add_series!(chart2, Dict(
		"name"=>"=Sheet1!\$C\$1",
		"categories"=>"=Sheet1!\$A\$2:\$A\$7",
		"values"=>"=Sheet1!\$C\$2:\$C\$7",
	))

	# Add a chart title and some axis
	 set_title!(chart2, Dict("name"=>"Stacked Chart"))
	 set_x_axis!(chart2, Dict("name"=>"Test number"))
	 set_y_axis!(chart2, Dict("name"=>"Sample length (mm)"))

	# Set an Excel chart
	 set_style!(chart2, 12)

	# Insert the chart into the worksheet (with an offset).
	 insert_chart!(ws, "D18", chart2, Dict("x_offset"=>25, "y_offset"=>10))

	#######################################################################
	#
	# Create a percentage stacked chart sub-
	#
	chart3 =  add_chart!(wb, Dict("type"=>"column", "subtype"=>"percent_stacked"))

	# Configure the first
	 add_series!(chart3, Dict(
		"name"=>"=Sheet1!\$B\$1",
		"categories"=>"=Sheet1!\$A\$2:\$A\$7",
		"values"=>"=Sheet1!\$B\$2:\$B\$7",
	))

	# Configure second
	 add_series!(chart3, Dict(
		"name"=>"=Sheet1!\$C\$1",
		"categories"=>"=Sheet1!\$A\$2:\$A\$7",
		"values"=>"=Sheet1!\$C\$2:\$C\$7",
	))

	# Add a chart title and some axis
	 set_title!(chart3, Dict("name"=>"Percent Stacked Chart"))
	 set_x_axis!(chart3, Dict("name"=>"Test number"))
	 set_y_axis!(chart3, Dict("name"=>"Sample length (mm)"))

	# Set an Excel chart
	 set_style!(chart3, 13)

	# Insert the chart into the worksheet (with an offset).
	 insert_chart!(ws, "D34", chart3, Dict("x_offset"=>25, "y_offset"=>10))

	 close(wb)
	 isfile("chart_column.xlsx")
end

test()


text:
using ParameterHandling: vec_to_tril, tril_to_vec

@testset "parameters_matrix.jl" begin
    @testset "orthogonal" begin
        is_almost_orthogonal(X::AbstractMatrix, tol) = norm(X'X - I) < tol

        @testset "nearest_orthogonal_matrix($T)" for T in [Float64, ComplexF64]
            X_orth = ParameterHandling.nearest_orthogonal_matrix(randn(T, 5, 4))
            @test is_almost_orthogonal(X_orth, 1e-9)
            X_orth_2 = ParameterHandling.nearest_orthogonal_matrix(X_orth)
            @test X_orth ≈ X_orth_2 # nearest_orthogonal_matrix is a projection.
        end

        X = orthogonal(randn(5, 4))
        @test X == X
        test_parameter_interface(X)
        @test is_almost_orthogonal(value(X), 1e-9)

        # We do not implement any custom rrules, so we only check that `Zygote` is able to
        # differentiate, and assume that the result is correct if it doesn't error.
        @testset "Zygote" begin
            _, pb = Zygote.pullback(X -> value(orthogonal(X)), randn(3, 2))
            @test only(pb(randn(3, 2))) isa Matrix{<:Real}
        end
    end

    @testset "positive_definite" begin
        @testset "vec_tril_conversion" begin
            X = tril!(rand(3, 3))
            @test vec_to_tril(tril_to_vec(X)) == X
            @test_throws ErrorException tril_to_vec(rand(4, 5))
        end
        X_mat = ParameterHandling.A_At(rand(3, 3)) # Create a positive definite object
        X = positive_definite(X_mat)
        @test X == X
        @test value(X) ≈ X_mat
        @test isposdef(value(X))
        @test vec_to_tril(X.L) ≈ cholesky(X_mat).L
        @test_throws ArgumentError positive_definite(rand(3, 3))
        test_parameter_interface(X)

        x, re = flatten(X)
        Δl = first(
            Zygote.gradient(x) do x
                X = re(x)
                return logdet(value(X))
            end,
        )
        ΔL = first(
            Zygote.gradient(vec_to_tril(X.L)) do L
                return logdet(L * L')
            end,
        )
        @test vec_to_tril(Δl) == tril(ΔL)
        ChainRulesTestUtils.test_rrule(vec_to_tril, x)
    end
end


text:
#
# This module represents an explicit, direct sparse 2D finite-difference Poisson solver for heterogeneous media,
# i.e. media having spatially varying (space-dependent) medium parameters.
# Current implementation assumes Neumann boundary conditions at all the boundaries.
#
# Consider the following Poisson experiment:
# ```math
# ∇⋅(σ(x,z)∇) ψ(t) = ∇⋅(Q(x,z)∇) p(t),
# ```
# ```math
# Q = k * Q_v / η.
# ```


text:
### Day 3 - Problem 2
## Squares With Three Sides
## Author: Thanasis Georgiou

# Read triangle data
input = readdlm("triangles.csv")

# Start counting
count = 0

# For each triangle
for i = 0:1:size(input, 1)-1
    # Calculate vertical index of triangle
    # i + 1: 1, 2, 3, 4, 5, 6, 7, 8, 9...
    # v    : 1, 1, 1, 4, 4, 4, 7, 7, 7...
    v = i - (i % 3) + 1

    # Grab triangle (verI + 2 cells below)
    triangle = sort(
        view(input, v:v + 2, (i % 3) + 1)
    )

    # Check if it's a valid triangle
    count += triangle[1] + triangle[2] > triangle[3]
end

# Print result
println("Possible Triangles: $(count)")


text:
# Julia/LLVM IR generation and transformation passes

function module_setup(mod::LLVM.Module)
    triple!(mod, "amdgcn-amd-amdhsa")

    # add debug info metadata
    push!(metadata(mod), "llvm.module.flags",
         MDNode([ConstantInt(Int32(1), JuliaContext()),    # llvm::Module::Error
                 MDString("Debug Info Version"),
                 ConstantInt(DEBUG_METADATA_VERSION(), JuliaContext())]))
end

# generate a pseudo-backtrace from a stack of methods being emitted
function backtrace(ctx::CompilerContext, method_stack::Vector{Core.MethodInstance})
    bt = StackTraces.StackFrame[]
    for method_instance in method_stack
        method = method_instance.def
        frame = StackTraces.StackFrame(method.name, method.file, method.line)
        pushfirst!(bt, frame)
    end
    bt
end

# NOTE: we use an exception to be able to display a stack trace using the logging framework
struct MethodSubstitutionWarning <: Exception
    original::Method
    substitute::Method
end
Base.showerror(io::IO, err::MethodSubstitutionWarning) =
    print(io, "You called $(err.original), maybe you intended to call $(err.substitute) instead?")

function compile_linfo(ctx::CompilerContext, linfo::Core.MethodInstance, world)
    # set-up the compiler interface
    function hook_module_setup(ref::Ptr{Cvoid})
        ref = convert(LLVM.API.LLVMModuleRef, ref)
        module_setup(LLVM.Module(ref))
    end
    dependencies = Vector{LLVM.Module}()
    function hook_module_activation(ref::Ptr{Cvoid})
        ref = convert(LLVM.API.LLVMModuleRef, ref)
        push!(dependencies, LLVM.Module(ref))
    end
    method_stack = Vector{Core.MethodInstance}()
    function hook_emit_function(method_instance, code, world)
        push!(method_stack, method_instance)

        # check for recursion
        if method_instance in method_stack[1:end-1]
            throw(KernelError(ctx, "recursion is currently not supported";
                              bt=backtrace(ctx, method_stack)))
        end

        # check for Base methods that exist in AMDGPUnative too
        # FIXME: this might be too coarse
        method = method_instance.def
        if Base.moduleroot(method.module) == Base &&
           isdefined(AMDGPUnative, method_instance.def.name)
            substitute_function = getfield(AMDGPUnative, method.name)
            tt = Tuple{method_instance.specTypes.parameters[2:end]...}
            if hasmethod(substitute_function, tt)
                method′ = which(substitute_function, tt)
                if Base.moduleroot(method′.module) == AMDGPUnative
                    @warn "calls to Base intrinsics might be GPU incompatible" exception=(MethodSubstitutionWarning(method, method′), backtrace(ctx, method_stack))
                end
            end
        end
    end
    function hook_emitted_function(method, code, world)
        @compiler_assert last(method_stack) == method ctx
        pop!(method_stack)
    end
    params = Base.CodegenParams(cached             = false,
                                track_allocations  = false,
                                code_coverage      = false,
                                static_alloc       = false,
                                prefer_specsig     = true,
                                module_setup       = hook_module_setup,
                                module_activation  = hook_module_activation,
                                emit_function      = hook_emit_function,
                                emitted_function   = hook_emitted_function)

    # get the code
    ref = ccall(:jl_get_llvmf_defn, LLVM.API.LLVMValueRef,
                (Any, UInt, Bool, Bool, Base.CodegenParams),
                linfo, world, #=wrapper=#false, #=optimize=#false, params)
    if ref == C_NULL
        throw(InternalCompilerError(ctx, "the Julia compiler could not generate LLVM IR"))
    end
    llvmf = LLVM.Function(ref)

    modules = [LLVM.parent(llvmf), dependencies...]
    return llvmf, modules
end

function irgen(ctx::CompilerContext)
    # get the method instance
    isa(ctx.f, Core.Builtin) && throw(KernelError(ctx, "function is not a generic function"))
    world = typemax(UInt)
    meth = which(ctx.f, ctx.tt)
    sig = Base.signature_type(ctx.f, ctx.tt)::Type
    (ti, env) = ccall(:jl_type_intersection_with_env, Any,
                      (Any, Any), sig, meth.sig)::Core.SimpleVector
    meth = Base.func_for_method_checked(meth, ti)
    linfo = ccall(:jl_specializations_get_linfo, Ref{Core.MethodInstance},
                  (Any, Any, Any, UInt), meth, ti, env, world)

    entry, modules = compile_linfo(ctx, linfo, world)

    # link in dependent modules
    mod = popfirst!(modules)
    for dep in modules
        link!(mod, dep)
    end

    # clean up incompatibilities
    for llvmf in functions(mod)
        llvmfn = LLVM.name(llvmf)

        # only occurs in debug builds
        delete!(function_attributes(llvmf), EnumAttribute("sspstrong", 0, JuliaContext()))

        # dependent modules might have brought in other jfptr wrappers, delete them
        if startswith(LLVM.name(llvmf), "jfptr_") && isempty(uses(llvmf))
            unsafe_delete!(mod, llvmf)
            continue
        end

        # llvmcall functions aren't to be called, so mark them internal (cleans up the IR)
        if startswith(llvmfn, "jl_llvmcall")
            linkage!(llvmf, LLVM.API.LLVMInternalLinkage)
            continue
        end

        # rename functions
        if !isdeclaration(llvmf)
            # Julia disambiguates local functions by prefixing with `#\d#`.
            # since we don't use a global function namespace, get rid of those tags.
            if occursin(r"^julia_#\d+#", llvmfn)
                llvmfn′ = replace(llvmfn, r"#\d+#"=>"")
                if !haskey(functions(mod), llvmfn′)
                    LLVM.name!(llvmf, llvmfn′)
                    llvmfn = llvmfn′
                end
            end

            # anonymous functions are just named `#\d`, make that somewhat more readable
            m = match(r"_#(\d+)_", llvmfn)
            if m !== nothing
                llvmfn′ = replace(llvmfn, m.match=>"_anonymous$(m.captures[1])_")
                LLVM.name!(llvmf, llvmfn′)
                llvmfn = llvmfn′
            end
        end
    end

    # rename the entry point
    llvmfn = replace(LLVM.name(entry), r"_\d+$"=>"")
    ## append a global unique counter
    global globalUnique
    globalUnique += 1
    llvmfn *= "_$globalUnique"
    LLVM.name!(entry, llvmfn)

    # minimal required optimization
    ModulePassManager() do pm
        global global_ctx
        global_ctx = ctx

        add!(pm, ModulePass("LowerThrow", lower_throw!))
        always_inliner!(pm)
        verifier!(pm)
        run!(pm, mod)
    end

    return mod, entry
end

# this pass lowers `jl_throw` and friends to GPU-compatible exceptions.
# this isn't strictly necessary, but has a couple of advantages:
# - we can kill off unused exception arguments that otherwise would allocate or invoke
# - we can fake debug information (lacking a stack unwinder)
#
# once we have thorough inference (ie. discarding `@nospecialize` and thus
# supporting exception arguments) and proper debug info to unwind the stack,
# this pass can go.
function lower_throw!(mod::LLVM.Module)
    ctx = global_ctx::CompilerContext
    changed = false

    throw_functions = Dict{String,String}(
        "jl_throw"                      => "exception",
        "jl_error"                      => "error",
        "jl_too_few_args"               => "too few arguments exception",
        "jl_too_many_args"              => "too many arguments exception",
        "jl_type_error_rt"              => "type error",
        "jl_undefined_var_error"        => "undefined variable error",
        "jl_bounds_error"               => "bounds error",
        "jl_bounds_error_v"             => "bounds error",
        "jl_bounds_error_int"           => "bounds error",
        "jl_bounds_error_tuple_int"     => "bounds error",
        "jl_bounds_error_unboxed_int"   => "bounds error",
        "jl_bounds_error_ints"          => "bounds error",
        "jl_eof_error"                  => "EOF error"
    )

    for (fn, name) in throw_functions
        if haskey(functions(mod), fn)
            f = functions(mod)[fn]

            for use in uses(f)
                call = user(use)::LLVM.CallInst

                # replace the throw with a GCN-compatible exception
                let builder = Builder(JuliaContext())
                    position!(builder, call)
                    emit_exception!(builder, name, call)
                    dispose(builder)
                end

                # remove the call
                call_args = collect(operands(call))[1:end-1] # last arg is function itself
                unsafe_delete!(LLVM.parent(call), call)

                # HACK: kill the exceptions' unused arguments
                for arg in call_args
                    # peek through casts
                    if isa(arg, LLVM.AddrSpaceCastInst)
                        cast = arg
                        arg = first(operands(cast))
                        isempty(uses(cast)) && unsafe_delete!(LLVM.parent(cast), cast)
                    end

                    if isa(arg, LLVM.Instruction) && isempty(uses(arg))
                        unsafe_delete!(LLVM.parent(arg), arg)
                    end
                end

                changed = true
            end

            @compiler_assert isempty(uses(f)) ctx
         end
     end

    return changed
end

# report an exception in a GPU-compatible manner
#
# the exact behavior depends on the debug level. in all cases, a `trap` will be emitted, On
# debug level 1, the exception name will be printed, and on debug level 2 the individual
# stack frames (as recovered from the LLVM debug information) will be printed as well.
function emit_exception!(builder, name, inst)
    bb = position(builder)
    fun = LLVM.parent(bb)
    mod = LLVM.parent(fun)

    # report the exception
    if Base.JLOptions().debug_level >= 1
        name = globalstring_ptr!(builder, name, "exception")
        if Base.JLOptions().debug_level == 1
            call!(builder, Runtime.get(:report_exception), [name])
        else
            call!(builder, Runtime.get(:report_exception_name), [name])
        end
    end

    # report each frame
    if Base.JLOptions().debug_level >= 2
        rt = Runtime.get(:report_exception_frame)
        bt = backtrace(inst)
        for (i,frame) in enumerate(bt)
            idx = ConstantInt(rt.llvm_types[1], i)
            func = globalstring_ptr!(builder, String(frame.func), "di_func")
            file = globalstring_ptr!(builder, String(frame.file), "di_file")
            line = ConstantInt(rt.llvm_types[4], frame.line)
            call!(builder, rt, [idx, func, file, line])
        end
    end

    trap = if haskey(functions(mod), "llvm.trap")
        functions(mod)["llvm.trap"]
    else
        LLVM.Function(mod, "llvm.trap", LLVM.FunctionType(LLVM.VoidType(JuliaContext())))
    end
    call!(builder, trap)
end


text:
module Abstractsmearing_module
    using LinearAlgebra
    import ..Wilsonloops_module:Wilson_loop_set,make_staples,Wilson_loop_set,
            make_cloverloops,Tensor_derivative_set, make_loops
    import ..AbstractGaugefields_module:AbstractGaugefields,Abstractfields,initialize_TA_Gaugefields,add_force!,
                        exp_aF_U!,clear_U!,add_U!,evaluate_wilson_loops!,exptU!,
                        Traceless_antihermitian_add!,set_wing_U!,Traceless_antihermitian,evaluate_gaugelinks!,
                        construct_Λmatrix_forSTOUT!,Traceless_antihermitian!,shift_U
    import Wilsonloop:Wilsonline,DwDU,make_loopforactions,make_Cμ,derive_U,derive_Udag,get_leftlinks,get_rightlinks
    import ..Verboseprint_mpi:Verbose_print,println_verbose_level1,println_verbose_level2,println_verbose_level3
    #import ..GaugeAction_module:GaugeAction

    abstract type Abstractsmearing end

    struct Nosmearing <: Abstractsmearing
    end

    abstract type CovLayer{Dim} end

    struct CovNeuralnet{Dim} <: Abstractsmearing
        #numlayers::Int64
        layers::Array{CovLayer{Dim},1}
    end

    function get_numlayers(c::CovNeuralnet)
        return length(c.layers)
    end

    function CovNeuralnet(;Dim=4)
        layers = CovLayer{Dim}[]
        return CovNeuralnet{Dim}(layers)
    end

    function Base.push!(c::CovNeuralnet,layer::CovLayer)
        #c.numlayers += 1
        push!(c.layers,layer)
    end

    function Base.getindex(c::CovNeuralnet,i)
        return c.layers[i]
    end

    function Base.show(c::CovNeuralnet)
        numlayers = get_numlayers(c)
        println("num. of layers: ",numlayers)
        for i=1:numlayers
            if i==1
                string = "st"
            elseif i==2
                string = "nd"
            elseif i==3
                string = "rd"
            else
                string = "th"
            end
            println("- $i-$string layer: $(get_name(c[i]))")
            show(c[i])
        end
    end

    function get_name(s::CovLayer)
        error("layer $s has no name")
    end

    include("./stout.jl")
    #include("./gradientflow.jl")

    function set_parameters(s::CovNeuralnet,i,v)
        set_parameters(s[i],v)
    end

    function construct_smearing(smearingparameters,loops_list,L,coefficients,numlayers)
        if smearingparameters == "nothing"
            smearing = Nosmearing()
        elseif smearingparameters == "stout"
            @assert loops_list != nothing "loops should be put if you want to use smearing schemes"
            loops = make_loops(loops_list,L)

            @assert coefficients != nothing "coefficients should be put if you want to use smearing schemes"
            println("stout smearing will be used")
            if numlayers == 1
                input_coefficients = [coefficients]
            else
                input_coefficients = coefficients
            end

            smearing = Stoutsmearing(loops,input_coefficients)
        elseif smearingparameters == "covnet_stout"
            if numlayers == 1
                input_coefficients = [coefficients]
            else
                input_coefficients = coefficients
            end
            println("covnet verion of the stout smearing will be used")
            smearing = CovNeuralnet_STOUT(loops_list,input_coefficients,L;Dim=length(L))
        else
            error("smearing = $smearing is not supported")
        end
        return smearing
    end

    function calc_smearedU(Uin::Array{T,1},smearing;calcdSdU = false,temps = nothing) where T<: AbstractGaugefields
        if smearing != nothing && typeof(smearing) != Nosmearing
            #println(smearing)
            #println(typeof(smearing))
            Uout_multi = apply_smearing_U(Uin,smearing)
            U = Uout_multi[end]

            #=
            if typeof(smearing) <: SmearingParam_single
                Uout_multi = nothing
                U = apply_smearing_U(Uin,smearing)
            elseif typeof(smearing) <: SmearingParam_multi
                Uout_multi = apply_smearing_U(Uin,smearing)
                U = Uout_multi[end]
            else
                error("something is wrong in calc_smearingU")
            end
            =#
            set_wing_U!(U)  #we want to remove this.
            if calcdSdU
                dSdU = [temps[end-3],temps[end-2],temps[end-1],temps[end]]
            else
                dSdU = nothing
            end
        else
            dSdU = nothing
            Uout_multi = nothing
            U = Uin
        end
        return U,Uout_multi,dSdU
    end

    function apply_smearing_U(Uin::Array{T,1},smearing) where T<: Abstractfields
        error("apply_smearing_U is not implemented in type $(typeof(Uin)) and smearing type $(typeof(smearing))")
    end

    function apply_smearing_U(Uin::Array{T,1},smearing::CovNeuralnet{Dim}) where {T<: Abstractfields,Dim}
        temp1  = similar(Uin[1])
        temp2  = similar(Uin[1])
        temp3  = similar(Uin[1])
        temp4  = similar(Uin[1])
        F0 = initialize_TA_Gaugefields(Uin[1])
        numlayers = get_numlayers(smearing)
        Uout_multi = Array{typeof(Uin),1}(undef,numlayers)
        for i=1:numlayers
            Uout_multi[i] = similar(Uin)
        end
        apply_neuralnet!(Uout_multi,smearing,Uin,[temp1,temp2,temp3,temp4],[F0])
        return Uout_multi

        error("apply_smearing_U is not implemented in type $(typeof(Uin)) and smearing type $(typeof(smearing))")
    end

    function back_prop(δL,net::CovNeuralnet{Dim},Uout_multi,Uin) where Dim
        temps = similar(Uout_multi[1])
        temps_F1 = initialize_TA_Gaugefields(temps[1])
        tempf = [temps_F1]

        layer = net.layers[get_numlayers(net)]
        δ_prev = similar(δL)
        δ_current = deepcopy(δL)
        set_wing_U!(δ_current)

        for i=get_numlayers(net):-1:2
            layer = net.layers[i]
            layer_pullback!(δ_prev,δ_current,layer,Uout_multi[i-1],temps,tempf)
            δ_current,δ_prev = δ_prev,δ_current
            #set_wing_U!(δ_current)
        end
        layer = net.layers[1]
        layer_pullback!(δ_prev,δ_current,layer,Uin,temps,tempf)
        δ_current,δ_prev = δ_prev,δ_current

        return δ_current
    end

    function get_parameter_derivatives(δL,net::CovNeuralnet{Dim},Uout_multi,Uin) where Dim
        temps = similar(Uout_multi[1])
        δs = get_δ_from_back_prop(δL,net,Uout_multi,Uin)
        numlayer = get_numlayers(net)
        dSdp = Array{Vector{Float64},1}(undef,numlayer)

        for i=1:numlayer
            layer = net.layers[i]
            U_current = Uout_multi[i]
            dSdp[i] = parameter_derivatives(δs[i],layer,U_current,temps)
        end
        return dSdp
    end

    function get_δ_from_back_prop(δL,net::CovNeuralnet{Dim},Uout_multi,Uin) where Dim
        temps = similar(Uout_multi[1])
        temps_F1 = initialize_TA_Gaugefields(temps[1])
        tempf = [temps_F1]

        layer = net.layers[get_numlayers(net)]
        δ_prev = similar(δL)
        δ_current = deepcopy(δL)
        set_wing_U!(δ_current)

        δs = Array{typeof(δL),1}(undef,get_numlayers(net))

        for i=get_numlayers(net):-1:2
            layer_pullback!(δ_prev,δ_current,layer,Uout_multi[i-1],temps,tempf)
            δs[i] = deepcopy(δ_prev)
            δ_current,δ_prev = δ_prev,δ_current
            #set_wing_U!(δ_current)
        end
        layer = net.layers[1]
        layer_pullback!(δ_prev,δ_current,layer,Uin,temps,tempf)
        δs[1] = deepcopy(δ_prev)
        δ_current,δ_prev = δ_prev,δ_current

        return δs
    end

    function apply_neuralnet!(Uout_multi,net::CovNeuralnet{Dim},Uin,temps,temps_F) where {Dim}
        layer = net.layers[1]
        println_verbose_level3(Uin[1],"apply_neuralnet! on 1st layer")
        apply_layer!(Uout_multi[1],layer,Uin,temps,temps_F)
        #set_wing_U!(Uout_multi[1])
        for i=2:get_numlayers(net)
            println_verbose_level3(Uin[1],"apply_neuralnet! on $i -th layer")
            layer = net.layers[i]
            apply_layer!(Uout_multi[i],layer,Uout_multi[i-1],temps,temps_F)
            #set_wing_U!(Uout_multi[i])
        end
    end

    function apply_layer!(Uout,layer::T,Uin,temps,temps_F) where T <: CovLayer
        error("apply_layer!(Uout,layer,Uin,temps,temps_F) is not implemented with type $(typeof(layer)) of layer.")
    end

    """
    layer_pullback!(δ_prev,δ_next,layer::T,Uprev,temps,tempf)
    This is a function for a back propagation
        δ_next,Uprev -> δ_prev
    """
    function layer_pullback!(δ_prev,δ_next,layer::T,Uprev,temps,tempf) where T <: CovLayer
        error("layer_pullback!(δ_prev,δ_next,layer::T,Uprev,temps,tempf) is not implemented with type $(typeof(layer)) of layer.")
    end

    function parameter_derivatives(δ_current,layer::T,U_current,temps) where T <: CovLayer
        error("parameter_derivatives(δ_current,layer::T,U_current,temps) is not implemented with type $(typeof(layer)) of layer.")
    end

    function apply_smearing_U(U::Array{T,1},smearing::Stoutsmearing) where T<: Abstractfields
        numlayer = length(smearing.ρs)
        Uout_multi = Array{typeof(U),1}(undef,numlayer)
        for i=1:numlayer
            Uout_multi[i] = similar(U)
        end
        #println("smearing.ρs ", smearing.ρs)
        #println("type U ",typeof(Uout_multi))
        #Uout = similar(U)
        calc_stout_multi!(Uout_multi,U,smearing.ρs,smearing.staples_for_stout)
        #Uout_multi = calc_stout_multi!(U,smearing.ρs,smearing.staples_for_stout)
        return Uout_multi
    end

    #=
    function calc_stout_multi(Uin::Array{<: AbstractGaugefields{NC,Dim},1},ρs::Array{Array{T,1},1},staples)  where {NC,Dim,T <: Number}
        numlayer = length(ρs)
        #println("numlayer = ",numlayer,"\t",ρs)
        Utype = eltype(Uin)
        Uout_multi = Array{Array{Utype,1}}(undef,numlayer)
        for i=1:numlayer
            Uout_multi[i] = similar(Uin)
        end
        calc_stout_multi!(Uout_multi,Uin,ρs,staples)

        return Uout_multi
    end
    =#

    function calc_stout_multi!(Uout_multi::Vector{<: Vector{<: AbstractGaugefields{NC,Dim}}},Uin::Array{<: AbstractGaugefields{NC,Dim},1},ρs::Array{Array{T,1},1},staples)  where {NC,Dim,T <: Number}
        numlayer = length(ρs)
        Utmp = similar(Uin)
        #Uout_multi = Array{Array{GaugeFields{SU{NC}},1}}(undef,numlayer)
        U = deepcopy(Uin)
        for i = 1:numlayer
            if i != numlayer
                apply_stout_smearing!(Utmp,U,ρs[i],staples)
                Uout_multi[i] = deepcopy(Utmp)
                Utmp,U = U,Utmp
            else
                apply_stout_smearing!(Uout_multi[i],U,ρs[i],staples)
            end
        end
    end

    function calc_stout_multi!(Uout::Array{<: AbstractGaugefields{NC,Dim},1},Uin::Array{<: AbstractGaugefields{NC,Dim},1},ρs::Array{Array{T,1},1},staples)  where {NC,Dim,T <: Number}
        numlayer = length(ρs)
        Utmp = similar(Uin)
        #Uout_multi = Array{Array{GaugeFields{SU{NC}},1}}(undef,numlayer)
        U = deepcopy(Uin)
        for i = 1:numlayer
            if i != numlayer
                apply_stout_smearing!(Utmp,U,ρs[i],staples)
                Utmp,U = U,Utmp
            else
                apply_stout_smearing!(Uout,U,ρs[i],staples)
            end
        end

    end

    function apply_stout_smearing!(Uout::Array{<: AbstractGaugefields{NC,Dim},1},
        U::Array{<: AbstractGaugefields{NC,Dim},1},ρs::Array{T,1},staples) where {NC,Dim,T <: Number}
        @assert Uout != U "input U and output U should not be same!"
        V  = similar(U[1])
        temp1  = similar(U[1])
        temp2  = similar(U[1])
        temp3  = similar(U[1])
        F0 = initialize_TA_Gaugefields(U[1])

        num = length(ρs)

        for μ=1:Dim
            clear_U!(V)
            for i=1:num
                loops = staples[i][μ]
                evaluate_wilson_loops!(temp3,loops,U,[temp1,temp2])
                add_U!(V,ρs[i],temp3)
            end
            mul!(temp1,V,U[μ]') #U U*V
            clear_U!(F0)
            Traceless_antihermitian_add!(F0,1,temp1)

            exptU!(temp3,1,F0,[temp1,temp2])

            mul!(Uout[μ],temp3,U[μ])
        end
        set_wing_U!(Uout)

        #error("ee")

        #error("not implemented")
    end

    #=

    abstract type SmearingParam end

    struct SmearingParam_nosmearing <: SmearingParam
    end

    abstract type SmearingParam_single <: SmearingParam
    end

    abstract type SmearingParam_multi <: SmearingParam
    end

    mutable struct SmearingParam_stout <: SmearingParam_single
        staples_for_stout::Array{Array{Wilson_loop_set,1},1}
        tensor_derivative::Array{Tensor_derivative_set,1}
        staples_for_stout_dag::Array{Array{Wilson_loop_set,1},1}
        tensor_derivative_dag::Array{Tensor_derivative_set,1}
        ρs::Array{Float64,1}
        #ρs::Array{Float64,1}
    end

    mutable struct SmearingParam_stout_multi <: SmearingParam_multi
        staples_for_stout::Array{Array{Wilson_loop_set,1},1}
        tensor_derivative::Array{Tensor_derivative_set,1}
        staples_for_stout_dag::Array{Array{Wilson_loop_set,1},1}
        tensor_derivative_dag::Array{Tensor_derivative_set,1}
        ρs::Array{Array{Float64,1},1}
        #ρs::Array{Float64,1}
    end

    const Nosmearing = SmearingParam_nosmearing
    const Stout = SmearingParam_stout

    =#
end


text:
@testset "Evaluation" begin

    f = x -> DependencyTrees.from_indices(x, form=2, label=3, head=4, id=1)

    gold = deptree(f, [
        (1, "The", "DT", 2),
        (2, "cat", "NSUBJ", 3),
        (3, "sleeps", "PRED", 0),
        (4, ".", ".", 3)])

    good = deptree(f, [
        (1, "The", "DT", 2),
        (2, "cat", "NSUBJ", 3),
        (3, "sleeps", "PRED", 0),
        (4, ".", ".", 3)])

    ok = deptree(f, [
        (1, "The", "DT", 2),
        (2, "cat", "NSUBJ", 3),
        (3, "sleeps", "PRED", 0),
        (4, ".", "WRONG", 3)])

    bad = deptree(f, [
        (1, "The", "DT", 2),
        (2, "cat", "NSUBJ", 3),
        (3, "sleeps", "PRED", 0),
        (4, ".", "WRONG", 2)])

    worse = deptree(f, [
        (1, "The", "DT", 2),
        (2, "cat", "NSUBJ", 3),
        (3, "sleeps", "PRED", 2),
        (4, ".", "WRONG", 0)])

    @test unlabeled_accuracy(gold, good) == 1
    @test unlabeled_accuracy(gold, ok)   == 1
    @test unlabeled_accuracy(gold, bad)  == 0.75
    @test unlabeled_accuracy(gold, worse)  == 0.5

    @test labeled_accuracy(gold, good) == 1
    @test labeled_accuracy(gold, ok)   == 0.75
    @test labeled_accuracy(gold, bad)  == 0.75
    @test labeled_accuracy(gold, worse)  == 0.5

    @test unlabeled_accuracy([gold], [good]) == 1
    @test unlabeled_accuracy([gold], [ok])   == 1
    @test unlabeled_accuracy([gold], [bad])  == 0.75
    @test unlabeled_accuracy([gold], [worse])  == 0.5

    @test labeled_accuracy([gold], [good]) == 1
    @test labeled_accuracy([gold], [ok])   == 0.75
    @test labeled_accuracy([gold], [bad])  == 0.75
    @test labeled_accuracy([gold], [worse])  == 0.5
end


text:
# Curves

struct NullCurve <: Curve end

mutable struct FittingCost <: CostFunction
  # value::Vector{Float64}
  # values::Vector{Float64}
  firstCashFlow::Vector{Int}
  curve::Curve
end

function FittingCost(size::Int, curve::Curve)
  # value = Vector{Float64}()
  # values = Vector{Float64}()
  firstCashFlow = zeros(Int, size)

  return FittingCost(firstCashFlow, curve)
end

# Interpolated Curve methods #
max_date(curve::InterpolatedCurve) = curve.dates[end]

function discount(curve::Curve, t::Float64)
  return discount_impl(curve, t)
end

function discount_impl(curve::InterpolatedCurve, t::Float64)
  calculate!(curve)
  if t <= curve.times[end]
    return QuantLib.Math.value(curve.interp, t)
  end

  # println("outside")
  # println(t)
  # println(curve.times)

  # do flat fwd extrapolation
  return 0.0
end

function perform_calculations!(curve::InterpolatedCurve)
  _calculate!(curve.boot, curve)
  return curve
end

function value(cf::CostFunction, x::Vector{T}) where {T}
  ref_date = cf.curve.referenceDate
  dc = cf.curve.dc
  squared_error = 0.0
  n = length(cf.curve.bonds)

  # for (i, bh) in enumerate(cf.curve.bonds)
  @inbounds @simd for i in eachindex(cf.curve.bonds)
    bond = cf.curve.bonds[i].bond
    bond_settlement = get_settlement_date(bond)
    model_price = -accrued_amount(bond, bond_settlement)
    leg = bond.cashflows
    for k = cf.firstCashFlow[i]:length(leg.coupons)
      # @inbounds df = discount_function(cf.curve.fittingMethod, x, year_fraction(dc, ref_date, date(leg.coupons[k])))
      model_price += amount(leg.coupons[k]) * discount_function(cf.curve.fittingMethod, x, year_fraction(dc, ref_date, date(leg.coupons[k])))
    end

    # redemption
    if leg.redemption != nothing
      @inbounds model_price += amount(leg.redemption) * discount_function(cf.curve.fittingMethod, x, year_fraction(dc, ref_date, date(leg.redemption)))
    end

    # adjust NPV for forward settlement
    if bond_settlement != ref_date
      model_price /= discount_function(cf.curve.fittingMethod, x, year_fraction(dc, ref_date, bond_settlement))
    end

    market_price = bond.faceAmount
    price_error = model_price - market_price
    weighted_error = cf.curve.fittingMethod.commons.weights[i] * price_error
    squared_error += weighted_error * weighted_error
  end

  return squared_error
end


text:
using ArgParse

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--scaling"
            default = 0.002
            arg_type = Float64
        "--nsim"
            default = 10
            arg_type = Int64
        "label"
            default = "nothing"
    end

    return parse_args(s)
end

parameters = parse_commandline()
println(parameters)

label = parameters["label"]
scaling = parameters["scaling"]
nsim = parameters["nsim"]


text:
export mis2

"""
Solving MIS problem with sofisticated branching algorithm.
"""
function mis2(eg::EliminateGraph)
    if nv(eg) == 0
        #@show "0" # CHECKED
        return 0
    elseif nv(eg) == 1
        return 1
    elseif nv(eg) == 2
        return 2 - (@inbounds isconnected(eg, eg.vertices[end-1], eg.vertices[end]))
    elseif nv(eg) == 3
        @inbounds a, b, c = eg.vertices[end-2:end]
        nedge = isconnected(eg, a, b) + isconnected(eg, a, c) + isconnected(eg, b, c)
        if nedge == 0
            return 3
        elseif nedge == 3
            return 1
        else
            return 2
        end
    else
        #@show "1" # CHECKED
        vmin, degmin = mindegree_vertex(eg)
        if degmin == 0  # DONE
            #@show "1.1(1)" # CHECKED
            return 1 + eliminate(mis2, eg, vmin)
        elseif degmin == 1  # DONE
            #@show "1.1(2)" # CHECKED
            return 1 + eliminate(mis2, eg, NeighborCover(vmin))
        elseif degmin == 2
            #@show "1.2" # CHECKED
            a, b = neighbors(eg, vmin)
            if isconnected(eg, a, b)
                #@show "1.2.1" # CHECKED
                return 1 + eliminate(mis2, eg, NeighborCover(vmin))
            else
                #@show "1.2.2" # CHECKED
                sn = neighbors2(eg, vmin)
                # NOTE: there is no degree one vertex!
                if length(sn) == 1
                    #@show "1.2.2.1" # CHECKED
                    w = sn[1]
                    #return max(2+eliminate(mis2, eg, NeighborCover(w) ∪ Neighbors{CLOSED,2}(vmin)),
                                #2+eliminate(mis2, eg, Neighbors{CLOSED,2}(vmin)),
                    # Note: it seems it must choose the latter. Gurantted if one removes one vertex, the MIS is non-increasing.
                    return 2+eliminate(mis2, eg, (vmin, w, a, b)
                    )
                else
                    #@show "1.2.2.2" # CHECKED
                    return max(1+eliminate(mis2, eg, NeighborCover(vmin)),
                                eliminate(mis2, eg, MirrorCover(vmin)))
                end
            end
        elseif degmin == 3 # DONE
            #@show "1.3" #CHECKED
            a, b, c = neighbors(eg, vmin)
            nedge = isconnected(eg, a, b) + isconnected(eg, a, c) + isconnected(eg, b, c)
            if nedge == 0
                #@show "1.3.1" #CHECKED
                ms = mirrorcover(eg, vmin)
                if length(ms) > 1
                    #@show "1.3.1.1" # CHECKED
                    return max(1+eliminate(mis2, eg, NeighborCover(vmin)),
                                eliminate(mis2, eg, ms))
                else
                    #@show "1.3.1.2" # CHECKED
                    return max(1+eliminate(mis2, eg, NeighborCover(vmin)),
                                2 + eliminate(mis2, eg, NeighborCover(a) ∪ NeighborCover(b)),
                                2 + eliminate(mis2, eg, NeighborCover(a) ∪ NeighborCover(c) ∪ Vertex(b)),
                                2 + eliminate(mis2, eg, NeighborCover(b) ∪ NeighborCover(c) ∪ Vertex(a)),
                                )
                end
            elseif nedge == 3
                #@show "1.3.2" # CHECKED
                return 1 + eliminate(mis2, eg, NeighborCover(vmin))
            else
                #@show "1.3.3" # CHECKED
                return max(1 + eliminate(mis2, eg, NeighborCover(vmin)),
                            eliminate(mis2, eg, MirrorCover(vmin)))
            end
        else # DONE
            #@show "1.4" # CHECKED
            vmax, degmax = maxdegree_vertex(eg)
            if degmax >= 6 # DONE
                #@show "1.4.1"
                return max(1+eliminate(mis2, eg, NeighborCover(vmax)),
                            eliminate(mis2, eg, vmax))
            elseif !isconnected(eg) # DONE
                #@show "1.4.2" # CHECKED
                cluster = find_cluster(eg, vmax)
                A = subgraph(eg, cluster)
                B = subgraph(eg, setdiff(vertices(eg), cluster))
                return mis2(A) + mis2(B)
            elseif degmin == degmax  # DONE
                #@show "1.4.3" # CHECKED
                return max(1+eliminate(mis2, eg, NeighborCover(vmax)),
                            eliminate(mis2, eg, MirrorCover(vmax)))
            else
                #@show "1.4.4" # CHECKED
                v4, v5 = adjacent45(eg)
                return max(1+eliminate(mis2, eg, NeighborCover(v5)),
                            1+eliminate(mis2, eg, MirrorCover(v5) ∪ NeighborCover(v4)),
                            eliminate(mis2, eg, MirrorCover(v5) ∪ Vertex(v4))
                            )
            end
        end
    end
end

"""find adjacent vertice with degree 4, 5."""
function adjacent45(eg::EliminateGraph)
    for vi in vertices(eg)
        deg = degree(eg, vi)
        for vj in neighbors(eg, vi)
            degree(eg, vj) != deg && return (vi,vj)
        end
    end
end


text:
tobig(x) = parse(BigInt, bytes2hex(reverse(x)), base=16)

function int2bytes(x::Integer)
    hex = string(x, base=16)
    if mod(length(hex), 2) != 0
        hex = string("0", hex)
    end

    return reverse(hex2bytes(hex))
end

style(x, n) = "\33[1;$(n)m$x\33[0m"

bitlength(::Type{T}) where T <: Integer = T.size * 8

function bitlength(p::Integer)

    bits = bitstring(p)
    start = findfirst(x -> x == '1', bits)
    N = length(bits) - start + 1

    return N
end

function bitlength(p::BigInt)

    bytes = int2bytes(p)
    bits = bitstring(bytes[end])
    start = findfirst(x -> x == '1', bits)
    N = length(bytes) * 8  - (start - 1)

    return N
end

bytelength(x::Integer) = div(bitlength(x) + 1, 8, RoundUp)

interpret(::Type{BigInt}, x::Vector{UInt8}) = tobig(reverse(x))

function interpret(::Type{T}, x::Vector{UInt8}) where T <: Integer

    L = bitlength(T) ÷ 8
    y = UInt8[zeros(UInt8, L - length(x))..., x...]

    r = reinterpret(T, reverse(y))[1]
    return r
end

interpret(::Type{Vector{UInt8}}, x::BigInt) = reverse(int2bytes(x))

interpret(::Type{Vector{UInt8}}, x::Integer) = reverse(reinterpret(UInt8, [x])) # Number of bytes are useful for construction for bytes.

function interpret(::Type{Vector{T}}, 𝐫::Vector{UInt8}, N::Int) where T <: Integer
    M = length(𝐫) ÷ N
    𝐮 = reshape(𝐫, (M, N))
    𝐭 = [interpret(T, 𝐮[:, i]) for i in 1:N]
    return 𝐭
end

import Base: &

mutable struct Report
    data::Vector{Tuple{String, Union{Bool, BitVector}}} # It could actually be made optional
    state::Bool # Whether all previous expressions have been
end

Report() = Report(Tuple{String, Union{Bool, BitVector}}[], true)

function (&)(report::Report, x::Tuple{String, Bool})
    push!(report.data, x)
    (name, value) = x
    report.state &= value
    return report
end

function (&)(report::Report, x::Tuple{String, BitVector})
    push!(report.data, x)
    (name, values) = x
    report.state &= prod(values)
    return report
end

(&)(report::Report, x::Tuple{String, Vector{Bool}}) = report & (x[1], BitVector(x[2]))

_entry(io::IO, key::String, value::Bool) = println(io, "$key: $value")

function _entry(io::IO, key::String, values::BitVector)

    t = prod(values)
    _entry(io, key, t)

    if t == false
        for (i, value) in enumerate(values)
            println(io, "    $i: $value")
        end
    end
end

function Base.println(io::IO, report::Report)

    println(io, "Report: $(isvalid(report))")

    (; data) = report

    for (key, value) in data
        print(io, "  ")
        _entry(io, key, value)
    end
end

Base.isvalid(report::Report) = report.state


text:
__precompile__()

module IncoherentScatter

const ℜ(z)=Base.real(z)
const ℑ(z)=Base.imag(z)
const ℂ(r, i) = Base.complex(r, i)

using FFTW,SpecialFunctions,Polynomials,PhysicalConstants.CODATA2018,Unitful
using CUDA
using LoopVectorization,Tullio

export Ion,Plasma,CollisionalPlasma,ScatterVolume,pwrspec,autocor,autocor!

friedconte(z) = im*√(π)*erfcx(-im*z)

struct Faddeeva
    n::Int
    l::Float64
    a::Polynomial{Float64}
    function Faddeeva(n=16)
        m  = 2*n
        l  = √(n/√(2))

        theta = (-m+1:m-1)*π/m
        t = l*tan.(theta/2)
        f = [0; exp.(-t.^2).*(l*l .+ t.^2)]
        a = real.(fft(fftshift(f)))[2:n-1]/(2*m)   # no flip because of convention in polyval
        new(n, l, Polynomial(a))
    end
end

const deff = Faddeeva()

"""
    faddeeva(z, f=deff)

   The Faddeeva function is defined as:

   w(z) = exp(-z^2) * erfc(-im*z)

   where erfc(x) is the complex complementary error function. Also,

   friedconte(z) = im * √(π) * faddeeva(z)

   where friedconte(z) is the plasma dispersion function as defined in [2].

   w = faddeeva(z, f) reuses previously calculated coefficients in struct f
   where f is initialized like f=Faddeeva(n) with n specifying the number of terms
   to truncate the expansion (see (13) in [1]). n = 16 is used as default, i.e. for f=Faddeeva().

   Example:
       x = linspace(-10,10,1001); [X,Y] = meshgrid(x,x);
       W = faddeeva(complex(X,Y));

   Reference:
   [1] J.A.C. Weideman, "Computation of the Complex Error Function," SIAM
       J. Numerical Analysis, pp. 1497-1518, No. 5, Vol. 31, Oct., 1994
       Available Online: http://www.jstor.org/stable/2158232

   [2] B.D. Fried, S.D. Conte, The plasma dispersion function. New York Academic Press, 1961.
"""
function faddeeva(z, f::Faddeeva=deff)
    if ℜ(z)==0      # purely imaginary input
        exp.(-z.^2).*erfc(ℑ(z))
    else               # complex valued input
        # make sure input in the upper half-plane (positive imag. values)
        zz = ℑ(z)<0 ? conj(z) : z
        yy = (f.l + im*zz)/(f.l - im*zz)
        w  = 2*polyval(f.a, yy)/(f.l - im*zz)^2 + (1/√(π))/(f.l - im*zz)
        # convert the upper half-plane results to the lower half-plane if necessary
        ℑ(z)<0 ? conj(2*exp(-zz^2) - w) : w
    end
end

const c  = ustrip(u"m/s", SpeedOfLightInVacuum)
const me = ustrip(u"kg",  ElectronMass)
const e  = ustrip(u"C",   ElementaryCharge)
const kB = ustrip(u"J/K", BoltzmannConstant)
const ϵ0 = ustrip(u"F/m", VacuumElectricPermittivity)
const amu= ustrip(u"kg",  AtomicMassConstant)

const fnormf = c*√(me/(2*kB))/2
const debyef = (4*π/(e*c))^2*ϵ0*kB

const mp = ustrip(u"kg",  ProtonMass)
const mn = ustrip(u"kg",  NeutronMass)
const sqrtmpme = √((mp+mn)/me/2)

"gfe(B in Tesla) -> electron gyrofrequency in Hz"
gfe(B) = e*B/me

"""
    ScatterVolume(tfeff, system, b, alpha)

has parameters of the scattering volume that are constant in a fit.
"""
mutable struct ScatterVolume{T<:Real}
    tfeff::T    # effective transmitter frequency in Hz = tf*cos(scattering angle)
    fnormc::T   # factor for frequency normalization f -> theta
    debyec::T   # Debye length correction factor

    system::T   # spectrum scale factor (depends on tx power, volume size, ... )

    b::T        # magnetic field strength in Tesla
    gfe::T      # electron gyrofrequency
    phiel::T    # normalized electron gyrofrequency
    alpha::T    # angle between b and scattering vector in rad
    sinfac::T   # correction factor for magnetic field
    cosal::T    #

    function ScatterVolume{T}(tfeff, system, b, alpha) where {T<:Real}
        fnormc = fnormf/tfeff
        debyec = debyef*tfeff^2
        gf     = gfe(b)
        phiel  = fnormc*gf
        sinfac = 0.5*(sin(alpha)/phiel)^2
        cosal  = cos(alpha)
        new(tfeff, fnormc, debyec, system, b, gf, phiel, alpha, sinfac, cosal)
    end
end
# ScatterVolume(tfeff::T, system, b, alpha) where {T<:Real} = ScatterVolume{T}(tfeff, system, b, alpha)
ScatterVolume(;tfeff::T=230e6, system=1e-6, b=50000e-9, alpha=deg2rad(30)) where {T<:Real} =
    ScatterVolume{T}(tfeff, system, b, alpha)

struct Ion{T<:Real}
    mass::T   # ion mass in amu
    srmr::T   # sqrt of mass ratio ion/electron
end
Ion(m=15.999) = Ion(m, √(m*amu/me))

"""
Plasma(ne, te, ti, ::ScatterVolume, [ions, fractn]

has parameters of the plasma that are constant over all frequencies of the spectrum.
"""
struct Plasma{T<:Real}
    ne::T
    te::T
    ti::T
    fractn1::T
    fractn::NTuple{2, T}    # up to two more ions
    ions::NTuple{3, Ion{T}}
    nion::Int

    spectrsf::T
    debyen::T   # normalized Debye length correction factor
    sqte::T
    sqti::T
    thetc::T
    tex::T
    tefac::T
    dpsidc::T

    # parameter for spectrum at zero frequency:
    cdtine::T
    ym::T

    function Plasma{T}(ne::T, te::T, ti::T, s::ScatterVolume{T}, fractn, ions, nion) where {T<:Real}
        sqti    = √(ti)
        thetc   = s.fnormc/sqti
        tex     = exp(-te*s.sinfac)
        tefac   = 0.5/te + s.sinfac/tex
        fractn1 = one(ne) - fractn[1] - fractn[2]
        dpsidc  = thetc*ions[1].srmr
        cdtine  = ti*s.debyec + ne
        ym      = ti*ne + te*cdtine
        new(ne, te, ti, fractn1, fractn, ions, nion,
            s.system*ne, s.debyec/ne, √(te), sqti, thetc,
            tex, tefac, dpsidc, cdtine, ym)
    end
end
Plasma(ne::T, te::T, ti::T,
       s::ScatterVolume{T}, nion=1, ions=(Ion(16*one(ne)), Ion(31*one(ne)), Ion(one(ne))),
       fractn=(zero(ne),zero(ne))) where {T<:Real} =
    Plasma{T}(ne, te, ti, s, fractn, ions, nion)
Plasma(s::ScatterVolume{T}; ne::T, te::T, ti::T,
       nion=1, ions=(Ion(16*one(ne)), Ion(31*one(ne)), Ion(one(ne))),
       fractn=(zero(ne),zero(ne))) where {T<:Real} =
    Plasma{T}(ne, te, ti, s, fractn, ions, nion)

struct CollisionalPlasma{T<:Real}
    p::Plasma{T}
    ν::T
end
CollisionalPlasma(p; ν) = CollisionalPlasma(p, ν)

"""
    ionadm(theta, z=friedconte(theta))
-> unmagnetized ion admittance for real normalized frequency
"""
ionadm(theta::T, z::Complex{T}=friedconte(theta)) where {T<:Real} =
    ℂ(theta*ℑ(z), 1 + theta*ℜ(z))

"""
    diadt1(::Plasma, ionadm::Complex)
-> derivative with respect to Ti of the unmagnetized ion admittance, 1st ion, zero frequency
"""
diadt1(p::Plasma{T}, ia::Complex{T}) where {T<:Real} =
    ℂ(0, p.fractn1*p.dpsidc*(√(π)*ℜ(ia) - 2*p.ions[1].srmr))

"""
    diadt1(θ, ::Plasma, ionadm::Complex)
-> derivative with respect to Ti of the unmagnetized ion admittance, 1st ion, real normalized frequency
"""
function diadt1(θ::T, p::Plasma{T}, ia::Complex{T}) where {T<:Real}
    r5ti = T(0.5)/p.ti
    p.fractn1*((θ^2/p.ti - r5ti)*ia + im*r5ti)
end

function diadc(theta::T, ia::Complex{T}=ionadm(theta)) where {T<:Real}
    th =1/theta
    ℂ(ℜ(ia)^2 + ℑ(ia)*(1 - ℑ(ia)))*th + ℑ(ia)*(th - 2*theta) - th,
      ℜ(ia)*(theta + (ℑ(ia) - 1)*th)
end

"""
    ionadm(freq, ::Plasma)
-> unmagnetized ion admittance for real unnormalized frequency and plasma parameters
"""
function ionadm(freq::T, p::Plasma{T}, s::ScatterVolume{T}) where {T<:Real}
    fnormt = s.fnormc/p.sqti
    thetai = fnormt*p.ions[1].srmr*freq
    yi     = ionadm(thetai)
    for k in 2:p.nion
        if k==2
            yi *= p.fractn1
        end
        thetai = fnormt*p.ions[k].srmr*freq
        yi    += p.fractn[k-1]*ionadm(thetai)
    end
    yi
end

"""
    ionadm(theta, z=friedconte(theta))
-> unmagnetized ion admittance for complex normalized frequency.
"""
function ionadm(theta::T, z::T=friedconte(theta)) where {T<:Complex}
    psizm = ℑ(theta)*(ℜ(z)^2 + ℑ(z)^2)
    yd    = 1 - ℑ(theta)*(2*ℑ(z) - psizm)
    ℂ(ℜ(theta)*(ℑ(z) - psizm)/yd, 1 + ℜ(theta)*ℜ(z)/yd)
end

"""
    ionadm(freq, ::CollisionalPlasma)
-> unmagnetized ion admittance for unnormalized frequency and collisional plasma parameters
"""
function ionadm(freq::T, coll::CollisionalPlasma{T}, s::ScatterVolume{T}) where {T<:Real}
    fnormt = s.fnormc/coll.p.sqti
    thetai = fnormt*coll.p.ions[1].srmr*ℂ(freq, coll.ν)
    yi     = ionadm(thetai)
    for k in 2:coll.p.nion
        if k==2
            yi *= coll.p.fractn1
        end
        thetai = fnormt*coll.p.ions[k].srmr*ℂ(freq, coll.ν)
        yi    += coll.p.fractn[k-1]*ionadm(thetai)
    end
    yi
end

"""
    elecadm(theta, tex, z)
-> magnetized electron admittance for real normalized frequency
"""
function elecadm(theta::T, tex::T, z::Complex{T}=friedconte(theta)) where {T<:Real}
    tt = tex*theta
    ℂ(tt*imag(z), 1 + tt*real(z)), imag(z)
end

# helper function, this part of the spectrum code is the same for non-collisional and collisional
"""
    espec(freq, ::Plasma, ::ScatterVolume, yi::Complex)
-> power spectral density for non-zero real frequency
"""
function espec(freq::T, p::Plasma{T}, s::ScatterVolume{T}, yi::Complex{T}) where {T<:Real}
    yiti   = yi/p.ti + im*p.debyen

    # electron admittance:
    thetae = s.fnormc/(p.sqte*s.cosal)*freq
    ye,    = elecadm(thetae, p.tex)
    yete   = ye/p.te

    the2   = thetae^2
    texthe = p.tex*thetae

    ym     = abs2(yete+yiti)
    yn     = abs2(yete)*real(yi) + abs2(yiti)*real(ye)

    (p.spectrsf/freq)*(yn/ym)
end

"""
    pwrspec(freq, ::Plasma, ::ScatterVolume)
-> power spectral density for non-zero real frequency
"""
function pwrspec(freq::T, p::Plasma{T}, s::ScatterVolume{T}) where {T<:Real}
    yi = ionadm(freq, p, s)
    espec(freq, p, s, yi)
end

"""
    pwrspec(freq, ::CollisionalPlasma, ::ScatterVolume)
-> power spectral density for non-zero frequency and collisional plasma
"""
function pwrspec(freq::T, coll::CollisionalPlasma{T}, s::ScatterVolume{T}) where {T<:Real}
    yi = ionadm(freq, coll, s)
    espec(freq, coll.p, s, yi)
end

# helper function
function yn0(p::Plasma{T}, s::ScatterVolume{T}, yir::T) where {T<:Real}
    tineq  = p.sqti*p.ne^2*yir*p.ti
    # sqtec  = p.cdtine*p.sqte*√(π)/s.cosal2 #  divide by colal, not colal*colal:
    sqtec  = p.cdtine*p.sqte*√(π)/s.cosal # zi(1) in fortran code = sqrt(π) for real zero frequency
    stecx  = sqtec*p.tex*p.te
    tineq + stecx*p.cdtine
end

"""
    pwrspec(::Plasma)
-> power spectral density for zero real frequency
"""
function pwrspec(p::Plasma, s::ScatterVolume)
    yir = sqrt(π)*p.ions[1].srmr
    for k in 2:p.nion
        if k==2
            yir *= p.fractn1
        end
        yir += p.fractn[k-1]*sqrt(π)*p.ions[k].srmr
    end
    yn = yn0(p, s, yir)
    s.fnormc*p.spectrsf*yn/p.ym^2
end

"""
    pwrspec(::Plasma, ::ScatterVolume, freq_itr)
-> power spectral density for iterable frequencies (with length)
"""
function pwrspec(p::Union{Plasma{T}, CollisionalPlasma{T}}, s::ScatterVolume{T}, freqs;
                 threads=true) where {T<:Real}
    if threads
        pwr = Vector{T}(undef, length(freqs))
        if freqs[1]==0
            rng = 2:length(freqs)
            pwr[1] = pwrspec(p, s)
        else
            rng = 1:length(freqs)
        end
        Threads.@threads for k ∈ rng
            pwr[k] = pwrspec(freqs[k], p, s)
        end
    else
        if freqs[1]==0
            pwr = vcat(pwrspec(p, s), @inbounds pwrspec.(freqs[2:end], [p], [s]))
        else
            pwr = @inbounds pwrspec.(freqs, [p], [s])
        end
    end
    return pwr
end

function yirc(coll::CollisionalPlasma, kion::Int)
    psiion = coll.p.dpsidc*coll.ν
    z      = friedconte(ℂ(zero(coll.p.ne), psiion))
    psiz   = 1 - psiion*ℑ(z)
    psizm  = coll.p.ions[kion].srmr/psiz
    psizm*ℑ(z)
end

"""
    pwrspec(::CollisionalPlasma)
-> power spectral density for zero frequency and collisional plasma
"""
function pwrspec(coll::CollisionalPlasma, s::ScatterVolume)
    yir = yirc(coll, 1)
    for k in 2:coll.p.nion
        if k==2
            yir *= coll.p.fractn1
        end
        z    = friedconte(ℂ(zero(coll.p.ne), p.thetc*coll.p.ions[k].srmr*coll.ν))
        yir += coll.p.fractn[k-1]*yirc(coll, k)
    end

    yn = yn0(coll.p, s, yir)
    s.fnormc*coll.p.spectrsf*yn/coll.p.ym^2
end

"""
    cosine(freqs, τ)
-> matrix for the transform spectrum->ACF, method by Wes Swartz, ISCATSPE.FOR
"""
function cosine(freqs, τ)
    c           = Matrix{eltype(freqs)}(undef, length(freqs), length(τ))
    c[    1, 1] = freqs[2]/2
    for k ∈ 2:length(freqs)-1
        c[k, 1] = (freqs[k+1]-freqs[k-1])/2
    end
    c[end, 1] = (freqs[end] - freqs[end-1])/2
    df        = diff(freqs)
    for l ∈ 2:length(τ)
        τ2    = 2*π*τ[l]
        c[1, l]       = (1 - cos(τ2*freqs[2]))/(τ2*τ2*freqs[2])
        c2            =  cos.(τ2*freqs[2:end-1])
        c[2:end-1, l] = ((c2 - cos.(τ2*freqs[1:end-2]))./df[1:end-1] +
                         (c2 - cos.(τ2*freqs[3:end  ]))./df[2:end  ])/(τ2*τ2)
        cend = cos(τ2*freqs[end])
        fend = 2*freqs[end]-freqs[end-1]
        c[end, l]     =  (2*cend - cos(τ2*freqs[end-1]) - cos(τ2*fend))/(df[end]*τ2*τ2)
    end
    return c
end

struct FT{T<:Real}
    cft::Matrix{T}
end
FT(freqs, τ) = FT(cosine(freqs, τ))

function pwr2acf(pwr::Vector{T}, ftm::Matrix{T}) where {T<:Real}
    @tullio acf[l] := ftm.cft[k, l]*pwr[k]
end

function autocor(cft::Matrix{T}, pwr::Vector{T}) where {T<:Real}
    @tullio acf[l] := cft[k, l]*pwr[k]
end

function autocor!(acf::Vector{T}, cft::Matrix{T}, pwr::Vector{T}) where {T<:Real}
    @tullio acf[l] = cft[k, l]*pwr[k]
end

greet() = print("Hello World!")

end   # module IncoherentScatter


text:
# https://adventofcode.com/2021/day/4
using AdventOfCode
input = begin
    inp=readlines("2021/data/day_4.txt")
    (draws=parse(inp[begin]),
    boards=replace(["[[";inp[begin+2:end].*";";"]]"],";"=>"],[") |>join|>parse)
end

function mark(draw, board)
    markidx=findfirst(x->x==draw,board)
    if isnothing(markidx) return(0) end
    board[markidx]=-1
    if all(<(0),(@view board[markidx[1],:])) || all(<(0),(@view board[:,markidx[2]])) 1 # win
        else 0 end
end
function part_1(input)
    boards=copy.(input.boards) # not to mutate the original input
    res=zeros(size(boards))
    score=0
    for draw in input.draws
        res.+=mark.(draw,boards)
        if any(>(0),res)
            score=sum(max.(boards[findfirst(>(0),res)],0))*draw
            break
        end
    end
    score
end

function part_2(input)
    boards=copy.(input.boards) # not to mutate the original input
    res=zeros(size(boards))
    score=0
    for draw in input.draws
        res.+=mark.(draw,boards)
        if all(>(0),res)
            score=sum(max.(boards[findmin(res)[2]],0))*draw
            break
        end
    end
    score
end


text:
"""
    set_scheduler_backup(
        filepath::AbstractString = "";
        migrate::Bool = false,
        delete_old::Bool = false,
        recover_settings::Bool = true,
        recover_queue::Bool = true
    )

Set the backup file of job scheduler.

If `filepath` was set to `""`, stop backup at exit.

If `filepath` was set to an existing file, `recover_settings` or `recover_queue` from `filepath` immediately.

If `filepath` was set to a new file, the backup file will be created at exit.

If `migrate=true` and the old `JobSchedulers.SCHEDULER_BACKUP_FILE` exists, the old backup file will be recovered before recovering from `filepath`.

"""
function set_scheduler_backup(
    filepath::AbstractString = "";
    migrate::Bool = false,
    delete_old::Bool = false,
    recover_settings::Bool = true,
    recover_queue::Bool = true
)
    global SCHEDULER_BACKUP_FILE
    if filepath == ""
        @info "Stop scheduler backup."
        delete_old && isfile(SCHEDULER_BACKUP_FILE) && rm(SCHEDULER_BACKUP_FILE, force=true)
        SCHEDULER_BACKUP_FILE = ""
        return SCHEDULER_BACKUP_FILE
    end

    filepath = abspath(filepath)

    if isdir(filepath) || isdirpath(filepath)
        @error "The file path is a directory." filepath
        return SCHEDULER_BACKUP_FILE
    end

    try
        mkpath(dirname(filepath))
    catch
        @error "Cannot set backup file: cannot create directory."
        return SCHEDULER_BACKUP_FILE
    end

    if migrate
        if isfile(SCHEDULER_BACKUP_FILE) # old one
            recover_backup(SCHEDULER_BACKUP_FILE; recover_settings = recover_settings, recover_queue = recover_settings)
        end
    end

    if isfile(filepath)
        recover_backup(filepath; recover_settings = recover_settings, recover_queue = recover_settings)
    end

    delete_old && isfile(SCHEDULER_BACKUP_FILE) && rm(SCHEDULER_BACKUP_FILE, force=true)
    SCHEDULER_BACKUP_FILE = filepath
    return SCHEDULER_BACKUP_FILE
end

"""
    backup()

Manually backup job scheduler settings and queues. The function is automatically triggered at exit.
"""
function backup()
    global SCHEDULER_BACKUP_FILE
    global SCHEDULER_MAX_CPU
    global SCHEDULER_MAX_MEM
    global SCHEDULER_UPDATE_SECOND
    global JOB_QUEUE_MAX_LENGTH
    global JOB_QUEUE
    global JOB_QUEUE_OK
    global QUEUING
    global RUNNING
    global CANCELLED

    # no backup when file is not exist
    SCHEDULER_BACKUP_FILE == "" && return

    scheduler_max_cpu = SCHEDULER_MAX_CPU
    scheduler_max_mem = SCHEDULER_MAX_MEM
    scheduler_update_second = SCHEDULER_UPDATE_SECOND
    job_queue_max_length = JOB_QUEUE_MAX_LENGTH

    wait_for_lock()
        # update state of jobs for the last time
        foreach(unsafe_update_state!, JOB_QUEUE)

        # job queue: set to cancelled
        job_queue = deepcopy(JOB_QUEUE)
        for job in job_queue
            job.task = nothing
            if job.state in (QUEUING, RUNNING)
                job.state = CANCELLED
            end
        end

        job_queue_ok = deepcopy(JOB_QUEUE_OK)
        for job in job_queue_ok
            job.task = nothing
        end
        append!(job_queue_ok, job_queue)
    release_lock()

    # clean old file
    rm(SCHEDULER_BACKUP_FILE, force=true)
    @save SCHEDULER_BACKUP_FILE scheduler_max_cpu scheduler_max_mem scheduler_update_second job_queue_max_length job_queue_ok
    @info "Scheduler backuped to $SCHEDULER_BACKUP_FILE"
end

atexit(backup)

"""
    recover_backup(filepath::AbstractString; recover_settings = true, recover_queue = true)

Recover job scheduler settings or job queues from file.
"""
function recover_backup(filepath::AbstractString; recover_settings::Bool = true, recover_queue::Bool = true)
    global SCHEDULER_BACKUP_FILE
    global SCHEDULER_MAX_CPU
    global SCHEDULER_MAX_MEM
    global SCHEDULER_UPDATE_SECOND
    global JOB_QUEUE_MAX_LENGTH
    global JOB_QUEUE
    global JOB_QUEUE_OK

    if !(recover_settings || recover_queue)
        return
    end

    if !is_valid_backup_file(filepath)
        @error "Cannot recover backup from an invalid file: $(filepath)"
        return
    end

    @load filepath scheduler_max_cpu scheduler_max_mem scheduler_update_second job_queue_max_length job_queue_ok

    if recover_settings
        @info "Settings recovered from the backup file ($filepath)"
        set_scheduler_max_cpu(scheduler_max_cpu)
        set_scheduler_max_mem(scheduler_max_mem)
        set_scheduler_update_second(scheduler_update_second)
        set_scheduler_max_job(job_queue_max_length)
    end

    if recover_queue
        wait_for_lock()
            current_jobs = Dict{Int64, Vector{Job}}()
            append_jobs_dict!(current_jobs, JOB_QUEUE)
            append_jobs_dict!(current_jobs, JOB_QUEUE_OK)

            job_indices_to_copy = Int[]

            for (ind, job) in enumerate(job_queue_ok)
                if !has_job_in(current_jobs, job)
                    push!(job_indices_to_copy, ind)
                end
            end
            @info "$(length(job_indices_to_copy)) jobs recovered from the backup file ($filepath)."
            pushfirst!(JOB_QUEUE_OK, view(job_queue_ok, job_indices_to_copy)...)
        release_lock()
    end

    nothing
end

function append_jobs_dict!(jobs_dict::Dict{Int64, Vector{Job}}, job_queue::Vector{Job})
    for job in job_queue
        job_vec = get(jobs_dict, job.id, nothing)
        if isnothing(job_vec)  # no key, create
            jobs_dict[job.id] = [job]
        elseif has_job_in(job_vec, job)  # has same job, skip
            continue
        else  # push
            push!(job_vec, job)
        end
    end
end

function is_same_job(a::Job, b::Job)
    a.id == b.id &&
    a.schedule_time == b.schedule_time &&
    a.create_time == b.create_time &&
    a.start_time == b.start_time &&
    a.stop_time == b.stop_time &&
    a.name == b.name &&
    a.user == b.user &&
    a.ncpu == b.ncpu &&
    a.mem == b.mem &&
    a.wall_time == b.wall_time &&
    a.priority == b.priority
end

function has_job_in(vec::Vector{Job}, job::Job)
    for j in vec
        if is_same_job(j, job)
            return true
        end
    end
    false
end

function has_job_in(jobs_dict::Dict{Int64, Vector{Job}}, job)
    job_vec = get(jobs_dict, job.id, nothing)
    if isnothing(job_vec)  # no key
        false
    elseif has_job_in(job_vec, job)  # has same job
        true
    else
        false
    end
end

function is_valid_backup_file(filepath::AbstractString)
    REQUIRED_FILE_HEADER = "HDF5-based Julia Data Format, version "
    io = open(filepath, "r")
    headermsg = String(read!(io, Vector{UInt8}(undef, length(REQUIRED_FILE_HEADER))))
    headermsg == REQUIRED_FILE_HEADER
end


text:
using Dates

mutable struct DiscretizedDiscountBond{L <: Lattice} <: DiscretizedAsset
  common::DiscretizedAssetCommon{L}
end

DiscretizedDiscountBond(lattice::L) where {L <: Lattice} = DiscretizedDiscountBond{L}(DiscretizedAssetCommon(lattice))

function reset!(dBond::DiscretizedDiscountBond, sz::Int)
  dBond.common.values = ones(sz)

  return dBond
end

mutable struct DiscretizedCallableFixedRateBond{L <: Lattice} <: DiscretizedAsset
  redemptionTime::Float64
  couponTimes::Vector{Float64}
  callabilityTimes::Vector{Float64}
  args::CallableBondArgs
  common::DiscretizedAssetCommon{L}
end

function DiscretizedCallableFixedRateBond(bond::CallableFixedRateBond, referenceDate::Date, dc::DayCount, lattice::L) where {L <: Lattice}
  args = CallableBondArgs(bond)
  redemptionTime = year_fraction(dc, referenceDate, args.redemptionDate)
  couponTimes = Vector{Float64}(undef, length(args.couponDates))

  # for i in eachindex(couponTimes)
  #   couponTimes[i] = year_fraction(dc, referenceDate, args.couponDates[i])
  # end
  map!(x -> year_fraction(dc, referenceDate, x), couponTimes, args.couponDates)

  callabilityTimes = Vector{Float64}(undef, length(args.callabilityDates))

  # for i in eachindex(callabilityTimes)
  #   callabilityTimes[i] = year_fraction(dc, referenceDate, args.callabilityDates[i])
  # end
  map!(x -> year_fraction(dc, referenceDate, x), callabilityTimes, args.callabilityDates)
  # similar to tree swaption engine, we collapse similar coupon and exercise dates to avoid mispricing.
  # Can remove if necessary

  for i in eachindex(callabilityTimes), j in eachindex(couponTimes)
    if within_next_week(callabilityTimes[i], couponTimes[j])
      couponTimes[j] = callabilityTimes[i]
    end
  end

  return DiscretizedCallableFixedRateBond{L}(redemptionTime, couponTimes, callabilityTimes, args, DiscretizedAssetCommon(lattice))
end

function reset!(dcb::DiscretizedCallableFixedRateBond, sz::Int)
  dcb.common.values = fill(dcb.args.redemption, sz)

  adjust_values!(dcb)
  return dcb
end

function mandatory_times(dcb::DiscretizedCallableFixedRateBond)
  times =  dcb.couponTimes[dcb.couponTimes .>= 0.0]
  if dcb.redemptionTime >= 0.0
    insert!(times, 1, dcb.redemptionTime)
  end
  callTimes = dcb.callabilityTimes[dcb.callabilityTimes .>= 0.0]
  times = vcat(times, callTimes)
  return times
end

function apply_callability!(dcb::DiscretizedCallableFixedRateBond, ::Call, i::Int)
  # for j in eachindex(dcb.common.values)
  #   dcb.common.values[j] = min(dcb.args.callabilityPrices[i], dcb.common.values[j])
  # end
  dcb.common.values = min.(dcb.args.callabilityPrices[i], dcb.common.values)
  return dcb
end

function apply_callability!(dcb::DiscretizedCallableFixedRateBond, ::Put, i::Int)
  # for j in eachindex(dcb.common.values)
  #   dcb.common.values[j] = max(dcb.common.values[j], dcb.args.callabilityPrices[i])
  # end
  dcb.common.values = max(dcb.args.callabilityPrices[i], dcb.common.values)
  return dcb
end

function add_coupon!(dcb::DiscretizedCallableFixedRateBond, i::Int)
  dcb.common.values .+= dcb.args.couponAmounts[i]

  return dcb
end

apply_callability!(dcb::DiscretizedCallableFixedRateBond, i::Int) = apply_callability!(dcb, dcb.args.putCallSchedule[i].optionType, i)

function post_adjust_values_impl!(dcb::DiscretizedCallableFixedRateBond)
  for i in eachindex(dcb.callabilityTimes)
    t = dcb.callabilityTimes[i]
    if t >= 0.0 && is_on_time(dcb, t)
      apply_callability!(dcb, i)
    end
  end

  for i in eachindex(dcb.couponTimes)
    t = dcb.couponTimes[i]
    if t >= 0.0 && is_on_time(dcb, t)
      add_coupon!(dcb, i)
    end
  end

  return dcb
end