text:
# TODO *Work in progress*, this is highly incomplete code

# Epi-composition, also known as infimal postcomposition.
# This is the dual operation to precomposition, see Rockafellar and Wets,
# "Variational Analisys", Theorem 11.23.
#
# Given a function f and a linear operator L, their epi-composition is:
#
#   g(y) = (Lf)(y) = inf_x { f(x) : Lx = y }.
#
# Plugging g directly in the definition of prox, one has:
#
#   prox_{\gamma g}(z) = argmin_y { (Lf)(y) + 1/(2\gamma)||y - z||^2 }
#     = argmin_y { inf_x { f(x) : Lx = y } + 1/(2\gamma)||y - z||^2 }
#     = L * argmin_x { f(x) + 1/(2\gamma)||Lx - z||^2 }.
#
# When L is such that L'*L = mu*Id, then this just requires prox_{\gamma f}.
# In some other cases the prox can be "easily" computed, such as when f is
# quadratic or extended-quadratic.

export Epicompose

mutable struct Epicompose{S <: AbstractMatrix, T <: ProximableFunction, F <: Factorization} <: ProximableFunction
  L::S
  f::T
  iter::Bool
  mu
  gamma
  fact::F
  function Epicompose{S, T, F}(L::S, f::T, iter::Bool) where {T <: ProximableFunction, S <: AbstractMatrix, F <: Factorization}
    new(L, f, iter, 0.0, -1.0)
  end
end

function Epicompose(L::S, f::T, iter::Bool=false) where {T <: ProximableFunction, S <: AbstractMatrix}
  # TODO this is incomplete
  Epicompose{S, T, Factorization}(L, f, iter)
end

function prox!(y, g::Epicompose{S, T}, x, gamma) where {T <: ProximableFunction, S <: AbstractMatrix}
  # this line here allocates stuff
  z = (g.L'*x)/f.mu
  p, v = prox(g.f, z, gamma/f.mu)
  A_mul_B!(y, f.L, p)
  return v
end

function prox!(y, g::Epicompose{S, T}, x, gamma) where {T <: Quadratic, S <: AbstractMatrix}
  if g.iter == false
    if gamma != g.gamma
      factor_step!(g, gamma)
    end
    # this line here allocates stuff
    p = g.fact\((g.L'*x)/gamma - g.f.q)
    fy = 0.5*vecdot(p, g.f.Q*p) + vecdot(p, g.f.q)
  else
    error("not implemented")
  end
  A_mul_B!(y, g.L, p)
  return fy
end

function factor_step!(g::Epicompose{S, T}, gamma::R) where {T <: Quadratic, S <: AbstractMatrix, R <: Real}
  g.gamma = gamma;
  g.fact = cholfact(g.f.Q + (g.L'*g.L)/gamma);
end

fun_expr(g::Epicompose) = "x ↦ inf { f(y) : Ly = x }"


text:
module CalcSimpleLoadBearing

include("./../../constraints/geometric/constants.jl")
using JuMP
using .Constants

export calc_simple_lb

function calc_simple_lb(m)

    @expression(m, accountable_mass[i=1:length(boxes), j=1:length(boxes); i!=j],
        m[:if_simple_lb_var][i, j] * masses[j]
    )

    @expression(m, sum_of_mass_on_top_of_i[i=1:length(boxes)],
        sum(m[:accountable_mass][i, j] for j in 1:length(boxes) if i!=j)
    )

    @constraint(m, c_sum_of_mass_on_top_of_i[i=1:length(boxes)-1],
        m[:sum_of_mass_on_top_of_i][i] <= mass_volume[i]
    )

    return m
end

end

ext: jl


text:
# Autogenerated wrapper script for louvain_jll for aarch64-linux-musl-cxx03
export convertlouvain, hierarchy, louvain

## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"
LIBPATH_default = ""

# Relative path to `convertlouvain`
const convertlouvain_splitpath = ["bin", "convert"]

# This will be filled out by __init__() for all products, as it must be done at runtime
convertlouvain_path = ""

# convertlouvain-specific global declaration
function convertlouvain(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(convertlouvain_path)
    end
end

# Relative path to `hierarchy`
const hierarchy_splitpath = ["bin", "hierarchy"]

# This will be filled out by __init__() for all products, as it must be done at runtime
hierarchy_path = ""

# hierarchy-specific global declaration
function hierarchy(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(hierarchy_path)
    end
end

# Relative path to `louvain`
const louvain_splitpath = ["bin", "louvain"]

# This will be filled out by __init__() for all products, as it must be done at runtime
louvain_path = ""

# louvain-specific global declaration
function louvain(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(louvain_path)
    end
end

# Inform that the wrapper is available for this platform
wrapper_available = true

"""
Open all libraries
"""
function __init__()
    # This either calls `@artifact_str()`, or returns a constant string if we're overridden.
    global artifact_dir = find_artifact_dir()

    global PATH_list, LIBPATH_list
    global convertlouvain_path = normpath(joinpath(artifact_dir, convertlouvain_splitpath...))

    push!(PATH_list, dirname(convertlouvain_path))
    global hierarchy_path = normpath(joinpath(artifact_dir, hierarchy_splitpath...))

    push!(PATH_list, dirname(hierarchy_path))
    global louvain_path = normpath(joinpath(artifact_dir, louvain_splitpath...))

    push!(PATH_list, dirname(louvain_path))
    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(vcat(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)]), ':')

end  # __init__()


text:
using LinearAlgebra
using IterativeSolvers
import LinearAlgebra: ldiv!

struct Projector{VT<:AbstractVector}
    v::VT
end

ldiv!(p::Projector, x) = x - p.v*(p.v'*x)
ldiv!(x, p::Projector) = x - p.v'*(x*p.v)

gs(A::AbstractMatrix) = (ei = eigen(A); (ei.values[1], ei.vectors[:,1]))
function gs_backward(A, e, v, de, dv)
    #x = gmres((I-v*v')*(e*I-A), (I-v*v')*dv)
    x = gmres((e*I-A)*(I-v*v'), dv)'
    #x = gmres(e*I-A, dv, Pl=Projector(v))
    #x = gmres(e*I-A, (I-v*v')*dv, Pl=Projector(v))
    #x = gmres(e*I-A, (I-v*v')*dv)
    #(de*v + (I-v*v')*pinv(e*I-A)*dv)*v'
    (de*v + x)*v'
end

function num_diff(A, de, dv)
    δ = 1e-5
    A = A |> copy
    dA = zeros(n, n)
    for i = 1:n
        for j = 1:n
            A[i,j]+=δ/2
            pos = gs(A)
            A[i,j]-=δ
            neg = gs(A)
            A[i,j]+=δ/2
            grad = (pos.-neg)./δ
            dA[i,j] = grad[1]*de + sum(grad[2].*dv)
        end
    end
    dA
end

n=3
A = randn(n,n); A += A'
e,v = gs(A)
de = randn()
dv = randn(n)
de = 0
dv = [1.0, 0,0]
dA = gs_backward(A, e, v, de, dv)

@show num_diff(A, e, v, de, dv)
@show dA

using Test
@test isapprox(dA .|> abs, num_diff() .|> abs)


text:
module ParamEstim

export
    C,
    V,
    observables,
    observables_index,
    Sim,
    Exp,
    Viz,
    param_values,
    initial_values,
    objective,
    get_search_index,
    get_search_region,
    decode_gene2val,
    encode_val2gene,
    encode_bestIndivVal2randGene,
    update_param,
    simulate_all,
    tpm2ival!,
    mul_ratio2mcf7!,
    ga_v1,
    ga_v1_continue,
    ga_v2,
    ga_v2_continue

using Printf
using LinearAlgebra
using Random
using StatsBase
using Statistics
using DelimitedFiles
using PyPlot

import Seaborn

include("model/name2idx/parameters.jl")
include("model/name2idx/species.jl")
include("model/set_model.jl")
include("model/observable.jl")
include("model/experimental_data.jl")
include("model/simulation.jl")
include("model/viz.jl")
include("model/fitness.jl")
include("model/set_search_param.jl")

include("ga/initial_population.jl")
include("ga/undxmgg.jl")
include("ga/converging.jl")
include("ga/local_search.jl")
include("ga/v1.jl")
include("ga/v2.jl")

include("dynamics/signaling_systems.jl")
include("dynamics/temporal_dynamics.jl")

include("load_csv.jl")

end  # module


text:
#! format: off
requires_initialization(::AbstractStorageFormulation) = false

get_variable_multiplier(_, ::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = NaN
get_expression_type_for_reserve(::ActivePowerReserveVariable, ::Type{<:PSY.Storage}, ::Type{<:PSY.Reserve{PSY.ReserveUp}}) = ReserveRangeExpressionUB
get_expression_type_for_reserve(::ActivePowerReserveVariable, ::Type{<:PSY.Storage}, ::Type{<:PSY.Reserve{PSY.ReserveDown}}) = ReserveRangeExpressionLB
########################### ActivePowerInVariable, Storage #################################
get_variable_binary(::ActivePowerInVariable, ::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = false
get_variable_lower_bound(::ActivePowerInVariable, d::PSY.Storage, ::AbstractStorageFormulation) = 0.0
get_variable_upper_bound(::ActivePowerInVariable, d::PSY.Storage, ::AbstractStorageFormulation) = PSY.get_input_active_power_limits(d).max
get_variable_multiplier(::ActivePowerInVariable, d::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = -1.0

########################### ActivePowerOutVariable, Storage #################################
get_variable_binary(::ActivePowerOutVariable, ::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = false
get_variable_lower_bound(::ActivePowerOutVariable, d::PSY.Storage, ::AbstractStorageFormulation) = 0.0
get_variable_upper_bound(::ActivePowerOutVariable, d::PSY.Storage, ::AbstractStorageFormulation) = PSY.get_output_active_power_limits(d).max
get_variable_multiplier(::ActivePowerOutVariable, d::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = 1.0

############## ReactivePowerVariable, Storage ####################
get_variable_multiplier(::ReactivePowerVariable, ::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = 1.0
get_variable_binary(::ReactivePowerVariable, ::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = false

############## EnergyVariable, Storage ####################
get_variable_binary(::EnergyVariable, ::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = false
get_variable_upper_bound(::EnergyVariable, d::PSY.Storage, ::AbstractStorageFormulation) = PSY.get_state_of_charge_limits(d).max
get_variable_lower_bound(::EnergyVariable, d::PSY.Storage, ::AbstractStorageFormulation) = PSY.get_state_of_charge_limits(d).min
get_variable_warm_start_value(::EnergyVariable, d::PSY.Storage, ::AbstractStorageFormulation) = PSY.get_initial_energy(d)

############## ReservationVariable, Storage ####################
get_variable_binary(::ReservationVariable, ::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = true
get_efficiency(v::T, var::Type{<:InitialConditionType}) where T <: PSY.Storage = PSY.get_efficiency(v)

############## EnergyShortageVariable, Storage ####################
get_variable_binary(::EnergyShortageVariable, ::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = false
get_variable_lower_bound(::EnergyShortageVariable, d::PSY.Storage, ::AbstractStorageFormulation) = 0.0
get_variable_upper_bound(::EnergyShortageVariable, d::PSY.Storage, ::AbstractStorageFormulation) = PSY.get_rating(d)

############## EnergySurplusVariable, Storage ####################
get_variable_binary(::EnergySurplusVariable, ::Type{<:PSY.Storage}, ::AbstractStorageFormulation) = false
get_variable_upper_bound(::EnergySurplusVariable, d::PSY.Storage, ::AbstractStorageFormulation) = 0.0
get_variable_lower_bound(::EnergySurplusVariable, d::PSY.Storage, ::AbstractStorageFormulation) = - PSY.get_rating(d)

#################### Initial Conditions for models ###############
initial_condition_default(::InitialEnergyLevel, d::PSY.Storage, ::AbstractStorageFormulation) = PSY.get_initial_energy(d)
initial_condition_variable(::InitialEnergyLevel, d::PSY.Storage, ::AbstractStorageFormulation) = EnergyVariable()

########################### Parameter related set functions ################################
get_parameter_multiplier(::VariableValueParameter, d::PSY.Storage, ::AbstractStorageFormulation) = 1.0
get_initial_parameter_value(::VariableValueParameter, d::PSY.Storage, ::AbstractStorageFormulation) = 1.0

########################Objective Function##################################################
objective_function_multiplier(::VariableType, ::AbstractStorageFormulation)=OBJECTIVE_FUNCTION_POSITIVE
objective_function_multiplier(::EnergySurplusVariable, ::EnergyTarget)=OBJECTIVE_FUNCTION_NEGATIVE
objective_function_multiplier(::EnergyShortageVariable, ::EnergyTarget)=OBJECTIVE_FUNCTION_POSITIVE

proportional_cost(cost::PSY.StorageManagementCost, ::EnergySurplusVariable, ::PSY.BatteryEMS, ::EnergyTarget)=PSY.get_energy_surplus_cost(cost)
proportional_cost(cost::PSY.StorageManagementCost, ::EnergyShortageVariable, ::PSY.BatteryEMS, ::EnergyTarget)=PSY.get_energy_shortage_cost(cost)

variable_cost(cost::PSY.StorageManagementCost, ::ActivePowerOutVariable, ::PSY.BatteryEMS, ::EnergyTarget)=PSY.get_variable(cost)

#! format: on

get_initial_conditions_device_model(
    ::OperationModel,
    ::DeviceModel{T, <:AbstractStorageFormulation},
) where {T <: PSY.Storage} = DeviceModel(T, BookKeeping)

get_multiplier_value(
    ::EnergyTargetTimeSeriesParameter,
    d::PSY.Storage,
    ::AbstractStorageFormulation,
) = PSY.get_rating(d)

function get_default_time_series_names(
    ::Type{D},
    ::Type{EnergyTarget},
) where {D <: PSY.Storage}
    return Dict{Type{<:TimeSeriesParameter}, String}(
        EnergyTargetTimeSeriesParameter => "storage_target",
    )
end

function get_default_time_series_names(
    ::Type{D},
    ::Type{<:Union{FixedOutput, AbstractStorageFormulation}},
) where {D <: PSY.Storage}
    return Dict{Type{<:TimeSeriesParameter}, String}()
end

function get_default_attributes(
    ::Type{D},
    ::Type{T},
) where {D <: PSY.Storage, T <: Union{FixedOutput, AbstractStorageFormulation}}
    return Dict{String, Any}("reservation" => true)
end

################################## output power constraints#################################

get_min_max_limits(
    device::PSY.Storage,
    ::Type{<:ReactivePowerVariableLimitsConstraint},
    ::Type{<:AbstractStorageFormulation},
) = PSY.get_reactive_power_limits(device)
get_min_max_limits(
    device::PSY.Storage,
    ::Type{<:InputActivePowerVariableLimitsConstraint},
    ::Type{<:AbstractStorageFormulation},
) = PSY.get_input_active_power_limits(device)
get_min_max_limits(
    device::PSY.Storage,
    ::Type{<:OutputActivePowerVariableLimitsConstraint},
    ::Type{<:AbstractStorageFormulation},
) = PSY.get_output_active_power_limits(device)
get_min_max_limits(
    device::PSY.Storage,
    ::Type{<:OutputActivePowerVariableLimitsConstraint},
    ::Type{BookKeeping},
) = PSY.get_output_active_power_limits(device)

function add_constraints!(
    container::OptimizationContainer,
    T::Type{<:PowerVariableLimitsConstraint},
    U::Type{<:Union{VariableType, ExpressionType}},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    X::Type{<:PM.AbstractPowerModel},
) where {V <: PSY.Storage, W <: AbstractStorageFormulation}
    if get_attribute(model, "reservation")
        add_reserve_range_constraints!(container, T, U, devices, model, X)
    else
        add_range_constraints!(container, T, U, devices, model, X)
    end
end

########################## Make initial Conditions for a Model #############################
function initial_conditions!(
    container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{St},
    formulation::AbstractStorageFormulation,
) where {St <: PSY.Storage}
    add_initial_condition!(container, devices, formulation, InitialEnergyLevel())
    return
end

############################ Energy Capacity Constraints####################################
"""
Min and max limits for Energy Capacity Constraint and AbstractStorageFormulation
"""
function get_min_max_limits(
    d,
    ::Type{EnergyCapacityConstraint},
    ::Type{<:AbstractStorageFormulation},
)
    return PSY.get_state_of_charge_limits(d)
end

"""
Add Energy Capacity Constraints for AbstractStorageFormulation
"""
function add_constraints!(
    container::OptimizationContainer,
    T::Type{EnergyCapacityConstraint},
    U::Type{<:VariableType},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    ::Type{X},
) where {V <: PSY.Storage, W <: AbstractStorageFormulation, X <: PM.AbstractPowerModel}
    add_range_constraints!(
        container,
        EnergyCapacityConstraint,
        EnergyVariable,
        devices,
        model,
        X,
    )
    return
end

############################ book keeping constraints ######################################

"""
Add Energy Balance Constraints for AbstractStorageFormulation
"""
function add_constraints!(
    container::OptimizationContainer,
    ::Type{EnergyBalanceConstraint},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    ::Type{X},
) where {V <: PSY.Storage, W <: AbstractStorageFormulation, X <: PM.AbstractPowerModel}
    time_steps = get_time_steps(container)
    resolution = get_resolution(container)
    fraction_of_hour = Dates.value(Dates.Minute(resolution)) / MINUTES_IN_HOUR
    names = [PSY.get_name(x) for x in devices]
    initial_conditions = get_initial_condition(container, InitialEnergyLevel(), V)
    energy_var = get_variable(container, EnergyVariable(), V)
    powerin_var = get_variable(container, ActivePowerInVariable(), V)
    powerout_var = get_variable(container, ActivePowerOutVariable(), V)

    constraint = add_constraints_container!(
        container,
        EnergyBalanceConstraint(),
        V,
        names,
        time_steps,
    )

    for ic in initial_conditions
        device = get_component(ic)
        efficiency = PSY.get_efficiency(device)
        name = PSY.get_name(device)
        constraint[name, 1] = JuMP.@constraint(
            container.JuMPmodel,
            energy_var[name, 1] ==
            get_value(ic) +
            (
                powerin_var[name, 1] * efficiency.in -
                (powerout_var[name, 1] / efficiency.out)
            ) * fraction_of_hour
        )

        for t in time_steps[2:end]
            constraint[name, t] = JuMP.@constraint(
                container.JuMPmodel,
                energy_var[name, t] ==
                energy_var[name, t - 1] +
                (
                    powerin_var[name, t] * efficiency.in -
                    (powerout_var[name, t] / efficiency.out)
                ) * fraction_of_hour
            )
        end
    end
    return
end

############################ reserve constraints ######################################
function add_constraints!(
    container::OptimizationContainer,
    ::Type{ReserveEnergyConstraint},
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, D},
    ::Type{<:PM.AbstractPowerModel},
) where {T <: PSY.Storage, D <: AbstractStorageFormulation}
    time_steps = get_time_steps(container)
    var_e = get_variable(container, EnergyVariable(), T)
    expr_up = get_expression(container, ReserveRangeExpressionUB(), T)
    expr_dn = get_expression(container, ReserveRangeExpressionLB(), T)
    names = [PSY.get_name(x) for x in devices]
    con_up = add_constraints_container!(
        container,
        ReserveEnergyConstraint(),
        T,
        names,
        time_steps,
        meta="up",
    )
    con_dn = add_constraints_container!(
        container,
        ReserveEnergyConstraint(),
        T,
        names,
        time_steps,
        meta="dn",
    )

    for d in devices, t in time_steps
        name = PSY.get_name(d)
        limits = PSY.get_state_of_charge_limits(d)
        efficiency = PSY.get_efficiency(d)
        con_up[name, t] = JuMP.@constraint(
            container.JuMPmodel,
            expr_up[name, t] <= (var_e[name, t] - limits.min) * efficiency.out
        )
        con_dn[name, t] = JuMP.@constraint(
            container.JuMPmodel,
            expr_dn[name, t] <= (limits.max - var_e[name, t]) / efficiency.in
        )
    end
    return
end

function add_constraints!(
    container::OptimizationContainer,
    ::Type{RangeLimitConstraint},
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, D},
    ::Type{<:PM.AbstractPowerModel},
) where {T <: PSY.Storage, D <: AbstractStorageFormulation}
    time_steps = get_time_steps(container)
    var_in = get_variable(container, ActivePowerInVariable(), T)
    var_out = get_variable(container, ActivePowerOutVariable(), T)
    expr_up = get_expression(container, ReserveRangeExpressionUB(), T)
    expr_dn = get_expression(container, ReserveRangeExpressionLB(), T)
    names = [PSY.get_name(x) for x in devices]
    con_up = add_constraints_container!(
        container,
        RangeLimitConstraint(),
        T,
        names,
        time_steps,
        meta="up",
    )
    con_dn = add_constraints_container!(
        container,
        RangeLimitConstraint(),
        T,
        names,
        time_steps,
        meta="dn",
    )

    for d in devices, t in time_steps
        name = PSY.get_name(d)
        out_limits = PSY.get_output_active_power_limits(d)
        in_limits = PSY.get_input_active_power_limits(d)
        con_up[name, t] = JuMP.@constraint(
            container.JuMPmodel,
            expr_up[name, t] <= var_in[name, t] + (out_limits.max - var_out[name, t])
        )
        con_dn[name, t] = JuMP.@constraint(
            container.JuMPmodel,
            expr_dn[name, t] <= var_out[name, t] + (in_limits.max - var_in[name, t])
        )
    end
    return
end

############################ Energy Management constraints ######################################
"""
Add Energy Target Constraints for EnergyTarget formulation
"""
function add_constraints!(
    container::OptimizationContainer,
    ::Type{EnergyTargetConstraint},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    ::Type{X},
) where {V <: PSY.Storage, W <: EnergyTarget, X <: PM.AbstractPowerModel}
    time_steps = get_time_steps(container)
    name_index = [PSY.get_name(d) for d in devices]
    energy_var = get_variable(container, EnergyVariable(), V)
    shortage_var = get_variable(container, EnergyShortageVariable(), V)
    surplus_var = get_variable(container, EnergySurplusVariable(), V)

    param = get_parameter_array(container, EnergyTargetTimeSeriesParameter(), V)
    multiplier =
        get_parameter_multiplier_array(container, EnergyTargetTimeSeriesParameter(), V)

    constraint = add_constraints_container!(
        container,
        EnergyTargetConstraint(),
        V,
        name_index,
        time_steps,
    )
    for d in devices
        name = PSY.get_name(d)
        shortage_cost = PSY.get_energy_shortage_cost(PSY.get_operation_cost(d))
        if shortage_cost == 0.0
            @warn(
                "Device $name has energy shortage cost set to 0.0, as a result the model will turnoff the EnergyShortageVariable to avoid infeasible/unbounded problem."
            )
            JuMP.delete_upper_bound.(shortage_var[name, :])
            JuMP.set_upper_bound.(shortage_var[name, :], 0.0)
        end
        for t in time_steps
            constraint[name, t] = JuMP.@constraint(
                container.JuMPmodel,
                energy_var[name, t] + shortage_var[name, t] + surplus_var[name, t] ==
                multiplier[name, t] * param[name, t]
            )
        end
    end
    return
end

##################################### Storage generation cost ############################
function objective_function!(
    container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    ::DeviceModel{T, U},
    ::Type{<:PM.AbstractPowerModel},
) where {T <: PSY.Storage, U <: AbstractStorageFormulation}
    add_proportional_cost!(container, ActivePowerOutVariable(), devices, U())
    add_proportional_cost!(container, ActivePowerInVariable(), devices, U())
    return
end

function objective_function!(
    container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{PSY.BatteryEMS},
    ::DeviceModel{PSY.BatteryEMS, T},
    ::Type{<:PM.AbstractPowerModel},
) where {T <: EnergyTarget}
    add_variable_cost!(container, ActivePowerOutVariable(), devices, T())
    add_proportional_cost!(container, EnergySurplusVariable(), devices, T())
    add_proportional_cost!(container, EnergyShortageVariable(), devices, T())
    return
end


text:
##############################################################################
# This implementation is based on https://github.com/julesjacobs/levenshtein
# and related blog post of Jules Jacobs
# http://julesjacobs.github.io/2015/06/17/disqus-levenshtein-simple-and-fast.html
struct String32{N} <: AbstractString
    s::NTuple{N, Char}
end

function String32(s::S) where {S <: AbstractString}
    return String32(Tuple(s))
end

String32(s::String32) = s

Base.:length(s32::String32) = length(s32.s)
Base.:(==)(s1::String32, s2::String32) = s1.s == s2.s

function Base.:iterate(s32::String32, state::Integer = 1)
    state > length(s32.s) && return nothing
    return (s32.s[state], state + 1)
end

Base.:ncodeunits(s32::String32) = length(s32)
Base.:isvalid(s::String32, i::Integer) = true
Base.:getindex(s32::String32, i::Integer) = s32.s[i]

"""
    LevenshteinAutomaton(term::String, n::Int)

Construct automaton for given `term` and distance `n`.
"""
struct LevenshteinAutomaton
    term::String32
    n::Int
end

"""
    start

Returns the start state.
"""
start(LA) = MVector{length(LA.term) + 1}(collect(0:length(LA.term)))

"""
    step

Returns the next state given a `state` and a character `c`.
"""
function step(LA, state, c)
    new_state = similar(state)
    new_state[1] = state[1] + 1
    for i in 1:length(state) - 1
        cost = LA.term[i] == c ? 0 : 1
        new_state[i + 1] = min(new_state[i] + 1, state[i] + cost, state[i + 1] + 1)
    end

    # operation of making finite space
    for i in 1:length(new_state)
        new_state[i] = new_state[i] <= LA.n+1 ? new_state[i] : LA.n+1
    end
    return new_state
end

"""
    is_match

Tells us whether a state is matching
"""
is_match(LA, state) = state[end] <= LA.n

"""
    can_match

tells us whether a `state` can become matching if we input more characters.
`can_match` is used to prune the search tree.

```julia
automaton = LevenshteinAutomaton("bannana", 1)

state0 = start(automaton)
state1 = step(automaton, state0, 'w')
can_amtch(automaton, state1) # `true`, 'w' can match "bannana" with distance 1
state2 = step(automaton, state1, 'o')
can_match(automaton, state2) # `false`, "wo" can't match "bannana" with distance 1
```
"""
can_match(LA, state) = minimum(state) <= LA.n

"""
    transitions

Computes the letters to try
"""
function transitions(LA, state)
    Set(c for (i, c) in enumerate(LA.term) if state[i] <= LA.n)
end

"""
    build_dfa

Generates graph of all possible states and transitions. This graph can be later converted
to DFA.
"""
function build_dfa(LA, state::S = start(LA)) where S
    states = Dict{S, Int}()
    matching = Int[]
    vtrans = Tuple{Int, Int, Char}[]

    explore(LA, state, states, matching, vtrans)

    # Here should be DFA output
    return states, matching, vtrans
end

function explore(LA, state, states, matching, vtrans)
    sleep(0.1)
    haskey(states, state) && return states[state]
    i = length(states) + 1
    states[state] = i
    is_match(LA, state) && push!(matching, i)
    for c in transitions(LA, state)
        newstate = step(LA, state, c)
        j = explore(LA, newstate, states, matching, vtrans)
        push!(vtrans, (i, j, c))
    end

    return i
end

########################################################

"""
    SparseLevenshteinAutomaton(term::String, n::Int)

Construct automaton for given `term` and distance `n`. Uses
different way to store intermidiate information
"""
struct SparseLevenshteinAutomaton
    term::String32
    n::Int
end

struct SparseLevenshteinAutomatonState
    indices::Vector{Int}
    values::Vector{Int}
end

"""
    start

Returns the start state.
"""
start(LA::SparseLevenshteinAutomaton) = SparseLevenshteinAutomatonState(collect(1:LA.n+1), collect(0:LA.n))

"""
    step

Returns the next state given a `state` and a character `c`.
"""
function step(LA::SparseLevenshteinAutomaton, state, c)
    if !isempty(state.indices) && state.indices[1] == 1 && state.values[1] < LA.n
        new_indices = Int[1]
        new_values = Int[state.values[1] + 1]
    else
        new_indices = Int[]
        new_values = Int[]
    end

    for (j, i) in enumerate(state.indices)
        i == length(LA.term) && break
        cost = LA.term[i] == c ? 0 : 1
        val = state.values[j] + cost

        if !isempty(new_values) && new_indices[end] == i
            val = min(val, new_values[end] + 1)
        end
        if j + 1 <= length(state.indices) && state.indices[j+1] == i+1
            val = min(val, state.values[j+1] + 1)
        end
        if val <= LA.n
            push!(new_indices, i+1)
            push!(new_values, val)
        end
    end
    return SparseLevenshteinAutomatonState(new_indices, new_values)
end

"""
    is_match

Tells us whether a state is matching
"""
is_match(LA::SparseLevenshteinAutomaton, state) = !isempty(state.indices) && state.indices[end] == length(LA.term)

"""
    can_match

tells us whether a `state` can become matching if we input more characters.
`can_match` is used to prune the search tree.

```julia
automaton = LevenshteinAutomaton("bannana", 1)

state0 = start(automaton)
state1 = step(automaton, state0, 'w')
can_amtch(automaton, state1) # `true`, 'w' can match "bannana" with distance 1
state2 = step(automaton, state1, 'o')
can_match(automaton, state2) # `false`, "wo" can't match "bannana" with distance 1
```
"""
can_match(LA::SparseLevenshteinAutomaton, state) = !isempty(state.indices)

"""
    transitions

Computes the letters to try
"""
function transitions(LA::SparseLevenshteinAutomaton, state)
    Set(LA.term[state.indices])
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -33.59531116254298 <= q <= 31.312008728451346)

                     add_NL_constraint(m, :(gelu(0.7464435003747418 + -0.7676428720221384*gelu(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4])) + -0.5937285886485864*gelu(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4])) + 0.70936536775501*gelu(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4])) + 0.9647291374973435*gelu(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4])) + 0.527810271934432*gelu(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))) + gelu(-0.35636357579153266 + 0.8665784119884141*gelu(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4])) + 0.03175667478129274*gelu(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4])) + -0.7799384758999897*gelu(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4])) + 0.4104070349373172*gelu(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4])) + -0.7416668360909755*gelu(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))) + gelu(0.5038932713674513 + 0.7223377417739489*gelu(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4])) + -0.17110213406275054*gelu(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4])) + 0.7726394930743403*gelu(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4])) + 0.6118790711423934*gelu(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4])) + 0.16250007841990932*gelu(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))) + gelu(-0.7190092552996492 + 0.6952008118484776*gelu(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4])) + 0.975536624348428*gelu(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4])) + 0.03367235210941821*gelu(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4])) + 0.8059867473566515*gelu(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4])) + 0.5874603423029647*gelu(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))) + gelu(-0.5715134646750415 + 0.6325466436773857*gelu(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4])) + -0.2633742432221431*gelu(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4])) + -0.3786600783719858*gelu(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4])) + 0.5923927477083897*gelu(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4])) + 0.9114166052727355*gelu(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
"""
$(SIGNATURES)

Prepare the initial state in proper type and dimension for ODE solvers. `type` specifies the dimension of the initial state: `:v` is 1-D state vector and `:m` is 2-D density matrix. `vectorize` indicate whether to vectorize the density matrix.
"""
function build_u0(raw_u0, type; vectorize::Bool = false)
    res = complex(raw_u0)
    if type == :v && ndims(res) != 1
        throw(ArgumentError("Cannot convert density matrix to state vector."))
    elseif type == :m && ndims(res) == 1
        res = res * res'
    elseif ndims(res) < 1 || ndims(res) > 2
        throw(ArgumentError("u0 can either be a vector or matrix."))
    end
    if vectorize == true
        res = res[:]
    end
    # use StaticArrays for initial state will cause some method to fail
    # if (ndims(res) == 1) && (length(res) <= 10)
    #     res = OpenQuantumBase.MVector{length(res)}(res)
    # elseif (ndims(res) == 2) && (size(res, 1) <= 10)
    #     res = OpenQuantumBase.MMatrix{size(res, 1),size(res, 2)}(res)
    # end
    res
end

function vectorize_cache(cache)
    # use StaticArray if the dimension is less than 3
    if size(cache, 1) <= 3
        OpenQuantumBase.@MMatrix zeros(eltype(cache), size(cache, 1)^2, size(cache, 2)^2)
    else
        one(cache) ⊗ cache
    end
end


text:
# Autogenerated wrapper script for PTSCOTCH_jll for x86_64-apple-darwin
export libptesmumps, libptscotch, libptscotcherr, libptscotcherrexit, libptscotchparmetis

using Zlib_jll
using MPICH_jll
using SCOTCH_jll
JLLWrappers.@generate_wrapper_header("PTSCOTCH")
JLLWrappers.@declare_library_product(libptesmumps, "@rpath/libesmumps.dylib")
JLLWrappers.@declare_library_product(libptscotch, "@rpath/libptscotch.dylib")
JLLWrappers.@declare_library_product(libptscotcherr, "@rpath/libptscotcherr.dylib")
JLLWrappers.@declare_library_product(libptscotcherrexit, "@rpath/libptscotcherrexit.dylib")
JLLWrappers.@declare_library_product(libptscotchparmetis, "@rpath/libptscotchparmetis.dylib")
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll, MPICH_jll, SCOTCH_jll)
    JLLWrappers.@init_library_product(
        libptesmumps,
        "lib/libptesmumps.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libptscotch,
        "lib/libptscotch.dylib",
        nothing,
    )

    JLLWrappers.@init_library_product(
        libptscotcherr,
        "lib/libptscotcherr.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libptscotcherrexit,
        "lib/libptscotcherrexit.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libptscotchparmetis,
        "lib/libptscotchparmetis.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
int main() {
  int[1] a;
  return 0;
}

max_stars_repo_name: fredefox/javalette

max_issues_repo_name: fredefox/javalette

max_forks_repo_name: fredefox/javalette


text:
using Printf

using BenchmarkTools
using CUDA
using DataFrames
using PrettyTables

using Oceananigans
using Oceananigans.Architectures
using JULES

using BenchmarkTools: prettytime, prettymemory

Archs = [CPU]
@hascuda Archs = [CPU, GPU]

Ns = [32, 128]
Tvars = [Energy, Entropy]
Gases = [DryEarth, DryEarth3]

tags = ["arch", "N", "gases", "tvar"]

suite = BenchmarkGroup(
    "cfl" => BenchmarkGroup(),
    "acoustic_cfl" => BenchmarkGroup()
)

_cfl(model) = cfl(model, 1)
_cfl(model::CompressibleModel{GPU}) = CUDA.@sync cfl(model, 1)
_acoustic_cfl(model) = acoustic_cfl(model, 1)
_acoustic_cfl(model::CompressibleModel{GPU}) = CUDA.@sync acoustic_cfl(model, 1)

for Arch in Archs, N in Ns, Gas in Gases, Tvar in Tvars
    @info "Running CFL benchmark [$Arch, N=$N, $Tvar, $Gas]..."

    grid = RegularCartesianGrid(size=(N, N, N), extent=(1, 1, 1))
    model = CompressibleModel(architecture=Arch(), grid=grid, thermodynamic_variable=Tvar(),
                              gases=Gas())

    # warmup
    _cfl(model)
    _acoustic_cfl(model)

    b_cfl = @benchmark _cfl($model) samples=10
    display(b_cfl)

    b_acfl = @benchmark _acoustic_cfl($model) samples=10
    display(b_acfl)

    key = (Arch, N, Gas, Tvar)
    suite["cfl"][key] = b_cfl
    suite["acoustic_cfl"][key] = b_acfl
end

function benchmarks_to_dataframe(suite)
    df = DataFrame(architecture=[], size=[], gases=[],
                   thermodynamic_variable=[], min=[], median=[],
                   mean=[], max=[], memory=[], allocs=[])

    for Arch in Archs, N in Ns, Gas in Gases, Tvar in Tvars
        b = suite[Arch, N, Gas, Tvar]

        d = Dict(
            "architecture" => Arch,
            "size" => "$(N)³",
            "gases" => Gas,
            "thermodynamic_variable" => Tvar,
            "min" => minimum(b.times) |> prettytime,
            "median" => median(b.times) |> prettytime,
            "mean" => mean(b.times) |> prettytime,
            "max" => maximum(b.times) |> prettytime,
            "memory" => prettymemory(b.memory),
            "allocs" => b.allocs
        )

        push!(df, d)
    end

    return df
end

header = ["Arch" "Size" "Gases" "ThermoVar" "min" "median" "mean" "max" "memory" "allocs"]

df = benchmarks_to_dataframe(suite["cfl"])
pretty_table(df, header, title="CFL benchmarks", nosubheader=true)

df = benchmarks_to_dataframe(suite["acoustic_cfl"])
pretty_table(df, header, title="Acoustic CFL benchmarks", nosubheader=true)

function gpu_speedups(suite)
    df = DataFrame(size=[], gases=[], thermodynamic_variable=[], speedup=[])

    for N in Ns, Gas in Gases, Tvar in Tvars
        b_cpu = suite[CPU, N, Gas, Tvar] |> median
        b_gpu = suite[GPU, N, Gas, Tvar] |> median
        b_ratio = ratio(b_cpu, b_gpu)

        d = Dict(
            "size" => "$(N)³",
            "gases" => Gas,
            "thermodynamic_variable" => Tvar,
            "speedup" => @sprintf("%.3fx", b_ratio.time)
        )

        push!(df, d)
    end

    return df
end

header = ["Size" "Gases" "ThermoVar" "speedup"]

df = gpu_speedups(suite["cfl"])
pretty_table(df, header, title="CFL speedups", nosubheader=true)

df = gpu_speedups(suite["acoustic_cfl"])
pretty_table(df, header, title="Acoustic CFL speedups", nosubheader=true)


text:
abstract type AbstractMessage end

struct NetworkEnvelope
    magic::UInt32
    command::String
    length::UInt32
    checksum::UInt32
    payload::Vector{UInt8}
    NetworkEnvelope(magic::UInt32, command::AbstractString, length::Integer, checksum::UInt32, payload::Vector{UInt8}) = new(magic, command, length, checksum, payload)
end

NetworkEnvelope(command::AbstractString, payload::Vector{UInt8}, testnet=false) = NetworkEnvelope(NETWORK_MAGIC[testnet], command, length(payload), reinterpret(UInt32, hash256(payload))[1], payload)
NetworkEnvelope(command::AbstractString, payload::AbstractMessage, testnet=false) = NetworkEnvelope(command, serialize(payload), testnet)

function show(io::IO, z::NetworkEnvelope)
    print(io, "Network Envelope - ", z.command, "\n", z.length, " bytes, checksum ", string(z.checksum, base=16), " on network ", string(z.magic, base=16), "\n--------\n", bytes2hex(z.payload))
end

"""
    IOBuffer, Bool -> Vector{NetworkEnvelope}

Takes a stream and creates a NetworkEnvelope
"""
function io2envelopes(s::IOBuffer, testnet::Bool=false)
    result = NetworkEnvelope[]
    while !(bytesavailable(s) < 24)
        try
            global magic = reinterpret(UInt32, read(s, 4))[1]
        catch BoundsError
            error("Connection reset!", bytes2hex(bin))
        end
        if magic != NETWORK_MAGIC[testnet]
            error("magic is not right ", magic, " vs ", NETWORK_MAGIC[testnet])
        end
        command = strip(String(read(s, 12)), '\0')
        payload_length = reinterpret(UInt32, read(s, 4))[1]
        checksum = reinterpret(UInt32, read(s, 4))[1]
        while bytesavailable(s) < payload_length
            println("Awaiting data from remote node...")
            sleep(0.02)
        end
        payload = read(s, payload_length)
        calculated_checksum = reinterpret(UInt32, hash256(payload)[1:4])[1]
        if calculated_checksum != checksum
            println("\n\nResult so far was : ", result,
                    "\nLast message had command ", command,
                    ", payload length ", payload_length,
                    " and payload ", bytes2hex(payload))
            error("Error parsing IO, calculated checksum does not match ", calculated_checksum, " vs ", checksum)
        end
        push!(result, NetworkEnvelope(magic, command, payload_length, checksum, payload))
    end
    result
end

function io2envelopes(bin::Vector{UInt8}, testnet::Bool=false)
    s = IOBuffer(bin)
    io2envelopes(s, testnet)
end

function io2envelope(s::IOBuffer, testnet::Bool=false)
    if !(bytesavailable(s) < 24)
        magic = reinterpret(UInt32, read(s, 4))[1]
        if magic != NETWORK_MAGIC[testnet]
            error("Network magic is not right!\nWe have received ", magic, " versus ", NETWORK_MAGIC[testnet])
        end
        command = strip(String(read(s, 12)), '\0')
        payload_length = reinterpret(UInt32, read(s, 4))[1]
        checksum = reinterpret(UInt32, read(s, 4))[1]
        while bytesavailable(s) < payload_length
            println("Buffer doesn't hold sufficient data to parse payoad, waiting for 5 s...")
            sleep(5)
        end
        payload = read(s, payload_length)
        calculated_checksum = reinterpret(UInt32, hash256(payload)[1:4])[1]
        if calculated_checksum != checksum
            println("\n\nResult so far was : ", result,
            "\nLast message had command ", command,
            ", payload length ", payload_length,
            " and payload ", bytes2hex(payload))
            error("Error parsing IO, calculated checksum does not match ", calculated_checksum, " vs ", checksum)
        end
        return NetworkEnvelope(magic, command, payload_length, checksum, payload)
    end
end

"""
Returns the byte serialization of the entire network message
"""
function serialize(envelope::NetworkEnvelope)
    result = Vector(reinterpret(UInt8, [htol(envelope.magic)]))
    append!(result, UInt8.(collect(envelope.command)))
    append!(result, fill(0x00, (12 - length(envelope.command))))
    append!(result, bytes(length(envelope.payload), len=4, little_endian=true))
    append!(result, bytes(envelope.checksum, len=4, little_endian=true))
    append!(result, envelope.payload)
end

struct Peer
    time::UInt32
    services::UInt64
    ip::IPAddr
    port::UInt16
    Peer(services::Integer, ip::IPAddr, port::Integer, time::Integer=0) = new(time, services, ip, port)
end

Peer(testnet::Bool=false) = Peer(DEFAULT["services"], DEFAULT["ip"], DEFAULT["port"][testnet])
Peer(ip::IPAddr, port::Integer) = Peer(DEFAULT["services"], ip, port)

function show(io::IO, z::Peer)
    print(io, z.services, "@", z.ip, ":", z.port)
end

"""
    IPv4 -> Vector{UInt8}

Return an 16 bytes UInt8 array representing the IP address

IPv6 address. Network byte order. The original client only supported IPv4 and
only read the last 4 bytes to get the IPv4 address. However, the IPv4 address
is written into the message as a 16 byte IPv4-mapped IPv6 address

(12 bytes 00 00 00 00 00 00 00 00 00 00 FF FF, followed by the 4 bytes of the
IPv4 address).
"""
function ip2bytes(ip::IPv4)
    result = fill(0x00, 10)
    append!(result, [0xff, 0xff])
    append!(result, Vector(reinterpret(UInt8, [hton(ip.host)])))
end

function bytes2ip(bin::Vector{UInt8})
    if bin[1:12] == IPV4_PREFIX
        IPv4(ntoh(reinterpret(UInt32, bin[13:16])[1]))
    else
        IPv6(ntoh(reinterpret(UInt128, bin)[1]))
    end
end

"""
    Peer -> Vector{UInt8}

Returns the serialization of a Peer
"""
function serialize(peer::Peer, versionmessage::Bool=false)
    result = UInt8[]
    if !versionmessage
        append!(result, Vector(reinterpret(UInt8, [htol(peer.time)])))
    end
    append!(result, Vector(reinterpret(UInt8, [htol(peer.services)])))
    append!(result, ip2bytes(peer.ip))
    append!(result, Vector(reinterpret(UInt8, [hton(peer.port)])))
end

"""
    parse_peer(payload::Vector{UInt8}) -> Peer

Parse Peer from bytes arrays
"""
function parse_peer(payload::Vector{UInt8}, versionmessage::Bool=false)
    io = IOBuffer(payload)
    if versionmessage
        time = ltoh(reinterpret(UInt32, read(io, 4))[1])
    else
        time = 0
    end
    services = ltoh(reinterpret(UInt64, read(io, 8))[1])
    ip = bytes2ip(read(io, 16))
    port = ntoh(reinterpret(UInt16, read(io, 8))[1])
    Peer(services, ip, port, time)
end

struct VersionMessage <: AbstractMessage
    command::String
    version::UInt32
    services::UInt64
    timestamp::UInt64
    receiver::Peer
    sender::Peer
    nonce::UInt64
    user_agent::String
    start_height::UInt32
    relay::Bool
    VersionMessage(version::Integer, services::Integer, timestamp::Integer, receiver::Peer, sender::Peer, nonce::Integer, user_agent::String, start_height::Integer, relay::Bool) = new("version", version, services, timestamp, receiver, sender, nonce, user_agent, start_height, relay)
end

VersionMessage(timestamp::UInt64, nonce::UInt64, testnet::Bool=false) =
    VersionMessage(DEFAULT["version"],
                   DEFAULT["services"],
                   timestamp,
                   Peer(testnet), Peer(testnet),
                   nonce,
                   USER_AGENT,
                   DEFAULT["start_height"],
                   DEFAULT["relay"])
VersionMessage(testnet::Bool=false) =
    VersionMessage(UInt64(round(time())), rand(UInt64), testnet)

function show(io::IO, z::VersionMessage)
    print(io, "Version Message\n--------\nVersion : ",
            z.version, " - User Agent : ", String(z.user_agent),
            "\nServices : ", z.services, " - Relay : ", z.relay,
            "\nTime Stamp : ", unix2datetime(z.timestamp),
            "\nLatest Block : ", string(z.start_height, base=16),
            "\nReceiver : ", z.receiver,
            "\nSender : ", z.sender)
end

"""
    serialize(version::VersionMessage) -> Vector{UInt8}

Serialize this message to send over the network
"""
function serialize(version::VersionMessage)
    result = Vector(reinterpret(UInt8, [htol(version.version)]))
    append!(result, Vector(reinterpret(UInt8, [htol(version.services)])))
    append!(result, Vector(reinterpret(UInt8, [htol(version.timestamp)])))
    append!(result, serialize(version.receiver, true))
    append!(result, serialize(version.sender, true))
    append!(result, Vector(reinterpret(UInt8, [htol(version.nonce)])))
    append!(result, serialize(VarString(version.user_agent)))
    append!(result, Vector(reinterpret(UInt8, [htol(version.start_height)])))
    version.relay ? append!(result, [0x01]) : append!(result, [0x00])
    return result
end

"""
    serialize(version::VersionMessage) -> Vector{UInt8}

"""
function payload2version(payload::Vector{UInt8})
    io = IOBuffer(payload)
    version = ltoh(reinterpret(UInt32, read(io, 4))[1])
    services = ltoh(reinterpret(UInt64, read(io, 8))[1])
    timestamp = ltoh(reinterpret(UInt64, read(io, 8))[1])
    sender = parse_peer(read(io, 26))
    receiver = parse_peer(read(io, 26))
    nonce = ltoh(reinterpret(UInt64, read(io, 8))[1])
    user_agent = io2varstring(io).str
    start_height = ltoh(reinterpret(UInt32, read(io, 4))[1])
    relay = Bool(read(io, 1)[1])
    VersionMessage(version, services, timestamp, receiver, sender, nonce, user_agent, start_height, relay)
end

struct VerAckMessage <: AbstractMessage
    command::String
    VerAckMessage() = new("verack")
end

payload2verack(payload::Vector{UInt8}) = VerAckMessage()

serialize(::VerAckMessage) = UInt8[]

struct PingMessage <: AbstractMessage
    command::String
    nonce::Vector{UInt8}
    PingMessage(nonce) = new("ping", nonce)
end

payload2ping(payload::Vector{UInt8}) = PingMessage(payload)

serialize(ping::PingMessage) = ping.nonce

struct PongMessage <: AbstractMessage
    command::String
    nonce::Vector{UInt8}
    PongMessage(nonce) = new("pong", nonce)
end

payload2pong(payload::Vector{UInt8}) = PongMessage(payload)

serialize(pong::PongMessage) = pong.nonce

struct GetHeadersMessage <: AbstractMessage
    command::String
    version::UInt32
    num_hashes::Integer
    start_block::Vector{UInt8}
    end_block::Vector{UInt8}
    GetHeadersMessage(version::Integer, num_hashes::Integer, start_block::Vector{UInt8}, end_block::Vector{UInt8}=fill(0x00, 32)) = new("getheaders", version, num_hashes, start_block, end_block)
end

GetHeadersMessage(start_block::Vector{UInt8}) = GetHeadersMessage(DEFAULT["version"], 1, start_block)

"""
Serialize this message to send over the network
    protocol version is 4 bytes little-endian
    number of hashes is a varint
    start block is in little-endian
    end block is also in little-endian
"""
function serialize(getheaders::GetHeadersMessage)
    result = bytes(getheaders.version, len=4, little_endian=true)
    append!(result, encode_varint(getheaders.num_hashes))
    append!(result, reverse!(copy(getheaders.start_block)))
    append!(result, reverse!(copy(getheaders.end_block)))
    return result
end

struct HeadersMessage <: AbstractMessage
    command::String
    headers::Vector{BlockHeader}
    HeadersMessage(headers::Vector{BlockHeader}) = new("headers", headers)
end

"""
    # number of headers is in a varint
    # initialize the headers array
    # loop through number of headers times
    # add a header to the headers array by parsing the stream
    # read the next varint (num_txs)
    # num_txs should be 0 or raise a RuntimeError
"""
function payload2headers(payload::Vector{UInt8})
    io = IOBuffer(payload)
    payload2headers(io)
end

function payload2headers(io::IOBuffer)
    num_headers = read_varint(io)
    headers = BlockHeader[]
    for i in 1:num_headers
        push!(headers, io2blockheader(io))
        num_txs = read_varint(io)
        if num_txs != 0
            error("number of txs not 0")
        end
    end
    return HeadersMessage(headers)
end

mutable struct GetDataMessage <: AbstractMessage
    command::String
    data::Vector{Tuple{Integer,Vector{UInt8}}}
    GetDataMessage(data::Vector{Tuple{Integer,Vector{UInt8}}}=Tuple{Integer,Vector{UInt8}}[]) = new("getdata", data)
end

import Base.append!

function append!(x::GetDataMessage, type::Integer, identifier::Vector{UInt8})
    push!(x.data, (type, identifier))
end

function serialize(x::GetDataMessage)
    result = encode_varint(length(x.data))
    for e in x.data
        append!(result, bytes(e[1], len=4, little_endian=true))
        append!(result, reverse!(copy(e[2])))
    end
    return result
end

mutable struct RejectMessage <: AbstractMessage
    command::String
    message::String
    ccode::Char
    reason::String
    data::Vector{UInt8}
    RejectMessage(message::String, ccode::Char, reason::String, data::Vector{UInt8}) = new("reject", message, ccode, reason, data)
end

function show(io::IO, z::RejectMessage)
    print(io, "Reject Message\n--------\nMessage : ",
            z.message, "\nccode : ", z.ccode,
            "\nReason : ", z.reason, "\n", bytes2hex(z.data))
end

"""
    payload2reject(payload::Vector{UInt8}) -> RejectMessage

Parse RejectMessage from NetworkEnvelope payload
"""
function payload2reject(payload::Vector{UInt8})
    io = IOBuffer(payload)
    message = io2varstring(io).str
    ccode = Char(read(io, 1)[1])
    reason = io2varstring(io).str
    data = read(io)
    RejectMessage(message, ccode, reason, data)
end

mutable struct MerkleBlockMessage <: AbstractMessage
    command::String
    header::BlockHeader
    tx_count::UInt32
    hash_count::Unsigned
    hashes::Vector{Vector{UInt8}}
    flag_byte_count::Unsigned
    flags::Vector{Bool}
    MerkleBlockMessage(header::BlockHeader, tx_count::Integer,
                  hash_count, hashes::Vector{Vector{UInt8}}, flag_byte_count,
                  flags::Vector{Bool}) = new("merkleblock",
                  header, tx_count, hash_count, hashes, flag_byte_count, flags)
end

function show(io::IO, z::MerkleBlockMessage)
    print(io, "Merkle Block Message\n--------\nMessage : ",
            z.header, "\ntx_count : ", z.tx_count,
            " flag_byte_count : ", z.flag_byte_count)
end

"""
    payload2merkleblock(payload::Vector{UInt8}) -> MerkleBlockMessage

Parse MerkleBlockMessage from NetworkEnvelope payload
"""
function payload2merkleblock(payload::Vector{UInt8})
    io = IOBuffer(payload)
    header = io2blockheader(io)
    tx_count = ltoh(reinterpret(UInt32, read(io, 4))[1])
    hash_count = read_varint(io)
    hashes = Vector{UInt8}[]
    for i in 1:hash_count
         push!(hashes, read(io, 32))
    end
    flag_byte_count = read_varint(io)
    flag_byte = read(io, flag_byte_count)
    flags = bytes2flags(flag_byte)
    MerkleBlockMessage(header, tx_count, hash_count, hashes, flag_byte_count, flags)
end

"""

Returns true if MerkleBlockMessage is valid
"""
function is_valid(mb::MerkleBlockMessage)
    tree = MerkleTree(mb.tx_count)
    populate!(tree, mb.flags, mb.hashes)
    root(tree) == mb.header.merkle_root
end

struct FilterAddMessage <: AbstractMessage
    command::String
    element_bytes::Unsigned
    element::Vector{UInt8}
    FilterAddMessage(element_bytes::Unsigned, element::Vector{UInt8}) = new("filteradd", element_bytes, element)
end

struct FilterClearMessage <: AbstractMessage
    command::String
    FilterClearMessage() = new("filterclear")
end

struct FilterLoadMessage <: AbstractMessage
    command::String
    n_filter_bytes::Unsigned
    filter::Vector{Bool}
    n_hash_funcs::UInt32
    n_tweak::UInt32
    n_flags::UInt8
    FilterLoadMessage(n_filter_bytes::Unsigned, filter::Vector{Bool}, n_hash_funcs::UInt32, n_tweak::UInt32, n_flags::UInt8) = new("filterload", n_filter_bytes, filter, n_hash_funcs, n_tweak, n_flags)
end

FilterLoadMessage(bf::BloomFilter, flag::UInt8=0x01) = FilterLoadMessage(bf.size, bf.bit_field, bf.function_count, bf.tweak, flag)

function serialize(msg::FilterLoadMessage)
    payload = encode_varint(msg.n_filter_bytes)
    append!(payload, flags2bytes(msg.filter))
    append!(payload, Vector(reinterpret(UInt8, [htol(msg.n_hash_funcs)])))
    append!(payload, Vector(reinterpret(UInt8, [htol(msg.n_tweak)])))
    append!(payload, msg.n_flags%UInt8)
end

struct SendHeadersMessage <: AbstractMessage
    command::String
    SendHeadersMessage() = new("sendheaders")
end

SendHeadersMessage(::Any) = SendHeadersMessage()

struct InventoryVector
    type::Integer
    hash::Vector{UInt8}
    InventoryVector(type, hash) = new(type, hash)
end

function show(io::IO, z::InventoryVector)
    print(io, DATA_MESSAGE_NAME[z.type],
          " :\n", bytes2hex(z.hash))
end

function payload2inventoryvector(payload::Vector{UInt8})
    type = ltoh(reinterpret(UInt32, payload[1:4])[1])
    hash = payload[5:26]
    InventoryVector(type, hash)
end

struct InvMessage <: AbstractMessage
    command::String
    count::Integer
    inventory::Vector{InventoryVector}
    InvMessage(inv_count::Integer, inventory::Vector{InventoryVector}) = new("inv", inv_count, inventory)
end

function payload2inv(payload::Vector{UInt8})
    io = IOBuffer(payload)
    inv_count = read_varint(io)
    inventory = InventoryVector[]
    for i ∈ 1:inv_count
        entry = payload2inventoryvector(read(io, 26))
        push!(inventory, entry)
    end
    InvMessage(inv_count, inventory)
end

struct SendCmpctMessage <: AbstractMessage
    command::String
    announce::Bool
    version::UInt64
    SendCmpctMessage(announce::Bool, version::UInt64) = new("sendcmpct", announce, version)
end

function payload2sendcmpct(payload::Vector{UInt8})
    annouce = Bool(payload[1])
    version = ltoh(reinterpret(UInt64, payload[2:9])[1])
    SendCmpctMessage(annouce, version)
end

struct FeeFilterMessage <: AbstractMessage
    command::String
    feerate::UInt64
    FeeFilterMessage(feerate::UInt64) = new("feefilter", feerate)
end

function payload2feefilter(payload::Vector{UInt8})
    feerate = ltoh(reinterpret(UInt64, payload[1:8])[1])
    SendCmpctMessage(feerate)
end

struct AddrMessage <: AbstractMessage
    n::Unsigned
    addr::Vector{IPAddr}
    FeeFilterMessage(n::Integer, addr::Vector{IPAddr}) = new("addr", n, addr)
end

function payload2addr(payload::Vector{UInt8})
    io = IOBuffer(payload)
    n = read_varint(io)
    addr = IPAddr[]
    for _ ∈ 1:n
        push!(addr, bytes2ip(read(io, 16)))
    end
    AddrMessage(n, addr)
end

PARSE_PAYLOAD = Dict([
    ("addr", payload2addr),
    ("feefilter", payload2feefilter),
    ("headers", payload2headers),
    ("inv", payload2inv),
    ("merkleblock", payload2merkleblock),
    ("ping", payload2ping),
    ("pong", payload2pong),
    ("reject", payload2reject),
    ("sendcmpct", payload2sendcmpct),
    ("sendheaders", SendHeadersMessage),
    ("tx", payload2tx),
    ("verack", payload2verack),
    ("version", payload2version)
])


text:
function Kprec( pw::PWGrid, psi::Array{ComplexF64,2} )

    Ngwx  = size(psi)[1]
    Nstates = size(psi)[2]
    G2    = pw.gvec.G2[pw.gvecw.idx_gw2r]
    Kpsi  = zeros( ComplexF64, Ngwx, Nstates )

    for ist = 1:Nstates
        for ig = 1:Ngwx
            Kpsi[ig,ist] = psi[ig,ist] / ( 1.0 + G2[ig] )
        end
    end
    return Kpsi
end

# For comparison with non-preconditioned
function Kprec(psi)
    return psi
end


text:
push!(LOAD_PATH, "../src/")
using Documenter, KiteConnect

makedocs(sitename="KiteConnect")


text:
using Stheno: CompositeMean, CompositeKernel, CompositeCrossKernel, ConstantMean,
    ZeroMean, LhsCross, RhsCross, OuterCross, OuterKernel, MeanFunction, Kernel,
    EmpiricalMean, pairwise

@testset "compose" begin

    # Test CompositeMean.
    let
        rng, N, D = MersenneTwister(123456), 100, 2
        μ1, μ2 = EmpiricalMean(randn(rng, N)), EmpiricalMean(randn(rng, N))
        X = eachindex(μ1)
        ν1, ν2 = CompositeMean(+, μ1, μ2), CompositeMean(*, μ1, μ2)
        @test ν1(X[1]) == μ1(X[1]) + μ2(X[1])
        @test ν2(X[3]) == μ1(X[3]) .* μ2(X[3])
        @test length(ν1) == N && length(ν2) == N
        @test eachindex(ν1) == eachindex(X)

        mean_function_tests(ν1, X)
        mean_function_tests(ν2, X)
        mean_function_tests(ν1, BlockData([X, X]))
        mean_function_tests(ν2, BlockData([X, X]))

        @test map(ν1, :) == μ1.μ .+ μ2.μ
        @test map(ν2, :) == μ1.μ .* μ2.μ
    end

    # Test CompositeKernel and CompositeCrossKernel.
    let
        rng, N, N′, D = MersenneTwister(123456), 5, 6, 2
        x0, x1, x2 = randn(rng, N), randn(rng, N), randn(rng, N′)
        X0, X1, X2 = ColsAreObs(randn(rng, D, N)), ColsAreObs(randn(rng, D, N)), ColsAreObs(randn(rng, D, N′))

        # Define base kernels and do composition.
        k1, k2 = EQ(), Linear(randn(rng))
        ν1, ν2 = CompositeKernel(+, k1, k2), CompositeCrossKernel(*, k1, k2)
        @test ν1(X0[1], X1[2]) == k1(X0[1], X1[2]) + k2(X0[1], X1[2])
        @test ν2(X0[1], X1[2]) == k1(X0[1], X1[2]) .* k2(X0[1], X1[2])
        @test size(ν1, 1) == Inf && size(ν1, 2) == Inf
        @test size(ν2, 1) == Inf && size(ν2, 2) == Inf
        @test !isstationary(ν1) && !isstationary(ν2)
        @test isstationary(CompositeKernel(+, k1, k1)) == true

        xb0, xb1, xb2 = BlockData([x0, x1]), BlockData([x1, x0]), BlockData([x2, x2])
        XB0, XB1, XB2 = BlockData([X0, X1]), BlockData([X1, X0]), BlockData([X2, X2])

        cross_kernel_tests(ν2, x0, x1, x2)
        cross_kernel_tests(ν2, X0, X1, X2)
        cross_kernel_tests(ν2, xb0, xb1, xb2)
        cross_kernel_tests(ν2, XB0, XB1, XB2)

        kernel_tests(ν1, x0, x1, x2)
        kernel_tests(ν1, X0, X1, X2)
        kernel_tests(ν1, xb0, xb1, xb2)
        kernel_tests(ν1, XB0, XB1, XB2)
    end

    # Test LhsCross functionality.
    let
        rng, N, N′, D = MersenneTwister(123456), 5, 6, 2
        x0, x1, x2 = randn(rng, N), randn(rng, N), randn(rng, N′)
        X0, X1, X2 = ColsAreObs(randn(rng, D, N)), ColsAreObs(randn(rng, D, N)), ColsAreObs(randn(rng, D, N′))

        f, k = x->sum(abs2, x), EQ()
        ν = LhsCross(f, k)

        @test ν(X0[1], X0[2]) == f(X0[1]) * k(X0[1], X0[2])
        @test size(ν, 1) == Inf && size(ν, 2) == Inf

        xb0, xb1, xb2 = BlockData([x0, x1]), BlockData([x1, x0]), BlockData([x2, x2])
        XB0, XB1, XB2 = BlockData([X0, X1]), BlockData([X1, X0]), BlockData([X2, X2])
        cross_kernel_tests(ν, x0, x1, x2)
        cross_kernel_tests(ν, X0, X1, X2)
        cross_kernel_tests(ν, xb0, xb1, xb2)
        cross_kernel_tests(ν, XB0, XB1, XB2)
    end

    # Test RhsCross functionality.
    let
        rng, N, N′, D = MersenneTwister(123456), 5, 6, 2
        x0, x1, x2 = randn(rng, N), randn(rng, N), randn(rng, N′)
        X0, X1, X2 = ColsAreObs(randn(rng, D, N)), ColsAreObs(randn(rng, D, N)), ColsAreObs(randn(rng, D, N′))

        f, k = x->sum(abs2, x), EQ()
        ν = RhsCross(k, f)

        @test ν(X0[1], X0[2]) == k(X0[1], X0[2]) * f(X0[2])
        @test size(ν, 1) == Inf && size(ν, 2) == Inf

        xb0, xb1, xb2 = BlockData([x0, x1]), BlockData([x1, x0]), BlockData([x2, x2])
        XB0, XB1, XB2 = BlockData([X0, X1]), BlockData([X1, X0]), BlockData([X2, X2])
        cross_kernel_tests(ν, x0, x1, x2)
        cross_kernel_tests(ν, X0, X1, X2)
        cross_kernel_tests(ν, xb0, xb1, xb2)
        cross_kernel_tests(ν, XB0, XB1, XB2)
    end

    # Test OuterCross.
    let
        rng, N, N′, D = MersenneTwister(123456), 5, 6, 2
        x0, x1, x2 = randn(rng, N), randn(rng, N), randn(rng, N′)
        X0, X1, X2 = ColsAreObs(randn(rng, D, N)), ColsAreObs(randn(rng, D, N)), ColsAreObs(randn(rng, D, N′))

        f, k = x->sum(abs2, x), EQ()
        ν = OuterCross(f, k)

        @test ν(X0[1], X0[2]) == f(X0[1]) * k(X0[1], X0[2]) * f(X0[2])
        @test size(ν, 1) == Inf && size(ν, 2) == Inf

        xb0, xb1, xb2 = BlockData([x0, x1]), BlockData([x1, x0]), BlockData([x2, x2])
        XB0, XB1, XB2 = BlockData([X0, X1]), BlockData([X1, X0]), BlockData([X2, X2])
        cross_kernel_tests(ν, x0, x1, x2)
        cross_kernel_tests(ν, X0, X1, X2)
        cross_kernel_tests(ν, xb0, xb1, xb2)
        cross_kernel_tests(ν, XB0, XB1, XB2)
    end

    # Test OuterKernel.
    let
        rng, N, N′, D = MersenneTwister(123456), 5, 6, 2
        x0, x1, x2 = randn(rng, N), randn(rng, N), randn(rng, N′)
        X0, X1, X2 = ColsAreObs(randn(rng, D, N)), ColsAreObs(randn(rng, D, N)), ColsAreObs(randn(rng, D, N′))

        f, k = x->sum(abs2, x), EQ()
        ν = OuterKernel(f, k)

        @test ν(X0[1], X0[2]) == f(X0[1]) * k(X0[1], X0[2]) * f(X0[2])
        @test size(ν, 1) == Inf && size(ν, 2) == Inf

        xb0, xb1, xb2 = BlockData([x0, x1]), BlockData([x1, x0]), BlockData([x2, x2])
        XB0, XB1, XB2 = BlockData([X0, X1]), BlockData([X1, X0]), BlockData([X2, X2])
        kernel_tests(ν, x0, x1, x2)
        kernel_tests(ν, X0, X1, X2)
        cross_kernel_tests(ν, xb0, xb1, xb2)
        cross_kernel_tests(ν, XB0, XB1, XB2)
    end

    # Test mean function composition.
    let
        rng, N, D = MersenneTwister(123456), 5, 2
        X, s = ColsAreObs(randn(rng, D, N)), randn(rng)
        μ, μ′ = ConstantMean(randn(rng)), CustomMean(x->sum(s .* x))

        # Test conversion and promotion.
        c = randn(rng)
        @test convert(MeanFunction, c) == ConstantMean(c)
        @test promote_rule(ConstantMean, Real) == MeanFunction
        @test promote(ConstantMean(c), c) == (ConstantMean(c), ConstantMean(c))

        # Test addition.
        @test μ + μ′ == CompositeMean(+, μ, μ′)
        @test map(μ + μ′, X) == map(μ, X) + map(μ′, X)
        @test map(μ + 5, X) == map(μ, X) .+ 5
        @test 2.34 + μ′ == CompositeMean(+, ConstantMean(2.34), μ′)

        # Test multiplication.
        @test μ * μ′ == CompositeMean(*, μ, μ′)
        @test map(μ * μ′, X) == map(μ, X) .* map(μ′, X)
        @test map(μ * 4.32, X) == 4.32 .* map(μ, X)
        @test 4.23 * μ′ == CompositeMean(*, ConstantMean(4.23), μ′)
    end

    # Test kernel composition.
    let
        rng, N, D = MersenneTwister(123456), 5, 2
        X, s = ColsAreObs(randn(rng, N, D)), randn(rng)
        k, k′ = EQ(), Linear(1)

        # Test conversion and promotion.
        @test convert(Kernel, s) == ConstantKernel(s)
        @test promote_rule(typeof(k), typeof(s)) == Kernel
        @test promote(k, s) == (k, ConstantKernel(s))
        @test promote(s, k) == (ConstantKernel(s), k)

        # Test addition.
        @test k + k′ == CompositeKernel(+, k, k′)
        @test k + s == CompositeKernel(+, k, ConstantKernel(s))
        @test s + k′ == CompositeKernel(+, ConstantKernel(s), k′)

        # Test multiplication.
        @test k * k′ == CompositeKernel(*, k, k′)
        @test k′ * s == CompositeKernel(*, k′, ConstantKernel(s))
        @test s * k == CompositeKernel(*, ConstantKernel(s), k)
    end
end


text:
############################################################
# Naive Vector&Dict implementation (as fast as manual loop),
############################################################
struct Orbit1{T, S} <: AbstractOrbit{T,S}
    elts::Vector{T}
    vals::Dict{T, S}
end
Orbit1(pt::T, val::S) where {T,S} = Orbit1{T,S}([pt], Dict(pt=>val))
Orbit1(pt::T) where T = Orbit1(pt, nothing)
# Orbit1(::Type{T}, ::Type{S}) where {T,S} = Orbit1(T[], Dict{T,S}())

function Orbit1(v::AbstractVector{T}) where T
    return Orbit1(Vector(v), Dict(i=>nothing for i in v))
end

@inline Base.in(x, orb::Orbit1) = haskey(orb.vals, x)
@inline function Base.push!(orb::Orbit1{T,S}, tup::Tuple{T,S}) where {T,S}
    push!(orb.elts, tup[1])
    orb.vals[tup[1]] = tup[2]
    return orb
end

@inline function Base.push!(orb::Orbit1{T,Nothing}, x::T) where {T}
    push!(orb.elts, x)
    orb.vals[x] = nothing
    return orb
end

@inline Base.getindex(orb::Orbit1, o) = orb.vals[o]

@inline Base.iterate(orb::Orbit1) = iterate(orb.elts)
@inline Base.iterate(orb::Orbit1, s) = iterate(orb.elts, s)
@inline Base.length(orb::Orbit1) = length(orb.elts)
@inline Base.:(==)(o1::Orbit1, o2::Orbit1) =
    o1.elts == o2.elts && o1.vals == o2.vals

@inline Base.first(orb::Orbit1) = orb.elts[1]
@inline Base.last(orb::Orbit1) = orb.elts[end]
@inline islast(orb::Orbit1, pt) = pt == last(orb)


text:
#=
Shows

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date: 2014-08-05

References
----------

http://quant-econ.net/career.html
=#
using PyPlot
using QuantEcon, QuantEcon.Models

# Quick meshgrid function
meshgrid(x::Vector, y::Vector) = (repmat(x, 1, length(y))',
                                  repmat(y, 1, length(x)))

srand(41)  # reproducible results
wp = CareerWorkerProblem()
v_init = fill(100.0, wp.N, wp.N)
func(x) = bellman_operator(wp, x)
v = compute_fixed_point(func, v_init, max_iter=500, verbose=false)

fig = figure(figsize=(8,6))
ax = fig[:add_subplot](111, projection="3d")
tg, eg = meshgrid(wp.theta, wp.epsilon)

ax[:plot_surface](tg, eg, v', rstride=2, cstride=2,
                cmap=ColorMap("jet"),
                alpha=0.5,
                linewidth=0.25)

ax[:set_zlim](150, 200)
ax[:set_xlabel]("theta", fontsize=14)
ax[:set_ylabel]("epsilon", fontsize=14)
plt.show()


text:
infile = "./src/Day16/input.txt"
input = readlines(infile)

function parse_input(infile)
    rules, myticket, nearby = open(infile) do f
        line = 1
        part = 1
        rules = Dict{String,Vector{UnitRange{Int64}}}()
        nearby = Vector{Vector{Int}}()
        while !eof(f)
            x = readline(f)
            if x == ""
                part += 1
                continue
            elseif x == "your ticket:"
                part = 2
                continue
            elseif x == "nearby tickets:"
                part = 3
                continue
            end
            if part == 1
                rulename, textranges = split(x, ": ")
                ranges = parserange.(split(textranges, " or "))
                rules[rulename] = ranges
            end
            if part == 2
                myticket = parse.(Int, split(x, ","))
            end
            if part == 3
                thisticket = parse.(Int, split(x, ","))
                push!(nearby, thisticket)
            end
        end
        return rules, myticket, nearby
    end
    return rules, myticket, nearby
end

function parserange(txt)
    left, right = parse.(Int, split(txt, "-"))
    return left:right
end

function find_invalid_values(tickets, rules)
    invalids = Vector{Int}()
    invalidtickets = Set{Int}()
    for (num, ticket) in enumerate(tickets)
        ticketvalid = false
        for (f,field) in enumerate(ticket)
            valid = false
            for rule in values(rules)
                if any(field .∈ rule)
                    valid = true
                    break
                end
            end
            if !valid
                push!(invalids, field)
                push!(invalidtickets, num)
                # println("Field $f on ticket $num invalid")
            end
        end
    end
    return invalids, invalidtickets
end

rules, myticket, nearby = parse_input(infile)

badfields, badtickets = find_invalid_values(nearby, rules)
validnearby = [nearby[x] for x in 1:length(nearby) if x ∉ badtickets]

# assign field to column
ruletocol = Dict{String, Set{Int}}()
coltorule = Dict{Int, Set{String}}()
for (rulename, ranges) in rules
    for col in 1:length(validnearby[1])
        valid = all([any(x[col] .∈ ranges) for x in validnearby])
        if valid
            if haskey(ruletocol, rulename)
                push!(ruletocol[rulename], col)
            else
                ruletocol[rulename] = Set(col)
            end
            if haskey(coltorule, col)
                push!(coltorule[col], rulename)
            else
                coltorule[col] = Set((rulename,))
            end
        end
    end
end

finalassignments = Dict{String, Int}()
# make logical assignments
# while any(length(values(coltorule)) .> 1)
for i in 1:20
    for (col, colrules) in coltorule
        if length(colrules) == 1
            # then this assignment is final and cannot be used elsewhere
            finalassignments[only(colrules)] = col
            # run through all other rules and remove this link
            for (othercol, otherrules) in coltorule
                if othercol != col
                    delete!(coltorule[othercol], only(colrules))
                end
            end
        end
    end
end

departurecols = [col for (name, col) in finalassignments if contains(name, "departure")]
departurevals = myticket[departurecols]
part2 = prod(departurevals)


text:
#using Plots
#gr()
module MyModule

    export plus
    function plus(x,y)
        z = x + y
        return z
    end

    function minus(x,y)
        z = x - y
        return z
    end

    function Print()
        for i=1:6
            println(Data[i])
        end
    end

    function OverWrite()
        Num = 1.
        println(Num)
    end

    function Pow(x,y)
        z = 1
        for i=1:y
            z = mul(z,x)
        end
        return z
    end

    function Pl(u,v)
        println(u)
        println(v)
        plot(u,v)
        savefig("test3.png")
        tles = read_tle("ISS_TLE.txt") #TLE読み込み
        println(tles.ω)
    end
end

max_stars_repo_name: yui-project/ACSim

max_issues_repo_name: yui-project/ACSim

max_forks_repo_name: yui-project/ACSim


text:
fname = "../externals/Core.Math.Samples/data/Pejcic_318.csv"
data = readdlm(fname, separator = ',')

data_read_02 = CSV.read("Pejcic_318.csv", header=true, delim=',', copycols = true)


text:
abstract type EigenVectorCalculation end
abstract type ShiftStrategy end

struct WithEigenVector <: EigenVectorCalculation end
struct WithoutEigenVector <: EigenVectorCalculation end
struct OnlyFirstRowEigenVector <: EigenVectorCalculation end

struct NoShift <: ShiftStrategy end
struct Overrelaxation <: ShiftStrategy end
struct CloseEigenValue <: ShiftStrategy end

struct TqrEigenDecomposition
  iter::Int
  d::Vector{Float64}
  ev::Matrix{Float64}
end

function TqrEigenDecomposition(dg::Vector{Float64}, subVec::Vector{Float64}, calc::EigenVectorCalculation, strategy::ShiftStrategy)
  iter = 0
  d = dg
  ev_cols = get_ev_cols(calc)
  n = length(d)
  ev = zeros(n, ev_cols)
  e = copy(subVec)
  # insert!(e, 1, 0.0)
  pushfirst!(e, 0.0)

  for i = 1:ev_cols
    ev[i, i] = 1.0
  end

  @inbounds for k = n:-1:2
    while ~off_diag_is_zero(k, e, d)
      l = k
      l -= 1
      while l > 1 && ~off_diag_is_zero(l, e, d)
        l -= 1
      end

      iter += 1

      q = d[l]

      t1 = sqrt(0.25 * (d[k] * d[k] + d[k - 1] * d[k - 1]) - 0.5 * d[k - 1] * d[k] + e[k] * e[k])
      t2 = 0.5 * (d[k] + d[k - 1])

      lambda = abs(t2 + t1 - d[k]) < abs(t2 - t1 - d[k]) ? t2 + t1 : t2 - t1

      q -= (k == n ? 1.25 : 1.0) * lambda

      # QR transformation
      sine = 1.0
      cosine = 1.0
      u = 0.0

      recoverUnderflow = false

      @inbounds for i = l+1:k
        if recoverUnderflow
          break
        end

        h = cosine * e[i]
        p = sine * e[i]

        e[i - 1] = sqrt(p * p + q * q)
        if e[i - 1] != 0.0
          sine = p / e[i - 1]
          cosine = q / e[i - 1]

          g = d[i - 1] - u
          t = (d[i] - g) * sine + 2.0 * cosine * h

          u = sine * t
          d[i - 1] = g + u
          q = cosine * t - h

          @inbounds for j = 1:ev_cols
            tmp = ev[i - 1, j]
            ev[i - 1, j] = sine * ev[i, j] + cosine * tmp
            ev[i, j] = cosine * ev[i, j] - sine * tmp
          end
        else
          # recover from underflow
          d[i - 1] -= u
          e[l] = 0.0
          recoverUnderflow = true
        end
      end

      if ~recoverUnderflow
        d[k] -= u
        e[k] = q
        e[l] = 0.0
      end
    end
  end

  # sort (eigenvalues, eigenvectors)
  ev = ev[sortperm(d, rev=true), :]
  sort!(d, rev=true)
  @simd for i = 1:n
    sgn = 1.0
    @inbounds if ev_cols > 0 && ev[i, 1] < 0.0
      sgn = -1.0
    end

    @inbounds for j = 1:ev_cols
      ev[i, j] = sgn * ev[i, j]
    end
  end

  return TqrEigenDecomposition(iter, d, ev)
end

get_ev_cols(::OnlyFirstRowEigenVector) = 1

off_diag_is_zero(k::Int, e::Vector{Float64}, d::Vector{Float64}) = abs(d[k - 1]) + abs(d[k]) == abs(d[k - 1]) + abs(d[k]) + abs(e[k])


text:
#
# Calculations related to CR3BP transfers
#

function closest_approach(man::Manifold{<:NormalizedSynodicCR3BPOrbit},
                          body_position::AbstractVector  = (secondary_synodic_position ∘ first ∘ first)(man);
                          max_velocity::Unitful.Velocity = 10u"km/s",
                          max_position::Unitful.Length   = 50_000u"km") where O <: NormalizedSynodicCR3BPOrbit

    @assert length(body_position) == 3 "The position vector provided must have 3 elements."
    orbit         = man |> first |> first
    max_velocity  = nondimensionalize_velocity(max_velocity, normalized_length_unit(orbit.system), normalized_time_unit(orbit.system))
    max_position  = nondimensionalize_length(max_position, normalized_length_unit(orbit))

    distance(orb)    = norm(position_vector(orb) - body_position)
    isrealistic(orb) = norm(position_vector(orb)) < max_position && norm(velocity_vector(orb)) < max_velocity

    orbits = filter(isrealistic, vcat(man...))

    !isempty(orbits) || throw(ArgumentError("No state within the manifold meets the `max_velocity` and `max_position` criteria! Use the keyword arguments to loosen the constraints."))
    return findmin(distance.(orbits))[2]

end

function optimal_approach(man::Manifold{<:NormalizedSynodicCR3BPOrbit},
                          body_position::AbstractVector = (secondary_synodic_position ∘ first ∘ first)(man),
                          cost_function::Function = (r,t)->√(t^2)) where O <: NormalizedSynodicCR3BPOrbit

    @assert length(body_position) == 3 "The position vector provided must have 3 elements."
    distance(orb)  = norm(position_vector(orb) - body_position)
    minimizer(orb) = cost_function(distance(orb), epoch(orb))
    max_distance = nondimensionalize_length(20_000u"km",
                                            normalized_length_unit(
                                                (first ∘ first)(man).system))
    orbits = filter(orbit->distance(orbit) ≤ max_distance, vcat(man...))
    getindex(orbits, findmin(minimizer.(orbits))[2])

end


text:
matmul2arg = [
(:*, :(A_mul_Bc(dy,x2)), :(Ac_mul_B(x1,dy))),
(:Ac_mul_B, :(A_mul_Bc(x2,dy)), :(x1*dy)),
(:A_mul_Bc, :(dy*x2), :(Ac_mul_B(dy,x1))),
(:Ac_mul_Bc,:(Ac_mul_Bc(x2,dy)), :(Ac_mul_Bc(dy,x1))),
(:At_mul_B, :(A_mul_Bt(x2,dy)), :(x1*dy)),
(:A_mul_Bt, :(dy*x2), :(At_mul_B(dy,x1))),
(:At_mul_Bt,:(At_mul_Bt(x2,dy)), :(At_mul_Bt(dy,x1))),
]

for (f,g1,g2) in matmul2arg
    @eval @primitive $f(x1,x2),dy,y  $g1  $g2
    addtest(f, rand(2,2), rand(2,2))
end

# Methods for multiplication:
# *(x::Float64, y::Float64) at float.jl:212  (same as .*)
# *(A::Number, B::AbstractArray{T,N}) at abstractarraymath.jl:54  (calls .*)
# *(A::AbstractArray{T,N}, B::Number) at abstractarraymath.jl:55  (calls .*)
# *(A::AbstractArray{T,2}, B::AbstractArray{T,1}) at linalg/matmul.jl:82
# *(A::AbstractArray{T,1}, B::AbstractArray{T,2}) at linalg/matmul.jl:89  (works only if size(B,1)==1)
# *(A::AbstractArray{T,2}, B::AbstractArray{T,2}) at linalg/matmul.jl:131
#
# The first three are handled by base/float.
# The final three implement matrix multiplication.

# TODO:

# eval
# arithtype: Not exported
# scale!
# scale
# vecdot
# dot
# Ac_mul_B
# At_mul_B
# *
# A_mul_B!
# At_mul_B!
# Ac_mul_B!
# A_mul_Bt
# A_mul_Bt!
# At_mul_Bt
# At_mul_Bt!
# A_mul_Bc
# A_mul_Bc!
# Ac_mul_Bc
# Ac_mul_Bc!
# Ac_mul_Bt: Not exported
# Ac_mul_Bt!: Not exported
# copytri!: Not exported
# gemv!: Not exported
# syrk_wrapper!: Not exported
# herk_wrapper!: Not exported
# gemm_wrapper: Not exported
# gemm_wrapper!: Not exported
# lapack_size: Not exported
# copy!
# copy_transpose!: Not exported
# generic_matvecmul!: Not exported
# generic_matmatmul: Not exported
# generic_matmatmul!: Not exported
# _generic_matmatmul!: Not exported
# matmul2x2: Not exported
# matmul2x2!: Not exported
# matmul3x3: Not exported
# matmul3x3!: Not exported


text:
#############################################################################
# backend
#

default_backend() = Backend(nothing)

#############################################################################
# vertex labels
#

function subscript_char(n::Integer)
  @assert 0 ≤ n ≤ 9
  return Char(0x2080 + n)
end

function subscript(n::Integer)
  ss = prod(Iterators.reverse((subscript_char(d) for d in digits(abs(n)))))
  if n < 0
    ss = "₋" * ss
  end
  return ss
end

default_vertex_labels_prefix(b::Backend, g) = "T"
function default_vertex_labels(
  b::Backend, g::AbstractGraph, vertex_labels_prefix=default_vertex_labels_prefix(b)
)
  return [string(vertex_labels_prefix, subscript(v)) for v in vertices(g)]
end

default_vertex_size(b::Backend, g) = 60
default_vertex_textsize(b::Backend, g) = 20

# TODO: customizable vertex marker
# nodeshapes="●", # ●, ▶, ◀, ■, █, ◩, ◪, ⧄, ⧅, ⦸, ⊘, ⬔, ⬕, ⬛, ⬤, 🔲, 🔳, 🔴, 🔵, ⚫
# edgeshapes="—", # ⇵, ⇶, ⇄, ⇅, ⇆, ⇇, ⇈, ⇉, ⇊, ⬱, —, –, ⟵, ⟶, ➖, −, ➡, ⬅, ⬆, ⬇

#############################################################################
# edge labels
#

default_edge_textsize(b::Backend) = 30

function default_edge_labels(b::Backend, g::AbstractGraph)
  return IndexLabels(b)
end

default_dims(b::Backend) = true
default_tags(b::Backend) = false
default_ids(b::Backend) = false
default_plevs(b::Backend) = false
default_qns(b::Backend) = false
default_newlines(b::Backend) = true

abstract type AbstractEdgeLabels end

(l::AbstractEdgeLabels)(g::AbstractGraph) = edge_labels(l, g)

struct IndexLabels <: AbstractEdgeLabels
  dims::Bool
  tags::Bool
  ids::Bool
  plevs::Bool
  qns::Bool
  newlines::Bool
end

IndexLabels(; kwargs...) = IndexLabels(Backend(); kwargs...)
IndexLabels(backend; kwargs...) = IndexLabels(Backend(backend); kwargs...)

function IndexLabels(
  b::Backend;
  dims=default_dims(b),
  tags=default_tags(b),
  ids=default_ids(b),
  plevs=default_plevs(b),
  qns=default_qns(b),
  newlines=default_newlines(b),
)
  return IndexLabels(dims, tags, ids, plevs, qns, newlines)
end

function edge_labels(b::Backend, l::IndexLabels, g::AbstractGraph)
  return edge_labels(l, g)
end

function edge_labels(l::IndexLabels, g::AbstractGraph)
  return String[edge_label(l, g, e) for e in edges(g)]
end

function edge_labels(b::Backend, params::NamedTuple, g::AbstractGraph)
  return IndexLabels(b; params...)(g)
end

function edge_label(l::IndexLabels, g::AbstractGraph, e)
  indsₑ = get_prop(g, e, :inds)
  return label_string(
    indsₑ;
    is_self_loop=is_self_loop(e),
    dims=l.dims,
    tags=l.tags,
    ids=l.ids,
    plevs=l.plevs,
    qns=l.qns,
    newlines=l.newlines,
  )
end

function _edge_label(l, g::Graph, e)
  return string(e)
end

edge_label(l::IndexLabels, g::Graph, e) = _edge_label(l, g, e)
edge_label(l, g::Graph, e) = _edge_label(l, g, e)

#function default_edge_labels(b::Backend, g; kwargs...)
#  return [edge_label(g, e; kwargs...) for e in edges(g)]
#end

plevstring(i::Index) = ITensors.primestring(plev(i))
idstring(i::Index) = string(id(i) % 1000)
tagsstring(i::Index) = string(tags(i))
qnstring(i::Index) = ""
function qnstring(i::QNIndex)
  str = "["
  for (n, qnblock) in pairs(space(i))
    str *= "$qnblock"
    if n ≠ lastindex(space(i))
      str *= ", "
    end
  end
  str *= "]"
  if dir(i) == ITensors.In
    str *= "†"
  end
  return str
end

function label_string(i::Index; dims, tags, plevs, ids, qns)
  str = ""
  if any((tags, plevs, ids, qns))
    str *= "("
  end
  if dims
    str *= string(dim(i))
  end
  if ids
    if dims
      str *= "|"
    end
    str *= idstring(i)
  end
  if tags
    if any((dims, ids))
      str *= "|"
    end
    str *= tagsstring(i)
  end
  if any((tags, plevs, ids, qns))
    str *= ")"
  end
  if plevs
    str *= plevstring(i)
  end
  if qns
    str *= qnstring(i)
  end
  return str
end

function label_string(is; is_self_loop=false, dims, tags, plevs, ids, qns, newlines)
  str = ""
  for n in eachindex(is)
    str *= label_string(is[n]; dims=dims, tags=tags, plevs=plevs, ids=ids, qns=qns)
    if n ≠ lastindex(is)
      if any((dims, tags, ids, qns))
        str *= "⊗"
      end
      if newlines && any((tags, ids, qns))
        str *= "\n"
      end
    end
  end
  return str
end

#############################################################################
# edge width
#

function width(inds)
  return log2(dim(inds)) + 1
end

function default_edge_widths(b::Backend, g::AbstractGraph)
  return Float64[width(get_prop(g, e, :inds)) for e in edges(g)]
end

function default_edge_widths(b::Backend, g::Graph)
  return [1.0 for e in edges(g)]
end

#############################################################################
# arrow
#

default_arrow_size(b::Backend, g) = 30

_hasqns(tn::Vector{ITensor}) = any(hasqns, tn)

function _hasqns(g::AbstractGraph)
  if iszero(ne(g))
    if has_prop(g, first(vertices(g)), :inds)
      return hasqns(get_prop(g, first(vertices(g)), :inds))
    else
      return hasqns(())
    end
  end
  return hasqns(get_prop(g, first(edges(g)), :inds))
end

_hasqns(g::Graph) = false

default_arrow_show(b::Backend, g) = _hasqns(g)

#############################################################################
# self-loop/siteinds direction
#

default_siteinds_direction(b::Backend, g) = Point2(0, -1)

#############################################################################
# dimensions
#

_ndims(::Any) = 2
_ndims(::NetworkLayout.AbstractLayout{N}) where {N} = N


text:
export tiny, huge
export TinyNumber, HugeNumber

import Base: show, maximum
import Base: convert, promote_rule

# Tiny number

"""
    TinyNumber <: Real

`TinyNumber` represents (wow!) tiny number that can be used in a various computations without unnecessary type promotions.

See also: [`HugeNumber`](@ref)
"""
struct TinyNumber <: Real end

Base.convert(::Type{ Float32 },  ::TinyNumber) = 1f-6
Base.convert(::Type{ Float64 },  ::TinyNumber) = 1e-12
Base.convert(::Type{ BigFloat }, ::TinyNumber) = big"1e-24"

Base.show(io::IO, ::TinyNumber) = print(io, "tiny")

Base.promote_rule(::Type{ TinyNumber }, ::Type{ I }) where { I <: Integer }  = promote_rule(TinyNumber, promote_type(I, Float64))
Base.promote_rule(::Type{ TinyNumber }, ::Type{ Float32 })                   = Float32
Base.promote_rule(::Type{ TinyNumber }, ::Type{ Float64 })                   = Float64
Base.promote_rule(::Type{ TinyNumber }, ::Type{ BigFloat })                  = BigFloat

# Huge number

"""
    HugeNumber <: Real

`HugeNumber` represents (wow!) huge number that can be used in a various computations without unnecessary type promotions.

See also: [`TinyNumber`](@ref)
"""
struct HugeNumber <: Real end

Base.convert(::Type{ Float32 },  ::HugeNumber) = 1f+6
Base.convert(::Type{ Float64 },  ::HugeNumber) = 1e+12
Base.convert(::Type{ BigFloat }, ::HugeNumber) = big"1e+24"

Base.show(io::IO, ::HugeNumber) = print(io, "huge")

Base.promote_rule(::Type{ HugeNumber }, ::Type{ I }) where { I <: Integer }  = promote_rule(HugeNumber, promote_type(I, Float64))
Base.promote_rule(::Type{ HugeNumber }, ::Type{ Float32 })                   = Float32
Base.promote_rule(::Type{ HugeNumber }, ::Type{ Float64 })                   = Float64
Base.promote_rule(::Type{ HugeNumber }, ::Type{ BigFloat })                  = BigFloat

##

Base.promote_type(::Type{T}, ::Type{ TinyNumber }, ::Type{ HugeNumber }) where T = promote_type(promote_type(T, TinyNumber), HugeNumber)
Base.promote_type(::Type{T}, ::Type{ HugeNumber }, ::Type{ TinyNumber }) where T = promote_type(promote_type(T, HugeNumber), TinyNumber)

#

"""
   tiny

An instance of a `TinyNumber`. Behaviour and actual value of the `tiny` number depends on the context.

# Example

```jldoctest
julia> tiny
tiny

julia> 1 + tiny
1.000000000001

julia> tiny + 1
1.000000000001

julia> 1f0 + tiny
1.000001f0

julia> big"1.0" + tiny
1.000000000000000000000001

julia> big"1" + tiny
1.000000000000000000000001
```

See also: [`huge`](@ref), [`TinyNumber`](@ref), [`HugeNumber`](@ref)
"""
const tiny = TinyNumber()

"""
   huge

An instance of a `HugeNumber`. Behaviour and actual value of the `huge` number depends on the context.

# Example

```jldoctest
julia> huge
huge

julia> 1 + huge
1.000000000001e12

julia> huge + 1
1.000000000001e12

julia> 1f0 + huge
1.000001f6

julia> big"1.0" + huge
1.000000000000000000000001e+24

julia> big"1" + huge
1.000000000000000000000001e+24
```

See also: [`tiny`](@ref), [`TinyNumber`](@ref), [`HugeNumber`](@ref)
"""
const huge = HugeNumber()

#

function softmax(v)
    ret = exp.(clamp.(v .- maximum(v), -100.0, 0.0))
    ret ./ sum(ret)
end

matrix_from_diagonal(diag::AbstractVector) = matrix_from_diagonal(eltype(diag), diag)

function matrix_from_diagonal(::Type{ T }, diag::AbstractVector) where T
    size   = length(diag)
    matrix = zeros(T, (size, size))
    for i in 1:size
        matrix[i, i] = diag[i]
    end
    return matrix
end


text:
function part1(input)
    joltages = sort(parse.(Int, readlines(input)))
    num_one_diffs = count(==(1), diff(vcat(0, joltages)))
    num_three_diffs = count(==(3), diff(vcat(0, joltages))) + 1
    return num_one_diffs * num_three_diffs
end

function part2(input)
    joltages = sort(parse.(Int, readlines(input)))
    num_combinations = Dict{Int, Int}()
    num_combinations[0] = 1
    for joltage in joltages
        num_combinations[joltage] = sum(get(num_combinations, joltage - n, 0)
                                        for n in 1:3)
    end
    return num_combinations[last(joltages)]
end