text:
using Distributions
using LinearAlgebra
using Random
using Test

using EnsembleKalmanProcesses.EnsembleKalmanProcessModule
using EnsembleKalmanProcesses.ParameterDistributionStorage
@testset "EnsembleKalmanProcessModule" begin

    # Seed for pseudo-random number generator
    rng_seed = 42
    Random.seed!(rng_seed)

    ### Generate data from a linear model: a regression problem with n_par parameters
    ### and 1 observation of G(u) = A \times u, where A : R^n_par -> R^n_obs
    n_obs = 10                  # dimension of synthetic observation from G(u)
    n_par = 2                  # Number of parameteres
    u_star = [-1.0, 2.0]          # True parameters
    noise_level = 0.05            # Defining the observation noise level (std)
    Γy = noise_level^2 * Matrix(I, n_obs, n_obs) # Independent noise for synthetic observations
    noise = MvNormal(zeros(n_obs), Γy)
    C = [1 -.9; -.9 1]          # Correlation structure for linear operator
    A = rand(MvNormal(zeros(2,), C), n_obs)'    # Linear operator in R^{n_par x n_obs}

    @test size(A) == (n_obs, n_par)

    #### Define linear model
    function G(u)
        A * u
    end

    y_star = G(u_star)
    y_obs = y_star .+ rand(noise)

    @test size(y_obs) == (n_obs,)

    # sum(y-G)^2 ~ n_obs*noise_level^2
    @test isapprox(norm(y_obs .- G(u_star))^2 - n_obs * noise_level^2, 0; atol = 0.05)

    #### Define prior information on parameters
    prior_distns = [Parameterized(Normal(0.0, 0.5)), Parameterized(Normal(3.0, 0.5))]
    constraints = [[no_constraint()], [no_constraint()]]
    prior_names = ["u1", "u2"]
    prior = ParameterDistribution(prior_distns, constraints, prior_names)

    prior_mean = get_mean(prior)

    # Assuming independence of u1 and u2
    prior_cov = get_cov(prior) #convert(Array, Diagonal([sqrt(2.), sqrt(2.)]))

    ###
    ###  Calibrate: Ensemble Kalman Sampler
    ###
    N_ens = 50 # number of ensemble members (NB for @test throws, make different to N_ens)
    N_iter = 20

    initial_ensemble = EnsembleKalmanProcessModule.construct_initial_ensemble(prior, N_ens; rng_seed = rng_seed)
    @test size(initial_ensemble) == (n_par, N_ens)

    eksobj =
        EnsembleKalmanProcessModule.EnsembleKalmanProcess(initial_ensemble, y_obs, Γy, Sampler(prior_mean, prior_cov))

    g_ens = G(get_u_final(eksobj))
    @test size(g_ens) == (n_obs, N_ens)
    # as the columns of g are the data, this should throw an error
    g_ens_t = permutedims(g_ens, (2, 1))

    # EKS iterations
    for i in 1:N_iter
        params_i = get_u_final(eksobj)
        g_ens = G(params_i)
        if i == 1
            g_ens_t = permutedims(g_ens, (2, 1))
            @test_throws DimensionMismatch EnsembleKalmanProcessModule.update_ensemble!(eksobj, g_ens_t)
        end

        EnsembleKalmanProcessModule.update_ensemble!(eksobj, g_ens)
    end
    # Plot evolution of the EKS particles
    eks_final_result = vec(mean(get_u_final(eksobj), dims = 2))

    if TEST_PLOT_OUTPUT
        gr()
        p = plot(
            get_u_prior(eksobj)[1, :],
            get_u_prior(eksobj)[2, :],
            seriestype = :scatter,
            label = "Initial ensemble",
        )
        plot!(get_u_final(eksobj)[1, :], get_u_final(eksobj)[2, :], seriestype = :scatter, label = "Final ensemble")
        plot!(
            [u_star[1]],
            xaxis = "u1",
            yaxis = "u2",
            seriestype = "vline",
            linestyle = :dash,
            linecolor = :red,
            label = :none,
        )
        plot!([u_star[2]], seriestype = "hline", linestyle = :dash, linecolor = :red, label = :none)
        savefig(p, "EKS_test.png")
    end

    ###
    ###  Calibrate: Ensemble Kalman Inversion
    ###

    #N_ens = 40 # number of ensemble members
    # We will try to compare the EKI and EKS results, therefore we need comparable timesteps and iterations
    #N_iter = 7 # number of EKI iterations
    Δ_vec = eksobj.Δt[2:end]

    ekiobj = EnsembleKalmanProcessModule.EnsembleKalmanProcess(initial_ensemble, y_obs, Γy, Inversion())

    # some checks
    g_ens = G(get_u_final(ekiobj))
    @test size(g_ens) == (n_obs, N_ens)
    # as the columns of g are the data, this should throw an error
    g_ens_t = permutedims(g_ens, (2, 1))
    @test_throws DimensionMismatch find_ekp_stepsize(ekiobj, g_ens_t)
    Δ = find_ekp_stepsize(ekiobj, g_ens)
    # huge collapse for linear problem so should find timestep should < 1
    @test Δ < 1
    # NOTE We don't use this info, this is just for the test.

    # EKI iterations
    params_i_vec = []
    g_ens_vec = []
    for i in 1:N_iter
        params_i = get_u_final(ekiobj)
        push!(params_i_vec, params_i)
        g_ens = G(params_i)
        push!(g_ens_vec, g_ens)
        if i == 1
            g_ens_t = permutedims(g_ens, (2, 1))
            @test_throws DimensionMismatch EnsembleKalmanProcessModule.update_ensemble!(ekiobj, g_ens_t)
        end
        Δ = Δ_vec[i] # use the same timesteps as EKS (EKI is more stable, so this is fine to do)
        EnsembleKalmanProcessModule.update_ensemble!(ekiobj, g_ens, Δt_new = Δ)

    end
    @test isapprox(ekiobj.Δt - eksobj.Δt, zeros(N_iter + 1))

    push!(params_i_vec, get_u_final(ekiobj))

    @test get_u_prior(ekiobj) == params_i_vec[1]
    @test get_u(ekiobj) == params_i_vec
    @test get_g(ekiobj) == g_ens_vec
    @test get_g_final(ekiobj) == g_ens_vec[end]
    @test get_error(ekiobj) == ekiobj.err

    # EKI results: Test if ensemble has collapsed toward the true parameter
    # values
    eki_final_result = vec(mean(get_u_final(ekiobj), dims = 2))

    # kind of arbitrary tolerance here,
    @test norm(u_star - eki_final_result) < 0.1

    # Plot evolution of the EKI particles
    #eki_final_result = vec(mean(get_u_final(ekiobj), dims = 2))
    if TEST_PLOT_OUTPUT
        gr()
        p = plot(
            get_u_prior(ekiobj)[1, :],
            get_u_prior(ekiobj)[2, :],
            seriestype = :scatter,
            label = "Initial ensemble",
        )
        plot!(get_u_final(ekiobj)[1, :], get_u_final(ekiobj)[2, :], seriestype = :scatter, label = "Final ensemble")
        plot!(
            [u_star[1]],
            xaxis = "u1",
            yaxis = "u2",
            seriestype = "vline",
            linestyle = :dash,
            linecolor = :red,
            label = :none,
        )
        plot!([u_star[2]], seriestype = "hline", linestyle = :dash, linecolor = :red, label = :none)
        savefig(p, "EKI_test.png")
    end

    posterior_cov_inv = (A' * (Γy \ A) + 1 * Matrix(I, n_par, n_par) / prior_cov)
    ols_mean = (A' * (Γy \ A)) \ (A' * (Γy \ y_obs))
    posterior_mean = posterior_cov_inv \ ((A' * (Γy \ A)) * ols_mean + (prior_cov \ prior_mean))

    #### This tests correspond to:
    # NOTE THESE ARE VERY SENSITIVE TO THE CORRESPONDING N_iter
    # EKI provides a solution closer to the ordinary Least Squares estimate
    @test norm(ols_mean - eki_final_result) < norm(ols_mean - eks_final_result)
    # EKS provides a solution closer to the posterior mean
    @test norm(posterior_mean - eks_final_result) < norm(posterior_mean - eki_final_result)

    ##### I expect this test to make sense:
    # In words: the ensemble covariance is still a bit ill-dispersed since the
    # algorithm employed still does not include the correction term for finite-sized
    # ensembles.
    @test abs(sum(diag(posterior_cov_inv \ cov(get_u_final(eksobj), dims = 2))) - n_par) > 1e-5

    ###
    ###  Calibrate: Unscented Kalman Inversion
    ###

    N_iter = 20 # number of UKI iterations
    α_reg = 1.0
    update_freq = 0
    process = Unscented(prior_mean, prior_cov; α_reg = α_reg, update_freq = update_freq)
    ukiobj = EnsembleKalmanProcessModule.EnsembleKalmanProcess(y_star, Γy, process)

    # UKI iterations
    params_i_vec = []
    g_ens_vec = []
    for i in 1:N_iter
        params_i = get_u_final(ukiobj)
        push!(params_i_vec, params_i)
        g_ens = G(params_i)
        push!(g_ens_vec, g_ens)
        if i == 1
            g_ens_t = permutedims(g_ens, (2, 1))
            @test_throws DimensionMismatch EnsembleKalmanProcessModule.update_ensemble!(ukiobj, g_ens_t)
        end
        EnsembleKalmanProcessModule.update_ensemble!(ukiobj, g_ens)
    end
    push!(params_i_vec, get_u_final(ukiobj))

    @test get_u_prior(ukiobj) == params_i_vec[1]
    @test get_u(ukiobj) == params_i_vec
    @test get_g(ukiobj) == g_ens_vec
    @test get_g_final(ukiobj) == g_ens_vec[end]
    @test get_error(ukiobj) == ukiobj.err

    # UKI results: Test if ensemble has collapsed toward the true parameter
    # values
    uki_final_result = get_u_mean_final(ukiobj)
    @test norm(u_star - uki_final_result) < 0.5

    if TEST_PLOT_OUTPUT
        gr()
        θ_mean_arr = hcat(ukiobj.process.u_mean...)
        N_θ = length(ukiobj.process.u_mean[1])
        θθ_std_arr = zeros(Float64, (N_θ, N_iter + 1))
        for i in 1:(N_iter + 1)
            for j in 1:N_θ
                θθ_std_arr[j, i] = sqrt(ukiobj.process.uu_cov[i][j, j])
            end
        end

        ites = Array(LinRange(1, N_iter + 1, N_iter + 1))
        p = plot(ites, grid = false, θ_mean_arr[1, :], yerror = 3.0 * θθ_std_arr[1, :], label = "u1")
        plot!(ites, fill(u_star[1], N_iter + 1), linestyle = :dash, linecolor = :grey, label = :none)
        plot!(ites, grid = false, θ_mean_arr[2, :], yerror = 3.0 * θθ_std_arr[2, :], label = "u2", xaxis = "Iterations")
        plot!(ites, fill(u_star[2], N_iter + 1), linestyle = :dash, linecolor = :grey, label = :none)
        savefig(p, "UKI_test.png")
    end

    ### Generate data from G₁(u) with a sparse u
    n_obs = 1                  # Number of synthetic observations from G₁(u)
    n_par = 2                  # Number of parameteres
    u_star = [1.0, 0.0]          # True parameters
    noise_level = 1e-3            # Defining the observation noise level
    Γy = noise_level * Matrix(I, n_obs, n_obs) # Independent noise for synthetic observations
    noise = MvNormal(zeros(n_obs), Γy)

    #### Define model
    function G₁(u)
        [sqrt((u[1] - u_star[1])^2 + (u[2] - u_star[2])^2)]
    end

    y_star = G₁(u_star)
    y_obs = y_star + rand(noise)

    @test size(y_star) == (n_obs,)

    #### Define prior information on parameters
    prior_distributions = [Parameterized(Normal(0, 2)), Parameterized(Normal(0, 2))]
    constraints = [[no_constraint()], [no_constraint()]]
    prior_names = ["u1", "u2"]
    prior = ParameterDistribution(prior_distns, constraints, prior_names)

    ###
    ###  Calibrate (4): Sparse Ensemble Kalman Inversion
    ###

    N_ens = 20 # number of ensemble members
    N_iter = 5 # number of EKI iterations
    initial_ensemble = EnsembleKalmanProcessModule.construct_initial_ensemble(prior, N_ens; rng_seed = rng_seed)
    @test size(initial_ensemble) == (n_par, N_ens)

    for (threshold_eki, threshold_value, test_name) in ((false, 1e-2, "test"), (true, 1e-2, "test_thresholded"))

        # Sparse EKI parameters
        γ = 1.0
        reg = 1e-4
        uc_idx = [1, 2]

        process = SparseInversion(γ, threshold_eki, threshold_value, reg, uc_idx)

        ekiobj = EnsembleKalmanProcessModule.EnsembleKalmanProcess(initial_ensemble, y_obs, Γy, process)

        # EKI iterations
        params_i_vec = []
        g_ens_vec = []
        for i in 1:N_iter
            params_i = get_u_final(ekiobj)
            push!(params_i_vec, params_i)
            g_ens = hcat([G₁(params_i[:, i]) for i in 1:N_ens]...)
            push!(g_ens_vec, g_ens)
            if i == 1
                g_ens_t = permutedims(g_ens, (2, 1))
                @test_throws DimensionMismatch EnsembleKalmanProcessModule.update_ensemble!(ekiobj, g_ens_t)
                EnsembleKalmanProcessModule.update_ensemble!(ekiobj, g_ens)
            else
                EnsembleKalmanProcessModule.update_ensemble!(ekiobj, g_ens, Δt_new = ekiobj.Δt[1])
            end
        end
        push!(params_i_vec, get_u_final(ekiobj))

        @test get_u_prior(ekiobj) == params_i_vec[1]
        @test get_u(ekiobj) == params_i_vec
        @test get_g(ekiobj) == g_ens_vec
        @test get_g_final(ekiobj) == g_ens_vec[end]
        @test get_error(ekiobj) == ekiobj.err

        # EKI results: Test if ensemble has collapsed toward the true parameter
        # values
        eki_final_result = vec(mean(get_u_final(ekiobj), dims = 2))
        @test norm(u_star - eki_final_result) < 0.1
        @test sum(eki_final_result .> 0.05) < size(eki_final_result)[1]

        # Plot evolution of the EKI particles
        eki_final_result = vec(mean(get_u_final(ekiobj), dims = 2))

        if TEST_PLOT_OUTPUT
            gr()
            p = plot(get_u_prior(ekiobj)[1, :], get_u_prior(ekiobj)[2, :], seriestype = :scatter)
            plot!(get_u_final(ekiobj)[1, :], get_u_final(ekiobj)[2, :], seriestype = :scatter)
            plot!([u_star[1]], xaxis = "u1", yaxis = "u2", seriestype = "vline", linestyle = :dash, linecolor = :red)
            plot!([u_star[2]], seriestype = "hline", linestyle = :dash, linecolor = :red)
            savefig(p, string("SparseEKI_", test_name, ".png"))
        end
    end
end


text:
function elastancefn(t::Float64,mparams::WK3.ModelParams,k::Float64)
    # two-Hill elastance
    tp = mod(t,sum(mparams.th));
    g1t = (tp/mparams.t1)^mparams.m1;
    g2t = (tp/mparams.t2)^mparams.m2;
    h1t = g1t/(1+g1t);
    h2t = 1/(1+g2t);
    E = k*h1t*h2t + mparams.Emin;

    # dE/dt
    g1p = mparams.m1/mparams.t1*(tp/mparams.t1)^(mparams.m1-1);
    g2p = mparams.m2/mparams.t2*(tp/mparams.t2)^(mparams.m2-1);
    h1p = g1p/((1+g1t)^2);
    h2p = -g2p/((1+g2t)^2);
    dE = k*(h1p*h2t + h2p*h1t);

    return E,dE
end

ext: jl


text:
# ps3000_set_ets.jl

export ps3000_set_ets, PS3000ETSMode

function ps3000_set_ets(handle::Int16, mode, ets_cycles, ets_interleave)
  effective_sample_time_ps = ccall((:ps3000_set_ets, ps3000driver), Int32,
    (Int16,Int16,Int16,Int16),handle,mode,ets_cycles,ets_interleave)
  if effective_sample_time_ps == 0
    error("ps3000 Error: ETS is disabled or one of the parameters is out of range")
  end
  return effective_sample_time_ps
end

baremodule PS3000ETSMode
  const OFF = 0 # disables ETS
  const FAST = 1 # enables ETS and provides ets_cycles cycles of data, which may contain data from previously returned cycles
  const SLOW = 2 # enables ETS and provides fresh data every ets_cycles cycles. PS3000_ETS_SLOW takes longer to provide each data set, but the data sets are more stable and unique
end  # module


text:
"""
    sequential_compare(f, vec)

Check that `f(vec[i], vec[i-1])` is true for each index of vec.
"""
function sequential_compare(f, vec)
    for i in Iterators.drop(eachindex(vec), 1)
        f(vec[i], vec[i - 1]) || return false
    end
    return true
end

"""
    is_increasing(vec; strict = false)
"""
function is_increasing(vec; strict=false)
    if strict
        return sequential_compare((x1, x0) -> x1 > x0, vec)
    else
        return sequential_compare((x1, x0) -> x1 >= x0, vec)
    end
end


text:
# MIT Mystery Hunt 2013 puzzle Wordplay
# http://www.mit.edu/~puzzle/2013/coinheist.com/get_smart/wordplay/index.html

@testset "wordplay" begin

    # Set 1
    @test best_feature(["ample", "adenoid", "music", "fifa"]).description == "is a hill word"

    # Set 2
    @test best_feature(["peeped", "isseis", "fee", "acacia", "salsas", "arrear"]).description == "is a pyramid word"

    # Set 3
    @test best_feature(["skort", "sporty", "yolks", "peccadillo", "unknot", "rosy"]).description == "is a valley word"

    # Set 4
    @test best_feature(["testset", "lol", "tenet", "malayalam"]).description == "is a palindrome"

    # Set 5
    @test best_feature(["hitchhiker", "kaashoek", "jellystone", "kierkegaard", "metallica", "maastrict", "menschheit"]).description == "contains 1 double letters"

    # Set 6
    @test best_feature(["aime", "eye", "eerie", "riaa", "oahu", "oeis"]).description == "has 1 unique consonants"
end


text:
using Documenter
using ConstraintSolver

makedocs(
    # See https://github.com/JuliaDocs/Documenter.jl/issues/868
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    strict = true,
    sitename = "ConstraintSolver",
    pages = [
        "Home" => "index.md",
        "Tutorial" => "tutorial.md",
        "How-To" => "how_to.md",
        "Explanation" => "explanation.md",
        "Reference" => "reference.md",
#        "Developer" => [],
#        "Library" => "library.md"
    ]
)

deploydocs(
    repo = "github.com/Wikunia/ConstraintSolver.jl.git",
)


text:
using VLModelParametersDB

# execute -
biophysical_dictionary = get_biophysical_parameters_dictionary_for_system(:HUMAN)

# get information -
key = :half_life_translation_capacity
parameter_value = value(biophysical_dictionary, key)
parameter_units = units(biophysical_dictionary, key)
parameter_reference = reference(biophysical_dictionary, key)


text:
# Chimera_states_v4.jl
# ====================   u.p. 13.9.21

using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot

function chimera_ODE(du,u,p,t)
    ω, R, alpha, nosc = p
    @inbounds for i = 1:nosc
        term = 0.0
        for j = i-R:i+R
            jj = 1 + mod(j-1,nosc)
            term = term + sin( u[i] - u[jj] + alpha )
        end
        du[i] = ω - term / (2*R )
    end
    return
end

#  ------------------------------  main program  --------------------------------
fig, axs = subplots(1,3; figsize=(figx,figy))

for axi in 1:2
if axi == 2
    # case a:  with chaotic transient
    nosc = 40 # 50#  40 # 100 # 40    # number of oscillators
    R = 14 # 15 # 14 # 35  # 14  # coupling range
    ish = 102  # select snapshot
    tsamp = 50.  # sampling time
    tspan = 25000.  # total lengthof time series
else
    # case b:  no chaotic transient (or transient end much later), but more oscillators
    nosc = 100 # 50#  40 # 100 # 40    # number of oscillators
    R = 35 # 15 # 14 # 35  # 14  # coupling range
    ish = 102  # select snapshot
    tsamp = 50.  # sampling time
    tspan = 25000.  # total lengthof time series
end

indx = LinRange(1,nosc,nosc)
ω = 0.
alpha = 1.46

p = ω, R, alpha, nosc

iseed = 7117 # 5117
rng = MersenneTwister(iseed) # random initial conditiosn in [0 ,2*pi]
u0 = 2 .* pi .* (rand(rng,nosc) .- 0.5)

prob = ODEProblem(chimera_ODE, u0, tspan, p)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8, saveat = tsamp)
uout = sol.u
tvec = sol.t
nts = length(tvec)

if axi == 1

axs[1].scatter(indx,mod.(uout[ish], 2*pi),color=COLORS[1])
axs[1].set_xlim(indx[1],indx[end])
axs[1].set_xticks(indx[1]:30:indx[end])
axs[1].set_ylim(0,2pi)
axs[1].set_xlabel(L"i"; labelpad = -20)
axs[1].set_ylabel(L"\phi"; labelpad=-10)
end

# plot evolution of phase velocity
# compute phase velocities (could perhaps be replaced by call of ODE function)
qmat = zeros(nosc,nts)
for its = 1:nts
    u = mod.(uout[its], 2*pi)
    for i = 1:nosc
        term = 0.0
        for j = i-R:i+R
            jj = 1 + mod(j-1,nosc)
            term = term + sin( u[i] - u[jj] + alpha )
        end
        qmat[i,its] = ω - term / (2*R )
    end
end

im1 = axs[1+axi].pcolormesh(tvec,indx,qmat,cmap = "inferno")
axs[1+axi].set_xlabel(L"t"; labelpad = -20)
axs[axi+1].set_ylabel(L"i"; labelpad = -20)
axs[axi+1].set_yticks([indx[1], indx[end]])
axs[axi+1].set_xticks([0, tvec[end]])

if axi == 2
    cbar = plt.colorbar(im1)  # , ticks=[-0.2, 0.2, 0.6, 1.], orientation="horizontal")
    # cbar.set_clim(1., 4.)
    cbar.set_label(L"d \phi /  dt")
end

end

# %%
add_identifiers!(fig, axs; xloc = 0.05)

fig.tight_layout(pad=0.3)

wsave(plotsdir("10", "chaimera"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
include("engine.jl")
include("parser.jl")

params = Dict()

params["file"] = ""
params["cells"] = 128
params["size"] = Int64
params["args"] = []
params["extended"] = false
params["build"] = false

i = 1

while i < length(ARGS)
    global i
    if ARGS[i] in ("--file", "-f")
        params["file"] = ARGS[i+=1]
    elseif ARGS[i] in ("--cells", "-c")
        try
            params["cells"] = Int128(ARGS[i+=1])
        catch cx
            exit()
        end
    elseif ARGS[i] in ("--size", "-s")
        try
            s = Int64(ARGS[i+=1])
            if s == 8
                params["size"] = Int8
            elseif s == 16
                params["size"] = Int16
            elseif s == 32
                params["size"] = Int32
            elseif s == 64
                params["size"] = Int64
            elseif s == 128
                params["size"] = Int128
            end
        catch cx
            exit()
        end
    elseif ARGS[i] in ("--args", "-a")
        while i < length(ARGS) && !(length(ARGS[i+1]) in ())
            err = false
            try
                a = parse(Int64, ARGS[i+1])
                if !err
                    push!(params["args"], a)
                end
            catch cx
                err = true
                push!(params["args"], Char(ARGS[i+1]))
            end
            i += 1
        end
    elseif ARGS[i] in ("--extended", "-e")
        params["extended"] = true
    elseif ARGS[i] in ("--build", "-b")
        params["build"] = true
    end
    i += 1
end

# TODO: Error Handling

source = open(params["file"]) do file
    read(file, String)
end

const ENGINE = create_engine(params["cells"], params["size"])

result = execute(ENGINE, source, params["args"])

for o = result
    print(typeof(o) == Char ? "$(o)" : "$(o) ")
end


text:
#=
MathOptInterface wrapper of Hypatia solver
=#

"""
$(TYPEDEF)

A MathOptInterface optimizer type for Hypatia.
"""
mutable struct Optimizer{T <: Real} <: MOI.AbstractOptimizer
    solver::Solvers.Solver{T} # Hypatia solver object

    # data for transforming certificates
    obj_sense::MOI.OptimizationSense
    zeros_idxs::Vector{UnitRange{Int}}
    moi_cones::Vector{MOI.AbstractVectorSet}
    moi_cone_idxs::Vector{UnitRange{Int}}

    function Optimizer{T}(; options...) where {T <: Real}
        opt = new{T}()
        opt.solver = Solvers.Solver{T}(; options...)
        return opt
    end
end

Optimizer(; options...) = Optimizer{Float64}(; options...) # default to Float64

MOI.is_empty(opt::Optimizer) = (opt.solver.status == Solvers.NotLoaded)

MOI.empty!(opt::Optimizer) = (opt.solver.status = Solvers.NotLoaded)

MOI.get(::Optimizer, ::MOI.SolverName) = "Hypatia"

MOI.get(opt::Optimizer, ::MOI.RawSolver) = opt.solver

MOI.supports(::Optimizer, ::MOI.Silent) = true

function MOI.set(opt::Optimizer, ::MOI.Silent, value::Bool)
    opt.solver.verbose = !value
    return
end

MOI.get(opt::Optimizer, ::MOI.Silent) = !opt.solver.verbose

MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true

function MOI.set(opt::Optimizer, ::MOI.TimeLimitSec, value::Union{Real, Nothing})
    opt.solver.time_limit = something(value, Inf)
    return
end

function MOI.get(opt::Optimizer, ::MOI.TimeLimitSec)
    if isfinite(opt.solver.time_limit)
        return opt.solver.time_limit
    end
    return
end

function MOI.get(opt::Optimizer, ::MOI.SolveTimeSec)
    if opt.solver.status == Solvers.NotLoaded
        error("solve has not been called")
    end
    return opt.solver.solve_time
end

MOI.get(opt::Optimizer, ::MOI.RawStatusString) = string(opt.solver.status)

MOI.get(opt::Optimizer, ::MOI.BarrierIterations) = opt.solver.num_iters

function MOI.set(opt::Optimizer, param::MOI.RawOptimizerAttribute, value)
    return setproperty!(opt.solver, Symbol(param.name), value)
end

function MOI.get(opt::Optimizer, param::MOI.RawOptimizerAttribute)
    return getproperty(opt.solver, Symbol(param.name))
end

MOI.supports(
    ::Optimizer{T},
    ::Union{MOI.ObjectiveSense, MOI.ObjectiveFunction{<:Union{VI, SAF{T}}}},
    ) where {T <: Real} = true

MOI.supports_constraint(
    ::Optimizer{T},
    ::Type{<:Union{VV, VAF{T}}},
    ::Type{<:Union{MOI.Zeros, SupportedCone{T}}},
    ) where {T <: Real} = true

function MOI.copy_to(
    opt::Optimizer{T},
    src::MOI.ModelLike,
    ) where {T <: Real}
    idx_map = MOI.Utilities.IndexMap()

    # variables
    n = MOI.get(src, MOI.NumberOfVariables()) # columns of A
    j = 0
    for vj in MOI.get(src, MOI.ListOfVariableIndices())
        j += 1
        idx_map[vj] = VI(j)
    end
    @assert j == n
    for attr in MOI.get(src, MOI.ListOfVariableAttributesSet())
        if attr == MOI.VariableName() || attr == MOI.VariablePrimalStart()
            continue
        end
        throw(MOI.UnsupportedAttribute(attr))
    end

    # objective
    opt.obj_sense = MOI.MIN_SENSE
    model_c = zeros(T, n)
    obj_offset = zero(T)
    for attr in MOI.get(src, MOI.ListOfModelAttributesSet())
        if attr == MOI.Name()
            continue
        elseif attr == MOI.ObjectiveSense()
            opt.obj_sense = MOI.get(src, MOI.ObjectiveSense())
        elseif attr isa MOI.ObjectiveFunction
            F = MOI.get(src, MOI.ObjectiveFunctionType())
            if !(F <: Union{VI, SAF{T}})
                error("objective function type $F not supported")
            end
            obj = convert(SAF{T}, MOI.get(src, MOI.ObjectiveFunction{F}()))
            for t in obj.terms
                model_c[idx_map[t.variable].value] += t.coefficient
            end
            obj_offset = obj.constant
        else
            throw(MOI.UnsupportedAttribute(attr))
        end
    end
    if opt.obj_sense == MOI.MAX_SENSE
        model_c .*= -1
        obj_offset *= -1
    end

    # constraints
    get_src_cons(F, S) = MOI.get(src, MOI.ListOfConstraintIndices{F, S}())
    get_con_fun(con_idx) = MOI.get(src, MOI.ConstraintFunction(), con_idx)
    get_con_set(con_idx) = MOI.get(src, MOI.ConstraintSet(), con_idx)

    # equality constraints
    (IA, JA, VA) = (Int[], Int[], T[])
    model_b = T[]
    opt.zeros_idxs = zeros_idxs = Vector{UnitRange{Int}}()
    for F in (VV, VAF{T}), ci in get_src_cons(F, MOI.Zeros)
        fi = get_con_fun(ci)
        si = get_con_set(ci)
        _con_IJV(IA, JA, VA, model_b, zeros_idxs, fi, si, idx_map)
        idx_map[ci] = MOI.ConstraintIndex{F, MOI.Zeros}(length(zeros_idxs))
    end
    model_A = dropzeros!(sparse(IA, JA, VA, length(model_b), n))

    # conic constraints
    (IG, JG, VG) = (Int[], Int[], T[])
    model_h = T[]
    moi_cones = MOI.AbstractVectorSet[]
    moi_cone_idxs = Vector{UnitRange{Int}}()
    cones = Cones.Cone{T}[]

    # build up one nonnegative cone
    for F in (VV, VAF{T}), ci in get_src_cons(F, MOI.Nonnegatives)
        fi = get_con_fun(ci)
        si = get_con_set(ci)
        _con_IJV(IG, JG, VG, model_h, moi_cone_idxs, fi, si, idx_map)
        push!(moi_cones, si)
        idx_map[ci] = MOI.ConstraintIndex{F, MOI.Nonnegatives}(length(moi_cones))
    end
    if !isempty(moi_cones)
        push!(cones, cone_from_moi(T, MOI.Nonnegatives(length(model_h))))
    end

    # other conic constraints
    for (F, S) in MOI.get(src, MOI.ListOfConstraintTypesPresent())
        if !MOI.supports_constraint(opt, F, S)
            throw(MOI.UnsupportedConstraint{F, S}())
        end
        for attr in MOI.get(src, MOI.ListOfConstraintAttributesSet{F, S}())
            if attr == MOI.ConstraintName() ||
               attr == MOI.ConstraintPrimalStart() ||
               attr == MOI.ConstraintDualStart()
                continue
            end
            throw(MOI.UnsupportedAttribute(attr))
        end
        if S == MOI.Zeros || S == MOI.Nonnegatives
            continue # already copied these constraints
        end

        for ci in get_src_cons(F, S)
            fi = get_con_fun(ci)
            si = get_con_set(ci)
            _con_IJV(IG, JG, VG, model_h, moi_cone_idxs, fi, si, idx_map)
            push!(cones, cone_from_moi(T, si))
            push!(moi_cones, si)
            idx_map[ci] = MOI.ConstraintIndex{F, S}(length(moi_cones))
        end
    end

    model_G = dropzeros!(sparse(IG, JG, VG, length(model_h), n))
    opt.moi_cone_idxs = moi_cone_idxs
    opt.moi_cones = moi_cones

    # finalize model and load into solver
    model = Models.Model{T}(model_c, model_A, model_b, model_G, model_h,
        cones; obj_offset = obj_offset)
    Solvers.load(opt.solver, model)

    return idx_map
end

MOI.optimize!(opt::Optimizer) = Solvers.solve(opt.solver)

function MOI.modify(
    opt::Optimizer{T},
    ::MOI.ObjectiveFunction{SAF{T}},
    chg::MOI.ScalarConstantChange{T},
    ) where {T}
    obj_offset = chg.new_constant
    if opt.obj_sense == MOI.MAX_SENSE
        obj_offset = -obj_offset
    end
    Solvers.modify_obj_offset(opt.solver, obj_offset)
    return
end

function MOI.modify(
    opt::Optimizer{T},
    ::MOI.ObjectiveFunction{SAF{T}},
    chg::MOI.ScalarCoefficientChange{T},
    ) where {T}
    new_c = chg.new_coefficient
    if opt.obj_sense == MOI.MAX_SENSE
        new_c = -new_c
    end
    Solvers.modify_c(opt.solver, [chg.variable.value], [new_c])
    return
end

function MOI.modify(
    opt::Optimizer{T},
    ci::MOI.ConstraintIndex{VAF{T}, MOI.Zeros},
    chg::MOI.VectorConstantChange{T},
    ) where {T}
    idxs = opt.zeros_idxs[ci.value]
    Solvers.modify_b(opt.solver, idxs, chg.new_constant)
    return
end

function MOI.modify(
    opt::Optimizer{T},
    ci::MOI.ConstraintIndex{VAF{T}, <:SupportedCone{T}},
    chg::MOI.VectorConstantChange{T},
    ) where {T}
    i = ci.value
    idxs = opt.moi_cone_idxs[i]
    set = opt.moi_cones[i]
    new_h = chg.new_constant
    if needs_permute(set)
        @assert !needs_rescale(set)
        new_h = permute_affine(set, new_h)
    end
    if needs_rescale(set)
        rescale_affine(set, new_h)
    end
    Solvers.modify_h(opt.solver, idxs, new_h)
    return
end

function MOI.get(opt::Optimizer, ::MOI.TerminationStatus)
    status = opt.solver.status
    if status in (Solvers.NotLoaded, Solvers.Loaded)
        return MOI.OPTIMIZE_NOT_CALLED
    elseif status == Solvers.Optimal
        return MOI.OPTIMAL
    elseif status == Solvers.NearOptimal
        return MOI.ALMOST_OPTIMAL
    elseif status == Solvers.PrimalInfeasible || status == Solvers.PrimalInconsistent
        return MOI.INFEASIBLE
    elseif status == Solvers.NearPrimalInfeasible
        return MOI.ALMOST_INFEASIBLE
    elseif status == Solvers.DualInfeasible || status == Solvers.DualInconsistent
        return MOI.DUAL_INFEASIBLE
    elseif status == Solvers.NearDualInfeasible
        return MOI.ALMOST_DUAL_INFEASIBLE
    elseif status == Solvers.SlowProgress
        return MOI.SLOW_PROGRESS
    elseif status == Solvers.IterationLimit
        return MOI.ITERATION_LIMIT
    elseif status == Solvers.TimeLimit
        return MOI.TIME_LIMIT
    elseif status == Solvers.NumericalFailure
        return MOI.NUMERICAL_ERROR
    elseif status in (Solvers.IllPosed, Solvers.NearIllPosed)
        return MOI.OTHER_LIMIT
    else
        @warn("Hypatia status $(opt.solver.status) not handled")
        return MOI.OTHER_ERROR
    end
end

function MOI.get(opt::Optimizer, attr::MOI.PrimalStatus)
    if attr.result_index != 1
        return MOI.NO_SOLUTION
    end
    status = opt.solver.status
    if status == Solvers.Optimal
        return MOI.FEASIBLE_POINT
    elseif status == Solvers.NearOptimal
        return MOI.NEARLY_FEASIBLE_POINT
    elseif status == Solvers.PrimalInfeasible
        return MOI.INFEASIBLE_POINT
    elseif status == Solvers.DualInfeasible
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif status == Solvers.NearDualInfeasible
        return MOI.NEARLY_INFEASIBILITY_CERTIFICATE
    elseif status in (Solvers.IllPosed, Solvers.NearIllPosed)
        return MOI.OTHER_RESULT_STATUS
    else
        return MOI.UNKNOWN_RESULT_STATUS
    end
end

function MOI.get(opt::Optimizer, attr::MOI.DualStatus)
    if attr.result_index != 1
        return MOI.NO_SOLUTION
    end
    status = opt.solver.status
    if status == Solvers.Optimal
        return MOI.FEASIBLE_POINT
    elseif status == Solvers.NearOptimal
        return MOI.NEARLY_FEASIBLE_POINT
    elseif status == Solvers.PrimalInfeasible
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif status == Solvers.NearPrimalInfeasible
        return MOI.NEARLY_INFEASIBILITY_CERTIFICATE
    elseif status == Solvers.DualInfeasible
        return MOI.INFEASIBLE_POINT
    elseif status in (Solvers.IllPosed, Solvers.NearIllPosed)
        return MOI.OTHER_RESULT_STATUS
    else
        return MOI.UNKNOWN_RESULT_STATUS
    end
end

_sense_val(sense::MOI.OptimizationSense) = (sense == MOI.MAX_SENSE ? -1 : 1)

function MOI.get(opt::Optimizer, attr::MOI.ObjectiveValue)
    MOI.check_result_index_bounds(opt, attr)
    return _sense_val(opt.obj_sense) * opt.solver.primal_obj
end

function MOI.get(opt::Optimizer, attr::MOI.DualObjectiveValue)
    MOI.check_result_index_bounds(opt, attr)
    return _sense_val(opt.obj_sense) * opt.solver.dual_obj
end

MOI.get(opt::Optimizer, ::MOI.ResultCount) = 1

function MOI.get(opt::Optimizer, attr::MOI.VariablePrimal, vi::VI)
    MOI.check_result_index_bounds(opt, attr)
    return opt.solver.result.x[vi.value]
end

function MOI.get(
    opt::Optimizer{T},
    attr::MOI.ConstraintDual,
    ci::MOI.ConstraintIndex{<:Union{VV, VAF{T}}, MOI.Zeros},
    ) where {T}
    MOI.check_result_index_bounds(opt, attr)
    return opt.solver.result.y[opt.zeros_idxs[ci.value]]
end

function MOI.get(
    opt::Optimizer{T},
    attr::MOI.ConstraintDual,
    ci::MOI.ConstraintIndex{<:Union{VV, VAF{T}}, <:SupportedCone{T}},
    ) where {T}
    MOI.check_result_index_bounds(opt, attr)
    i = ci.value
    z_i = opt.solver.result.z[opt.moi_cone_idxs[i]]
    return untransform_affine(opt.moi_cones[i], z_i)
end

function MOI.get(
    opt::Optimizer{T},
    attr::MOI.ConstraintPrimal,
    ci::MOI.ConstraintIndex{<:Union{VV, VAF{T}}, <:SupportedCone{T}},
    ) where {T}
    MOI.check_result_index_bounds(opt, attr)
    i = ci.value
    s_i = opt.solver.result.s[opt.moi_cone_idxs[i]]
    return untransform_affine(opt.moi_cones[i], s_i)
end

function _con_IJV(
    IM::Vector{Int},
    JM::Vector{Int},
    VM::Vector{T},
    vect::Vector{T},
    idxs_vect::Vector{UnitRange{Int}},
    func::VV,
    set::MOI.AbstractVectorSet,
    idx_map::MOI.IndexMap,
    ) where {T <: Real}
    dim = MOI.output_dimension(func)
    start = length(vect)
    idxs = start .+ (1:dim)
    push!(idxs_vect, idxs)
    append!(vect, zero(T) for _ in 1:dim)
    if needs_permute(set)
        perm_idxs = permute_affine(set, 1:dim)
        perm_idxs .+= start
        append!(IM, perm_idxs)
    else
        append!(IM, idxs)
    end
    append!(JM, idx_map[vi].value for vi in func.variables)
    append!(VM, -one(T) for _ in 1:dim)
    if needs_rescale(set)
        @views rescale_affine(set, VM[(end - dim + 1):end])
    end
    return
end

function _con_IJV(
    IM::Vector{Int},
    JM::Vector{Int},
    VM::Vector{T},
    vect::Vector{T},
    idxs_vect::Vector{UnitRange{Int}},
    func::VAF{T},
    set::MOI.AbstractVectorSet,
    idx_map::MOI.IndexMap,
    ) where {T <: Real}
    dim = MOI.output_dimension(func)
    start = length(vect)
    idxs = start .+ (1:dim)
    push!(idxs_vect, idxs)
    if needs_permute(set)
        @assert !needs_rescale(set)
        append!(vect, permute_affine(set, func.constants))
        perm_idxs = permute_affine(set, func)
        perm_idxs .+= start
        append!(IM, perm_idxs)
    else
        append!(vect, func.constants)
        append!(IM, start + vt.output_index for vt in func.terms)
    end
    append!(JM, idx_map[vt.scalar_term.variable].value for vt in func.terms)
    append!(VM, -vt.scalar_term.coefficient for vt in func.terms)
    if needs_rescale(set)
        @views vm = VM[(end - length(func.terms) + 1):end]
        rescale_affine(set, func, vm)
        @views rescale_affine(set, vect[idxs])
    end
    return
end


text:
function deferred_acceptance(m_prefs, f_prefs)

    m_num = length(m_prefs[1, :])#4
    f_num = length(f_prefs[1, :])#3 男４女３のケース
    singlemale = [i for i in 1:m_num] #[1,2,3,4]番がsingle
    m_unmatched = 0#自分とマッチの状態と定義
    f_unmatched = 0

    m_matched = [0 for i in 1:m_num]#[0,0,0,0]
    f_matched = [0 for i in 1:f_num]#[0,0,0] マッチしているのは自分

    target = [1 for i in 1:m_num]#[1,1,1,1] 各人が自分の第何希望の相手にプロポーズするか　はじめはみんな第一希望にプロポーズ

    while length(singlemale)>0#singleの男性がいなくなれば終了
        x = shift!(singlemale)#singlemaleの初項 singleの要素は削除されていく　男の番号
        y = m_prefs[:, x][target[x]]#男[x]の次に告白する相手＝女性y

        if y == 0#f_matched[y]を使っているのでyが0だとindexが取り出せない　告白相手が自分
            m_matched[x]=0

        else#告白相手が自分以外
            if f_matched[y] == 0#女性yの相手が自分自身の時
                if findfirst(f_prefs[:, y],f_matched[y])<findfirst(f_prefs[:, y],x)#現在のマッチ＞告白してきたx
                    push!(singlemale, x)#xをsinglemaleに戻す
                    target[x]+=1#次に好きな相手にtargetを変更
                else##今の相手<新たな男  ＝新たな男xとマッチするほうがいい f_matched[y] == 0なので捨てられる男はいない
                    f_matched[y]=x
                    m_matched[x]=y#新しいマッチに変更
                end

            else#女性に相手誰かしらのマッチ相手がいる場合
                if findfirst(f_prefs[:, y],f_matched[y])<findfirst(f_prefs[:, y],x)#現在のマッチ＞告白してきたx
                    push!(singlemale, x)#xをsinglemaleに戻す
                    target[x]+=1#次に好きな相手にtargetを変更
                else#今の相手<新たな男  ＝新たな男xとマッチするほうがいい
                    push!(singlemale, f_matched[y])#捨てられる男をsingleに戻す
                    m_matched[f_matched[y]]=0#捨てられた男は初期状態の自分とマッチに戻す
                    f_matched[y]=x
                    m_matched[x]=y#新しいマッチに変更
                end
            end
        end
    end
    return m_matched, f_matched
end


text:
function dagify(doc::Document)
    g = Dict{Symbol,Set{Symbol}}()

    node_ids = Dict{Symbol,Int}()
    sorted_definitions = Any[]

    # list of user defined types
    for (i, d) in enumerate(doc)
        !(d isa TypeDefinition) && continue

        typ = d.type
        if typ isa ObjectTypeDefinition ||
           typ isa InputObjectTypeDefinition ||
           typ isa UnionTypeDefinition
            n = node_id(d)
            g[n] = Set{Symbol}()
            node_ids[n] = i

            edges = collect_edges(d)
            for e in edges
                (e in BUILTIN_GRAPHQL_TYPES) && continue
                push!(g[n], e)
            end
        elseif typ isa ScalarTypeDefinition || typ isa EnumTypeDefinition
            # toplevel
            push!(sorted_definitions, d)
        end
    end

    sorted_nodes, revisited_graph = topological_sort(g)
    for n in sorted_nodes
        id = get(node_ids, n, nothing)
        !isnothing(id) && push!(sorted_definitions, doc[id])
    end

    return Document(sorted_definitions), revisited_graph
end

function topological_sort(graph::Dict{Symbol,Set{Symbol}})
    nodes = collect(keys(graph))
    visited_nodes = Set{Symbol}()
    temporary_nodes = Set{Symbol}()
    sorted_nodes = Symbol[]
    revisited_graph = Dict{Symbol,Set{Symbol}}()

    while !isempty(nodes)
        n = pop!(nodes)
        visit!(
            (nothing, n),
            graph,
            sorted_nodes,
            visited_nodes,
            temporary_nodes,
            revisited_graph,
        )
    end

    return sorted_nodes, revisited_graph
end

function visit!(
    edge::Tuple{Union{Symbol,Nothing},Symbol},
    graph::Dict{Symbol,Set{Symbol}},
    sorted_nodes::Vector{Symbol},
    visited_nodes::Set{Symbol},
    temporary_nodes::Set{Symbol},
    revisited_graph::Dict{Symbol,Set{Symbol}},
)
    from, node = edge[1], edge[2]

    node in visited_nodes && return

    if node in temporary_nodes
        node == from && return

        if haskey(revisited_graph, from)
            push!(revisited_graph[from], node)
        else
            revisited_graph[from] = Set([node])
        end
        return
    end

    push!(temporary_nodes, node)

    for n in get(graph, node, Symbol[])
        visit!(
            (node, n),
            graph,
            sorted_nodes,
            visited_nodes,
            temporary_nodes,
            revisited_graph,
        )
    end

    delete!(temporary_nodes, node)
    push!(visited_nodes, node)
    push!(sorted_nodes, node)

    return
end

collect_edges(t) = Symbol[]
collect_edges(t::TypeDefinition) = collect_edges(t.type)

function collect_edges(t::ObjectTypeDefinition)
    map(t.fields_definition) do fd
        node_id(fd)
    end
end

function collect_edges(t::InputObjectTypeDefinition)
    map(t.fields) do fd
        node_id(fd)
    end
end

function collect_edges(t::UnionTypeDefinition)
    map(t.types) do fd
        node_id(fd)
    end
end

node_id(x::TypeDefinition) = node_id(x.type)
node_id(x::ObjectTypeDefinition) = node_id(x.name)
node_id(x::InputObjectTypeDefinition) = node_id(x.name)
node_id(x::UnionTypeDefinition) = node_id(x.name)
node_id(x::ScalarTypeDefinition) = node_id(x.name)
node_id(x::EnumTypeDefinition) = node_id(x.name)
node_id(x::InputValueDefinition) = node_id(x.type)
node_id(x::FieldDefinition) = node_id(x.type)
node_id(x::NamedType) = node_id(x.type)
node_id(x::ListType) = node_id(x.type)
node_id(x::RBNF.Token) = Symbol(x.str)


text:
@testset "ComponentArrays" begin
    x_seed_1 = [0.0]
    x_seed_2 = [0.0]
    x_seed = ComponentArray(x_seed_1=x_seed_1, x_seed_2=x_seed_2)
    g_seed_1 = [0.0]
    g_seed_2 = [0.0]
    g_seed = ComponentArray(g_seed_1=g_seed_1, g_seed_2=g_seed_2)
    f_x_seed = 8157.682077608529

    nd = NonDifferentiable(exponential, x_seed)
    @test nd.f == exponential
    @test value(nd) == 0.0
    @test nd.f_calls == [0]
    od = OnceDifferentiable(exponential, exponential_gradient!, nothing, x_seed, 0.0, g_seed)
    @test od.f == exponential
    @test od.df == exponential_gradient!
    @test value(od) == 0.0
    @test od.f_calls == [0]
    @test od.df_calls == [0]
    @test typeof(od.DF) <: ComponentArray
    @test typeof(od.x_f) <: ComponentArray
    @test typeof(od.x_df) <: ComponentArray
end
@testset "Matrix OnceDifferentiable" begin
    x_seed = fill(0.0, 1, 2)
    g_seed = fill(0.0, 1, 2)
    f_x_seed = 8157.682077608529

    nd = NonDifferentiable(exponential, x_seed)
    @test nd.f == exponential
    @test value(nd) == 0.0
    @test nd.f_calls == [0]
    od = OnceDifferentiable(exponential, exponential_gradient!, nothing, x_seed, 0.0, g_seed)
    @test od.f == exponential
    @test od.df == exponential_gradient!
    @test value(od) == 0.0
    @test od.f_calls == [0]
    @test od.df_calls == [0]
    @test typeof(od.DF) <: Matrix
    @test typeof(od.x_f) <: Matrix
    @test typeof(od.x_df) <: Matrix
end
@testset "RecursiveArrays" begin
    x_seed_1 = [0.0]
    x_seed_2 = [0.0]
    x_seed = ArrayPartition(x_seed_1, x_seed_2)
    g_seed_1 = [0.0]
    g_seed_2 = [0.0]
    g_seed = ArrayPartition(g_seed_1, g_seed_2)
    f_x_seed = 8157.682077608529

    nd = NonDifferentiable(exponential, x_seed)
    @test nd.f == exponential
    @test value(nd) == 0.0
    @test nd.f_calls == [0]
    od = OnceDifferentiable(exponential, exponential_gradient!, nothing, x_seed, 0.0, g_seed)
    @test od.f == exponential
    @test od.df == exponential_gradient!
    @test value(od) == 0.0
    @test od.f_calls == [0]
    @test od.df_calls == [0]
    @test typeof(od.DF) <: ArrayPartition
    @test typeof(od.x_f) <: ArrayPartition
    @test typeof(od.x_df) <: ArrayPartition
end


text:
#!/usr/bin/env julia

include("read.jl")

# Function to call variation based on given alt reads and fraction cut-off

function recall(paths = ARGS[1:end], min_alt_read::Int=2, min_alt_frac=0.05)
	for file in paths
		vcf = read(file)

		for line in vcf
			if ismatch(r"^#|CHROM", line)
				print(line)
			else
				cols = split(line, "\t")
				alt_reads  = []
				tot_reads = []
				for gt in cols[5:end]
					push!(alt_reads, parse(Int, split(gt, ':')[1]))
					push!(tot_reads, parse(Int, split(gt, ':')[2]))
				end
				gt = maximum(alt_reads)
				af = maximum(alt_reads ./ tot_reads)
				if gt >= min_alt_read && af >= min_alt_frac
					print(line)
				end
			end
		end
	end
end

recall()


text:
using Test
import GDAL
import ArchGDAL; const AG = ArchGDAL

@testset "Tests for field defn" begin
    AG.createfielddefn("fieldname",GDAL.OFTInteger) do fd
        @test sprint(print, fd) == "fieldname (OFTInteger)"
        @test AG.getname(fd) == "fieldname"
        AG.setname!(fd, "newname")
        @test AG.getname(fd) == "newname"
        @test AG.gettype(fd) == GDAL.OFTInteger
        AG.settype!(fd, GDAL.OFTDate)
        @test AG.gettype(fd) == GDAL.OFTDate
        AG.settype!(fd, GDAL.OFTInteger)
        @test AG.getsubtype(fd) == GDAL.OFSTNone
        AG.setsubtype!(fd, GDAL.OFSTInt16)
        @test AG.getsubtype(fd) == GDAL.OFSTInt16
        AG.setsubtype!(fd, GDAL.OFSTBoolean)
        @test AG.getsubtype(fd) == GDAL.OFSTBoolean
        AG.setsubtype!(fd, GDAL.OFSTNone)
        @test AG.getjustify(fd) == GDAL.OJUndefined
        AG.setjustify!(fd, GDAL.OJLeft)
        @test AG.getjustify(fd) == GDAL.OJLeft
        @test AG.getwidth(fd) == 0
        AG.setwidth!(fd, 10)
        @test AG.getwidth(fd) == 10
        @test AG.getprecision(fd) == 0
        AG.setprecision!(fd, 20)
        @test AG.getprecision(fd) == 20
        AG.setparams!(fd, "finalname", GDAL.OFTDate, nwidth=5, nprecision=2,
                      justify=GDAL.OJRight)
        @test AG.gettype(fd) == GDAL.OFTDate
        @test AG.getname(fd) == "finalname"
        @test AG.getsubtype(fd) == GDAL.OFSTNone
        @test AG.getjustify(fd) == GDAL.OJRight
        @test AG.getwidth(fd) == 5
        @test AG.getprecision(fd) == 2
        @test AG.isignored(fd) == false
        AG.setignored!(fd, true)
        @test AG.isignored(fd) == true
        AG.setignored!(fd, false)
        @test AG.isignored(fd) == false

        @test AG.isnullable(fd) == true
        AG.setnullable!(fd, false)
        @test AG.isnullable(fd) == false
        AG.setnullable!(fd, true)
        @test AG.isnullable(fd) == true

        @test AG.getdefault(fd) == ""
        AG.setdefault!(fd, "0001/01/01 00:00:00")
        @test AG.getdefault(fd) == "0001/01/01 00:00:00"
        @test AG.isdefaultdriverspecific(fd) == true
    end
end

@testset "Tests for Geom Field Defn" begin
    AG.creategeomdefn("geomname", GDAL.wkbPolygon) do gfd
        @test AG.getname(gfd) == "geomname"
        AG.setname!(gfd, "my name!")
        @test AG.getname(gfd) == "my name!"

        @test AG.gettype(gfd) == GDAL.wkbPolygon
        AG.settype!(gfd, GDAL.wkbPolyhedralSurface)
        @test AG.gettype(gfd) == GDAL.wkbPolyhedralSurface

        @test sprint(print, AG.getspatialref(gfd)) == "NULL Spatial Reference System"
        AG.getspatialref(gfd) do spref
            @test sprint(print, spref) == "NULL Spatial Reference System"
        end
        AG.setspatialref!(gfd, AG.importEPSG(4326))
        @test sprint(print, AG.getspatialref(gfd)) == "Spatial Reference System: +proj=longlat +datum=WGS84 +no_defs"
        AG.getspatialref(gfd) do spref
            if !Sys.iswindows()
                @test sprint(print, spref) == "Spatial Reference System: +proj=longlat +datum=WGS84 +no_defs"
                # NOTE(yeesian): we get the following error on x86 in Appveyor
                # Expression: sprint(print, spref) == "Spatial Reference System: +proj=longlat +datum=WGS84 +no_defs "
                #   GDALError (CE_Failure, code 6):
                #     No translation for an empty SRS to PROJ.4 format is known.
            end
        end

        @test AG.isignored(gfd) == false
        AG.setignored!(gfd, true)
        @test AG.isignored(gfd) == true
        AG.setignored!(gfd, false)
        @test AG.isignored(gfd) == false

        @test AG.isnullable(gfd) == true
        AG.setnullable!(gfd, false)
        @test AG.isnullable(gfd) == false
        AG.setnullable!(gfd, true)
        @test AG.isnullable(gfd) == true
    end
end

@testset "Tests for Feature Defn" begin
    AG.createfeaturedefn("new_feature") do fd
        @test AG.nreference(fd) == 0
        AG.reference(fd)
        @test AG.nreference(fd) == 1
        AG.reference(fd)
        @test AG.nreference(fd) == 2
        AG.release(fd)
        @test AG.nreference(fd) == 1
        AG.dereference(fd)
        @test AG.nreference(fd) == 0
        AG.createfielddefn("fieldname",GDAL.OFTInteger) do fielddef
            @test AG.nfield(fd) == 0
            AG.addfielddefn!(fd, fielddef)
            @test AG.nfield(fd) == 1
            AG.addfielddefn!(fd, fielddef)
            @test AG.nfield(fd) == 2
            AG.addfielddefn!(fd, fielddef)
            @test AG.nfield(fd) == 3
            AG.createfielddefn("newfield",GDAL.OFTInteger) do fielddef2
                AG.addfielddefn!(fd, fielddef2)
                @test AG.nfield(fd) == 4
                AG.getname(AG.getfielddefn(fd,0)) == "fieldname"
                AG.getname(AG.getfielddefn(fd,1)) == "fieldname"
                AG.getname(AG.getfielddefn(fd,2)) == "fieldname"
                AG.getname(AG.getfielddefn(fd,3)) == "newfield"
            end
        end
        AG.deletefielddefn!(fd, 0)
        @test AG.nfield(fd) == 3
        AG.getname(AG.getfielddefn(fd,0)) == "fieldname"
        AG.getname(AG.getfielddefn(fd,1)) == "fieldname"
        AG.getname(AG.getfielddefn(fd,2)) == "newfield"

        AG.reorderfielddefns!(fd, Cint[2,1,0])
        @test AG.nfield(fd) == 3
        AG.getname(AG.getfielddefn(fd,0)) == "newfield"
        AG.getname(AG.getfielddefn(fd,1)) == "fieldname"
        AG.getname(AG.getfielddefn(fd,2)) == "fieldname"

        @test AG.ngeom(fd) == 1
        @test AG.getgeomtype(fd) == GDAL.wkbUnknown
        AG.setgeomtype!(fd, GDAL.wkbPolygon)
        @test AG.getgeomtype(fd) == GDAL.wkbPolygon
        @test AG.ngeom(fd) == 1

        @test AG.isgeomignored(fd) == false
        AG.setgeomignored!(fd, true)
        @test AG.isgeomignored(fd) == true
        AG.setgeomignored!(fd, false)
        @test AG.isgeomignored(fd) == false

        @test AG.isstyleignored(fd) == false
        AG.setstyleignored!(fd, true)
        @test AG.isstyleignored(fd) == true
        AG.setstyleignored!(fd, false)
        @test AG.isstyleignored(fd) == false

        @test AG.findgeomindex(fd) == 0
        gfd0 = AG.getgeomdefn(fd, 0)
        @test AG.ngeom(fd) == 1
        AG.addgeomdefn!(fd, gfd0)
        @test AG.ngeom(fd) == 2
        gfd1 = AG.getgeomdefn(fd, 1)
        AG.setname!(gfd0, "name0")
        AG.setname!(gfd1, "name1")
        @test AG.findgeomindex(fd, "") == -1
        @test AG.findgeomindex(fd, "name0") == 0
        @test AG.findgeomindex(fd, "name1") == 1
        AG.deletegeomdefn!(fd, 0)
        @test AG.ngeom(fd) == 1
        @test AG.findgeomindex(fd, "") == -1
        @test AG.findgeomindex(fd, "name0") == -1
        @test AG.findgeomindex(fd, "name1") == 0
        @test AG.nreference(fd) == 0
        AG.createfeature(fd) do f
            @test AG.nreference(fd) == 2 # artificially inflated
            @test AG.issame(AG.getfeaturedefn(f), fd) == true
        end
        @test AG.nreference(fd) == 0
    end
end


text:
@memoize function make_plan(size::Tuple)
    plan_fft(zeros(ComplexF64, size), flags=FFTW.MEASURE)
end

@memoize function make_plan!(size::Tuple)
    plan_fft!(zeros(ComplexF64, size), flags=FFTW.MEASURE)
end

function planned_fft(x)
    plan = make_plan(size(x))
    plan * x
end

function planned_ifft(x)
    plan = make_plan(size(x))
    plan \ x
end

function planned_fft!(x)
    plan = make_plan!(size(x))
    plan * x
end

function planned_ifft!(x)
    plan = make_plan!(size(x))
    plan \ x
end

function rrule(::typeof(planned_fft), x)
    N = prod(size(x))
    function pullback(∂x)
        NoTangent(), N * planned_ifft(∂x)
    end
    return planned_fft(x), pullback
end

function rrule(::typeof(planned_ifft), x)
    N = prod(size(x))
    function pullback(∂x)
        NoTangent(), planned_fft(∂x) / N
    end
    return planned_ifft(x), pullback
end


text:
module webcrawler

using HTTP
using Gumbo
using Cascadia
using .Cascadia: nodeText
import Cascadia: matchFirst
import Cascadia: nodeText

export fetchrandom, fetchpage, articleinfo

const PAGE_URL = "https://en.wikipedia.org/wiki/Julia_(programming_language)"
const LINKS = String[]

function fetchpage(url)
  url = startswith(url, "/") ? buildurl(url) : url

  response = HTTP.get(url)

  if response.status == 200 && length(response.body) > 0
    String(response.body)
  else
    ""
  end
end

function extractlinks(elem)
  if  isa(elem, HTMLElement) &&  tag(elem) == :a && in("href", collect(keys(attrs(elem))))
      url = getattr(elem, "href")
      startswith(url, "/wiki/") && push!(LINKS, url)
  end

  for child in children(elem)
    extractlinks(child)
  end
end

function extracttitle(elem)
  matchFirst(Selector("#section_0"), elem) |> nodeText
end

function extractimage(elem)
  e = matchFirst(Selector(".content a.image img"), elem)
  isa(e, Nothing) ? "" : e.attributes["src"]
end

function fetchrandom()
  fetchpage(RANDOM_PAGE_URL)
end

function articledom(content)
  if ! isempty(content)
    return Gumbo.parsehtml(content)
  end

  error("Article content can not be parsed into DOM")
end

function articleinfo(content)
  dom = articledom(content)

  Dict( :content => content,
        :links => extractlinks(dom.root),
        :title => extracttitle(dom.root),
        :image => extractimage(dom.root)
  )
end

function buildurl(article_url)
  PROTOCOL * DOMAIN_NAME * article_url
end

end


text:
# Copyright 2021 The FirstOrderLp Authors
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

# A script for calling SCS for benchmark comparisons with FirstOrderLp. Run
# with --help for a pretty description of the arguments. See the comments on
# solve_instance_and_output for a description of the output formats.

# Built-in Julia packages.
import SparseArrays

# Third-party Julia packages.
import ArgParse
import HiGHS
import JSON3
import JuMP
import MathOptInterface
import SCS

import FirstOrderLp

const MOI = MathOptInterface

function write_vector_to_file(filename, vector)
  open(filename, "w") do io
    for x in vector
      println(io, x)
    end
  end
end

"""
Creates a JuMP model based on a QuadraticProgrammingProblem. Returns the model,
the decision variables, and the references to the constraints.
"""
function instance_to_model(problem::FirstOrderLp.QuadraticProgrammingProblem)
  if !iszero(problem.objective_matrix)
    error("QP support is not yet implemented.")
  end

  model = JuMP.Model()
  num_variables = length(problem.variable_lower_bound)

  JuMP.@variable(model, x[i = 1:num_variables])
  for i in 1:num_variables
    if isfinite(problem.variable_lower_bound[i])
      JuMP.set_lower_bound(x[i], problem.variable_lower_bound[i])
    end
    if isfinite(problem.variable_upper_bound[i])
      JuMP.set_upper_bound(x[i], problem.variable_upper_bound[i])
    end
  end
  JuMP.@objective(model, Min, problem.objective_vector' * x)

  equality_range = FirstOrderLp.equality_range(problem)
  equalities = JuMP.@constraint(
    model,
    problem.constraint_matrix[equality_range, :] * x .==
    problem.right_hand_side[equality_range]
  )

  inequality_range = FirstOrderLp.inequality_range(problem)
  inequalities = JuMP.@constraint(
    model,
    problem.constraint_matrix[inequality_range, :] * x .>=
    problem.right_hand_side[inequality_range]
  )

  return model, x, [equalities; inequalities]
end

function extract_termination_reason(jump_model)
  termination_status = JuMP.termination_status(jump_model)
  if termination_status == MOI.OPTIMAL
    return FirstOrderLp.TERMINATION_REASON_OPTIMAL
  elseif termination_status == MOI.INFEASIBLE
    return FirstOrderLp.TERMINATION_REASON_PRIMAL_INFEASIBLE
  elseif termination_status == MOI.DUAL_INFEASIBLE
    return FirstOrderLp.TERMINATION_REASON_DUAL_INFEASIBLE
  elseif termination_status == MOI.TIME_LIMIT
    return FirstOrderLp.TERMINATION_REASON_TIME_LIMIT
  elseif termination_status == MOI.ITERATION_LIMIT
    return FirstOrderLp.TERMINATION_REASON_ITERATION_LIMIT
  elseif termination_status == MOI.NUMERICAL_ERROR
    return FirstOrderLp.TERMINATION_REASON_NUMERICAL_ERROR
  else
    return FirstOrderLp.TERMINATION_REASON_OTHER
  end
end

"""
Solves a linear programming problem using the provided solver. Takes a path to
an instance. The instance must have the extension .mps, or .mps.gz. Creates
`instance_summary.json` with a SolveLog serialized in JSON format,
`instance_primal.txt` with the primal solution, `instance_dual.txt` with the
dual solution, and `instance_stderr.txt`, `instance_stdout.txt`.
"""
function solve_instance_and_output(
  optimizer_with_attributes::MathOptInterface.OptimizerWithAttributes,
  output_dir::String,
  instance_path::String,
  print_stdout::Bool,
)
  if !isdir(output_dir)
    mkpath(output_dir)
  end

  if endswith(instance_path, ".mps") || endswith(instance_path, ".mps.gz")
    lp = FirstOrderLp.qps_reader_to_standard_form(instance_path)
  else
    error("Instance has unrecognized file extension: ", basename(instance_path))
  end

  presolve_info = FirstOrderLp.presolve(lp)

  jump_model, jump_variables, jump_constraints = instance_to_model(lp)

  JuMP.set_optimizer(jump_model, optimizer_with_attributes)

  instance_name =
    replace(replace(basename(instance_path), ".mps.gz" => ""), ".mps" => "")

  println("Instance: ", instance_name)

  stdout_path = joinpath(output_dir, instance_name * "_stdout.txt")
  stderr_path = joinpath(output_dir, instance_name * "_stderr.txt")

  open(stderr_path, "w") do stderr_io
    redirect_stderr(stderr_io) do
      # stdout is handled differently from stderr because we need to read the
      # solver output to populate the SolveLog.
      old_stdout = stdout
      stdout_read, stdout_write = redirect_stdout()

      JuMP.optimize!(jump_model)

      Base.Libc.flush_cstdio()
      redirect_stdout(old_stdout)
      close(stdout_write)

      solve_log = read(stdout_read, String)
      write(stdout_path, solve_log)
      if print_stdout
        println(solve_log)
      end

      running_time = JuMP.solve_time(jump_model)
      println("Recorded solve time: $running_time sec")

      log = FirstOrderLp.SolveLog()
      log.instance_name = instance_name
      log.command_line_invocation = join([PROGRAM_FILE; ARGS...], " ")
      log.termination_reason = extract_termination_reason(jump_model)
      log.termination_string = JuMP.raw_status(jump_model)
      log.solve_time_sec = running_time
      cumulative_kkt_matrix_passes = NaN
      if optimizer_with_attributes.optimizer_constructor == SCS.Optimizer
        log.iteration_count = Int32(MOI.get(jump_model, SCS.ADMMIterations()))
        # We read the log to extract "avg # CG iterations". This is printed only
        # in SCS's indirect mode.
        for line in split(solve_log, '\n')
          if occursin("avg # CG iterations", line)
            avg_cg_iters = parse(Float64, rstrip(split(line)[6], ','))
            # This formula was given by Brendan O'Donoghue.
            cumulative_kkt_matrix_passes =
              log.iteration_count + log.iteration_count * avg_cg_iters
            break
          end
        end
      end

      if JuMP.result_count(jump_model) > 0
        primal_sol = JuMP.value.(jump_variables)
        dual_sol = JuMP.dual.(jump_constraints)
        # eps_optimal_absolute and eps_optimal_relative are set to 1 because the
        # the result depends only on their ratio, which is fixed in our
        # experiments.
        last_iteration_stats = FirstOrderLp.compute_iteration_stats(
          lp,
          FirstOrderLp.cached_quadratic_program_info(lp),
          primal_sol,
          dual_sol,
          primal_sol,  # primal_ray_estimate
          dual_sol,  # dual_ray_estimate
          log.iteration_count,
          cumulative_kkt_matrix_passes,
          running_time,
          1.0,  # eps_optimal_absolute
          1.0,  # eps_optimal_relative
          0.0,  # step_size
          0.0,  # primal_weight
          FirstOrderLp.POINT_TYPE_CURRENT_ITERATE,
        )
        log.solution_stats = last_iteration_stats
        log.solution_type = FirstOrderLp.POINT_TYPE_CURRENT_ITERATE
        primal_sol, dual_sol =
          FirstOrderLp.undo_presolve(presolve_info, primal_sol, dual_sol)

        primal_output_path = joinpath(output_dir, instance_name * "_primal.txt")
        write_vector_to_file(primal_output_path, primal_sol)

        dual_output_path = joinpath(output_dir, instance_name * "_dual.txt")
        write_vector_to_file(dual_output_path, dual_sol)
      else
        log.solution_type = FirstOrderLp.POINT_TYPE_NONE
        log.solution_stats = FirstOrderLp.IterationStats()
        log.solution_stats.cumulative_time_sec = running_time
      end

      # Complete iteration stats are not available, so we write only the
      # summary.
      summary_output_path =
        joinpath(output_dir, instance_name * "_summary.json")
      open(summary_output_path, "w") do io
        write(io, JSON3.write(log, allow_inf = true))
      end
    end
  end

  return
end

"""
Defines parses and args.

# Returns
A dictionary with the values of the command-line arguments.
"""
function parse_command_line()
  arg_parse = ArgParse.ArgParseSettings()

  help_instance_path = "The path to the instance to solve in .mps.gz, or .mps format."

  help_print_stdout =
    "In addition to redirecting stdout, also print it at " *
    " the end for convenience."

  ArgParse.@add_arg_table! arg_parse begin
    "--solver"
    help =
      "The solver to use. May be 'scs-direct', 'scs-indirect', " *
      "'highs-simplex', or 'highs-ipm'."
    arg_type = String
    required = true

    "--output_dir"
    help = "The directory for output files."
    arg_type = String
    required = true

    "--instance_path"
    help = help_instance_path
    arg_type = String
    required = true

    "--verbose"
    help = "Enable verbose solver output."
    arg_type = Bool
    default = true

    "--print_stdout"
    help = help_print_stdout
    arg_type = Bool
    default = true

    "--tolerance"
    help =
      "For SCS, the value to set for 'eps', a convergence tolerance." *
      "For HiGHS, the primal and dual feasibility tolerance. NOTE:" *
      "Tolerances are interpreted differently by each solver!"
    arg_type = Float64
    required = true

    # SCS parameters are defined at
    # https://github.com/cvxgrp/scs/blob/e5bb794ac014b7a86d127ac03651d2c8a12ecba8/include/scs.h#L44
    # with default values at
    # https://github.com/cvxgrp/scs/blob/e5bb794ac014b7a86d127ac03651d2c8a12ecba8/include/glbopts.h#L30.
    "--scs-normalize"
    help = "For SCS only, apply SCS's internal rescaling heuristic."
    arg_type = Bool
    default = false

    "--scs-acceleration_lookback"
    help = "For SCS only, the number of iterations used for Anderson Acceleration"
    arg_type = Int64
    default = 10

    "--scs-cg_rate"
    help = "For scs-indirect only, the rate at which the CG convergence tolerance decreases as a function of the iteration number."
    arg_type = Float64
    default = 2.0

    "--scs-scale"
    help = "For SCS only, when normalize is set, the iterates are additionally scaled by this factor."
    arg_type = Float64
    default = 1.0

    "--iteration_limit"
    help = "Maximum number of iterations to run."
    arg_type = Int64
    default = typemax(Int64)

    # NOTE: This flag is required for compatibility with experiment scripts.
    "--redirect_stdio"
    help = "Redirect stdout and stderr to files (for batch runs). The false value is not supported."
    arg_type = Bool
    default = true

  end

  return ArgParse.parse_args(arg_parse)
end

function main()
  parsed_args = parse_command_line()

  if !parsed_args["redirect_stdio"]
    error("The current implementation doesn't support redirect_stdio=false.")
  end

  iteration_limit = parsed_args["iteration_limit"]

  if parsed_args["solver"] == "scs-indirect" ||
     parsed_args["solver"] == "scs-direct"
    parameters = []
    push!(
      parameters,
      ("acceleration_lookback", parsed_args["scs-acceleration_lookback"]),
    )
    # This is the over-relaxation parameter. 1.0 is the best value for LP
    # according to SCS's author.
    push!(parameters, ("alpha", 1.0))
    # Disable rescaling to remove confounding factors.
    if !parsed_args["scs-normalize"]
      push!(parameters, ("normalize", 0))
    end
    push!(parameters, ("scale", parsed_args["scs-scale"]))
    push!(parameters, ("eps", parsed_args["tolerance"]))
    if iteration_limit != typemax(Int64)
      push!(parameters, ("max_iters", iteration_limit))
    end
    if parsed_args["verbose"]
      push!(parameters, ("verbose", 1))
    else
      push!(parameters, ("verbose", 0))
    end
    if parsed_args["solver"] == "scs-indirect"
      push!(parameters, ("linear_solver", SCS.IndirectSolver))
      push!(parameters, ("cg_rate", parsed_args["scs-cg_rate"]))
    else
      push!(parameters, ("linear_solver", SCS.DirectSolver))
    end
    optimizer_with_attributes = MathOptInterface.OptimizerWithAttributes(
      SCS.Optimizer,
      [MOI.RawParameter(p) => v for (p, v) in parameters],
    )
  elseif parsed_args["solver"] == "highs-simplex" ||
         parsed_args["solver"] == "highs-ipm"
    # See https://www.maths.ed.ac.uk/hall/HiGHS/HighsOptions.set for the
    # possible options.
    parameters = []
    if parsed_args["solver"] == "highs-simplex"
      push!(parameters, ("solver", "simplex"))
      if iteration_limit != typemax(Int64)
        push!(parameters, ("simplex_iteration_limit", iteration_limit))
      end
    else
      push!(parameters, ("solver", "ipm"))
      if iteration_limit != typemax(Int64)
        push!(parameters, ("ipm_iteration_limit", iteration_limit))
      end
      # TODO: Consider disabling crossover using the "run_crossover" option.
      # This option appears to be incompatible with presolving. See:
      # https://github.com/ERGO-Code/HiGHS/blob/cfe064e60fdca9dc8008fd956fff04d7b562e210/src/lp_data/Highs.cpp#L507
    end
    push!(parameters, ("parallel", "off"))
    push!(
      parameters,
      ("primal_feasibility_tolerance", parsed_args["tolerance"]),
    )
    push!(parameters, ("dual_feasibility_tolerance", parsed_args["tolerance"]))
    optimizer_with_attributes = MathOptInterface.OptimizerWithAttributes(
      HiGHS.Optimizer,
      [MOI.RawParameter(p) => v for (p, v) in parameters],
    )
  else
    error("Unrecognized solver $(parsed_args["solver"]).")
  end

  solve_instance_and_output(
    optimizer_with_attributes,
    parsed_args["output_dir"],
    parsed_args["instance_path"],
    parsed_args["print_stdout"],
  )
end

main()


text:
using Test
using ConstraintSolver
using JSON
using MathOptInterface, JuMP

const MOI = MathOptInterface
const CS = ConstraintSolver

include("fct_tests.jl")
include("moi.jl")

include("sudoku_fcts.jl")
include("small_special_tests.jl")
include("sudoku_tests.jl")
include("killer_sudoku_tests.jl")
include("graph_color_tests.jl")


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "tectonic"
version = v"0.1.15"

# Collection of sources required to build tar
sources = [
    ArchiveSource(
        "https://github.com/tectonic-typesetting/tectonic/archive/tectonic@$(version).tar.gz",
        "0e55188eafc1b58f3660a303fcdd6adc071051b9eb728119837fbeed2309914f"
    )
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/tectonic-*/
cargo build --release -j${nproc}
cp target/${rust_target}/release/tectonic${exeext} ${bindir}/
"""

# Some platforms disabled for now due issues with rust and musl cross compilation. See #1673.
platforms = [
    Platform("x86_64", "freebsd"),
    Platform("aarch64", "linux"; libc="glibc"),
    # Platform("aarch64", "linux"; libc="musl"),
    Platform("armv7l", "linux"; libc="glibc"),
    # Platform("armv7l", "linux"; libc="musl"),
    Platform("i686", "linux"; libc="glibc"),
    # Platform("i686", "linux"; libc="musl"),
    Platform("powerpc64le", "linux"; libc="glibc"),
    Platform("x86_64", "linux"; libc="glibc"),
    # Platform("x86_64", "linux"; libc="musl"),
    Platform("x86_64", "macos"),
    # Platform("i686", "windows"),
    Platform("x86_64", "windows"),
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("tectonic", :tectonic),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Fontconfig_jll"),
    Dependency("FreeType2_jll"),
    Dependency("Graphite2_jll"),
    Dependency("HarfBuzz_jll"),
    Dependency(PackageSpec(; name="ICU_jll", version=v"67.1.0")),
    Dependency("OpenSSL_jll"),
    Dependency("Zlib_jll"),
    Dependency("libpng_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; compilers=[:c, :rust], preferred_gcc_version=v"7", lock_microarchitecture=false)

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil


text:
@isdefined(postproc) || (postproc = false)
@isdefined(plotresults) || (plotresults = false)
plotresults && (postproc = true)

if postproc

    Φ, Θ = [0.0], range(0,stop=π,length=100)
    pts = [point(cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)) for ϕ in Φ for θ in Θ]
    ffd = potential(MWFarField3D(wavenumber=κ), pts, u, X)

    fcr, geo = facecurrents(u, X)

    nx, nz = 50, 100
    xs, zs = range(-2,stop=2,length=nx), range(-4,stop=4,length=nz)
    gridpoints = [point(x,0,z) for x in xs, z in zs]
    nfd = potential(MWSingleLayerField3D(wavenumber = κ), gridpoints, u, X)
    nfd = reshape(nfd, (nx,nz))
    nfd .-= E.(gridpoints)
end


text:
module LTI

include("_1_lti_primitives.jl")
include("_2_lti_systems.jl")
include("_3_lti_pprint.jl")
include("_4_lti_kernel.jl")
include("_5_lti_freq_resp.jl")
include("_6_lti_utils.jl")

end # module


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

@recipe function f(data::Data, vars=nothing)
  # valid variables
  validvars = Tables.schema(values(data)).names

  # plot all variables by default
  isnothing(vars) && (vars = validvars)
  @assert vars ⊆ validvars "invalid variable names"

  # underlying domain where variables live
  dom = domain(data)

  # shared plot specs
  layout --> length(vars)

  for (i, var) in enumerate(vars)
    vals = data[var]

    # handle missing values
    vals = replace(vals, missing => NaN)

    @series begin
      subplot := i
      title --> string(var)
      legend --> false
      dom, vals
    end
  end
end


text:
import StatsBase.mean

"""
    standard_te(driver, response;
        estimator = :tetogrid,
        E = nothing, v = nothing,
        ϵ = nothing, n_ϵ = 5, η = 1, τ = 1, dim = 3,
        k1 = 2, k2 = 3, distance_metric = Chebyshev(),
        which_is_surr = :none,
        surr_func = aaft,
        min_numbins = nothing,
        max_numbins = nothing)

Compute transfer entropy from a `driver` time series to a `response` time
series.

## Method
`method` sets the transfer entropy estimator to use. There are two types of
estimators: grid-based approaches, and nearest-neighbor based approaches.

## Embedding instructions
Specify embedding and instructions for how to compute marginal
entropies are with `E::AbstractEmbedding` and `v::TEVars`. If either
is missing, the data is embedded using the provided forward prediction lag `η`,
with embedding dimension `dim` and embedding lag `η`.

## Grid-based estimators

- `:transferentropy_transferoperator_grid`, or `:tetogrid`. Grid-based transfer entropy
    estimator based on the transfer operator.
- `::transferentropy_visitfreq`, or `:tefreq`. Simple visitation frequency based estimator.

For the grid based TE estimators, `ϵ` sets the bin size. If `ϵ == nothing`, then
the algorithm uses a bin size corresponding to a number of subdivisions `N`
along each axis so that `N ≦ npoints^(1/(dim+1))`. Transfer entropy is then
computed as an average over `n_ϵ` different bin sizes corresponding to
`N-1` to `N` subdivisions along each axis. If `ϵ` is specified and consists
of multiple bin sizes, an average
TE estimate over the bin sizes is returened. If `ϵ` is a single value,
then TE is estimated for that bin size only.

## Nearest neighbor based estimator.
- `:transferentropy_kraskov`, `:tekraskov`, or `:tekNN`. A nearest-neighbor
    based estimator that computes transfer entropy as the sum of two mutual
    information terms. `k1` and `k2` sets the
    number of nearest neighbors used to estimate the mutual information terms.
    `distance_metric` sets the distance metric (must be an instance
    of a valid metric from `Distances.jl`).

## Surrogate analysis
A surrogate analysis will be run if `which_is_surr` is set to
    `:both`, `:driver`, `:response`. Default is `:none`.
- `which_is_surr = :both` will replace both time series with a surrogate.
- `which_is_surr = :driver` will replace the driver time series with a surrogate.
- `which_is_surr = :response` will replace the response time series with a surrogate.

The type of surrogate must be either `randomshuffle`, `randomphases`,
`randomamplitudes`, `aaft` or `iaaft`.
"""
function standard_te(driver::AbstractVector{T}, response::AbstractVector{T};
        dim = 3,
        τ = 1,
        η = 1,
        estimator = :tetogrid,
        E = nothing, v = nothing,
        ϵ = nothing, n_ϵ = 5,
        k1 = 2, k2 = 3, distance_metric = Chebyshev(),
        which_is_surr = :none,
        surr_func = aaft,
        nbin_range = 2,
        min_numbins = nothing,
        max_numbins = nothing) where T

    estims_te_grid = [:transferentropy_transferoperator_grid, :tetogrid]
    estims_te_kNN = [:transferentropy_kraskov, :tekraskov, :tekNN]
    estims_te_freq = [:transferentropy_visitfreq, :tefreq]
    valid_estimators = estims_te_grid ∪ estims_te_kNN ∪ estims_te_freq

    if dim < 3
        error("Dimension $dim is too low. Must be at least 3.")
    end

    if !(estimator ∈ valid_estimators)
        error("Transfer entropy method $estimator not valid.")
    end

    valid_surr_funcs = [TimeseriesSurrogates.randomshuffle,
                        TimeseriesSurrogates.randomphases,
                        TimeseriesSurrogates.randomamplitudes,
                        TimeseriesSurrogates.aaft,
                        TimeseriesSurrogates.iaaft]

    if !(surr_func ∈ valid_surr_funcs)
        error("Surrogate type $surr_func is not valid.")
    end
    if which_is_surr == :driver
        driver[:] = surr_func(driver)
    elseif which_is_surr == :response
        response[:] = surr_func(response)
    elseif which_is_surr == :both
        driver[:] = surr_func(driver)
        response[:] = surr_func(response)
    elseif !(which_is_surr == :none)
        error("$which_is_surr is not a valid surrogate choice." )
    end
    D = Dataset(driver, response)
    # Compute transfer entropy analogously to how Krakovskà et al. (2018)
    # computes conditional mutual information.
    if !(typeof(E) == Embeddings.AbstractEmbedding) || !(typeof(v) == TEVars)
        if dim == 3
            E = customembed(D, Positions(2, 2, 1), Lags(η, 0, 0))
            v = TEVars(Tf = [1], Tpp = [2], Spp = [3])
        elseif dim == 4
            E = customembed(D, Positions(2, 2, 2, 1), Lags(η, 0, -τ, 0))
            v = TEVars(Tf = [1], Tpp = [2, 3], Spp = [4])
        elseif dim == 5
            E = customembed(D, Positions(2, 2, 2, 2, 1), Lags(η, 0, -τ, -2*τ, 0))
            v = TEVars(Tf = [1], Tpp = [2, 3, 4], Spp = [5])
        elseif dim == 6
            E = customembed(D, Positions(2, 2, 2, 2, 2, 1), Lags(η, 0, -τ, -2*τ, -3*τ, 0))
            v = TEVars(Tf = [1], Tpp = [2, 3, 4, 5], Spp = [6])
        elseif dim == 7
            E = customembed(D, Positions(2, 2, 2, 2, 2, 2, 1), Lags(η, 0, -τ, -2*τ, -3*τ, -4*τ, 0))
            v = TEVars(Tf = [1], Tpp = [2, 3, 4, 5, 6], Spp = [7])
        elseif dim == 8
            E = customembed(D, Positions(2, 2, 2, 2, 2, 2, 2, 1), Lags(η, 0, -τ, -2*τ, -3*τ, -4*τ, -5*τ, 0))
            v = TEVars(Tf = [1], Tpp = [2, 3, 4, 5, 6, 7], Spp = [8])
        end
    end

    if ϵ == nothing
        ϵ = zeros(Float64, n_ϵ, dim)
        # We a range of bin sizes along each axis.
        for d in 1:dim
            dmin = Base.minimum(E.reconstructed_pts[:, d])
            dmax = Base.maximum(E.reconstructed_pts[:, d])
            maxrange = dmax - dmin

            # If minimum and maximum number of bins are not provided,
            # use the Palus recommendation of N^(1/dim) as a max.
            if !(typeof(min_numbins) == Int && typeof(max_numbins) == Int)
                max_numbins = ceil(Int, length(E.reconstructed_pts)^(1/(dim+1)))
                min_numbins = max_numbins - nbin_range
            end

            max_bin_size = maxrange/max_numbins
            min_bin_size = maxrange/min_numbins

            ϵ[:, d] = LinRange(log(10, min_bin_size), log(10, max_bin_size), n_ϵ)
        end
        binsizes = [ϵ[i, :] for i = 1:n_ϵ]
        ϵs = binsizes
    else
        if length(ϵ) == 1
            ϵs = [ϵ]
        else
            ϵs = ϵ
        end
    end

    tes = zeros(Float64, length(ϵs))

    # Loop over the bin sizes
    for i = 1:length(ϵs)
        if estimator ∈ (:transferentropy_transferoperator_grid,
                        :tetogrid)
            #tes[i] = transferentropy_transferoperator_grid(E, ϵs[i], v)
            tes[i] = transferentropy(E.reconstructed_pts, v, RectangularBinning(ϵs[i]), TransferOperatorGrid())

        elseif estimator ∈ (:transferentropy_kraskov,
                            :tekraskov,
                            :tekNN)
            tes[i] = transferentropy_kraskov(E, k1, k2, v, metric = distance_metric)

        elseif estimator ∈ (:transferentropy_visitfreq,
                            :tefreq)
            #tes[i] = transferentropy_visitfreq(E, ϵs[i], v)
            tes[i] = transferentropy(E.reconstructed_pts, v, RectangularBinning(ϵs[i]), VisitationFreq())

        elseif estimator == :transferoperator_triang
                warn("high level wrapper for $estimator is not implemented yet.")
        end

    end

    # If we only use one binsize, we can't take the integral. This is always
    # the case for the kNN based estimators, because they don't use a bin size.
    if estimator ∈ estims_te_kNN
        return tes[1]
    else
        return mean(tes)
    end
end

"""
    te_reg(source::AbstractArray{<:Real, 1},
        response::AbstractArray{<:Real, 1},
        k::Int, l::Int, m::Int;
        η = 1, τ = 1,
        estimator = VisitationFrequency(),
        n_subdivs = 1,
        b = 2)

## TE estimation with default discretization scheme(s)

Calculate transfer entropy from `source` to `response` using the provided
`estimator` on a rectangular discretization of a `k + l + m`-dimensional
delay embedding of the input data, using an embedding delay of `τ` across
all embedding components. `η` is the prediction lag.

## Arguments
- **`source`**: The source data series.
- **`target`**: The target data series.
- **`k`**: The dimension of the ``T_{f}`` component of the embedding.
- **`l`**: The dimension of the ``T_{pp}`` component of the embedding.
- **`m`**: The dimension of the ``S_{pp}`` component of the embedding.

## Keyword arguments

- **`τ`**: The embedding lag. Default is `τ = 1`.
- **`η`**: The prediction lag. Default is `η = 1`.
- **`estimator`**: The transfer entropy estimator to use. The default
    is `VisitationFrequency()`.
- **`n_subdivs`**: The number of different partitions of varying coarseness
    to compute TE over, as described below. Default is `n_subdivs = 2`.
    (this way, TE is computed over two separate partitions). Unless `n_subdivs = 0`,
    make sure that `n_subdivs` is the same across analyses if they
    are to be compared. T TE
- **`b`**: Base of the logarithm. The default (`b = 2`) gives the TE in bits.
## More about the embedding

To compute transfer entropy, we need an appropriate delay embedding
of `source` (``S``) and `target` (``T``). For convenience, define

```math
\\begin{align}
T_f^{(k)} &= \\{(T(t+\\eta_k), \\ldots, T(t+\\eta_2), T(t+\\eta_1)) \\} \\\\
T_{pp}^{(l)} &= \\{ (T(t), T(t-\\tau_1), T(t-\\tau_2), \\ldots, T(t - \\tau_{l - 1})) \\} \\\\
S_{pp}^{(m)} &= \\{ (S(t), S(t-\\tau_1), S(t-\\tau_2), \\ldots, S(t-\\tau_{m - 1})) \\}
\\end{align}
```

where ``T_f`` denotes the `k`-dimensional set of vectors furnishing the future states of ``T``,
``T_{pp}`` denotes the `l`-dimensional set of vectors furnishing the past and present states of ``T``,
and ``S_{pp}`` denotes the `m`-dimensional set of vectors furnishing the past and present of ``S``.
``\\eta`` is the prediction lag. This convenience function uses ``\\tau_1`` = `τ`,
``\\tau_2`` = `2*τ`, ``\\tau_3`` = `3*τ`, and so on.

Combined, we get the generalised embedding ``\\mathbb{E} = (T_f^{(k)}, T_{pp}^{(l)}, S_{pp}^{(m)})``,
which is discretized as described below.

## More about discretization

To compute TE, we coarse-grain the `k+l+m`-dimensional generalised embedding
``\\mathbb{E} = (T_f^{(k)}, T_{pp}^{(l)}, S_{pp}^{(m)})`` into hyperrectangular boxes.
The magnitude of the TE may be biased by the particular choice of binning scheme,
so we compute TE across a number of different box sizes, determined as follows.

Let ``L`` be the number of observations in ``S`` (and ``T``). The coarsest box size is
determined by subdiving the i-th coordinate axis into
``N = ceiling(L^\\frac{1}{k + l + m + 1})`` intervals of equal lengths,
resulting in a box size of ``|max(dim_{i}) - min(dim_{i})|/N``. The next box size
is given by ``|max(dim_{i}) - min(dim_{i})|/(N+1)``, then
``|max(dim_{i}) - min(dim_{i})|/(N+2)``, and so on, until the finest box
size which is given by ``|max(dim_{i}) - min(dim_{i})|/(N+N_{subdivs})``.

## Transfer entropy computation

TE is then computed as

```math
\\begin{align}
TE_{S \\rightarrow T} = \\int_{\\mathbb{E}} P(T_f, T_{pp}, S_{pp}) \\log_{b}{\\left(\\frac{P(T_f | T_{pp}, S_{pp})}{P(T_f | T_{pp})}\\right)}
\\end{align}
```
using the provided `estimator` (default = `VisitationFrequency`) for
each of the discretizations. A vector of the TE estimates for each discretization
is returned.
"""
function te_reg(source::AbstractArray{<:Real, 1},
        response::AbstractArray{<:Real, 1},
        k::Int, l::Int, m::Int;
        η = 1, τ = 1,
        estimator = VisitationFrequency(),
        n_subdivs = 2,
        b = 2)

    k + l + m >= 3 || throw(ArgumentError("`dim = k + l + m` must be 3 or higher for regular TE"))

    pts, vars = te_embed(source, response, k, l, m, η = η, τ = τ)

    # Determine appropriate binnings from time series length (roughly according to
    # Krakovska et al. (2018)'s recommendations)
    L = length(source)

    # one subdivision coarser than Krakovska as the coarest partition
    nbins_coarsest = ceil(Int, L^(1/(k + l + m + 1))) - 1

    # with the default n_subdivs = 2, we cover both the Krakovska recommendation,
    # as well as one finer bin size and one coarser bin size.
    nbins_finest = nbins_coarsest + n_subdivs
    binning_scheme = map(n-> RectangularBinning(n), nbins_coarsest:nbins_finest)
    # Compute TE over different partitions
    # ====================================
    tes = map(binscheme -> transferentropy(pts, vars, binscheme, estimator, b = b), binning_scheme)

    return tes
end

"""
    te_reg(source::AbstractArray{<:Real, 1},
        response::AbstractArray{<:Real, 1},
        k::Int, l::Int, m::Int,
        binning_scheme::Vector{RectangularBinning};
        η = 1, τ = 1,
        estimator = VisitationFrequency(),
        b = 2)

# TE with user-provided discretization scheme(s)

Calculate transfer entropy from `source` to `response` using the provided
`estimator` on discretizations constructed by the provided `binning_scheme`(s)
over a `k + l + m`-dimensional delay embedding of the input data,
using an embedding delay of `τ` across all embedding components.
`η` is the prediction lag.

## Arguments
- **`source`**: The source data series.
- **`target`**: The target data series.
- **`k`**: The dimension of the ``T_{f}`` component of the embedding.
- **`l`**: The dimension of the ``T_{pp}`` component of the embedding.
- **`m`**: The dimension of the ``S_{pp}`` component of the embedding.
- **`binning_scheme`**: The binning scheme(s) used to construct the partitions
    over which TE is computed. Must be either one or several instances of
    `RectangularBinning`s (provided as a vector). TE is computed for each
    of the resulting partitions.

## Keyword arguments

- **`τ`**: The embedding lag. Default is `τ = 1`.
- **`η`**: The prediction lag. Default is `η = 1`.
- **`estimator`**: The transfer entropy estimator to use. The default
    is `VisitationFrequency()`.
- **`b`**: Base of the logarithm. The default (`b = 2`) gives the TE in bits.

## More about the embedding

To compute transfer entropy, we need an appropriate delay embedding
of `source` (``S``) and `target` (``T``). For convenience, define

```math
\\begin{align}
T_f^{(k)} &= \\{(T(t+\\eta_k), \\ldots, T(t+\\eta_2), T(t+\\eta_1)) \\} \\\\
T_{pp}^{(l)} &= \\{ (T(t), T(t-\\tau_1), T(t-\\tau_2), \\ldots, T(t - \\tau_{l - 1})) \\} \\\\
S_{pp}^{(m)} &= \\{ (S(t), S(t-\\tau_1), S(t-\\tau_2), \\ldots, S(t-\\tau_{m - 1})) \\}
\\end{align}
```

where ``T_f`` denotes the `k`-dimensional set of vectors furnishing the future states of ``T``,
``T_{pp}`` denotes the `l`-dimensional set of vectors furnishing the past and present states of ``T``,
and ``S_{pp}`` denotes the `m`-dimensional set of vectors furnishing the past and present of ``S``.
``\\eta`` is the prediction lag. This convenience function uses ``\\tau_1`` = `τ`,
``\\tau_2`` = `2*τ`, ``\\tau_3`` = `3*τ`, and so on.

Combined, we get the generalised embedding ``\\mathbb{E} = (T_f^{(k)}, T_{pp}^{(l)}, S_{pp}^{(m)})``,
which is discretized as described below.

## More about discretization

The discretization scheme must be either a single `RectangularBinning` instance, or a vector of
`RectangularBinning` instances. Run `?RectangularBinning` after loading `CausalityTools` for
details.

## Transfer entropy computation

TE is then computed as

```math
\\begin{align}
TE_{S \\rightarrow T} = \\int_{\\mathbb{E}} P(T_f, T_{pp}, S_{pp}) \\log_{b}{\\left(\\frac{P(T_f | T_{pp}, S_{pp})}{P(T_f | T_{pp})}\\right)}
\\end{align}
```
using the provided `estimator` (default = `VisitationFrequency`) for
each of the discretizations. A vector of the TE estimates for each discretization
is returned.
"""
function te_reg(source::AbstractArray{<:Real, 1},
        response::AbstractArray{<:Real, 1},
        k::Int, l::Int, m::Int,
        binning_scheme::Vector{RectangularBinning},
        η = 1, τ = 1,
        estimator::TransferEntropyEstimator = VisitationFrequency(),
        b = 2)

    k + l + m >= 3 || throw(ArgumentError("`dim = k + l + m` must be 3 or higher for regular TE"))

    pts, vars = te_embed(source, response, k, l, m, η = η, τ = τ)

    # Compute TE over the partitions constructed from the provided binning schemes
    # ====================================
    tes = [transferentropy(pts, vars, bs, estimator, b = b) for bs in binning_scheme]

    return tes
end

function te_reg(source::AbstractArray{<:Real, 1},
        response::AbstractArray{<:Real, 1},
        k::Int, l::Int, m::Int,
        binning_scheme::RectangularBinning,
        η = 1, τ = 1,
        estimator::TransferEntropyEstimator = VisitationFrequency(),
        b = 2)

    k + l + m >= 3 || throw(ArgumentError("`dim = k + l + m` must be 3 or higher for regular TE"))

    pts, vars = te_embed(source, response, k, l, m, η = η, τ = τ)

    # Compute TE over the partitions constructed from the provided binning schemes
    # ====================================
    te = transferentropy(pts, vars, binning_scheme, estimator, b = b)

    return te
end

export te_reg


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

@doc doc"""
    abstract type MCMCAlgorithm end

!!! note

    The details of the `MCMCIterator` and `MCMCAlgorithm` API (see below)
    currently do not form part of the BAT-API, and may change without notice.

The following methods must be defined for subtypes (e.g.
for `SomeAlgorithm<:MCMCAlgorithm`):

```julia
MCMCIterator(
    rng::AbstractRNG,
    algorithm::SomeAlgorithm,
    density::AbstractDensity,
    chainid::Integer,
    [startpos::AbstractVector{<:Real}]
)
```

To implement a new MCMC algorithm, subtypes of both `MCMCAlgorithm` and
[`MCMCIterator`](@ref) are required.
"""
abstract type MCMCAlgorithm end
export MCMCAlgorithm

function get_mcmc_tuning end

@doc doc"""
abstract type MCMCInitAlgorithm end

Abstract type for MCMC initialization algorithms.
"""
abstract type MCMCInitAlgorithm end
export MCMCInitAlgorithm

"""
MCMCTuningAlgorithm

Abstract super-type for MCMC tuning algorithms.
"""
abstract type MCMCTuningAlgorithm end
export MCMCTuningAlgorithm

"""
MCMCBurninAlgorithm

Abstract super-type for MCMC burn-in algorithms.
"""
abstract type MCMCBurninAlgorithm end
export MCMCBurninAlgorithm

"""
MCMCBurninAlgorithm

Abstract super-type for MCMC convergence tests.
"""
abstract type MCMCConvergenceTest end
export MCMCConvergenceTest

@with_kw struct MCMCIteratorInfo
    id::Int32
    cycle::Int32
    tuned::Bool
    converged::Bool
end

@doc doc"""
    abstract type MCMCIterator end

Represents the current state of a MCMC chain.

!!! note

    The details of the `MCMCIterator` and `MCMCAlgorithm` API (see below)
    currently do not form part of the BAT-API, and may change without notice.

To implement a new MCMC algorithm, subtypes of both [`MCMCAlgorithm`](@ref)
and `MCMCIterator` are required.

The following methods must be defined for subtypes of `MCMCIterator` (e.g.
`SomeMCMCIter<:MCMCIterator`):

```julia
BAT.getalgorithm(chain::SomeMCMCIter)::MCMCAlgorithm

BAT.getdensity(chain::SomeMCMCIter)::AbstractDensity

BAT.getrng(chain::SomeMCMCIter)::AbstractRNG

BAT.mcmc_info(chain::SomeMCMCIter)::MCMCIteratorInfo

BAT.nsteps(chain::SomeMCMCIter)::Int

BAT.nsamples(chain::SomeMCMCIter)::Int

BAT.current_sample(chain::SomeMCMCIter)::DensitySample

BAT.sample_type(chain::SomeMCMCIter)::Type{<:DensitySample}

BAT.samples_available(chain::SomeMCMCIter, nonzero_weights::Bool = false)::Bool

BAT.get_samples!(samples::DensitySampleVector, chain::SomeMCMCIter, nonzero_weights::Bool)::typeof(samples)

BAT.next_cycle!(chain::SomeMCMCIter)::SomeMCMCIter

BAT.mcmc_step!(
    chain::SomeMCMCIter
    callback::Function,
)::nothing
```

The following methods are implemented by default:

```julia
getalgorithm(chain::MCMCIterator)
getdensity(chain::MCMCIterator)
DensitySampleVector(chain::MCMCIterator)
mcmc_iterate!(chain::MCMCIterator, ...)
mcmc_iterate!(chains::AbstractVector{<:MCMCIterator}, ...)
isvalidchain(chain::MCMCIterator)
isviablechain(chain::MCMCIterator)
```
"""
abstract type MCMCIterator end
export MCMCIterator

function getalgorithm end

function getdensity end

function getrng end

function mcmc_info end

function nsteps end

function nsamples end

function current_sample end

function sample_type end

function samples_available end

function get_samples! end

function next_cycle! end

function mcmc_step! end

DensitySampleVector(chain::MCMCIterator) = DensitySampleVector(sample_type(chain), totalndof(getdensity(chain)))

function mcmc_iterate! end

function mcmc_iterate!(
    output::Union{DensitySampleVector,Nothing},
    chain::MCMCIterator;
    max_nsteps::Integer = 1,
    max_time::Real = Inf,
    nonzero_weights::Bool = true,
    callback::Function = nop_func
)
    @debug "Starting iteration over MCMC chain $(chain.info.id), max_nsteps = $max_nsteps, max_time = $max_time"

    start_time = time()
    start_nsteps = nsteps(chain)
    start_nsamples = nsamples(chain)

    while (
        (nsteps(chain) - start_nsteps) < max_nsteps &&
        (time() - start_time) < max_time
    )
        mcmc_step!(chain)
        callback(Val(:mcmc_step), chain)
        if !isnothing(output)
            get_samples!(output, chain, nonzero_weights)
        end
    end

    end_time = time()
    elapsed_time = end_time - start_time

    @debug "Finished iteration over MCMC chain $(chain.info.id), nsteps = $(nsteps(chain)), time = $(Float32(elapsed_time))"

    nothing
end

function mcmc_iterate!(
    outputs::Union{AbstractVector{<:DensitySampleVector},Nothing},
    chains::AbstractVector{<:MCMCIterator};
    kwargs...
)
    if isempty(chains)
        @debug "No MCMC chain(s) to iterate over."
        return chains
    else
        @debug "Starting iteration over $(length(chains)) MCMC chain(s)"
    end

    if !isnothing(outputs)
        @sync for i in eachindex(outputs, chains)
            Base.Threads.@spawn mcmc_iterate!(outputs[i], chains[i]; kwargs...)
        end
    else
        @sync for i in eachindex(chains)
            Base.Threads.@spawn mcmc_iterate!(nothing, chains[i]; kwargs...)
        end
    end

    nothing
end

isvalidchain(chain::MCMCIterator) = current_sample(chain).logd > -Inf

isviablechain(chain::MCMCIterator) = nsamples(chain) >= 2

"""
    BAT.MCMCSampleGenerator

*BAT-internal, not part of stable public API.*

MCMC sample generator.

Constructors:

```julia
MCMCSampleGenerator(chain::AbstractVector{<:MCMCIterator})
```
"""
struct MCMCSampleGenerator{T<:AbstractVector{<:MCMCIterator}} <: AbstractSampleGenerator
    chains::T
end

getalgorithm(sg::MCMCSampleGenerator) = sg.chains[1].spec.algorithm

abstract type AbstractMCMCTunerInstance end

function tuning_init! end

function mcmc_init! end

function mcmc_burnin! end

function isvalidchain end

function isviablechain end


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

using Markdown
using InteractiveUtils

# ╔═╡ d2ac9e5e-8bb6-4fa1-9cb3-957d238d7355
md"""
# Despacho múltiple / multiple dispatch
"""

# ╔═╡ 60cc9b65-62bd-4795-a9d9-30ff6a5fdc48
md"""
En este notebook vamos a explorar **multiple dispatch**, un concepto fundamental en Julia.

Multiple dispatch permite software:
- rápido
- extendible
- programable
- divertido para experimentar
"""

# ╔═╡ dda25ade-80d1-429a-a285-dbc477305966
md"""
Para entender el despacho múltiple en Julia, observemos el operador `+`.

Si llamamos `methods()` sobre `+`, podemos ver todas las definiciones de `+`
"""

# ╔═╡ d232a85b-62c8-4de5-8798-637af323b492
md"""
Podemos usar el macro de `@which` para saber qué método en particular estamos usando de `+`.

Distintos métodos se usan en cada uno de estos ejemplos.
"""

# ╔═╡ 495ee98c-e2b8-475f-8c36-b870c816d33b
md"""
Aún más, pues podemos definir nuevos métodos de `+`.

Primero tenemos que importar `+` de Base.
"""

# ╔═╡ 37b52b24-d323-47e4-87db-fcbe204ca96e
import Base: +

# ╔═╡ 19b8a372-4ea5-4ff4-8a6a-f48aa34430f4
methods(+);

# ╔═╡ 6426c7fa-c3ed-4d06-be67-355ac10ebfea
@which 3 + 3

# ╔═╡ b297c837-e937-4da0-acd8-933df7fee791
@which 3.0 + 3.0

# ╔═╡ 9896765a-8906-4877-b8e5-b5c0fc4b09f8
@which 3 + 3.0

# ╔═╡ 76151c1b-139f-40c6-b111-0b0634473e17
md"""
Digamos que queremos concatenar elementos con  `+`. Sin extender el método, no funciona
"""

# ╔═╡ 4f27b734-90ce-4502-9df2-b7d69c230266
"hello " + "world!"

# ╔═╡ 317c59a1-fe9d-47bc-acd6-d25f8a2bcbc7
@which "hello " + "world!"

# ╔═╡ 74997f65-4f75-4aa8-88f1-49a96a724e82
md"""
Entonces agregamos a `+` un método que toma dos cadenas y las concatena
"""

# ╔═╡ 1df76000-3434-4be6-b49d-6839436f1c6e
Base.:+(x::String, y::String) = string(x, y)

# ╔═╡ 37d06d98-4d74-4d29-9ec9-39fb0d31d31a
"hello " + "world!"

# ╔═╡ 1dbf2eb6-87ff-425e-b9bb-4039538c54ce
md"""
¡Funciona! Y si queremos más, podemos comprobarnos que Julia ha despachado sobre los tipos de "hello" y "world!", sobre el método que acabamos de definir
"""

# ╔═╡ 66dccffc-d0f1-4974-a673-4941cdc61429
@which "hello " + "world!"

# ╔═╡ b87d4755-802e-45a5-9478-2aa69af92a62
md"""
Vamos por un ejemplo más
"""

# ╔═╡ c5397d75-f538-4346-95e5-0c30f4d203e0
begin
	foo(x, y) = println("duck-typed foo!")
	foo(x::Int, y::Float64) = println("foo con entero y flotante!")
	foo(x::Float64, y::Float64) = println("foo con dos flotantes!")
	foo(x::Int, y::Int) = println("foo con dos enteros!")
end

# ╔═╡ b16e8931-5887-4335-bafa-8a06969d679f
foo(1, 1)

# ╔═╡ 33e60ea0-8788-477e-a45d-dd5c92f36929
foo(1., 1.)

# ╔═╡ a007b5cf-3ebf-473d-b640-32c1d9182fe9
foo(1, 1.0)

# ╔═╡ c135586a-4c72-4489-b143-c18226000bcc
foo(true, false)

# ╔═╡ 4a1d67f1-ee84-49e5-94ec-84a8c1ec965e
md"""
Nota que el último ejemplo aplica por default el caso de 'duck-typed foo' porque no había un método definido exclusivamente para dos booleanos.
"""

# ╔═╡ 976d280a-1474-47aa-9856-eb5e4ee6ab49
md"""
### Ejercicios

#### 9.1

Agrega un método para `+` que aplique un cifrado de César a una cadena (cómo en el notebook 6) tal que

```julia
"hello" + 4 == "lipps"
```
"""

# ╔═╡ ddd66730-c1d2-4a28-9b9e-c4aa66057e4e

# ╔═╡ 9cfd983b-0b80-4712-a394-d95b33e79ca8
md"""
#### 9.2

Checa que has extendido propiamente `+` recorriendo la próxima cadena para atrás por 3 letras:

"Gr#qrw#phggoh#lq#wkh#diidluv#ri#gudjrqv#iru#|rx#duh#fuxqfk|#dqg#wdvwh#jrrg#zlwk#nhwfkxs1"
"""

# ╔═╡ 09d662a9-9ce6-4718-9740-a79cf6ac18dd

# ╔═╡ Cell order:
# ╟─d2ac9e5e-8bb6-4fa1-9cb3-957d238d7355
# ╟─60cc9b65-62bd-4795-a9d9-30ff6a5fdc48
# ╟─dda25ade-80d1-429a-a285-dbc477305966
# ╠═19b8a372-4ea5-4ff4-8a6a-f48aa34430f4
# ╟─d232a85b-62c8-4de5-8798-637af323b492
# ╠═6426c7fa-c3ed-4d06-be67-355ac10ebfea
# ╠═b297c837-e937-4da0-acd8-933df7fee791
# ╠═9896765a-8906-4877-b8e5-b5c0fc4b09f8
# ╟─495ee98c-e2b8-475f-8c36-b870c816d33b
# ╠═37b52b24-d323-47e4-87db-fcbe204ca96e
# ╟─76151c1b-139f-40c6-b111-0b0634473e17
# ╠═4f27b734-90ce-4502-9df2-b7d69c230266
# ╠═317c59a1-fe9d-47bc-acd6-d25f8a2bcbc7
# ╟─74997f65-4f75-4aa8-88f1-49a96a724e82
# ╠═1df76000-3434-4be6-b49d-6839436f1c6e
# ╠═37d06d98-4d74-4d29-9ec9-39fb0d31d31a
# ╟─1dbf2eb6-87ff-425e-b9bb-4039538c54ce
# ╠═66dccffc-d0f1-4974-a673-4941cdc61429
# ╟─b87d4755-802e-45a5-9478-2aa69af92a62
# ╠═c5397d75-f538-4346-95e5-0c30f4d203e0
# ╠═b16e8931-5887-4335-bafa-8a06969d679f
# ╠═33e60ea0-8788-477e-a45d-dd5c92f36929
# ╠═a007b5cf-3ebf-473d-b640-32c1d9182fe9
# ╠═c135586a-4c72-4489-b143-c18226000bcc
# ╟─4a1d67f1-ee84-49e5-94ec-84a8c1ec965e
# ╟─976d280a-1474-47aa-9856-eb5e4ee6ab49
# ╠═ddd66730-c1d2-4a28-9b9e-c4aa66057e4e
# ╟─9cfd983b-0b80-4712-a394-d95b33e79ca8
# ╠═09d662a9-9ce6-4718-9740-a79cf6ac18dd

max_stars_repo_name: miguelraz/JuliaTutorials

max_issues_repo_name: miguelraz/JuliaTutorials

max_forks_repo_name: miguelraz/JuliaTutorials


text:
## 10.jl
# 行数のカウント

f = open("./dat/hightemp.txt")
println(length(readlines(f)))
close(f)


text:
using GPUArrays, Test

include("testsuite.jl")

@testset "JLArray" begin
    using GPUArrays.JLArrays
    TestSuite.test(JLArray)
end