text:
export sample_priors
sample_priors(planet::Planet) = rand.(ComponentArray(planet.priors.priors))
sample_priors(planet::Planet,N) = [sample_priors(planet) for _ in 1:N]

function sample_priors(system::System)
    sampled = ComponentVector(
        merge(NamedTuple(rand.(ComponentArray(system.priors.priors))),
        # (;planets=[sample_priors(planet) for planet in system.planets])
        (;planets=namedtuple(collect(keys(system.planets)), [
            ComponentArray(NamedTuple([k=>v for (k,v) in pairs(NamedTuple(sample_priors(planet)))]))
            for planet in system.planets
        ]))
    ))
    return getdata(sampled)
end

sample_priors(system::System,N) = [sample_priors(system) for _ in 1:N]

# Instead of just calling mean for the distributions, we sample and then take the mean of the data.
# This does add a little jitter, but some distributions do not directly define the mean function!
# Specifically, truncated(InverseGamma()) does not work, and this is very useful.
# mean_priors(planet::Planet) = Statistics.mean.(Statistics.rand.(planet.priors.priors,1000))
function mean_priors(system::System)
    priors_all = ComponentVector(;
        NamedTuple(system.priors.priors)...,
        planets=[planet.priors.priors for planet in system.planets]
    )
    # return Statistics.mean.(Statistics.rand.(priors_all,1000))
    return Statistics.mean.(priors_all)
end

function guess_starting_position(system, N=500_000)

    @info "Guessing a good starting location by sampling from priors" N
    # TODO: this shouldn't have to allocate anything, we can just loop keeping the best.
    θ = sample_priors(system, N)
    arr2nt = DirectDetections.make_arr2nt(system)
    θr = arr2nt.(θ)

    posts = zeros(N)
    ln_prior = make_ln_prior(system)
    arr2nt = DirectDetections.make_arr2nt(system)
    Threads.@threads for i in eachindex(posts)
        θ_res = arr2nt(θ[i])
        posts[i] = ln_prior(θ[i]) + ln_like(θ_res, system)
    end
    # posts = map(eachrow(A)) do c
    #     DirectDetections.ln_post(ComponentVector(c, ax), system)
    # end
    mapv,mapi = findmax(posts)
    best = θ[mapi]

    # @info "Found good location" mapv best=NamedTuple(best)

    return best
end

using Optim
function optimize_starting_position(ℓπ, initial_θ_t)

    @info "Optimizing starting location "
    results = optimize(θ_t-> -ℓπ(θ_t), initial_θ_t, LBFGS(), autodiff = :forward, Optim.Options(iterations=100_000, time_limit=5, allow_f_increases=true))
    # results = optimize(θ_t-> -ℓπ(θ_t), initial_θ_t, NelderMead(), Optim.Options(iterations=100_000, time_limit=5))
    # results = optimize(θ_t-> -ℓπ(θ_t), initial_θ_t, ParticleSwarm(), Optim.Options(iterations=100_000, time_limit=5))

    intial_objective = ℓπ(initial_θ_t)

    minimizer = Optim.minimizer(results)

    @info "Starting location improved" logimprovement=(ℓπ(minimizer) - intial_objective)
    return minimizer
end

"""
    construct_elements(θ_system, θ_planet)

Given a named tuple for of parameters from a System (θ_system) and Planet (θ_planet),
return a `KeplerianElements` from DirectOrbits.jl.
"""
function construct_elements(θ_system, θ_planet)
    return KeplerianElements((;
        θ_system.μ,
        θ_system.plx,
        θ_planet.i,
        θ_planet.Ω,
        θ_planet.ω,
        θ_planet.e,
        θ_planet.τ,
        θ_planet.a,
    ))
end

"""
    construct_elements(chains, :b, 4)

Given a Chains object, a symbol matching the name of a planet, and an index,
construct a `KeplerianElements` from DirectOrbits of that planet from that
index of the chains.
"""
function construct_elements(chain::Chains, planet_key::Union{String,Symbol}, i::Union{Integer,CartesianIndex})
    pk = string(planet_key)
    return KeplerianElements((;
        μ=chain["μ"][i],
        plx=chain["plx"][i],
        i=chain[pk*"[i]"][i],
        Ω=chain[pk*"[Ω]"][i],
        ω=chain[pk*"[ω]"][i],
        e=chain[pk*"[e]"][i],
        τ=chain[pk*"[τ]"][i],
        a=chain[pk*"[a]"][i],
    ))
end

"""
    construct_elements(chains, :b, [4,5,10])

Given a Chains object, a symbol matching the name of a planet, and an array of indices,
construct a `KeplerianElements` from DirectOrbits of that planet from those indices
of the chains.
"""
function construct_elements(chain::Chains, planet_key::Union{String,Symbol}, ii::AbstractArray{<:Union{Integer,CartesianIndex}})
    pk = string(planet_key)
    μs=chain["μ"]
    plxs=chain["plx"]
    is=chain[pk*"[i]"]
    Ωs=chain[pk*"[Ω]"]
    ωs=chain[pk*"[ω]"]
    es=chain[pk*"[e]"]
    τs=chain[pk*"[τ]"]
    as=chain[pk*"[a]"]
    return map(ii) do i
        KeplerianElements((;
            μ=μs[i],
            plx=plxs[i],
            i=is[i],
            Ω=Ωs[i],
            ω=ωs[i],
            e=es[i],
            τ=τs[i],
            a=as[i],
        ))
    end
end

# Fallback when no random number generator is provided (as is usually the case)
function hmc(system::System, target_accept::Number=0.8, ensemble::AbstractMCMC.AbstractMCMCEnsemble=MCMCSerial(); kwargs...)
    return hmc(Random.default_rng(), system, target_accept, ensemble; kwargs...)
end

function hmc(
    rng::Random.AbstractRNG,
    system::System, target_accept::Number=0.8,
    ensemble::AbstractMCMC.AbstractMCMCEnsemble=MCMCSerial();
    num_chains=1,
    adaptation,
    iterations,
    discard_initial=adaptation,
    tree_depth=10,
    initial_samples=50_000,
    initial_parameters=nothing,
    step_size=nothing,
    progress=true
)

    # Choose parameter dimensionality and initial parameter value
    initial_θ_0 = sample_priors(system)
    D = length(initial_θ_0)

    ln_prior_transformed = make_ln_prior_transformed(system)
    arr2nt = DirectDetections.make_arr2nt(system)

    priors_vec = _list_priors(system)
    Bijector_invlinkvec = make_Bijector_invlinkvec(priors_vec)

    # Capture these variables in a let binding to improve performance
    ℓπ = let system=system, ln_prior_transformed=ln_prior_transformed, arr2nt=arr2nt, Bijector_invlinkvec=Bijector_invlinkvec
        function (θ_t)
            # Transform back from the unconstrained support to constrained support for the likelihood function
            θ = Bijector_invlinkvec(θ_t)
            θ_res = arr2nt(θ)
            ll = ln_prior_transformed(θ) + ln_like(θ_res, system)
            return ll
        end
    end

    if isnothing(initial_parameters)
        initial_θ = guess_starting_position(system,initial_samples)
        # Transform from constrained support to unconstrained support
        initial_θ_t = Bijectors.link.(priors_vec, initial_θ)

        # initial_θ_guess = guess_starting_position(system,initial_samples)
        # # Transform from constrained support to unconstrained support
        # initial_θ_guess_t = Bijectors.link.(priors_vec, initial_θ_guess)

        # # @show initial_θ_guess
        # initial_θ_t = optimize_starting_position(ℓπ, initial_θ_guess_t)

        # # Just for display
        # initial_θ = Bijectors.invlink.(priors_vec, initial_θ_t)

        # # @show initial_θ_guess initial_θ
        # # @show priors_vec

    else
        initial_θ = initial_parameters
        # Transform from constrained support to unconstrained support
        initial_θ_t = Bijectors.link.(priors_vec, initial_θ)
    end

    # Define a Hamiltonian system
    metric = DenseEuclideanMetric(D)
    hamiltonian = Hamiltonian(metric, ℓπ, ForwardDiff)

    if !isnothing(step_size)
        initial_ϵ = step_size
    else
        initial_ϵ = find_good_stepsize(hamiltonian, initial_θ_t)
        @info "Found initial stepsize" initial_ϵ
    end

    integrator = Leapfrog(initial_ϵ)
    # integrator = TemperedLeapfrog(initial_ϵ, 1.05)

    # # We have change some parameters when running with image data
    if !isnothing(system.images) && target_accept > 0.6 && isnothing(step_size)
        @warn "Sampling from images with target accept greater than 0.6. This may lead to insufficient exploration."
    end

    mma = MassMatrixAdaptor(metric)
    if isnothing(step_size)
        @info "Will adapt step size and mass matrix"
        ssa = StepSizeAdaptor(target_accept, integrator)
        adaptor = StanHMCAdaptor(mma, ssa)
    else
        @info "Will adapt mass matrix only"
        adaptor = MassMatrixAdaptor(metric)
    end

    model = AdvancedHMC.DifferentiableDensityModel(ℓπ, ForwardDiff)

    # κ = NUTS{MultinomialTS,GeneralisedNoUTurn}(integrator, max_depth=tree_depth)
    # κ = NUTS{SliceTS, StrictGeneralisedNoUTurn}(integrator, max_depth=tree_depth)

    # Had some good results with this one:
    # κ = NUTS{MultinomialTS, StrictGeneralisedNoUTurn}(integrator, max_depth=tree_depth)

    κ = NUTS(integrator, max_depth=tree_depth)
    sampler = AdvancedHMC.HMCSampler(κ, metric, adaptor)

    AbstractMCMC.setprogress!(true)
    start_time = time()

    # Neat: it's possible to return a live iterator
    # We could use this to build e.g. live plotting while the code is running
    # once the analysis code is ported to Makie.
    # return  AbstractMCMC.steps(
    #     rng,
    #     model,
    #     sampler,
    #     nadapts = adaptation,
    #     init_params = initial_θ_t,
    #     discard_initial = adaptation,
    #     progress=progress,
    #     verbose=false
    # )

    # function callback(rng, model, sampler, transition, state, iteration; kwargs...)
    #     # nadapts init_params
    #     @show propertynames(kwargs)
    #     @show iteration
    # end

    mc_samples_all_chains = sample(
        rng,
        model,
        sampler,
        ensemble,
        iterations+adaptation, num_chains;
        nadapts = adaptation,
        init_params = initial_θ_t,
        discard_initial,
        progress=progress,
        # callback
    )
    stop_time = time()

    @info "Sampling compete. Building chains."
    # Go through each chain and repackage results
    chains = MCMCChains.Chains[]
    logposts = Vector{Float64}[]
    for (i,mc_samples) in enumerate(mc_samples_all_chains)
        stat = map(s->s.stat, mc_samples)
        logpost = map(s->s.z.ℓπ.value, mc_samples)

        mean_accept = mean(getproperty.(stat, :acceptance_rate))
        num_err_frac = mean(getproperty.(stat, :numerical_error))
        mean_tree_depth = mean(getproperty.(stat, :tree_depth))
        max_tree_depth_frac = mean(getproperty.(stat, :tree_depth) .== tree_depth)

        println("""\
        Sampling report for chain $i:
        mean_accept =         $mean_accept
        num_err_frac =        $num_err_frac
        mean_tree_depth =     $mean_tree_depth
        max_tree_depth_frac = $max_tree_depth_frac
        """)

        # Report some warnings if sampling did not work well
        if num_err_frac == 1.0
            @error "Numerical errors encountered in ALL iterations. Check model and priors."
        elseif num_err_frac > 0.1
            @warn "Numerical errors encountered in more than 10% of iterations" num_err_frac
        end
        if max_tree_depth_frac > 0.1
            @warn "Maximum tree depth hit in more than 10% of iterations (reduced efficiency)" max_tree_depth_frac
        end

        logpost = map(s->s.z.ℓπ.value, mc_samples)

        # Transform samples back to constrained support
        samples = map(mc_samples) do s
            θ_t = s.z.θ
            θ = Bijectors.invlink.(priors_vec, θ_t)
            return θ
        end
        chain_res = arr2nt.(samples)
        push!(chains, DirectDetections.result2mcmcchain(system, chain_res))
        push!(logposts, logpost)
    end

    # Concatenate the independent chains now that we have remapped / resolved the variables.
    mcmcchains = AbstractMCMC.chainscat(chains...)

    # Concatenate the log posteriors and make them the same shape as the chains (N_iters,N_vars,N_chains)
    logposts_mat = reduce(hcat, logposts)
    mcmcchains_with_info = MCMCChains.setinfo(
        mcmcchains,
        (;
            start_time,
            stop_time,
            model=system,
            adaptor,
            mc_samples=mc_samples_all_chains,
            sampler,
            logpost=logposts_mat,
            initial_ϵ
        )
    )
    return mcmcchains_with_info
end

include("tempered-sampling.jl")

"""
Convert a vector of component arrays returned from sampling into an MCMCChains.Chains
object.
"""
function result2mcmcchain(system, chains_in_0)
    chains_in = ComponentArray.(chains_in_0)
    # `system` not currently used, but a more efficient/robust mapping in future might require it.

    # There is a specific column name convention used by MCMCChains to indicate
    # that multiple parameters form a group. Instead of planets.X.a, we adapt our to X[a]
    # accordingly
    flattened_labels = replace.(labels(first(chains_in)), r"planets\.([^\.]+).([^\.]+)" => s"\1[\2]")
    c = Chains(
        reduce(vcat, getdata.(chains_in)'),
        flattened_labels
    )
    return c
end


text:
# NOTE: This is an example of how to do Mv Adaptive Metropolis for a GMM.
# Updates are done iteratively for each set of parameters, instead of jointly
# for all parameters.

import Pkg; Pkg.activate("../")
using MCMC
using BangBang
using Bijectors
using Distributions
import Random

# Set random seed for reproducibility
Random.seed!(0)

# Define model.
struct Gmm <: MCMC.Model
  y::Vector{Float64}
  K::Int
end

# Create initial state maker.
MCMC.make_init_state(m::Gmm) = (mu=randn(m.K), sigma=rand(m.K), eta=fill(1/m.K, m.K))

# Define loglikelihood and full conditionals
loglike(m, s) = sum(logpdf.(UnivariateGMM(s.mu, s.sigma, Categorical(s.eta)), m.y))
function sigma_cond(m, s, x)
  return sum(logpdf.(LogNormal(0, 1), x)) + loglike(m, setproperty!!(s, :sigma, x))
end
function mu_cond(m, s, x)
  return sum(logpdf.(Normal(0, 1), x)) + loglike(m, setproperty!!(s, :mu, x))
end
function eta_cond(m, s, x)
  return logpdf(Dirichlet(m.K, 1/m.K), x) + loglike(m, setproperty!!(s, :eta, x))
end

# Define true model.
true_model = MixtureModel(Normal.([-2, 0, 2], [.2, .2, .3]), [.3, .4, .3])

# Simulate data.
y = rand(true_model, 500)

# Instantiate model.
model = Gmm(y, 3)

# Initialize model state.
init = MCMC.make_init_state(model)

# Create sampler.
r1 = RWM(:mu, mu_cond, mvarwm(init.mu))
r2 = RWM(:sigma, sigma_cond, mvarwm(init.sigma), bijector=Bijectors.Log{1}())
r3 = RWM(:eta, eta_cond, mvarwm(init.eta), bijector=Bijectors.SimplexBijector{true}())
spl = MCMC.Gibbs(model, r1, r2, r3)

# Sample via Adaptive Metropolis.
chain = mcmc(spl, 2000, nburn=10000, thin=2).chain

# Check model fit.
ord = sortperm(mean(getindex.(chain, :mu)))
@assert isapprox(mean(getindex.(chain, :mu))[ord], mean.(true_model.components), atol=0.1)
@assert isapprox(mean(getindex.(chain, :sigma))[ord], std.(true_model.components), atol=0.1)
@assert isapprox(mean(getindex.(chain, :eta))[ord], true_model.prior.p, atol=0.05)

ext: jl


text:
# Autogenerated wrapper script for JpegTurbo_jll for x86_64-linux-musl
export cjpeg, djpeg, jpegtran, libjpeg, libturbojpeg

JLLWrappers.@generate_wrapper_header("JpegTurbo")
JLLWrappers.@declare_executable_product(cjpeg)
JLLWrappers.@declare_executable_product(djpeg)
JLLWrappers.@declare_executable_product(jpegtran)
JLLWrappers.@declare_library_product(libjpeg, "libjpeg.so.62")
JLLWrappers.@declare_library_product(libturbojpeg, "libturbojpeg.so.0")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        cjpeg,
        "bin/cjpeg",
    )

    JLLWrappers.@init_executable_product(
        djpeg,
        "bin/djpeg",
    )

    JLLWrappers.@init_executable_product(
        jpegtran,
        "bin/jpegtran",
    )

    JLLWrappers.@init_library_product(
        libjpeg,
        "lib/libjpeg.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libturbojpeg,
        "lib/libturbojpeg.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# test nullhypo with n-dim partial

using Test
using IncrementalInference

##

@testset "test n-dimensional partial" begin

##

fg = initfg()
addVariable!(fg, :x0, ContinuousEuclid{3})

addFactor!(fg, [:x0;], PartialPrior(MvNormal(zeros(2), ones(2)), (2,3)) )

addVariable!(fg, :x1, ContinuousEuclid{3})
addFactor!(fg, [:x1;], PartialPrior(Normal(10,1),(1,)))
addFactor!(fg, [:x0; :x1], LinearRelative(MvNormal([10;0;0.0], ones(3))) )

##

initAll!(fg)

##

destlbl = :x0

dens = Vector{ManifoldKernelDensity}()
factors = getFactor.(fg, ls(fg, destlbl))
inferdim = IIF.proposalbeliefs!(fg, destlbl, factors, dens )

oldBel = getBelief(fg, destlbl)
oldpts = getPoints(oldBel)

varType = getVariableType(fg, destlbl)
pGM = getPoints( AMP.manifoldProduct(dens, getManifold(varType), N=100, oldPoints=oldpts), false )
# pGM = AMP.productbelief(oldpts, getManifold(varType), dens, 100, asPartial=false )

##

densPts, inferdim = predictbelief(fg, :x0, :, needFreshMeasurements=true )

##

solveTree!(fg);

##

@warn "WIP on testPartialNH.jl during transition to Manifolds.jl"
@test isapprox( getPPE(fg, :x0).suggested, [0;0;0], atol=1)
@test isapprox( getPPE(fg, :x1).suggested, [10;0;0], atol=1)

##

end

@testset "test n-dimensional partial with nullhypo" begin

##

fg = initfg()
addVariable!(fg, :x0, ContinuousEuclid{3})

addFactor!(fg, [:x0;], PartialPrior(MvNormal(zeros(2), ones(2)),(2,3)) , nullhypo=0.2)

addVariable!(fg, :x1, ContinuousEuclid{3})
addFactor!(fg, [:x1;], PartialPrior(Normal(10,1),(1,)))
addFactor!(fg, [:x0; :x1], LinearRelative(MvNormal([10;0;0.0], ones(3))) , nullhypo=0.2)

##

solveTree!(fg);

##

@warn "WIP testPartialNH.jl during transition to Manifolds.jl"
@test isapprox( getPPE(fg, :x0).suggested, [0;0;0], atol=1)
@test isapprox( getPPE(fg, :x1).suggested, [10;0;0], atol=1)

##

end

#


text:
"""
    psMF(lnν)

Press-Schechter halo multiplicity function.

*Reference*: Press & Schechter, 187, 425 (1974)
- See also Bond, Cole, Efstathiou & Kaiser, ApJ, 379, 440 (1991)
- See Sheth & Tormen, MNRAS, 329, 61 (2002) for the notation used here. `psMF(lnν)` corresponds to νf(ν) in this reference.

# Arguments
- `lnν::Real`: natural logarithm of a threshold, ν, i.e., `lnν` = log(ν), defined by ν ≡ [δc/σ(R)]^2. Here, δc = 1.6865 and σ(R) is the r.m.s. mass fluctuation within a top-hat smoothing of scale R.

Press & Schechter's multiplicity function is normalized, i.e,

``∫_-∞^∞ dlnν psMF(lnν) = 1``

This function is based on [Cosmology Routine Library (CRL)](https://wwwmpa.mpa-garching.mpg.de/~komatsu/crl/).
"""
function psMF(lnν::Real)
    ν = exp(lnν)
    mf = exp(-ν / 2) * √(ν / 2π)
end

"""
    stMF(lnν)

Sheth-Tormen halo multiplicity function.

*Reference*: Equation (2) of Sheth & Tormen, MNRAS, 329, 61 (2002)
- See also Sheth, Mo & Tormen, MNRAS, 323, 1 (2001)
- `stMF(lnν)` corresponds to νf(ν) in these references

# Arguments
- `lnν::Real`: natural logarithm of a threshold, ν, i.e., `lnν` = log(ν), defined by ν ≡ [δc/σ(R)]^2. Here, δc = 1.6865 and σ(R) is the r.m.s. mass fluctuation within a top-hat smoothing of scale R.

Sheth & Tormen's multiplicity function is normalized, i.e,

``∫_-∞^∞ dlnν stMF(lnν) = 1``

This function is based on [Cosmology Routine Library (CRL)](https://wwwmpa.mpa-garching.mpg.de/~komatsu/crl/).
"""
function stMF(lnν::Real)
    A, p, q = 0.3222, 0.3, 0.707
    ν = exp(lnν)
    mf = A * (1 + (q * ν)^-p) * exp(-q * ν / 2) * √(q * ν / 2π)
end

"""
    jenkinsMF(lnν)

Jenkins et al.'s halo multiplicity function.

*Reference*: Equation (B3) of Jenkins et al., MNRAS, 321, 372 (2001)
- `jenkinsMF(lnν)` corresponds to f(σ)/2 in this reference, where σ=√(δc/ν).
- We use the fitting function for Δm = 180. Here, Δm is an overdensity within a spherical region of radius R, whose mean density is equal to Δm times the mean **mass** density of the Universe.

# Arguments
- `lnν::Real`: natural logarithm of a threshold, ν, i.e., `lnν` = log(ν), defined by ν ≡ [δc/σ(R)]^2. Here, δc = 1.6865 and σ(R) is the r.m.s. mass fluctuation within a top-hat smoothing of scale R.

Jenkins et al.'s multiplicity function is NOT normalized, i.e,

``∫_-∞^∞ dlnν jenkinsMF(lnν) ≠ 1``

This function is based on [Cosmology Routine Library (CRL)](https://wwwmpa.mpa-garching.mpg.de/~komatsu/crl/).
"""
function jenkinsMF(lnν::Real)
    ν = exp(lnν)
    σ = 1.6865 / √ν
    mf = 0.5 * 0.301 * exp(-abs(log(1 / σ) + 0.64)^3.82)
end


text:
function get_thermal_standard_uc_template()
    template = OperationsProblemTemplate(CopperPlatePowerModel)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    set_device_model!(template, ThermalStandard, ThermalStandardUnitCommitment)
    return template
end

function get_thermal_dispatch_template_network(network = CopperPlatePowerModel)
    template = OperationsProblemTemplate(network)
    set_device_model!(template, ThermalStandard, ThermalDispatch)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    set_device_model!(template, MonitoredLine, StaticBranchBounds)
    set_device_model!(template, Line, StaticBranch)
    set_device_model!(template, Transformer2W, StaticBranch)
    set_device_model!(template, TapTransformer, StaticBranch)
    set_device_model!(template, HVDCLine, HVDCLossless)
    return template
end

function get_template_basic_uc_simulation()
    template = OperationsProblemTemplate(CopperPlatePowerModel)
    set_device_model!(template, ThermalStandard, ThermalBasicUnitCommitment)
    set_device_model!(template, RenewableDispatch, RenewableFullDispatch)
    set_device_model!(template, PowerLoad, StaticPowerLoad),
    set_device_model!(template, InterruptibleLoad, StaticPowerLoad)
    set_device_model!(template, HydroEnergyReservoir, HydroDispatchRunOfRiver)
    return template
end

function get_template_standard_uc_simulation()
    template = get_template_basic_uc_simulation()
    set_device_model!(template, ThermalStandard, ThermalStandardUnitCommitment)
    return template
end

function get_template_nomin_ed_simulation(network = CopperPlatePowerModel)
    template = OperationsProblemTemplate(network)
    set_device_model!(template, ThermalStandard, ThermalDispatchNoMin)
    set_device_model!(template, RenewableDispatch, RenewableFullDispatch)
    set_device_model!(template, PowerLoad, StaticPowerLoad),
    set_device_model!(template, InterruptibleLoad, DispatchablePowerLoad)
    set_device_model!(template, HydroEnergyReservoir, HydroDispatchRunOfRiver)
    return template
end

function get_template_hydro_st_uc(network = CopperPlatePowerModel)
    template = OperationsProblemTemplate(network)
    set_device_model!(template, ThermalStandard, ThermalStandardUnitCommitment),
    set_device_model!(template, RenewableDispatch, RenewableFullDispatch),
    set_device_model!(template, PowerLoad, StaticPowerLoad),
    set_device_model!(template, InterruptibleLoad, DispatchablePowerLoad),
    set_device_model!(template, HydroEnergyReservoir, HydroDispatchReservoirStorage),
    return template
end

function get_template_hydro_st_ed(network = CopperPlatePowerModel)
    template = OperationsProblemTemplate(network)
    set_device_model!(template, ThermalStandard, ThermalDispatch)
    set_device_model!(template, RenewableDispatch, RenewableFullDispatch)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    set_device_model!(template, InterruptibleLoad, DispatchablePowerLoad)
    set_device_model!(template, HydroEnergyReservoir, HydroDispatchReservoirStorage)
    return template
end

function get_template_dispatch_with_network(network = StandardPTDFModel)
    template = OperationsProblemTemplate(network)
    set_device_model!(template, PowerLoad, StaticPowerLoad)
    set_device_model!(template, ThermalStandard, ThermalDispatch)
    set_device_model!(template, Line, StaticBranch)
    set_device_model!(template, Transformer2W, StaticBranch)
    set_device_model!(template, TapTransformer, StaticBranch)
    set_device_model!(template, HVDCLine, HVDCLossless)
    return template
end


text:
                #subfuncs to handle sums of log probabilities that may include -Inf (ie p=0), returning -Inf in this case rather than NaNs
                function lps(adjuvants)
                    prob = sum(adjuvants) ; isnan(prob) ? - Inf : prob
                end

                function lps(base, adjuvants...)
                    prob = base+sum(adjuvants) ; isnan(prob) ? -Inf : prob
                end


text:
using Documenter, Trace

makedocs(
    modules=[Trace],
    format=:html,
    pages=[
        "Home" => "index.md",
        "API" => "api.md",
    ],
    repo="https://github.com/invenia/Trace.jl/blob/{commit}{path}#L{line}",
    sitename="Trace.jl",
    authors="Invenia Technical Computing Corporation",
    assets=["assets/invenia.css"],
)

deploydocs(
    repo = "github.com/invenia/Trace.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
)


text:
using Test
using RendezvousHashing

using DataStructures: DefaultDict
using HypothesisTests: ChisqTest, pvalue
using Random: randstring

NODE_ID_WEIGHT_LIST = [
    ("node1", 100.0),
    ("node2", 200.0),
    ("node3", 300.0),
    ("node4", 400.0),
]

@testset "Smoke" begin
    coordinators = [
        WrhCoordinator(),
        WdhtCoordinator(),
    ]

    for coord in coordinators
        println("----------------")
        for (id, w) in NODE_ID_WEIGHT_LIST
            insert_node!(coord, id, w)
        end
        remove_node!(coord, "node4")

        println(choose_nodes(coord, "foo"))
        println(choose_nodes(coord, "bar", 2))
        println(choose_nodes(coord, "hello", 1))
    end
    println()
end

function compute_id_len(n::Int)::Int
    Int(ceil(log(n, 62))) + 1
end

function generate_nodes(n::Int)::Dict{String, Float64}
    id_len = compute_id_len(n)
    d = Dict{String, Float64}()
    while length(d) < n
        w = rand()
        d[randstring(id_len)] = w
    end
    w_sum = sum(values(d))
    for k in keys(d)
        d[k] /= w_sum
    end
    @assert sum(values(d)) ≈ 1.0
    d
end

@testset "Weighted Rendezvous Hashing" begin
    # `WdhtCoordinator` usually does NOT pass the Chi-squared test. When `N_TESTS` is large, both
    # will not pass the Chi-squared test unless `reseed!` is called regularly. The reason is
    # unknown. I have scrutinized my code but haven't found any obvious bugs.
    N_TESTS = 10000
    N_NODES = 5
    dataset = generate_nodes(N_NODES)
    coordinators = [
        WrhCoordinator(),
        WdhtCoordinator(),
    ]

    for coord in coordinators
        println("----------------")
        for (id, w) in dataset
            insert_node!(coord, id, w)
        end

        ID_LEN_TEST = compute_id_len(N_TESTS)
        counter = DefaultDict{String, Int}(0)
        for i in 1:N_TESTS
            data = randstring(ID_LEN_TEST)
            chosen = choose_nodes(coord, data, 1)[1]
            counter[chosen] += 1

            if i % 50 == 0
                reseed!(coord)
            end
        end

        for (id, freq) in counter
            println("Expected = $(dataset[id]), actual = $(freq / N_TESTS)")
        end

        observed = [counter[k] for k in keys(dataset)]
        expected = [dataset[k] for k in keys(dataset)]

        p = pvalue(ChisqTest(observed, expected))
        if p < 0.05
            println("(!) The implementation of $(typeof(coord)) is probably wrong, p = $p")
        end
    end
end

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
using Documenter
using Equate

makedocs(
    sitename = "Equate.jl",
    format = Documenter.HTML(
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://takuizum.github.io/Equate.jl/stable",
        assets=String[],
    ),
    modules = [Equate],
    pages = ["Home" => "index.md"]
)

deploydocs(;
    repo="github.com/takuizum/Equate.jl",
)


text:
@testset "504.base-7.jl" begin
    @test convert_to_base7(100) == "202"
    @test convert_to_base7(-7) == "-10"
    @test convert_to_base7(0) == "0"
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 9404760f-4015-46b1-8d2b-1bec68024dbb
using Colors, Manopt, Manifolds, PlutoUI

# ╔═╡ 0fabb176-c3e8-11ec-1304-f3d3346be36f
md"""
# Illustration of Jacobi Fields

This tutorial illustrates the usage of Jacobi Fields within
`Manopt.jl`.
For this tutorial you should be familiar with the basic terminology on a
manifold like the exponential and logarithmic map as well as
[shortest geodesic](https://juliamanifolds.github.io/Manifolds.jl/stable/interface.html#ManifoldsBase.shortest_geodesic-Tuple{Manifold,Any,Any})s.

We first initialize the packages we need
"""

# ╔═╡ cde64039-7fd7-4897-bc31-2e45b7ef18cd
md"and we define some colors from [Paul Tol](https://personal.sron.nl/~pault/)"

# ╔═╡ c8590a06-4825-4e4d-b3fa-7a82a8982063
begin
    black = RGBA{Float64}(colorant"#000000")
    TolVibrantMagenta = RGBA{Float64}(colorant"#EE3377")
    TolVibrantOrange = RGBA{Float64}(colorant"#EE7733")
    TolVibrantCyan = RGBA{Float64}(colorant"#33BBEE")
    TolVibrantTeal = RGBA{Float64}(colorant"#009988")
end;

# ╔═╡ 2028d1ab-3dcb-4f17-a02f-f9ef8c2dfacb
md"And setup our graphics paths"

# ╔═╡ 4e1fc0af-2837-40cc-b6ff-3d8d6d641b6b
begin
    localpath = join(splitpath(@__FILE__)[1:(end - 1)], "/") # files folder
    image_prefix = localpath * "/jacobi_fields"
    @info image_prefix
    render_asy = false # on CI or when you do not have asymptote, this should be false
end;

# ╔═╡ 17803b9a-5b91-4d16-ae52-06a83653dea2
md"""
Assume we have two points on the equator of the [Sphere](https://juliamanifolds.github.io/Manifolds.jl/stable/manifolds/sphere.html) $\mathcal M = \mathbb S^2$
"""

# ╔═╡ bab2db18-0596-4716-a6ce-8bfaf2b8ca24
M = Sphere(2)

# ╔═╡ 9e0e4098-d6c4-49dd-9564-961c73be88d5
p, q = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0]]

# ╔═╡ ee9a8c21-431e-4448-b758-01e552b1df8c
md"""
their connecting [shortest geodesic](https://juliamanifolds.github.io/Manifolds.jl/stable/interface.html#ManifoldsBase.shortest_geodesic-Tuple{Manifold,Any,Any}) (sampled at `100` points)
"""

# ╔═╡ 8a0d1039-5cfe-4e4d-9a97-89c3b5705b36
geodesic_curve = shortest_geodesic(M, p, q, [0:0.1:1.0...]);

# ╔═╡ 5e0d7686-7129-41f6-87d7-7298ddb131a5
md"Is just a curve along the equator"

# ╔═╡ b75c32b0-e322-42e5-a8fa-f74893655e18
render_asy && asymptote_export_S2_signals(
    image_prefix * "/jacobi_geodesic.asy";
    curves=[geodesic_curve],
    points=[[p, q]],
    colors=Dict(:curves => [black], :points => [TolVibrantOrange]),
    dot_size=3.5,
    line_width=0.75,
    camera_position=(1.0, 1.0, 0.5),
);

# ╔═╡ 7e3d4b24-ad82-4146-93c4-3bd4011ab151
render_asy && render_asymptote(image_prefix * "/jacobi_geodesic.asy"; render=2);

# ╔═╡ 34622df6-2979-4f05-8be9-224a2ec413f9
PlutoUI.LocalResource(image_prefix * "/jacobi_geodesic.png")

# ╔═╡ f0118459-5343-45dd-aa62-494795019044
md"""
where ``p`` is on the left and ``q`` on the right.
We know solve the following task:

Given a direction $X ∈ T_p\mathcal M$, for example
"""

# ╔═╡ a7b379ac-293e-454e-8ede-6f81bb7caa6e
X = [0.0, 0.4, 0.5]

# ╔═╡ f338d7d0-eef5-4264-af90-33e5bc6ecc62
md"""
we move the start point ``x`` into, how does any point on the geodesic move?

Or mathematically: Compute ``D_p g(t; p,q)`` for some fixed ``t∈[0,1]``
and a given direction ``X_p``.
Of course two cases are quite easy: For ``t=0`` we are in ``x`` and how ``x`` “moves”
is already known, so ``D_x g(0;p,q) = X``. On the other side, for ``t=1``,
``g(1; p,q) = q`` which is fixed, so ``D_p g(1; p,q)`` is the zero tangent vector
(in ``T_q\mathcal M``).

For all other cases we employ a [`jacobi_field`](https://manoptjl.org/stable/functions/Jacobi_fields.html), which is a (tangent)
vector field along the [shortest geodesic](https://juliamanifolds.github.io/Manifolds.jl/stable/interface.html#ManifoldsBase.shortest_geodesic-Tuple{Manifold,Any,Any}) given as follows: The _geodesic variation_
``\Gamma_{g,X}(s,t)`` is defined for some ``\varepsilon > 0`` as
````math
\Gamma_{g,X}(s,t):=\exp{\gamma_{p,X}(s)}[t\log_{g(s;p,X)}p],\qquad s∈(-\varepsilon,\varepsilon),\ t∈[0,1].
````
Intuitively we make a small step ``s`` into direction ``ξ`` using the geodesic
``g(⋅; p,X)`` and from ``r=g(s; p,X)`` we follow (in ``t``) the geodesic
``g(⋅; r,q)``. The corresponding Jacobi field ``J_{g,X}``
along ``g(⋅; p,q)`` is given by

````math
J_{g,X}(t):=\frac{D}{\partial s}\Gamma_{g,X}(s,t)\Bigl\rvert_{s=0}
````

which is an ODE and we know the boundary conditions ``J_{g,X}(0)=X`` and
``J_{g,X}(t) = 0``. In symmetric spaces we can compute the solution, since the
system of ODEs decouples, see for example [^doCarmo1992],
Chapter 4.2. Within `Manopt.jl` this is implemented as
`jacobi_field(M,p,q,t,X[,β])`, where the optional parameter (function)
``β`` specifies, which Jacobi field we want to evaluate and the one used here is
the default.

We can hence evaluate that on the points on the geodesic at
"""

# ╔═╡ b457c7c8-7f5c-4b92-a064-0f36b188f7a1
T = [0:0.1:1.0...];

# ╔═╡ c2ba61e5-7f74-4cb7-b68d-9057ab896859
md"namely"

# ╔═╡ 3e044a29-dff6-4174-9adf-f1228672b91d
r = shortest_geodesic(M, p, q, T);

# ╔═╡ c17abc34-6e5d-41d9-b428-b984bc8f7cd9
md"The tangent vector now moves as a differential along the geodesic as"

# ╔═╡ d7308356-100f-4766-9afe-588960b925dd
W = jacobi_field.(Ref(M), Ref(p), Ref(q), T, Ref(X));

# ╔═╡ fab5b16b-e05e-4a8a-b001-236422d9e245
md"""Which can also be called using `differential_geodesic_startpoint`.
We can add to the image above by creating extended tangent vectors
the include their base points"""

# ╔═╡ 0713167c-80c6-4907-b1e3-e104d598a804
V = [Tuple([a, b]) for (a, b) in zip(r, W)];

# ╔═╡ e26bb3fb-5f7e-470a-a2a1-7be37eb1bcb9
md"to add the vectors as one further set to the Asymptote export."

# ╔═╡ aa621823-57e3-4f53-8b65-1afbb80943bc
render_asy && asymptote_export_S2_signals(
    image_prefix * "/jacobi_geodesic_diff_start.asy";
    curves=[geodesic_curve],
    points=[[p, q], r],
    tangent_vectors=[V],
    colors=Dict(
        :curves => [black],
        :points => [TolVibrantOrange, TolVibrantCyan],
        :tvectors => [TolVibrantCyan],
    ),
    dot_sizes=[3.5, 2.0],
    line_width=0.75,
    camera_position=(1.0, 1.0, 0.5),
);

# ╔═╡ d02cff4d-8e88-4979-b536-6f367af8ab32
render_asy && render_asymptote(image_prefix * "/jacobi_geodesic_diff_start.asy"; render=2);

# ╔═╡ a17f81a9-066f-435c-b5ba-ae4809ac835b
PlutoUI.LocalResource(image_prefix * "/jacobi_geodesic_diff_start.png")

# ╔═╡ b2ac06aa-4b9f-4869-8b4c-c6f30661b8fe
md"""
The interpretation is as follows: If an infinitesimal change of the start point in direction ``X`` would happen, the infinitesimal change of a point along the line would change as indicated.

Note that each new vector is a tangent vector to its point (up to a small numerical tolerance), so the blue vectors are not just “shifted and scaled versions” of ``X``.
"""

# ╔═╡ fa6e15b8-2201-4095-8ebd-ddcd749fe11d
all([is_vector(M, a[1], a[2]; atol=1e-15) for a in V])

# ╔═╡ 99526567-14aa-4526-9a8d-48e2a3ba15a6
md"If we further move the end point, too, we can derive that Differential in direction"

# ╔═╡ 84090442-ea5f-4a1e-9c66-b0e828463853
begin
    Y = [0.2, 0.0, -0.5]
    W2 = differential_geodesic_endpoint.(Ref(M), Ref(p), Ref(q), T, Ref(Y))
    V2 = [Tuple([a, b]) for (a, b) in zip(r, W2)]
end;

# ╔═╡ bf3205fa-5348-42f9-8c1d-16e3fbcbc06c
md"and we can combine both vector fields we obtained"

# ╔═╡ c319b655-b3e5-4190-bf3c-d5af55da42ba
V3 = [Tuple([a, b]) for (a, b) in zip(r, W2 + W)];

# ╔═╡ 63b0d8ae-9a1b-46ff-978c-753f1cbdefe9
render_asy && asymptote_export_S2_signals(
    image_prefix * "/jacobi_geodesic_complete.asy";
    curves=[geodesic_curve],
    points=[[p, q], r],
    tangent_vectors=[V, V2, V3],
    colors=Dict(
        :curves => [black],
        :points => [TolVibrantOrange, TolVibrantCyan],
        :tvectors => [TolVibrantCyan, TolVibrantMagenta, TolVibrantTeal],
    ),
    dot_sizes=[3.5, 2.0],
    line_width=0.75,
    camera_position=(1.0, 1.0, 0.0),
);

# ╔═╡ 18278235-a049-4bd5-91dd-bddd31bcb8fb
render_asy && render_asymptote(image_prefix * "/jacobi_geodesic_complete.asy"; render=2);

# ╔═╡ cd1fa853-16ea-4975-b91b-5e6c1943523d
PlutoUI.LocalResource(image_prefix * "/jacobi_geodesic_complete.png")

# ╔═╡ ee803867-bc00-43b3-85c5-c4471370b6a8
md"""
Here the first vector field is still in blue, the second is in magenta, and their combined effect is in teal. Sure as a differential this does not make much sense, maybe as an infinitesimal movement of both start and end point cmobined.
"""

# ╔═╡ 83d85bac-85c7-462b-bf1c-097e27a3ca9a
md"""
## Literature

[^doCarmo1992]:
    > do Carmo, Manfredo __Riemannian Geometry__,  Birkhäuser Boston, 1992, ISBN: 0-8176-3490-8.
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
Manifolds = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
Manopt = "0fc0a36d-df90-57f3-8f93-d78a9fc72bb5"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
Colors = "~0.12.8"
Manifolds = "~0.7.8"
Manopt = "~0.3.23"
PlutoUI = "~0.7.38"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.2"
manifest_format = "2.0"

[[deps.AbstractFFTs]]
deps = ["ChainRulesCore", "LinearAlgebra"]
git-tree-sha1 = "6f1d9bc1c08f9f4a8fa92e3ea3cb50153a1b40d4"
uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
version = "1.1.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

[[deps.Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "af92965fb30777147966f58acb05da51c5616b5f"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.3"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.ArnoldiMethod]]
deps = ["LinearAlgebra", "Random", "StaticArrays"]
git-tree-sha1 = "62e51b39331de8911e4a7ff6f5aaf38a5f4cc0ae"
uuid = "ec485272-7323-5ecc-a04f-4719b315124d"
version = "0.2.0"

[[deps.ArrayInterface]]
deps = ["Compat", "IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"]
git-tree-sha1 = "c933ce606f6535a7c7b98e1d86d5d1014f730596"
uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
version = "5.0.7"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "9950387274246d08af38f6eef8cb5480862a435f"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.14.0"

[[deps.ChangesOfVariables]]
deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
git-tree-sha1 = "bf98fa45a0a4cee295de98d4c1462be26345b9a1"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.2"

[[deps.ColorSchemes]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random"]
git-tree-sha1 = "12fc73e5e0af68ad3137b886e3f7c1eacfca2640"
uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
version = "3.17.1"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "b153278a25dd42c65abbf4e62344f9d22e59191b"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.43.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.CovarianceEstimation]]
deps = ["LinearAlgebra", "Statistics", "StatsBase"]
git-tree-sha1 = "a3e070133acab996660d31dcf479ea42849e368f"
uuid = "587fd27a-f159-11e8-2dae-1979310e6154"
version = "0.2.7"

[[deps.DataAPI]]
git-tree-sha1 = "fb5f5316dd3fd4c5e7c30a24d50643b73e37cd40"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.10.0"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.11"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.DensityInterface]]
deps = ["InverseFunctions", "Test"]
git-tree-sha1 = "80c3e8639e3353e5d2912fb3a1916b8455e2494b"
uuid = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
version = "0.4.0"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.Distributions]]
deps = ["ChainRulesCore", "DensityInterface", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns", "Test"]
git-tree-sha1 = "5a4168170ede913a2cd679e53c2123cb4b889795"
uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
version = "0.25.53"

[[deps.DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.Einsum]]
deps = ["Compat"]
git-tree-sha1 = "4a6b3eee0161c89700b6c1949feae8b851da5494"
uuid = "b7d42ee7-0b51-5a75-98ca-779d3107e4c0"
version = "0.4.1"

[[deps.FillArrays]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
git-tree-sha1 = "246621d23d1f43e3b9c368bf3b72b2331a27c286"
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
version = "0.13.2"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Graphs]]
deps = ["ArnoldiMethod", "Compat", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"]
git-tree-sha1 = "57c021de207e234108a6f1454003120a1bf350c4"
uuid = "86223c79-3864-5bf0-83f7-82e725a168b6"
version = "1.6.0"

[[deps.HybridArrays]]
deps = ["LinearAlgebra", "Requires", "StaticArrays"]
git-tree-sha1 = "eb6b23460f5544c5d09efae0818b86736cefcd3d"
uuid = "1baab800-613f-4b0a-84e4-9cd3431bfbb9"
version = "0.4.10"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.IfElse]]
git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1"
uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
version = "0.1.1"

[[deps.Inflate]]
git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c"
uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
version = "0.1.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "91b5dcf362c5add98049e6c29ee756910b03051d"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.3"

[[deps.IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

[[deps.JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.4.1"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.3"

[[deps.Kronecker]]
deps = ["LinearAlgebra", "NamedDims", "SparseArrays", "StatsBase"]
git-tree-sha1 = "a51f46415c844dee694cb8b20a3fcbe6dba342c2"
uuid = "2c470bb0-bcc8-11e8-3dad-c9649493f05e"
version = "0.5.1"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.LogExpFunctions]]
deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "a970d55c2ad8084ca317a4658ba6ce99b7523571"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.12"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.9"

[[deps.Manifolds]]
deps = ["Colors", "Distributions", "Einsum", "Graphs", "HybridArrays", "Kronecker", "LinearAlgebra", "ManifoldsBase", "Markdown", "Random", "RecipesBase", "RecursiveArrayTools", "Requires", "SimpleWeightedGraphs", "SpecialFunctions", "StaticArrays", "Statistics", "StatsBase"]
git-tree-sha1 = "cc5e34ecd061bf9bedd62a8b8a28abe96c3b7c6e"
uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
version = "0.7.8"

[[deps.ManifoldsBase]]
deps = ["LinearAlgebra", "Markdown"]
git-tree-sha1 = "b8518cbf0f06f1375c95479759b266bd75d89ad9"
uuid = "3362f125-f0bb-47a3-aa74-596ffd7ef2fb"
version = "0.12.13"

[[deps.Manopt]]
deps = ["ColorSchemes", "ColorTypes", "Colors", "DataStructures", "Dates", "LinearAlgebra", "ManifoldsBase", "Markdown", "Printf", "Random", "Requires", "SparseArrays", "StaticArrays", "Statistics", "Test"]
git-tree-sha1 = "a8a292a422b0a6924f46d4e022fc8263721120c5"
uuid = "0fc0a36d-df90-57f3-8f93-d78a9fc72bb5"
version = "0.3.23"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NamedDims]]
deps = ["AbstractFFTs", "ChainRulesCore", "CovarianceEstimation", "LinearAlgebra", "Pkg", "Requires", "Statistics"]
git-tree-sha1 = "0856b62716585eb90cc1dada226ac9eab5f69aa5"
uuid = "356022a1-0364-5f58-8944-0da4b18d706f"
version = "0.2.47"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"

[[deps.OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.PDMats]]
deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
git-tree-sha1 = "e8185b83b9fc56eb6456200e873ce598ebc7f262"
uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
version = "0.11.7"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "3b429f37de37f1fc603cc1de4a799dc7fbe4c0b6"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.3.0"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "670e559e5c8e191ded66fa9ea89c97f10376bb4c"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.38"

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "47e5f437cc0e7ef2ce8406ce1e7e24d44915f88d"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.3.0"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.QuadGK]]
deps = ["DataStructures", "LinearAlgebra"]
git-tree-sha1 = "78aadffb3efd2155af139781b8a8df1ef279ea39"
uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
version = "2.4.2"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.RecipesBase]]
git-tree-sha1 = "6bf3f380ff52ce0832ddd3a2a7b9538ed1bcca7d"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.2.1"

[[deps.RecursiveArrayTools]]
deps = ["Adapt", "ArrayInterface", "ChainRulesCore", "DocStringExtensions", "FillArrays", "LinearAlgebra", "RecipesBase", "Requires", "StaticArrays", "Statistics", "ZygoteRules"]
git-tree-sha1 = "bfe14f127f3e7def02a6c2b1940b39d0dabaa3ef"
uuid = "731186ca-8d62-57ce-b412-fbd966d074cd"
version = "2.26.3"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.3.0"

[[deps.Rmath]]
deps = ["Random", "Rmath_jll"]
git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f"
uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
version = "0.7.0"

[[deps.Rmath_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7"
uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
version = "0.3.0+0"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.SimpleTraits]]
deps = ["InteractiveUtils", "MacroTools"]
git-tree-sha1 = "5d7e3f4e11935503d3ecaf7186eac40602e7d231"
uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
version = "0.9.4"

[[deps.SimpleWeightedGraphs]]
deps = ["Graphs", "LinearAlgebra", "Markdown", "SparseArrays", "Test"]
git-tree-sha1 = "a6f404cc44d3d3b28c793ec0eb59af709d827e4e"
uuid = "47aef6b3-ad0c-573a-a1e2-d07658019622"
version = "1.2.1"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.SpecialFunctions]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "5ba658aeecaaf96923dce0da9e703bd1fe7666f9"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "2.1.4"

[[deps.Static]]
deps = ["IfElse"]
git-tree-sha1 = "87e9954dfa33fd145694e42337bdd3d5b07021a6"
uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
version = "0.6.0"

[[deps.StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "cd56bf18ed715e8b09f06ef8c6b781e6cdc49911"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.4.4"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.StatsAPI]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "c82aaa13b44ea00134f8c9c89819477bd3986ecd"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.3.0"

[[deps.StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "8977b17906b0a1cc74ab2e3a05faa16cf08a8291"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.16"

[[deps.StatsFuns]]
deps = ["ChainRulesCore", "InverseFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
git-tree-sha1 = "5950925ff997ed6fb3e985dcce8eb1ba42a0bbe7"
uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
version = "0.9.18"

[[deps.SuiteSparse]]
deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.ZygoteRules]]
deps = ["MacroTools"]
git-tree-sha1 = "8c1a8e4dfacb1fd631745552c8db35d0deb09ea0"
uuid = "700de1a5-db45-46bc-99cf-38207098b444"
version = "0.2.2"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─0fabb176-c3e8-11ec-1304-f3d3346be36f
# ╠═9404760f-4015-46b1-8d2b-1bec68024dbb
# ╟─cde64039-7fd7-4897-bc31-2e45b7ef18cd
# ╠═c8590a06-4825-4e4d-b3fa-7a82a8982063
# ╟─2028d1ab-3dcb-4f17-a02f-f9ef8c2dfacb
# ╠═4e1fc0af-2837-40cc-b6ff-3d8d6d641b6b
# ╟─17803b9a-5b91-4d16-ae52-06a83653dea2
# ╠═bab2db18-0596-4716-a6ce-8bfaf2b8ca24
# ╠═9e0e4098-d6c4-49dd-9564-961c73be88d5
# ╟─ee9a8c21-431e-4448-b758-01e552b1df8c
# ╠═8a0d1039-5cfe-4e4d-9a97-89c3b5705b36
# ╟─5e0d7686-7129-41f6-87d7-7298ddb131a5
# ╠═b75c32b0-e322-42e5-a8fa-f74893655e18
# ╠═7e3d4b24-ad82-4146-93c4-3bd4011ab151
# ╠═34622df6-2979-4f05-8be9-224a2ec413f9
# ╟─f0118459-5343-45dd-aa62-494795019044
# ╟─a7b379ac-293e-454e-8ede-6f81bb7caa6e
# ╟─f338d7d0-eef5-4264-af90-33e5bc6ecc62
# ╠═b457c7c8-7f5c-4b92-a064-0f36b188f7a1
# ╟─c2ba61e5-7f74-4cb7-b68d-9057ab896859
# ╠═3e044a29-dff6-4174-9adf-f1228672b91d
# ╟─c17abc34-6e5d-41d9-b428-b984bc8f7cd9
# ╠═d7308356-100f-4766-9afe-588960b925dd
# ╟─fab5b16b-e05e-4a8a-b001-236422d9e245
# ╠═0713167c-80c6-4907-b1e3-e104d598a804
# ╟─e26bb3fb-5f7e-470a-a2a1-7be37eb1bcb9
# ╠═aa621823-57e3-4f53-8b65-1afbb80943bc
# ╠═d02cff4d-8e88-4979-b536-6f367af8ab32
# ╟─a17f81a9-066f-435c-b5ba-ae4809ac835b
# ╟─b2ac06aa-4b9f-4869-8b4c-c6f30661b8fe
# ╠═fa6e15b8-2201-4095-8ebd-ddcd749fe11d
# ╟─99526567-14aa-4526-9a8d-48e2a3ba15a6
# ╠═84090442-ea5f-4a1e-9c66-b0e828463853
# ╠═bf3205fa-5348-42f9-8c1d-16e3fbcbc06c
# ╠═c319b655-b3e5-4190-bf3c-d5af55da42ba
# ╠═63b0d8ae-9a1b-46ff-978c-753f1cbdefe9
# ╠═18278235-a049-4bd5-91dd-bddd31bcb8fb
# ╠═cd1fa853-16ea-4975-b91b-5e6c1943523d
# ╟─ee803867-bc00-43b3-85c5-c4471370b6a8
# ╟─83d85bac-85c7-462b-bf1c-097e27a3ca9a
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
################################################################################
#
# RBC perturbation.jl solves the RBC model using Schmitt-Grohé & Uribe (2004) method
#
# Shane McMiken
# June 2021
# Boston College
#
################################################################################

#--------------------------------#
#         PREAMBLE
#--------------------------------#

using Plots
using ForwardDiff
using LinearAlgebra # schur(A::StridedMatrix, B::StridedMatrix) -> F::GeneralizedSchur

#--------------------------------#
#         RBC MODEL
#--------------------------------#

# 1 - β [ (c^η/cp^η) (1 + rp - δ) ]
# w - l^ψ c^η
# (1-δ) k + i - kp
# c + i - y
# y - a * k^α * l^(1-α)
# w - (1-α)*(y/l)
# r - α * (y/k)
# log(ap) - ρ*log(a)

# parameters
η   = 0.8   # relative risk aversion
ψ   = 0.2   # inverse Frisch elasticity
β   = 0.95  # discount rate
α   = 0.3   # capital share of income
δ   = 0.05  # depreciation rate
ρ   = 0.9   # persistence parameter
σ   = 0.02  # Standard deviation of TFP shocks

# steady state - levels
abar     = 1.0
rbar     = β^(-1)-(1-δ)
klbar    = (rbar/α)^(1/(α-1))
clbar    = abar*klbar^α-δ*klbar
lbar     = ((1-α)*abar*klbar^α/(abar*klbar^α-δ*klbar)^η)^(1/(ψ+η))
kbar     = klbar*lbar
cbar     = clbar*lbar
wbar     = (1-α)*abar*klbar^α
ybar     = kbar^α*lbar^(1-α)
ibar     = δ*kbar

# Define vectors of controls y, states x
Yss = [ybar, cbar, lbar, wbar, rbar, ibar]
Xss = [kbar, abar]

ny  = length(Yss) # Number of control variables
nx  = length(Xss) # Number of state variables
nyp = length(Yss)
nxp = length(Xss)

# take logs
lYss = log.(Yss)
lXss = log.(Xss)

xs  = vcat(lYss,lXss,lYss,lXss)

#--------------------------------#
#       MODEL EQUATIONS
#--------------------------------#

# 8 equations & 8 unknowns
#
#
# y c l w r i k a
# 1 2 3 4 5 6 7 8
#
# yp cp lp wp rp ip kp ap
# 9  10 11 12 13 14 15 16
#
# 1. c^-η - β [ cp^-η (1 + rp - δ) ]
# 2. l^ψ c^η - w
# 3. (1-δ) k + i - kp
# 4. c + i - y
# 5. y - a * k^α * l^(1-α)
# 6. w - (1-α)*(k/l)^α
# 7. r - α * (l/k)^(1-α)
# 8. log(ap) - ρ*log(a)

f1(x::AbstractArray) = exp(x[2])^(-η) - β*(exp(x[10])^(-η)*(1 + exp(x[13]) - δ))
f2(x::AbstractArray) = exp(x[4]) - exp(x[3])^ψ * exp(x[2])^η
f3(x::AbstractArray) = (1-δ)*exp(x[7]) + exp(x[6]) - exp(x[15])
f4(x::AbstractArray) = exp(x[2]) + exp(x[6]) - exp(x[1])
f5(x::AbstractArray) = exp(x[1]) - exp(x[8]) * exp(x[7])^α * exp(x[3])^(1-α)
f6(x::AbstractArray) = exp(x[4]) - (1-α)*(exp(x[7])/exp(x[3]))^α
f7(x::AbstractArray) = exp(x[5]) - α*(exp(x[3])/exp(x[7]))^(1-α)
f8(x::AbstractArray) = log(exp(x[16])) - ρ*log(exp(x[8]))

f = [f1;f2;f3;f4;f5;f6;f7;f8]

#--------------------------------#
#   GET NUMERICAL DERIVATIVES
#--------------------------------#

function Numerical_derivatives(f::Array{Function},y::Array,x::Array,yp::Array,xp::Array; ny=ny, nx=nx, nyp=nyp, nxp=nxp, xs=xs)

    # compute jacobian
    J   = zeros(length(f),length(xs))
    for (ii,func) in enumerate(f)
        J[ii,:] = ForwardDiff.gradient(func,xs)
    end

    # compute first order matices
    fy  = J[:,1:ny]
    fx  = J[:,(1+ny):(nx+ny)]
    fyp = J[:,(1+nx+ny):(nx+ny+nyp)]
    fxp = J[:,(1+nx+ny+nyp):(nx+ny+nyp+nxp)]

    return fy,fx,fyp,fxp
end

@time fy,fx,fyp,fxp = Numerical_derivatives(f,lYss,lXss,lYss,lXss)

#--------------------------------#
#   FIRST ORDER APPROXIMATION
#--------------------------------#

# compute the gx, hx matrix
function gx_hx(fy::Array,fx::Array,fyp::Array,fxp::Array; ny=ny, nx=nx, nyp=nyp, nxp=nxp)

    # Create system matrices A,B
    A   = [-fxp -fyp]
    B   = [fx fy]
    NK  = size(fx,2)

    # Complex Schur decomposition
    F = schur(A,B)

    # Pick non-explosive (stable) eigenvalues
    slt = abs.(diag(F.T)) .< abs.(diag(F.S))
    nk = sum(slt)

    # Reorder the system with stable eigs in upper-left
    S,T,Q,Z = ordschur(F,slt)

    # Split up results appropriately
    Z21 = Z[nk+1:end,1:nk]
    Z11 = Z[1:nk,1:nk]
    S11 = S[1:nk,1:nk]
    T11 = T[1:nk,1:nk]

    # Identify cases with no or multiple solutions
    if nk > NK
        print("Warning gx_hx: the equilibrium is locally indeterminate")
    elseif nk < NK
        print("Warning gx_hx: no local equilibrium exists")
    end

    # Check invertibility
    if rank(Z11)<nk
        print("Invertibility condition violated")
    end

    # Compute solution
    Z11i    = Z11\I(nk)
    gx      = real(Z21*Z11i)
    hx      = real(Z11*(S11\T11)*Z11i)

    return gx, hx
end

@time gx, hx = gx_hx(fy,fx,fyp,fxp)

#--------------------------------#
#   IMPULSE RESPONSES
#--------------------------------#

T   = 50
eta = [0,σ]

function ir(gx::Array,hx::Array,x0,T::Int)

  #x0=x0(:)
  pd=size(x0,1)
  MX=[gx;I(pd)]
  IR=zeros(T,size(MX,1))
  x=x0

  for t=1:T
  IR[t,:]=(MX*x)'
  x = hx * x
  end

  iry = IR[:,1:ny]
  irx = IR[:,ny+1:end]

  return iry,irx
end

irfy,irfx = ir(gx,hx,eta,T)

# Output
iry1 = plot(1:T,irfy[:,1], title="Output", label="IRF")
# Consumption
iry2 = plot(1:T,irfy[:,2], title="Consumption", label="IRF")
# Hours
iry3 = plot(1:T,irfy[:,3], title="Hours", label="IRF")
# Wages
iry4 = plot(1:T,irfy[:,4], title="Wages", label="IRF")
# Investment
iry5 = plot(1:T,irfy[:,5], title="Investment", label="IRF")
# Real Rate
iry6 = plot(1:T,irfy[:,6], title="Real Rate", label="IRF")

# Capital
irx1 = plot(1:T,irfx[:,1], title="Capital", label="IRF")
# TFP
irx2 = plot(1:T,irfx[:,2], title="Technology", label="IRF")

impulseresponses = plot(iry1,iry2,iry3,iry4,iry5,iry6,irx1,irx2, layout = (4,2), legend = false)
display(impulseresponses)
savefig(impulseresponses,"./impulseresponses.png")


text:
# Since it takes so long time to load the Embeddings we will instead cache
# calculated word distances so we only need to add the embeddings lazily, i.e.
# when the distance is needed to a new word.
vocabulary2indices(embtable) = Dict(word=>ii for (ii,word) in enumerate(embtable.vocab))

const DefaultWordCacheFilename = ".word_distances_cache.json"
const DefaultEmbeddingSystem = :word2vec
const DefaultEmbeddingFilenumber = 1

# We use our own mapping from embedding system names so users need not
# learn about the insides of the Embeddings package.
const GloVeLangFiles = language_files(GloVe{:en})
gloveindex(filename) = findfirst(f -> f == filename, GloVeLangFiles)

const EmbeddingName2Embedding = Dict(
    :word2vec      => (Embeddings.Word2Vec, 1),
    :glove         => (Embeddings.GloVe{:en}, gloveindex("glove.6B/glove.6B.300d.txt")),

    :glove6b50d    => (Embeddings.GloVe{:en}, gloveindex("glove.6B/glove.6B.50d.txt")),
    :glove6b100d   => (Embeddings.GloVe{:en}, gloveindex("glove.6B/glove.6B.100d.txt")),
    :glove6b200d   => (Embeddings.GloVe{:en}, gloveindex("glove.6B/glove.6B.200d.txt")),
    :glove6b300d   => (Embeddings.GloVe{:en}, gloveindex("glove.6B/glove.6B.300d.txt")),

    :glove42b300d  => (Embeddings.GloVe{:en}, gloveindex("glove.42B.300d/glove.42B.300d.txt")),
    :glove840b300d => (Embeddings.GloVe{:en}, gloveindex("glove.840B.300d/glove.840B.300d.txt")),
)

embeddingname2embedding(name) = lookup_embeddingname(Symbol(lowercase(string(name))))

current_embeddings() = collect(keys(EmbeddingName2Embedding))

function lookup_embeddingname(name::Symbol)
    try
        return EmbeddingName2Embedding[name]
    catch err
        allembeddings = join(map(string, collect(keys(EmbeddingName2Embedding))), ", ")
        error("Cannot find an embedding for $name.\nAvailable embeddings: $allembeddings")
    end
end

# We cache embeddings and their word indices also so we need not load
# and calc them multiple times.
const EmbeddingsCache = Dict()

function cache_load_embeddings(embsys, embfilenumber)
    if !haskey(EmbeddingsCache, (embsys, embfilenumber))
        println("Loading the embeddings for ($embsys, $embfilenumber). This often can take 20-60 seconds...")
        embtable = load_embeddings(embsys, embfilenumber)
        wordindex = vocabulary2indices(embtable)
        EmbeddingsCache[(embsys, embfilenumber)] = (embtable, wordindex)
    end
    EmbeddingsCache[(embsys, embfilenumber)]
end

mutable struct WordDistanceCache
    dirname::String
    cachefile::String
    embeddingsystem::Symbol
    embtable::Union{Embeddings.EmbeddingTable, Nothing}
    wordindex::Dict{String, Int}
    worddistances::Dict{Tuple{String,String}, Float64}

    WordDistanceCache(dir::String, embeddingsystem::Symbol;
        cachefile = DefaultWordCacheFilename) = begin
        new(dir, cachefile, embeddingsystem,
            nothing,
            Dict{String, Int}(),
            Dict{Tuple{String,String}, Float64}())
    end
end

function WordDistanceCache(cachefilepath::String = joinpath(".", DefaultWordCacheFilename))
    if isfile(cachefilepath)
        load_word_distance_cache(cachefilepath)
    else
        WordDistanceCache(dirname(cachefilepath), DefaultEmbeddingSystem)
    end
end

function hasembedding(wdc::WordDistanceCache, word::AbstractString)
    wdc.embtable != nothing || load_embeddings!(wdc)
    haskey(wdc.wordindex, string(word))
end

function load_embeddings!(wdc::WordDistanceCache)
    if wdc.embtable == nothing
        embsys, embfilenum = embeddingname2embedding(wdc.embeddingsystem)
        wdc.embtable, wdc.wordindex = cache_load_embeddings(embsys, embfilenum)
    end
end

function worddistance(wdc::WordDistanceCache, w1::String, w2::String)
    get!(wdc.worddistances, (w1, w2)) do
        calcworddistance(wdc, w1, w2)
    end
end

const EmbeddingDistance = Euclidean()

function calcworddistance(wdc::WordDistanceCache, w1::String, w2::String)
    evaluate(EmbeddingDistance, embedding(wdc, w1), embedding(wdc, w2))
end

function embedding(wdc::WordDistanceCache, w::String)
    wdc.embtable != nothing || load_embeddings!(wdc)
    try
        ind = wdc.wordindex[w]
        return wdc.embtable.embeddings[:,ind]
    catch err
        # Very strange but "and" is not available in Word2Vec while "And" is. Maybe
        # we want to fix this?
        error("No embedding found for word $w")
    end
end

# This is for when we load from json only...
function add_distance!(wdc::WordDistanceCache, w1::String, w2::String, d::Number)
    wdc.worddistances[(w1, w2)] = Float64(d)
end

# When saving to json we must output distances in an array so we don't
# lose its structure. JSON cannot handle Julia tuples.
function distances_as_array_for_json(wdc::WordDistanceCache)
    Any[Any[w1, w2, d] for ((w1, w2), d) in wdc.worddistances]
end

import JSON.json
function json(wdc::WordDistanceCache)
    metadict = Dict(
        "embeddingsystem" => string(wdc.embeddingsystem),
        "lastsave" => Libc.strftime("%Y-%m-%d %H:%M.%S", time())
    )
    js = json(Dict("meta" => metadict, "distances" => distances_as_array_for_json(wdc)), 2)
    return js
end

function save(wdc::WordDistanceCache)
    fn = joinpath(wdc.dirname, wdc.cachefile)
    open(fn, "w") do fh
        println(fh, json(wdc))
    end
    println("Cached word distances to file $fn")
end

function load_word_distance_cache(cachefilepath::String)
    d = JSON.parse(read(cachefilepath, String))
    println("Reading word distance cache that was last saved: ", d["meta"]["lastsave"])

    cachefile = last(split(cachefilepath, "/"))
    embeddingsystem = Symbol(get(d["meta"], "embeddingsystem", "Word2Vec"))
    wdc = WordDistanceCache(dirname(cachefilepath), embeddingsystem;
        cachefile = cachefile)

    # In json the distances are saved as arrays [word1, word2, dist] so we add them
    # manually based on that.
    for entry in d["distances"]
        add_distance!(wdc, entry[1], entry[2], entry[3])
    end

    return wdc
end


text:
# Important, for broadcast joins we cannot assume c and d have same number of columns:
# c could have more columns than d
@generated function rowcmp(c::StructVector, i, d::StructVector{C, D}, j) where {C, D}
    N = fieldcount(D)
    ex = :(cmp(getfield(fieldarrays(c),$N)[i], getfield(fieldarrays(d),$N)[j]))
    for n in N-1:-1:1
        ex = quote
            let k = rowcmp(getfield(fieldarrays(c),$n), i, getfield(fieldarrays(d),$n), j)
                (k == 0) ? ($ex) : k
            end
        end
    end
    ex
end

@inline function rowcmp(c::AbstractVector, i, d::AbstractVector, j)
    cmp(c[i], d[j])
end

struct GroupJoinPerm{LP<:GroupPerm, RP<:GroupPerm}
    left::LP
    right::RP
end

GroupJoinPerm(lkeys::AbstractVector, rkeys::AbstractVector, lperm=sortperm(lkeys), rperm=sortperm(rkeys)) =
    GroupJoinPerm(GroupPerm(lkeys, lperm), GroupPerm(rkeys, rperm))

function _pick(s, a, b)
    if a === nothing && b === nothing
        return nothing
    elseif a === nothing
        return (1:0, b[1]), (1, a, b)
    elseif b === nothing
        return (a[1], 1:0), (-1, a, b)
    else
        lp = sortperm(s.left)
        rp = sortperm(s.right)
        cmp = rowcmp(parent(s.left), lp[first(a[1])], parent(s.right), rp[first(b[1])])
        if cmp < 0
            return (a[1], 1:0), (-1, a, b)
        elseif cmp == 0
            return (a[1], b[1]), (0, a, b)
        else
            return (1:0, b[1]), (1, a, b)
        end
    end
end

function Base.iterate(s::GroupJoinPerm)
    l = iterate(s.left)
    r = iterate(s.right)
    _pick(s, l, r)
end

function Base.iterate(s::GroupJoinPerm, (select, l, r))
    (select <= 0) && (l = iterate(s.left, l[2]))
    (select >= 0) && (r = iterate(s.right, r[2]))
    _pick(s, l, r)
end

Base.IteratorSize(::Type{<:GroupJoinPerm}) = Base.SizeUnknown()
Base.IteratorEltype(::Type{<:GroupJoinPerm}) = Base.HasEltype()
Base.eltype(::Type{<:GroupJoinPerm}) = Tuple{UnitRange{Int}, UnitRange{Int}}


text:
using Weave
using Test

cleanup = true

VER = "$(VERSION.major).$(VERSION.minor)"

weave("documents/chunk_options.noweb")
result =  read("documents/chunk_options.md", String)
ref =  read("documents/chunk_options_ref.md", String)
@test result == ref
cleanup && rm("documents/chunk_options.md")

tangle("documents/chunk_options.noweb", out_path = "documents/tangle")
result =  read("documents/tangle/chunk_options.jl", String)
ref =  read("documents/tangle/chunk_options.jl.ref", String)
@test ref == result
cleanup && rm("documents/tangle/chunk_options.jl")


text:
#=
    Mesh Types to reproduce the structure of an OpenFOAM (poly) mesh
=#

abstract type AbstractOfMesh end

struct OfPoints <: AbstractOfMesh
    no::Int64
    x::Vector{Float64}
    y::Vector{Float64}
    z::Vector{Float64}
end

struct OfFace <: AbstractOfMesh
    no_points::Int64
    point_list::Vector{Int64}
end

struct OfFaces <: AbstractOfMesh
    no::Int64
    face_list::Vector{OfFace}
end

struct OfPatch <: AbstractOfMesh
    no_faces::Int64
    start_face::Int64
    type::String # Wall, Patch, Empty
end

mutable struct OfMeshInfo <: AbstractOfMesh
    version::Float64
    format::String
    class::String
    location::String
    casepath::String
end

struct OfMesh <: AbstractOfMesh
    info::OfMeshInfo
    points::OfPoints
    faces::OfFaces
    patches::Vector{OfPatch}
    face_owner_cell::Vector{Int64}
    face_neighbor_cell::Vector{Int64}
end

max_stars_repo_name: c-schubert/linearMeshIO

max_issues_repo_name: c-schubert/linearMeshIO

max_forks_repo_name: c-schubert/linearMeshIO


text:
module TraceGraph

using Cassette
using LightGraphs
import LightGraphs: add_vertex!, add_edge!
using DataStructures
using InteractiveUtils

"""
    struct TGraph
Object that represents the traced graph of the function call.
fields:
- g::SimpleDiGraph  : Object of type SimpleDiGraph
- nodelabel::Vector{String} : List of names corresponding to the nodes
- nodeop::Vector{String} : List of names of operations of the nodes
- nodeval::Vector{Any} : List of values of ndoes
"""
struct TGraph
    g::SimpleDiGraph
    nodelabel::Vector{String}
    nodeop::Vector{String}
    nodevalue::Vector{Any}
end

include("NoRecurseList.jl")
include("GenerateGraph.jl")

export TGraph, tracegraph
export add_norecurse, rm_norecurse, show_norecurse

end # module


text:
import FinanceLib.FixedIncomes.Bonds as Bd
import FinanceLib.FixedIncomes.Bonds.Rates as Rt
import FinanceLib as Fl

@testset "FinanceLib.FixedIncomes.Bonds                               " begin

    # @test Bd.priceCouponBond(Fl.PeriodSeries([(1, 0.09), (2, 0.10)]), 0.06) * 1000 == 931.0789294108725
    # @test Bd.priceCouponBond(Fl.PeriodSeries([(1, 0.09), (2, 0.10), (3, 0.11)]), 0.05) * 100 == 85.49448240486119

    @test Bd.priceCouponBond(0.03, Bd.CouponBond(0.05, 2, 3)) == 1.0569718716547531

    @test Bd.priceCouponBond(Fl.RateCurve{Fl.NomRate}([0.0016, 0.0021, 0.0027, 0.0033, 0.0037, 0.0041],2), Bd.CouponBond(0.05, 2, 3)) == 1.1369147941993403

    @test Bd.ytmCouponBond(Bd.CouponBond(0.05, 2, 3), 1.1369147941993403) == 0.004038639185260929

    # rS = Fl.PeriodSeries([(1, 0.05), (2, 0.06), (3, 0.07), (4, 0.08), (5, 0.09)])
    # bd = Bd.priceCouponBond(rS, 0.1)*100
    # @test bd == 105.42950371364294
    # @test Bd.ytmCouponBond(rS, 0.1) == 0.0861792633877846

    @testset "FinanceLib.FixedIncomes.Bonds.Rates                     " begin

      @test Rt.parToSpotRates(Fl.RateCurve{Fl.NomRate}([0.05, 0.0597, 0.0691, 0.0781], 1)).rate[4] ≈ 0.08002666689660898

      er = Rt.parToSpotRates(Fl.RateCurve{Fl.NomRate}([0.020000, 0.024000, 0.027600, 0.030840, 0.033756, 0.036380], 2)).rate

      er[1] ≈ 0.02
      er[4] ≈ 0.030974
      er[5] ≈ 0.033975
      er[6] ≈ 0.036701

      et = Rt.spotToParRates(Fl.RateCurve{Fl.NomRate}(er, 2)).rate

      @test et[3] ≈ 0.027600
      @test et[4] ≈ 0.030840
      @test et[6] ≈ 0.036380

    end

end


text:
using Gadfly, DataArrays, RDatasets

plot(dataset("car", "UN"), x="GDP", y="InfantMortality",
     Geom.histogram2d, Scale.x_log10, Scale.y_log10)


text:
include(find_in_path("PatternDispatch.jl"))

module TestImmutable
using PatternDispatch.Immutable

@immutable type T
    x::Int
end
@immutable type S
    x
end
@immutable type XY
    x
    y
end

@assert !(T(7)::T === T(42)::T)
@assert T(42)::T === T(42)::T

@assert !(S(7)::S === S(42)::S)
@assert S(42)::S === S(42)::S
@assert S(42)::S === S(42.0)::S # current semantics, subject to change...

@assert !(S(42)::S === T(42)::T)

@assert !(XY(1,2)::XY === XY(3,4)::XY)
@assert !(XY(1,2)::XY === XY(1,4)::XY)
@assert !(XY(1,2)::XY === XY(3,2)::XY)
@assert (XY(1,2)::XY === XY(1,2)::XY)
@assert (XY(3,4)::XY === XY(3,4)::XY)

end


text:
using Documenter

makedocs(
    format = Documenter.HTML(),
    sitename = "JuMPTutorials.jl",
    authors  = "JuMP-dev",
    pages = [
        "Home" => "index.md",
        "Adding a new tutorial" => "new.md",
        "Notes" => "notes.md",
        "Future work" => "future.md",
    ],
)

Documenter.deploydocs(repo = "github.com/jump-dev/JuMPTutorials.jl.git")

const converter_path = joinpath(@__DIR__, "../converter")

# notebook deployment phase
using Pkg
Pkg.activate(converter_path)
Pkg.instantiate()

include(joinpath(converter_path, "convert_pages.jl"))

Documenter.deploydocs(
    repo = "github.com/jump-dev/JuMPTutorials.jl.git",
    target = "../notebook",
    branch = "notebook",
)


text:
#-----------------------------------------------------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------------------------------------------------
#
#   CALCULATIONS OF BAND STRUCTURE
#
#-----------------------------------------------------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------------------------------------------------

#-----------------------------------------------------------------------------------------------------------------------------
#
#   METHODS FOR CALCULATING THE BAND STRUCTURE OF A MATRIX IN K SPACE
#
#-----------------------------------------------------------------------------------------------------------------------------

#-----------------------------------------------------------------------------------------------------------------------------
#
#   BAND STRUCTURE ALONG A PATH IN K SPACE
#
#   path has to be of following format
#   Array[
#       ["name1", [coordinates1]],
#       ["name2", [coordinates2]],
#       ...
#   ]
#   path does not close! to close, insert first point again
#
#   Parameters (necessary):
#   - lattice: The lattice object of which the band structure should be calculated
#   - path: The path along which the band structure should be calculated (in format given above)
#
#   Parameters (optional):
#   - reduceLattice: If the lattice should be reduced to a 1x1(x1) lattice of the original unitcell (i.e. for purposes of Luttinger Tisza, etc.)
#   - percentages: Either array of floats giving the individual percentages of the segments or the String "EQUAL" for equal length in the plot
#   - resolution: How many k-points to calculate in total
#   - enforce_hermitian: If the matrix should be enforced to be hermitian
#   - limits_energy: The y-axis (energy axis) limits of the plot
#   - plot_title: A title for the whole plot ("AUTO" for automatic title)
#   - plot_color: The color of bands
#   - figsize: the figure size as given to PyPlot
#
#-----------------------------------------------------------------------------------------------------------------------------
# BAND STRUCTURE ALONG PATH
function calculateBandStructureAlongPath(
        interaction_matrix::Array{Complex,2},
        path;
        percentages="EQUAL",
        resolution=1000,
        enforce_hermitian=false,
        limits_energy="AUTO",
        plot_title="",
        plot_color="b",
        figsize=(6,4),
        showPlot=true
            )

    # normalize percentages
    if percentages == "EQUAL"
        percentages = ones(size(path,1)-1)
    end
    percentages = percentages ./ sum(percentages)
    # build up the segment list
    segments = Array[]
    for i in 1:size(path,1)-1
        segment = [i, i+1, percentages[i]]
        push!(segments, segment)
    end
    # if LT is checked
    #if check_LT
    #    LT_k = Array[]
    #end
    # segment data, i.e. the bandstructure over the segments
    segments_data = Array[]
    resolution_actual = 0
    hlines = []
    # iterate over all segments
    for segment in segments
        # get the grid in between two points
        resolution_local = Int(floor(segment[3]*resolution))
        multipliers = linspace(0, 1, resolution_local)
        resolution_actual += resolution_local
        push!(hlines, resolution_actual+0)
        #println(segment)
        k1 = convert(Array{Float64,1}, path[Int(segment[1])][2:end])
        k2 = convert(Array{Float64,1}, path[Int(segment[2])][2:end])
        #println(k1)
        #println(k2)
        # insert bands
        bands = Array[]
        for b in 1:size(interaction_matrix,1)
            push!(bands, zeros(resolution_local))
        end
        # calculate all energies
        for i in 1:resolution_local
            # get the current k
            k = k2 .* multipliers[i] .+ k1 .* (1-multipliers[i])
            # if LT is checked, push current k
            #if check_LT
            #    push!(LT_k, k)
            #end
            # get the interaction matrix for this k
            matrix = copy(interaction_matrix)
            # diagonalize the matrix
            eigenvalues = eigvals(matrix)
            # save all the eigenvalues
            for b in 1:size(bands, 1)
                if imag(eigenvalues[b]) > 0
                    if imag(eigenvalues[b]) > 1e-15
                        println(imag(eigenvalues[b]))
                        println(matrix)
                        bands[b][i] = eigenvalues[b]
                    else
                        bands[b][i] = real(eigenvalues[b])
                    end
                else
                    bands[b][i] = eigenvalues[b]
                end
            end
        end
        # push the obtained back structure into the data array
        push!(segments_data, bands)
    end
    # generate the complete band structure
    bandstructure = Array[zeros(resolution_actual) for b in segments_data[1]]
    index = 1
    for i in 1:size(segments_data,1)
        segment = segments[i]
        data = segments_data[i]
        for b in 1:size(bandstructure,1)
            bandstructure[b][index:hlines[i]] = data[b]
        end
        index = hlines[i]+1
    end
    # if LT is checked, give the results
    #if check_LT
    #    LT_v = checkLuttingerTisza(lattice, LT_k, only_GS=false)
    #    println("$(100.0*sum(LT_v)/length(LT_v)) % of all eigenvalues are valid in LT")
    #end
    # plot the eigenvalues
    rc("font", family="serif")
    fig = figure(figsize=figsize)
    if plot_title == "AUTO"
        title("energy spectrum along path of interaction matrix")
    elseif plot_title == ""
        # do nothing title related
    else
        title(plot_title)
    end
    for l in hlines[1:end-1]
        axvline(l,color=[0.6, 0.6, 0.6], linestyle="--")
    end
    xlabel("momentum")
    ylabel("energy")
    for b in bandstructure
        plot(collect(1:resolution_actual), b, "-$(plot_color)")
    end
    ax = gca()
    axx = ax[:get_xaxis]()
    xtpos = []
    push!(xtpos, 0)
    for h in hlines
        push!(xtpos, h)
    end
    xtlabs = [p[1] for p in path]
    xticks(xtpos, xtlabs)
    #axx[:set_ticks]([])
    axx[:set_tick_params](which="both", direction="out")
    axx[:set_tick_params](which="top", color="none")
    axy = ax[:get_yaxis]()
    axy[:set_tick_params](which="both", direction="out")
    # check if specific boundaries are desired
    if !(limits_energy == "AUTO")
        ylim(limits_energy[1], limits_energy[2])
    end
    # tighten the layout
    tight_layout()
    # save the plot
    figurename = "interaction_matrix"
    figurename1 = "$(FOLDER_SPECTRA)bandstructure_path_$(figurename[1:end-4]).pdf"
    figurename2 = "$(FOLDER_SPECTRA)bandstructure_path_$(figurename[1:end-4]).png"
    buildFolderSpectra()
    savefig(figurename1)
    savefig(figurename2)
    if showPlot
        show()
        #print("Continue? ")
        #readline()
    end
    return fig
end
function calculateBandStructureAlongPath(
        lattice::Lattice,
        path;
        reduceLattice=true,
        percentages="EQUAL",
        resolution=1000,
        enforce_hermitian=false,
        limits_energy="AUTO",
        plot_title="",
        plot_color="b",
        figsize=(6,4),
        showPlot=true
            )

    # check if to reduce the lattice
    if reduceLattice && lattice.unitcell.filename != UNITCELL_DUMMY_FILENAME
        lattice = getLatticePeriodic(lattice.unitcell, ones(Int64, size(lattice.unitcell.lattice_vectors,1)), save=false)
    end
    # normalize percentages
    if percentages == "EQUAL"
        percentages = ones(size(path,1)-1)
    end
    percentages = percentages ./ sum(percentages)
    # build up the segment list
    segments = Array[]
    for i in 1:size(path,1)-1
        segment = [i, i+1, percentages[i]]
        push!(segments, segment)
    end
    # if LT is checked
    #if check_LT
    #    LT_k = Array[]
    #end
    # segment data, i.e. the bandstructure over the segments
    segments_data = Array[]
    resolution_actual = 0
    hlines = []
    # iterate over all segments
    for segment in segments
        # get the grid in between two points
        resolution_local = Int(floor(segment[3]*resolution))
        multipliers = linspace(0, 1, resolution_local)
        resolution_actual += resolution_local
        push!(hlines, resolution_actual+0)
        #println(segment)
        k1 = convert(Array{Float64,1}, path[Int(segment[1])][2:end])
        k2 = convert(Array{Float64,1}, path[Int(segment[2])][2:end])
        #println(k1)
        #println(k2)
        # insert bands
        bands = Array[]
        for b in 1:size(lattice.positions,1)
            push!(bands, zeros(resolution_local))
        end
        # calculate all energies
        for i in 1:resolution_local
            # get the current k
            k = k2 .* multipliers[i] .+ k1 .* (1-multipliers[i])
            # if LT is checked, push current k
            #if check_LT
            #    push!(LT_k, k)
            #end
            # get the interaction matrix for this k
            matrix = getInteractionMatrixKSpace(lattice, k, enforce_hermitian=enforce_hermitian)
            # diagonalize the matrix
            eigenvalues = eigvals(matrix)
            # save all the eigenvalues
            for b in 1:size(bands, 1)
                if imag(eigenvalues[b]) > 0
                    if imag(eigenvalues[b]) > 1e-15
                        println(imag(eigenvalues[b]))
                        println(matrix)
                        bands[b][i] = eigenvalues[b]
                    else
                        bands[b][i] = real(eigenvalues[b])
                    end
                else
                    bands[b][i] = eigenvalues[b]
                end
            end
        end
        # push the obtained back structure into the data array
        push!(segments_data, bands)
    end
    # generate the complete band structure
    bandstructure = Array[zeros(resolution_actual) for b in segments_data[1]]
    index = 1
    for i in 1:size(segments_data,1)
        segment = segments[i]
        data = segments_data[i]
        for b in 1:size(bandstructure,1)
            bandstructure[b][index:hlines[i]] = data[b]
        end
        index = hlines[i]+1
    end
    # if LT is checked, give the results
    #if check_LT
    #    LT_v = checkLuttingerTisza(lattice, LT_k, only_GS=false)
    #    println("$(100.0*sum(LT_v)/length(LT_v)) % of all eigenvalues are valid in LT")
    #end
    # plot the eigenvalues
    rc("font", family="serif")
    fig = figure(figsize=figsize)
    if plot_title == "AUTO"
        title("energy spectrum along path of lattice \"$(lattice.filename)\"")
    elseif plot_title == ""
        # do nothing title related
    else
        title(plot_title)
    end
    for l in hlines[1:end-1]
        axvline(l,color=[0.6, 0.6, 0.6], linestyle="--")
    end
    xlabel("momentum")
    ylabel("energy")
    for b in bandstructure
        plot(collect(1:resolution_actual), b, "-$(plot_color)")
    end
    ax = gca()
    axx = ax[:get_xaxis]()
    xtpos = []
    push!(xtpos, 0)
    for h in hlines
        push!(xtpos, h)
    end
    xtlabs = [p[1] for p in path]
    xticks(xtpos, xtlabs)
    #axx[:set_ticks]([])
    axx[:set_tick_params](which="both", direction="out")
    axx[:set_tick_params](which="top", color="none")
    axy = ax[:get_yaxis]()
    axy[:set_tick_params](which="both", direction="out")
    # check if specific boundaries are desired
    if !(limits_energy == "AUTO")
        ylim(limits_energy[1], limits_energy[2])
    end
    # tighten the layout
    tight_layout()
    # save the plot
    figurename = split(lattice.filename, FOLDER_SPECTRA[end])[end]
    figurename1 = "$(FOLDER_SPECTRA)bandstructure_path_$(figurename[1:end-4]).pdf"
    figurename2 = "$(FOLDER_SPECTRA)bandstructure_path_$(figurename[1:end-4]).png"
    #buildFolderSpectra()
    savefig(figurename1)
    savefig(figurename2)
    if showPlot
        show()
        #print("Continue? ")
        #readline()
    end
    return fig
end
function calculateBandStructureAlongPath(
        unitcell::Unitcell,
        path;
        percentages="EQUAL",
        resolution=1000,
        enforce_hermitian=false,
        limits_energy="AUTO",
        plot_title="",
        plot_color="b",
        figsize=(6,4),
        showPlot=true
            )

    # make a lattice from the unitcell
    lattice = getLatticePeriodic(unitcell, ones(Int64, size(unitcell.lattice_vectors,1)), save=false)
    lattice.filename = replace(identity(unitcell.filename), FOLDER_UNITCELLS, FOLDER_LATTICES)
    # push to lattice based method and return the result
    return calculateBandStructureAlongPath(
        lattice,
        path;
        reduceLattice=false,
        percentages=percentages,
        resolution=resolution,
        enforce_hermitian=enforce_hermitian,
        limits_energy=limits_energy,
        plot_title=plot_title,
        plot_color=plot_color,
        figsize=figsize,
        showPlot=showPlot
            )
end
export calculateBandStructureAlongPath

# BAND STRUCTURE OF STRIP (1D periodic)
function calculateBandStructureOfStrip(
        unitcell::Unitcell,
        periodic_direction::Int64,
        finite_direction_N::Int64;
        percentages="EQUAL",
        resolution=1000,
        enforce_hermitian=false,
        limits_energy="AUTO",
        plot_title="",
        plot_color="b",
        figsize=(6,4),
        showPlot=true,
        majorana=false
    )

    # obtain the lattice
    dimensions = ones(Int64, size(unitcell.lattice_vectors, 1)) .* finite_direction_N
    dimensions[periodic_direction] = -1
    lattice = getLattice(unitcell, dimensions)

    # obtain the path
    last_vector = lattice.lattice_vectors[1]
    last_vector = last_vector
    a = sqrt(sum(last_vector.*last_vector))
    gamma_zero  = last_vector .* 0.0
    gamma_minus = last_vector .* (-2*pi) / (a*a)
    gamma_plus  = last_vector .* ( 2*pi) / (a*a)
    K_minus     = last_vector .* (-1*pi) / (a*a)
    K_plus      = last_vector .* ( 1*pi) / (a*a)

    path = Array[
        ["Gamma (-1)";  gamma_minus],
        ["-pi";         K_minus],
        ["Gamma (0)";   gamma_zero],
        ["pi";          K_plus],
        ["Gamma (+1)";  gamma_plus]
    ]

    # calculate the band structure
    calculateBandStructureAlongPath(
        lattice,
        path;
        reduceLattice=false,
        percentages=percentages,
        resolution=resolution,
        enforce_hermitian=enforce_hermitian,
        limits_energy=limits_energy,
        plot_title=plot_title,
        plot_color=plot_color,
        figsize=figsize,
        showPlot=showPlot,
        majorana=majorana
    )
end
export calculateBandStructureOfStrip

# FULL 2D Bandstructure
function calculateBandStructure2D(
        unitcell::Unitcell,
        kx, ky;
        enforce_hermitian=false,
        limits_energy="AUTO",
        plot_title="",
        plot_color="b",
        figsize=(6,4),
        showPlot=true
            )

    # insert bands
    bands = Array[]
    kx_vals = zeros(length(kx), length(ky))
    ky_vals = zeros(length(kx), length(ky))
    for b in 1:size(unitcell.basis,1)
        push!(bands, zeros(length(kx), length(ky)))
    end
    # calculate all energies
    for i in 1:length(kx)
    for j in 1:length(ky)
        # get the current k
        k = [kx[i], ky[j]]
        kx_vals[i,j] = k[1]
        ky_vals[i,j] = k[2]
        # get the interaction matrix for this k
        matrix = getInteractionMatrixKSpace(unitcell, k, enforce_hermitian=enforce_hermitian)
        # diagonalize the matrix
        eigenvalues = eigvals(matrix)
        # save all the eigenvalues
        for b in 1:size(bands, 1)
            if imag(eigenvalues[b]) > 0
                if imag(eigenvalues[b]) > 1e-15
                    println(imag(eigenvalues[b]))
                    println(matrix)
                    bands[b][i,j] = eigenvalues[b]
                else
                    bands[b][i,j] = real(eigenvalues[b])
                end
            else
                bands[b][i,j] = eigenvalues[b]
            end
        end
    end
    end
    # generate the complete band structure
    bandstructure = bands
    # if LT is checked, give the results
    #if check_LT
    #    LT_v = checkLuttingerTisza(lattice, LT_k, only_GS=false)
    #    println("$(100.0*sum(LT_v)/length(LT_v)) % of all eigenvalues are valid in LT")
    #end
    # plot the eigenvalues
    rc("font", family="serif")
    fig = figure(figsize=figsize)
    if plot_title == "AUTO"
        if majorana
            title("majorana energy spectrum along path of unitcell \"$(unitcell.filename)\"")
        else
            title("energy spectrum along path of unitcell \"$(unitcell.filename)\"")
        end
    elseif plot_title == ""
        # do nothing title related
    else
        title(plot_title)
    end
    xlabel("momentum kx")
    ylabel("momentum ky")
    zlabel("energy")
    for b in bandstructure
        plot_surface(kx_vals,ky_vals,b, rstride=1, cstride=1, cmap="coolwarm", linewidth=0)
    end
    # tighten the layout
    tight_layout()
    # save the plot
    figurename = split(unitcell.filename, FOLDER_SPECTRA[end])[end]
    figurename1 = "$(FOLDER_SPECTRA)bandstructure_path_$(figurename[1:end-4]).pdf"
    figurename2 = "$(FOLDER_SPECTRA)bandstructure_path_$(figurename[1:end-4]).png"
    buildFolderSpectra()
    savefig(figurename1)
    savefig(figurename2)
    if showPlot
        show()
        #print("Continue? ")
        #readline()
    end
    return fig
end
export calculateBandStructure2D

# SOME DEFAULT PATHS
DEFAULT_PATH_FCC = Array[
    ["gamma"; [0,0,0]],
    ["X"; [2*pi, 0, 0]],
    ["W"; [2*pi, pi, 0]],
    ["L"; [pi, pi, pi]],
    ["gamma"; [0,0,0]],
    ["K"; [3*pi/2, 3*pi/2, 0]],
    ["X"; [2*pi, 0, 0]]
]
export DEFAULT_PATH_FCC

DEFAULT_PATH_TRIANGULAR = Array[
    ["gamma"; [0,0]],
    ["K"; [2*pi/sqrt(3.0), 2*pi/3]],
    ["M"; [2*pi/sqrt(3.0), 0]],
    ["gamma"; [0,0]]
]
export DEFAULT_PATH_TRIANGULAR

DEFAULT_PATH_SQUAREOCTAGON_2 = Array[
    ["gamma"; [0,0]],
    ["K"; [2,  0].*(pi / (1.0 + 1.0/sqrt(2.0)))],
    ["M"; [1, -1].*(pi / (1.0 + 1.0/sqrt(2.0)))],
    ["gamma"; [0,0]]
]
export DEFAULT_PATH_SQUAREOCTAGON_2

DEFAULT_PATH_SQUARE_LONG = Array[
    ["M";     [pi,   0.0]],
    ["Gamma"; [0.0,  0.0]],
    ["K'";    [pi,   -pi]],
    ["M";     [pi,   0.0]],
    ["M'";    [0.0,   pi]],
    ["K";     [pi,    pi]],
    ["Gamma"; [0.0,  0.0]]
]
export DEFAULT_PATH_SQUARE_LONG

DEFAULT_PATH_SQUARE_SHORT = Array[
    ["Gamma"; [0.0,  0.0]],
    ["M";     [pi,   0.0]],
    ["K";     [pi,    pi]],
    ["Gamma"; [0.0,  0.0]]
]
export DEFAULT_PATH_SQUARE_SHORT

DEFAULT_PATH_SQUARE = DEFAULT_PATH_SQUARE_SHORT
export DEFAULT_PATH_SQUARE


text:
@testset "5.4.2 Exponentials of inverse cotangent" begin
    (a, c, n, x, ) = @variables a c n x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*involving*exponentials*of*inverse*cotangent=#

    #= ::Section::Closed:: =#
    #=Integrands*involving*exponentials*of*inverse*cotangents=#

    @test_int [ℯ^acot(x), x, 2, (4/5 + (8*im)/5)*((-im + x)/x)^(1 + im/2)*((im + x)/x)^(-1 - im/2)*HypergeometricFunctions._₂F₁(1 + im/2, 2, 2 + im/2, (1 - im/x)/(1 + im/x))]

    @test_int [ℯ^acot(x)/(a + a*x^2), x, 1, -(ℯ^acot(x)/a)]
    @test_int [ℯ^acot(x)/(a + a*x^2)^2, x, 2, -((2*ℯ^acot(x))/(5*a^2)) - (ℯ^acot(x)*(1 - 2*x))/(5*a^2*(1 + x^2))]
    @test_int [ℯ^acot(x)/(a + a*x^2)^3, x, 3, -((24*ℯ^acot(x))/(85*a^3)) - (ℯ^acot(x)*(1 - 4*x))/(17*a^3*(1 + x^2)^2) - (12*ℯ^acot(x)*(1 - 2*x))/(85*a^3*(1 + x^2))]

    @test_int [ℯ^acot(x)/(a + a*x^2)^(3/2), x, 1, -((ℯ^acot(x)*(1 - x))/(2*a*sqrt(a + a*x^2)))]
    @test_int [ℯ^acot(x)/(a + a*x^2)^(5/2), x, 2, -((ℯ^acot(x)*(1 - 3*x))/(10*a*(a + a*x^2)^(3/2))) - (3*ℯ^acot(x)*(1 - x))/(10*a^2*sqrt(a + a*x^2))]
    @test_int [ℯ^acot(x)/(a + a*x^2)^(7/2), x, 3, -((ℯ^acot(x)*(1 - 5*x))/(26*a*(a + a*x^2)^(5/2))) - (ℯ^acot(x)*(1 - 3*x))/(13*a^2*(a + a*x^2)^(3/2)) - (3*ℯ^acot(x)*(1 - x))/(13*a^3*sqrt(a + a*x^2))]

    @test_int [ℯ^(n*acot(a*x))/(c + a^2*c*x^2)^(1/3), x, 3, (3*(1 + 1/(a^2*x^2))^(1/3)*((a - im/x)/(a + im/x))^((1/6)*(2 - 3*im*n))*(1 - im/(a*x))^((1/6)*(-2 + 3*im*n))*(1 + im/(a*x))^((1/6)*(4 - 3*im*n))*x*HypergeometricFunctions._₂F₁(-(1/3), (1/6)*(2 - 3*im*n), 2/3, (2*im)/((a + im/x)*x)))/(c + a^2*c*x^2)^(1/3)]
    @test_int [ℯ^(n*acot(a*x))/(c + a^2*c*x^2)^(2/3), x, 3, -((3*(1 + 1/(a^2*x^2))^(2/3)*((a - im/x)/(a + im/x))^((1/6)*(4 - 3*im*n))*(1 - im/(a*x))^((1/6)*(-4 + 3*im*n))*(1 + im/(a*x))^((1/6)*(2 - 3*im*n))*x*HypergeometricFunctions._₂F₁(1/3, (1/6)*(4 - 3*im*n), 4/3, (2*im)/((a + im/x)*x)))/(c + a^2*c*x^2)^(2/3))]
    @test_int [ℯ^(n*acot(a*x))/(c + a^2*c*x^2)^(4/3), x, 4, -((3*ℯ^(n*acot(a*x))*(3*n - 2*a*x))/(a*c*(4 + 9*n^2)*(c + a^2*c*x^2)^(1/3))) - (6*(1 + 1/(a^2*x^2))^(1/3)*((a - im/x)/(a + im/x))^((1/6)*(2 - 3*im*n))*(1 - im/(a*x))^((1/6)*(-2 + 3*im*n))*(1 + im/(a*x))^((1/6)*(4 - 3*im*n))*x*HypergeometricFunctions._₂F₁(-(1/3), (1/6)*(2 - 3*im*n), 2/3, (2*im)/((a + im/x)*x)))/(c*(4 + 9*n^2)*(c + a^2*c*x^2)^(1/3))]
    @test_int [ℯ^(n*acot(a*x))/(c + a^2*c*x^2)^(5/3), x, 4, -((3*ℯ^(n*acot(a*x))*(3*n - 4*a*x))/(a*c*(16 + 9*n^2)*(c + a^2*c*x^2)^(2/3))) - (12*(1 + 1/(a^2*x^2))^(2/3)*((a - im/x)/(a + im/x))^((1/6)*(4 - 3*im*n))*(1 - im/(a*x))^((1/6)*(-4 + 3*im*n))*(1 + im/(a*x))^((1/6)*(2 - 3*im*n))*x*HypergeometricFunctions._₂F₁(1/3, (1/6)*(4 - 3*im*n), 4/3, (2*im)/((a + im/x)*x)))/(c*(16 + 9*n^2)*(c + a^2*c*x^2)^(2/3))]
    @test_int [ℯ^(n*acot(a*x))/(c + a^2*c*x^2)^(7/3), x, 5, If(13>=8, -((3*ℯ^(n*acot(a*x))*(3*n - 8*a*x))/(a*c*(64 + 9*n^2)*(c + a^2*c*x^2)^(4/3))) - (120*ℯ^(n*acot(a*x))*(3*n - 2*a*x))/(a*c^2*(4 + 9*n^2)*(64 + 9*n^2)*(c + a^2*c*x^2)^(1/3)) - (240*(1 + 1/(a^2*x^2))^(1/3)*((a - im/x)/(a + im/x))^((1/6)*(2 - 3*im*n))*(1 - im/(a*x))^((1/6)*(-2 + 3*im*n))*(1 + im/(a*x))^((1/6)*(4 - 3*im*n))*x*HypergeometricFunctions._₂F₁(-(1/3), (1/6)*(2 - 3*im*n), 2/3, (2*im)/((a + im/x)*x)))/(c^2*(4 + 9*n^2)*(64 + 9*n^2)*(c + a^2*c*x^2)^(1/3)), -((3*ℯ^(n*acot(a*x))*(3*n - 8*a*x))/(a*c*(64 + 9*n^2)*(c + a^2*c*x^2)^(4/3))) - (120*ℯ^(n*acot(a*x))*(3*n - 2*a*x))/(a*c^2*(256 + 612*n^2 + 81*n^4)*(c + a^2*c*x^2)^(1/3)) - (240*(1 + 1/(a^2*x^2))^(1/3)*((a - im/x)/(a + im/x))^((1/6)*(2 - 3*im*n))*(1 - im/(a*x))^((1/6)*(-2 + 3*im*n))*(1 + im/(a*x))^((1/6)*(4 - 3*im*n))*x*HypergeometricFunctions._₂F₁(-(1/3), (1/6)*(2 - 3*im*n), 2/3, (2*im)/((a + im/x)*x)))/(c^2*(256 + 612*n^2 + 81*n^4)*(c + a^2*c*x^2)^(1/3)))]
    end

max_stars_repo_path: test/5 Inverse trig functions/5.4 Inverse cotangent/5.4.2 Exponentials of inverse cotangent.jl


text:
"""

ex = num

sum = e("@ex + @ex")
prod = e("@ex * @ex")
div = e("@ex / @ex")
num = digit
digit = 0:9

digit excluding zero = "1" | "2" | "3"

"""

@grammar begin
  start = block
  command = turnleft | turnright | forward | if_statement | loop
  turnleft = E("turn_left(ant)") # -> expr
  turnright = E("turn_right(ant)") # -> expr
  forward = E("forward(world, ant)") # -> expr
  if_statement = E("if @condition; @block else @block") # expr
  if_statement = Expr(:if, condition, block, block)
  condition = food_ahead
  food_ahead = Expr(:call, :food_ahead, :world, :ant)
  loop = Expr(:for, Expr(:(=), :i, Expr(:(:), 1, digit)), block) # this is nasty
  block[make_block] = (command)^(1:5)
  digit = 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9
end