text:
# Multiple supresion using the parabolic Radon transform. Real data example.

using PyPlot, Seismic

close("all")

# 1- Download data set from web site (1 NMO-corrected CMP gather from the GOM)
download("http://seismic.physics.ualberta.ca/data/gom_cdp_nmo.su",
         "gom_cdp_nmo.su")

# 2- Convert SU to internal SeismicJulia format
SegyToSeis("gom_cdp_nmo.su", "gom_cdp_nmo.seis", format="su", input_type="ieee",
           swap_bytes=true)

# 3- Read data d and trace headers h
d, h = SeisRead("gom_cdp_nmo.seis")
i1 = 801; i2 = 1201
d = d[i1:i2, :]       # Use subset of the data
imute = findin(d, 0)  # Extract indexes of muted data
offset = Seismic.ExtractHeader(h, "h")
offset = offset*0.3048  # I like SI units
href = maximum(abs(offset))
nt = size(d, 1)
dt = h[1].d1

# 4- Model data in the Radom domain via inersion with parabolic Radon transform
p  = collect(linspace(-0.3, 0.9, 200))
param = Dict(:order=>"parab", :dt=>dt, :p=>p, :h=>offset, :href=>href,
             :flow=>0.0, :fhigh=>125.0)
m = SeisRadonFreqInv(d; mu=0.01, param...)

# 5- Filter primaries and keep multiples in Radon gather
mf = copy(m)
pcut = 0.1
icut = find(p .< pcut)
mf[:, icut] = 0

# 6- Forward model multiples from Radon panel
d_mult = SeisRadonFreqFor(mf, nt; param...)

# 7- Compute primaries = input data - modelled multiples
d_prim = d - d_mult
d_prim[imute] = 0

# 8- Plot input data, Radon panel, and estimated primaries
figure(1, figsize=(9,6));
t1 = (i1-1)*dt
subplot(131)
SeisPlot(d, title="Input CMP", xlabel="Offset [m]", ylabel="Time [s]",
         cmap="PuOr", ox=offset[1], dx=offset[2]-offset[1], oy=t1, dy=dt,
         fignum=1)
subplot(132)
SeisPlot(m[1:nt, :], title="Radon panel", xlabel="Residual moveout [s]",
         cmap="PuOr", ox=p[1], dx=p[2]-p[1], oy=t1, dy=dt, fignum=1)
subplot(133)
SeisPlot(d_prim, title="Estimated primaries", xlabel="Offset [m]", cmap="PuOr",
         ox=offset[1], dx=offset[2]-offset[1], oy=t1, dy=dt, fignum=1)
tight_layout()


text:
function generate_available(domain::Domain, state::State,
                            domain_type::Symbol, state_type::Symbol,
                            action_name::Symbol, action_type::Symbol)
    action = get_actions(domain)[action_name]
    varmap = Dict{Var,Any}(a => :(args[$i].name) for (i, a) in
                           enumerate(get_argvars(action)))
    precond = generate_check_expr(domain, state, get_precond(action), varmap)
    if domain isa AbstractedDomain
        available_def = quote
            function available(domain::$domain_type, state::$state_type,
                               action::$action_type, args)
                precond = @inbounds $precond
                return precond == true || precond == both
            end
        end
    else
        available_def = quote
            function available(domain::$domain_type, state::$state_type,
                               action::$action_type, args)
                return @inbounds $precond
            end
        end
    end
    return available_def
end

function generate_available(domain::Domain, state::State,
                            domain_type::Symbol, state_type::Symbol)
    available_def = quote
        function available(domain::$domain_type, state::$state_type, term::Term)
            available(domain, state, get_actions(domain)[term.name], term.args)
        end
        function available(domain::$domain_type, state::$state_type)
            f = named_action -> begin
                name, action = named_action
                grounded_args = groundargs(domain, state, action)
                return (Compound(name, collect(args)) for args in grounded_args
                        if available(domain, state, action, args))
            end
            actions = Base.Generator(f, pairs(get_actions(domain)))
            return Iterators.flatten(actions)
        end
    end
    return available_def
end

ext: jl


text:
include("../src/InsarTimeseries.jl")
using .InsarTimeseries
using Test

@time begin
    @testset "Common" begin
        include("test_common.jl")
    end
    @testset "Stack averaging" begin
        include("test_stackavg.jl")
    end
    @testset "SBAS inversion" begin
        include("test_sbas.jl")
    end
end


text:
using Fakerjl
using Test
using SafeTestsets

@safetestset "Fakerjl Test" begin include("persontest.jl") end

Fakerjl.init()

println(Fakerjl.firstname())
println(Fakerjl.firstmalename())
println(Fakerjl.firstfemalename())
println(Fakerjl.phone())
println(Fakerjl.email())


text:
function prepare_metnet(proj, base_model; inf_medium = false)

    base_model = deepcopy(base_model)

    ## ---------------------------------------------------------------------
    println("FBOut: ", MetLP.fba!(base_model, HG.HUMAN_BIOMASS_IDER).obj_val)

    ## ---------------------------------------------------------------------
    ## Description

    # This script prepare the original HumanGEM model for modeling. It include some modifications to it (see comment in code) extracted from several data sources.

    ## ---------------------------------------------------------------------
    make_supp_files = (proj == HG)
    make_supp_files && @info("Creating support files")

    ## ---------------------------------------------------------------------
    M, N = size(base_model)
    println("Base Model: ", (M, N))

    ## ---------------------------------------------------------------------
    # mets_readable_ids
    # we build a more human readable met ids
    if make_supp_files
        println("\nMet readable ids")
        mat_model = HG.load_humangem_mat_model()
        met_readable_ids = Dict()
        for i in 1:M
            readable_id = mat_model[:metNames][i] * "[" * mat_model[:comps][mat_model[:metComps][i]] * "]";
            met_readable_ids[base_model.mets[i]] = readable_id
            met_readable_ids[readable_id] = base_model.mets[i]
        end

        # Saving
        sdat(HG, met_readable_ids,
            "met_readable_ids", ".jls";
            verbose = true
        )
    else
        met_readable_ids = HG.load_met_readable_ids()
    end

    ## ---------------------------------------------------------------------
    # Exchanges
    # bkwd and fwd splitted reactions are troublemakers for EP, but they are necessary to model enzymatic costs.
    # So, we leave as least as possible. We unified the exchanges (make them a unique rxn),
    # and let them in a semi-open state (intake bloked, outtake open)

    # I will delete all the Boundary (x) (see comps in the matmodel) metabilites,
    # leaving only the Extracellular (s) metabolites in the exchange reactions.
    # Why? they are not required
    println("\nDeleting Boundary (x) metabolites ")
    to_del = [met for met in base_model.mets if endswith(met, "x")];
    println("todel: ", length(to_del))
    println("before: ", size(base_model))
    MetNets.del_met!(base_model, to_del);
    MetNets.del_void_rxns_mets!(base_model)
    base_model = MetNets.compacted_model(base_model)
    println("after: ", size(base_model))

    # redefining the Chemostat exchange criteria
    exch_subsys_hint = "Exchange/demand"
    function is_exchange(model::MetNets.MetNet, ider::MetNets.IDER_TYPE)
        idx = MetNets.rxnindex(model, ider)
        subsys = model.subSystems[idx]
        return occursin(exch_subsys_hint, string(subsys))
    end

    ## ---------------------------------------------------------------------
    println("FBOut: ", MetLP.fba!(base_model, HG.HUMAN_BIOMASS_IDER).obj_val)

    ## ---------------------------------------------------------------------
    # Exchanges
    exchs = []
    for exch in filter((ider) -> is_exchange(base_model, ider), base_model.rxns)
        exch_i = MetNets.rxnindex(base_model, exch)

        # First close it, later if it is what I want, open the outtake
        MetNets.lb!(base_model, exch_i, 0.0)
        MetNets.ub!(base_model, exch_i, 0.0)

        mets = MetNets.rxn_mets(base_model, exch_i)
        reacts = MetNets.rxn_reacts(base_model, exch_i)

        # I want only the forward monomoleculars
        !(length(mets) == length(reacts) == 1) && continue
        # I want only the exchanges 's'
        isempty(reacts) && !endswith(base_model.mets[first(reacts)], "s") && continue

        # Because, this reactions are forward unbalanced (A <-> nothing)
        # positibe (+) bounds limit the outtake of the cell and
        # negative (-) bounds limit the intake.
        # Because in the Chemostat the intakes is
        # controlled by the medium, we'll close all intakes to handle with them later
        # We'll open all outtakes
        MetNets.lb!(base_model, exch_i, 0.0)
        MetNets.ub!(base_model, exch_i, HG.ABS_MAX_BOUND)

        push!(exchs, exch_i)

    end
    println("\nExchanges: ", exchs |> length)

    ## ---------------------------------------------------------------------
    println("FBOut: ", MetLP.fba!(base_model, HG.HUMAN_BIOMASS_IDER).obj_val)

    ## ---------------------------------------------------------------------
    # Exch Met map
    # A fast way to get the exch reaction from the metabolite and viceversa
    if make_supp_files
        exch_met_map = Dict()
        for rxn in exchs
            mets = MetNets.rxn_mets(base_model, rxn)
            if length(mets) == 1 # only the monomoleculars
                met = base_model.mets[mets[1]]
                rxn = base_model.rxns[rxn]

                exch_met_map[met] = rxn
                exch_met_map[rxn] = met
            end
        end
        println("\nExchanges met map: ", exch_met_map |> length)

        # Saving
        sdat(HG, exch_met_map,
            "exch_met_map", ".jls";
            verbose = true
        )
    else
        exch_met_map = HG.load_exch_met_map()
    end

    ## ---------------------------------------------------------------------
    println("FBOut: ", MetLP.fba!(base_model, HG.HUMAN_BIOMASS_IDER).obj_val)

    ## ---------------------------------------------------------------------
    # Base intake info
    # The Base model will have a medium (expressed as open intake fluxes) that resamble the cultivation
    # at xi = 1 using the set up in Rath 2017 exp A. So the lb of the intakes will be directly the (negative)
    # concentration in 42_MAX_UB standard medium (see Cossio's paper). Also, a few intakes, not justified in
    # the standard medium will be add based in the intakes of the original model FBA analysis.
    base_intake_info = HG.load_base_intake_info(; inf_medium)

    ## ---------------------------------------------------------------------
    println("FBOut: ", MetLP.fba!(base_model, HG.HUMAN_BIOMASS_IDER).obj_val)

    ## ---------------------------------------------------------------------
    # Apply medium
    # The steady state assumption in the context of the Chemostat culture impose a constraint over the intakes dependent of xi and c
    ξ = 1
    println("\nApplaying chemostat bound, xi: ", ξ)
    Chemostat.apply_bound!(base_model, ξ, base_intake_info;
        emptyfirst = true, ignore_miss = true
    );

    ## ---------------------------------------------------------------------
    println("FBOut: ", MetLP.fba!(base_model, HG.HUMAN_BIOMASS_IDER).obj_val)

    ## ---------------------------------------------------------------------
    # Niklas Biomasss

    # I will modified the biomass equation of MODEL1105100000 model with data
    # derived from Niklas (2013): https://doi.org/10.1016/j.ymben.2013.01.002. Table1. (see README)
    # I do not touch the energetic part of the equation, atp + h20 -> adp + h2 + pi
    println("\nApplying Niklas Biomass")

    biomass_idx = MetNets.rxnindex(base_model, HG.HUMAN_BIOMASS_IDER)
    base_model.S[:, biomass_idx] .= zeros(size(base_model, 1))
    for (met, y) in HG.load_Niklas_biomass()
        MetNets.S!(base_model, met, biomass_idx, y)
    end

    ## ---------------------------------------------------------------------
    println("FBOut: ", MetLP.fba!(base_model, HG.HUMAN_BIOMASS_IDER).obj_val)

    ## ---------------------------------------------------------------------
    # ATPM demand

    # +
    # Cell density ρ = 0.25 pgDW/μm³ from Niklas(2011) https://doi.org/10.1007/s00449-010-0502-y.
    # pgDW/μm³ * 1e9 = pgDW/μL
    # pgDW/μL * 1e6 = pgDW/L
    # pgDW/L * 1e-12 = gDW/L
    # atpm_flux = 0.3 mol ATP/ L hr from Fernandez-de-Cossio-Diaz,
    # Jorge, and Alexei Vazquez. https://doi.org/10.1038/s41598-018-26724-7.
    atpm_flux = 0.3 / (0.25 * 1e9 * 1e6 * 1e-12) * 1e3  # mmol/gWD hr

    # This reaction is foward defined with respect to atp
    # atp + more_reacts <-> adp + more_products
    # so we limit the lower bounds as the minimum atp demand
    # that the cell metabolism must fullfill
    if HG.HUMAN_ATPM_IDER in base_model.rxns
        MetNets.lb!(base_model, HG.HUMAN_ATPM_IDER, atpm_flux)

        println("\nATP demand")
        println(MetNets.rxn_str(base_model, HG.HUMAN_ATPM_IDER), " ", MetNets.bounds(base_model, HG.HUMAN_ATPM_IDER))
    end

    ## ---------------------------------------------------------------------
    println("FBOut: ", MetLP.fba!(base_model, HG.HUMAN_BIOMASS_IDER).obj_val)

    ## ---------------------------------------------------------------------
    # dimentions
    base_model = MetNets.force_dims(base_model)

    ## ---------------------------------------------------------------------
    println("FBOut: ", MetLP.fba!(base_model, HG.HUMAN_BIOMASS_IDER).obj_val)

    ## ---------------------------------------------------------------------
    # summary
    println("\nBase Model summary")
    MetNets.summary(base_model)

    ## ---------------------------------------------------------------------
    #  FBA Test
    fbaout = MetLP.fba!(base_model, HG.HUMAN_BIOMASS_IDER)
    println("\nFBAout summary")
    MetNets.summary(base_model, fbaout)

    ## ---------------------------------------------------------------------
    println("\nComparing with experiments")
    compare_with_experimets(base_model)

    return base_model
end


text:
"""
    point foot quadruped
    q = (p, r, f1, f2, f3, f4)
        p - body position
        r - body orientation (modified Rodriques parameters)
        f1 - foot 1 position
        f2 - foot 2 position
        f3 - foot 3 position
        f4 - foot 4 position
"""
mutable struct PointFootQuadruped120{T} <: Model{T}
    # dimensions
	nq::Int # generalized coordinates
    nu::Int # controls
    nw::Int # parameters
    nc::Int # contact points

	# parameters
	body_height::T
	foot_x::T
	foot_y::T

	# parameters
    mass_body::T
    inertia_body::Matrix{T}
    mass_foot::T

	# environment
	μ_world::T
	joint_friction::Vector{T}
	spring_stiffness_joint::Vector{T}
	orientation_friction::Vector{T}
	spring_stiffness_orientation::Vector{T}
	g::T

	# fast methods
	base
	dyn
end

# function L_mult(x)
#     [x[1] -transpose(x[2:4]);
#      x[2:4] x[1] * I(3) + skew(x[2:4])]
# end
#
# # right quaternion multiply as matrix
# function R_mult(x)
#     [x[1] -transpose(x[2:4]); x[2:4] x[1] * I(3) - skew(x[2:4])]
# end
#
# # rotation matrix
# function quat_rot_mat(q)
#     H = [zeros(1, 3); I(3)]
#     transpose(H) * L_mult(q) * transpose(R_mult(q)) * H
# end
#
# function quat_from_mrp(p)
#     """Quaternion (scalar first) from MRP"""
#     return (1.0 / (1.0 + dot(p, p))) * [(1 - dot(p, p)); 2.0 * p]
# end
#
# function mrp_rot_mat(x)
#     quat_rot_mat(quat_from_mrp(x))
# end

# Kinematics
function kinematics(model::PointFootQuadruped120, q)
	q[6 .+ (1:12)]
end

lagrangian(model::PointFootQuadruped120, q, q̇) = 0.0

function M_func(model::PointFootQuadruped120, q)
    cat(
        model.mass_body * Diagonal(ones(3)),     # body position
        model.inertia_body,                      # body orienation
        model.mass_foot * Diagonal(ones(3 * 4)), # feet position
        dims=(1, 2)
        )
end

function C_func(model::PointFootQuadruped120, q, q̇)
	offsets = [
		[+model.foot_x, +model.foot_y, -model.body_height],
		[+model.foot_x, -model.foot_y, -model.body_height],
		[-model.foot_x, +model.foot_y, -model.body_height],
		[-model.foot_x, -model.foot_y, -model.body_height],
		]

	joint_spring = model.spring_stiffness_joint
	orientation_spring = model.spring_stiffness_orientation
	[
	    model.mass_body * [0,0,model.g] + joint_spring .* (+sum(offsets) + 4 * q[1:3] - q[7:9] - q[10:12] - q[13:15] - q[16:18]);     # body position
	    skew(q̇[4:6]) * model.inertia_body * q̇[4:6] .+ orientation_spring .* q[4:6]; # body orientation
	    model.mass_foot * [0,0,model.g] + joint_spring .* (-offsets[1] + q[7:9] - q[1:3]);
	    model.mass_foot * [0,0,model.g] + joint_spring .* (-offsets[2] + q[10:12] - q[1:3]);
	    model.mass_foot * [0,0,model.g] + joint_spring .* (-offsets[3] + q[13:15] - q[1:3]);
	    model.mass_foot * [0,0,model.g] + joint_spring .* (-offsets[4] + q[16:18] - q[1:3]);
	]
end

function ϕ_func(model::PointFootQuadruped120, env::Environment, q)

    position_foot1 = q[6 .+ (1:3)]
    position_foot2 = q[9 .+ (1:3)]
    position_foot3 = q[12 .+ (1:3)]
	position_foot4 = q[15 .+ (1:3)]

	return [position_foot1[3]; position_foot2[3]; position_foot3[3]; position_foot4[3]]
end

function B_func(model::PointFootQuadruped120, q)
    position_body = q[1:3]
    orientation_body = q[3 .+ (1:3)]
	# R = mrp_rot_mat(orientation_body)
	R = euler_rotation_matrix(orientation_body)

	# kinematics in world frame
	r1 = q[6 .+ (1:3)] - position_body
	r2 = q[9 .+ (1:3)] - position_body
	r3 = q[12 .+ (1:3)] - position_body
	r4 = q[15 .+ (1:3)] - position_body

	z3 = zeros(3, 3)

	transpose([
        I(3) I(3) I(3) I(3);
		transpose(R) * skew(r1) transpose(R) * skew(r2) transpose(R) * skew(r3) transpose(R) * skew(r4);
        -I(3)    z3    z3   z3;
        z3    -I(3)    z3   z3;
        z3       z3 -I(3)   z3;
        z3       z3    z3 -I(3)
    ])
end

function A_func(model::PointFootQuadruped120, q)
    @SMatrix [1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
              0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
			  0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]
end

function J_func(model::PointFootQuadruped120, env::Environment, q)
    z3 = zeros(3, 3)

    [
        z3   z3 I(3)   z3   z3   z3;
        z3   z3   z3 I(3)   z3   z3;
        z3   z3   z3   z3 I(3)   z3;
        z3   z3   z3   z3   z3 I(3);
    ]
end

function contact_forces(model::PointFootQuadruped120, env::Environment{<:World, LinearizedCone}, γ1, b1, q2, k)
	m = friction_mapping(env)

	SVector{12}([
		m * b1[1:4]; γ1[1];
		m * b1[5:8]; γ1[2];
		m * b1[9:12]; γ1[3];
		m * b1[13:16]; γ1[4];
		])
end

function velocity_stack(model::PointFootQuadruped120, env::Environment{<:World, LinearizedCone}, q1, q2, k, h)
	v = J_func(model, env, q2) * (q2 - q1) / h[1]
	SVector{16}([
		transpose(friction_mapping(env)) * v[1:2];
		transpose(friction_mapping(env)) * v[4:5];
		transpose(friction_mapping(env)) * v[7:8];
		transpose(friction_mapping(env)) * v[10:11];
	])
end

function friction_coefficients(model::PointFootQuadruped120)
	return [model.μ_world]
end

function initialize_z!(z, model::PointFootQuadruped120, idx::RoboDojo.IndicesZ, q)
    z .= 1.0
    z[idx.q] .= q
end

# function relative_state_cost(qbody, qorientation, qfoot)
# 	# cost function on state: 1/2 * qbody'*Qbody*qbody
# 		# 1/2 * qbody'*Qbody*qbody
# 		# 1/2 * qorientation'*Qorientation*qorientation
# 		# 1/2 * (qfoot-qbody)'*Qfoot*(qfoot-qbody)
# 	Q = zeros(18,18)
# 	Q[1:3,1:3] = Diagonal(qbody)
# 	Q[4:6,4:6] = Diagonal(qorientation)
# 	for i = 1:4
# 		Q[1:3,1:3] += Diagonal(qfoot)
# 		Q[3+3i .+ (1:3), 3+3i .+ (1:3)] += Diagonal(qfoot)
# 		Q[1:3, 3+3i .+ (1:3)] += -Diagonal(qfoot)
# 		Q[3+3i .+ (1:3), 1:3] += -Diagonal(qfoot)
# 	end
# 	return Q
# end

function nominal_configuration(model::PointFootQuadruped120)
    [
        0.0; 0.0; model.body_height;
        0.0; 0.0; 0.0;
        model.foot_x; model.foot_y; 0.0;
        model.foot_x;-model.foot_y; 0.0;
       -model.foot_x; model.foot_y; 0.0;
       -model.foot_x;-model.foot_y; 0.0;
    ]
end

function nominal_state(model::PointFootQuadruped120)
	[nominal_configuration(model); zeros(model.nq)]
end

function dynamics(model::PointFootQuadruped120, h, q0, q1, u1, w1, Λ1, q2)
	# evalutate at midpoint
	qm1 = 0.5 * (q0 + q1)
    vm1 = (q1 - q0) / h[1]
    qm2 = 0.5 * (q1 + q2)
    vm2 = (q2 - q1) / h[1]

	D1L1, D2L1 = lagrangian_derivatives(model, qm1, vm1)
	D1L2, D2L2 = lagrangian_derivatives(model, qm2, vm2)

	d = 0.5 * h[1] * D1L1 + D2L1 + 0.5 * h[1] * D1L2 - D2L2 # variational integrator (midpoint)
	d .+= transpose(B_fast(model, qm2)) * u1        # control inputs
	d .+= transpose(A_fast(model, qm2)) * w1        # control inputs
	d .+= Λ1                                        # contact impulses

	d[4:6] .-= h[1] * model.orientation_friction .* vm2[4:6] # orientation friction
	d[1:3] .-= h[1] * model.joint_friction .* (vm2[1:3] - vm2[7:9]) # joint friction
	d[1:3] .-= h[1] * model.joint_friction .* (vm2[1:3] - vm2[10:12]) # joint friction
	d[1:3] .-= h[1] * model.joint_friction .* (vm2[1:3] - vm2[13:15]) # joint friction
	d[1:3] .-= h[1] * model.joint_friction .* (vm2[1:3] - vm2[16:18]) # joint friction

	d[7:9]   .-= h[1] * model.joint_friction .* (vm2[7:9] - vm2[1:3]) # joint friction
	d[10:12] .-= h[1] * model.joint_friction .* (vm2[10:12] - vm2[1:3]) # joint friction
	d[13:15] .-= h[1] * model.joint_friction .* (vm2[13:15] - vm2[1:3]) # joint friction
	d[16:18] .-= h[1] * model.joint_friction .* (vm2[16:18] - vm2[1:3]) # joint friction
	return d
end

# dimensions
nq = 3 + 3 + 3 * 4       # generalized coordinates
nu = 3 * 4               # controls
nw = 3                   # parameters
nc = 4                   # contact points

# parameters
body_height = 0.3
foot_x = 0.17
foot_y = 0.15

# parameters
g = 9.81                 # gravity
μ_world = 0.3
# V0
joint_friction = 30 * ones(3)
spring_stiffness_joint = 0.0 * ones(3)
orientation_friction = 0.0 * ones(3)
spring_stiffness_orientation = 0.0 * ones(3)

# V1
joint_friction = 00.0 * ones(3)
spring_stiffness_joint = 00.0 * ones(3)
orientation_friction = 5 * ones(3)
spring_stiffness_orientation = 0.5 * ones(3)
# inertial properties
mass_body = 13.5
i_xx = 0.0178533
i_xy = 0.0
i_xz = 0.0
i_yz = 0.0
i_yy = 0.0377999
i_zz = 0.0456542
inertia_body = Array(Diagonal([i_xx, i_yy, i_zz]))
mass_foot = 0.2

point_foot_quadruped = PointFootQuadruped120(nq, nu, nw, nc,
				body_height,
				foot_x,
				foot_y,
				mass_body,
                inertia_body,
                mass_foot,
				μ_world,
				joint_friction,
				spring_stiffness_joint,
				orientation_friction,
                spring_stiffness_orientation,
				g,
				BaseMethods(),
				DynamicsMethods(),
                )


text:
@testset "O(logn)" begin
    function latex_ground(args...)
        translate(-200, -100)
        background("white")
        sethue("black")
        fontsize(50)
    end

    function latex_blend(latex_string, fsize)
        orangeblue = blend(Point(50, 70), Point(50, 120), "orange", "blue")
        setblend(orangeblue)
        fontsize(fsize)
        latex(latex_string, Point(50, 70))
    end

    video = Video(400, 200)
    javis(
        video,
        [
            BackgroundAction(1:1, latex_ground),
            Action((args...) -> latex(L"\mathcal{O}(\log{n})")), # default fontsize 50
            Action(
                (args...) -> latex_blend(L"\mathcal{O}\left(\frac{\log{x}}{2}\right)", 20),
            ),
        ],
        tempdirectory = "images",
        pathname = "",
    )
    @test_reference "refs/ologn.png" load("images/0000000001.png")
    rm("images/0000000001.png")
end

@testset "latex pos in function" begin
    function latex_ground(args...)
        background("white")
        sethue("black")
        fontsize(30)
    end

    function foreground(latex_string)
        latex(latex_string, Point(50, 40))
        circle(O, 20, :fill) # should be in the center and not affected by latex
    end

    video = Video(400, 200)
    javis(
        video,
        [
            BackgroundAction(1:1, latex_ground),
            Action((args...) -> foreground(L"\mathcal{O}(\log{n})")),
        ],
        tempdirectory = "images",
        pathname = "",
    )
    @test_reference "refs/ologn_circ.png" load("images/0000000001.png")
    rm("images/0000000001.png")
end

@testset "latex pos in function DEPRECATED" begin
    function latex_ground(args...)
        background("white")
        sethue("black")
        fontsize(30)
    end

    function foreground(latex_string)
        translate(50, 40)
        fsize = get_fontsize()
        latex(latex_string, fsize) # is automatically changed to fontsize(30); latex(latex_string)
        translate(-50, -40)
        circle(O, 20, :fill) # should be in the center and not affected by latex
    end

    video = Video(400, 200)
    javis(
        video,
        [
            BackgroundAction(1:1, latex_ground),
            Action((args...) -> foreground(L"\mathcal{O}(\log{n})")),
        ],
        tempdirectory = "images",
        pathname = "",
    )
    @test_reference "refs/ologn_circ.png" load("images/0000000001.png")
    rm("images/0000000001.png")
end


text:
module LightField

include("params.jl")
include("psf.jl")
include("psfsize.jl")

include("projection.jl")

using .psf, .psfsize, .params, .projection
using Pkg.TOML

function setup(filename::String)

    config = TOML.parsefile(filename)
    mladict = pop!(config, "mla")
    pitch = pop!(mladict, "pitch")
    fml = pop!(mladict, "fml")

    opt = MicOptics(pop!(config, "optics"), fml)
    sim = SimulationParams(pop!(config, "simparams"), pitch)
    con = Constants(opt)

    xobj::Vector{Float64} = (sim.pixelpitch / opt.M) .* collect(-fld(sim.vpix,2):1:fld(sim.vpix,2))
    yobj::Vector{Float64} = (sim.pixelpitch / opt.M) .* collect(-fld(sim.vpix,2):1:fld(sim.vpix,2))
    zobj::Vector{Float64} = sim.zmin:sim.zstep:sim.zmax

    xml::Vector{Float64} = collect(-fld((sim.subvpix - 1), 2):1:fld((sim.subvpix - 1), 2)) .* sim.subpixelpitch
    yml::Vector{Float64} = collect(-fld((sim.subvpix - 1), 2):1:fld((sim.subvpix - 1), 2)) .* sim.subpixelpitch
    obj = Space(xobj, yobj, zobj)
    mlaspace = Space(xml, yml)
    par = ParameterSet(opt, sim, con)
    img = calcsize(par,obj)
    mla = MicroLensArray(pitch, fml, mlaspace, img, par)
    lf = LightFieldSimulation(par, mla, img, obj)
    return lf
end

export setup, propagate

end # module


text:
using DataFrames, DataFramesMeta, Statistics
import FileIO:load

# Helper functions to merge files of simulations/optimization runs
# that were ran separately

# Get original df files (all p_cs with 10^5 steps)
# and merge them with the results for p_c = 1e-4 with 10^6 steps. (replace for this p_c value)
function df_edit_lowpc(;
                    filepath = "/path/to/jldfiles_GAUSSIAN/",
                    filepath_lowpc = "/path/to/jldfiles_GAUSSIA_10to6/")
    savepath = homedir()
    jldfiles = searchdirendings(filepath, ".jld2")

    for i in 1:length(jldfiles)
        file1 = load(joinpath(filepath, jldfiles[i]))
        df1 = file1["df"]

        file_lowpc = load(joinpath(filepath_lowpc, jldfiles[i]))
        df_lowpc = file_lowpc["df"]

        if in(:nposteriors, names(df1))
            deletecols!(df1, :nposteriors)
        end
        if in(:nposteriors, names(df_lowpc))
            deletecols!(df_lowpc, :nposteriors)
        end
        if in(:seedparticlefilter, names(df1))
            rename!(df1, :seedparticlefilter => :seedlearner)
        end
        if in(:seedparticlefilter, names(df_lowpc))
            rename!(df_lowpc, :seedparticlefilter => :seedlearner)
        end
        # Remove old results of pc=1e-4
        df1 = df1[df1[:, :changeprobability] .!= 1e-4, :]
        # Marge the two dfs
        masterdf = vcat(df1, df_lowpc)
        # Save
        save(joinpath(savepath, jldfiles[i]),
            Dict("df" => masterdf,
                "specs" => file1["specs"],
                "specs_lowpc" => file_lowpc["specs"]))
    end

end

# Get original df files (all p_cs with 10^5 steps)
# and merge them with the results for p_c = 1e-3 with 10^6 steps
# and for p_c = 1e-4 with 10^6 steps. (replace for these p_c value)
function df_edit_midlowlowpc(;
                    filepath = "/path/to/jldfiles_GAUSSIAN/",
                    filepath_midlowpc = "/path/to/jldfiles_GAUSSIAN_1e-3_10to6/",
                    filepath_lowpc = "/path/to/jldfiles_GAUSSIAN_1e-4_10to6/")
    savepath = homedir()
    jldfiles = searchdirendings(filepath, ".jld2")

    for i in 1:length(jldfiles)

        file1 = load(joinpath(filepath, jldfiles[i]))
        df1 = file1["df"]
        file_midlowpc = load(joinpath(filepath_midlowpc, jldfiles[i]))
        df_midlowpc = file_midlowpc["df"]
        file_lowpc = load(joinpath(filepath_lowpc, jldfiles[i]))
        df_lowpc = file_lowpc["df"]

        if in(:nposteriors, names(df1))
            deletecols!(df1, :nposteriors)
        end
        if in(:nposteriors, names(df_midlowpc))
            deletecols!(df_midlowpc, :nposteriors)
        end
        if in(:nposteriors, names(df_lowpc))
            deletecols!(df_lowpc, :nposteriors)
        end

        if in(:seedparticlefilter, names(df1))
            rename!(df1, :seedparticlefilter => :seedlearner)
        end
        if in(:seedparticlefilter, names(df_midlowpc))
            rename!(df_midlowpc, :seedparticlefilter => :seedlearner)
        end
        if in(:seedparticlefilter, names(df_lowpc))
            rename!(df_lowpc, :seedparticlefilter => :seedlearner)
        end
        # Remove old results of pc=1e-3
        df1 = df1[df1[:, :changeprobability] .!= 1e-3, :]
        df1 = df1[df1[:, :changeprobability] .!= 1e-4, :]
        # Merge the two dfs
        masterdf = vcat(df1, df_midlowpc)
        masterdf = vcat(masterdf, df_lowpc)

        @show unique(masterdf[:, :changeprobability])
        #@show unique(masterdf[:, :stochasticity])
        @show unique(masterdf[:, :sigma])
        # Save
        save(joinpath(savepath, jldfiles[i]),
            Dict("df" => masterdf,
                "specsoriginal" => file1["specs"],
                "sigmasmerged" => [file1["specs"].sigmas,
                                                file_midlowpc["specs"].sigmas,
                                                file_lowpc["specs"].sigmas],
                # "stochasticitiesmerged" => [file1["specs"].stochasticities,
                #                             file_midlowpc["specs"].stochasticities,
                #                             file_lowpc["specs"].stochasticities],
                "changeprobabilitiesmerged" => [file1["specs"].changeprobabilities,
                                                file_midlowpc["specs"].changeprobabilities,
                                                file_lowpc["specs"].changeprobabilities],
                "nstepsmerged" => [file1["specs"].nsteps,
                                file_midlowpc["specs"].nsteps,
                                file_lowpc["specs"].nsteps],
                "seedsmerged" => [file1["specs"].seedenvpolicy,
                                file_midlowpc["specs"].seedenvpolicy,
                                file_lowpc["specs"].seedenvpolicy]
                )
            )
    end
end

function df_merge_highlow(;
                    filepath_highpc = "/path/to/jldfiles_GAUSSIAN_high/",
                    filepath_lowpc = "/path/to/jldfiles_GAUSSIAN_low/"
                    )

    savepath = homedir()
    jldfiles = searchdirendings(filepath_highpc, ".jld2")

    for i in 1:length(jldfiles)

        file_high = load(joinpath(filepath_highpc, jldfiles[i]))
        df_high = file_high["df"]

        file_low = load(joinpath(filepath_lowpc, jldfiles[i]))
        df_low = file_low["df"]

        if in(:nposteriors, names(df_high))
            deletecols!(df_high, :nposteriors)
        end
        if in(:nposteriors, names(df_low))
            deletecols!(df_low, :nposteriors)
        end

        if in(:seedparticlefilter, names(df_high))
            rename!(df_high, :seedparticlefilter => :seedlearner)
        end
        if in(:seedparticlefilter, names(df_low))
            rename!(df_low, :seedparticlefilter => :seedlearner)
        end

        # Merge the two dfs
        masterdf = vcat(df_high, df_low)

        @show unique(masterdf[:, :changeprobability])
        if in(:sigma, names(masterdf))
            @show unique(masterdf[:, :sigma])
            save(joinpath(savepath, jldfiles[i]),
                Dict("df" => masterdf,
                    "specshigh" => file_high["specs"],
                    "sigmasmerged" => [file_high["specs"].sigmas,
                                        file_low["specs"].sigmas],
                    "changeprobabilitiesmerged" => [file_high["specs"].changeprobabilities,
                                                    file_low["specs"].changeprobabilities],
                    "nstepsmerged" => [file_high["specs"].nsteps,
                                    file_low["specs"].nsteps],
                    "seedsmerged" => [file_high["specs"].seedenvpolicy,
                                    file_low["specs"].seedenvpolicy]
                    )
                )
        else
            @show unique(masterdf[:, :stochasticity])
            save(joinpath(savepath, jldfiles[i]),
                Dict("df" => masterdf,
                    "specshigh" => file_high["specs"],
                    "stochasticitiesmerged" => [file_high["specs"].stochasticities,
                                                file_low["specs"].stochasticities],
                    "changeprobabilitiesmerged" => [file_high["specs"].changeprobabilities,
                                                    file_low["specs"].changeprobabilities],
                    "nstepsmerged" => [file_high["specs"].nsteps,
                                    file_low["specs"].nsteps],
                    "seedsmerged" => [file_high["specs"].seedenvpolicy,
                                    file_low["specs"].seedenvpolicy]
                    )
                )
        end

    end
end


text:
import DPFEHM
import GaussianRandomFields
import Optim
import PyPlot
import Random
import Zygote

Random.seed!(0)

function getobsnodes(coords, obslocs)
	obsnodes = Array{Int}(undef, length(obslocs))
	for i = 1:length(obslocs)
		obsnodes[i] = findmin(map(j->sum((obslocs[i] .- coords[[1, 3], j]) .^ 2), 1:size(coords, 2)))[2]
	end
	return obsnodes
end

mins = [0, 0]#meters
maxs = [100, 10]#meters
#ns = [11, 11]
#ns = [21, 21]
#ns = [51, 51]
#ns = [101, 101]
#ns = [201, 201]
#ns = [401, 401]
ns = [801, 801]
num_eigenvectors = 200
x_true = randn(num_eigenvectors)
x0 = zeros(num_eigenvectors)
sqrtnumobs = 16
obslocs_x = range(mins[1], maxs[1]; length=sqrtnumobs + 2)[2:end - 1]
obslocs_z = range(mins[2], maxs[2]; length=sqrtnumobs + 2)[2:end - 1]
obslocs = collect(Iterators.product(obslocs_x, obslocs_z))[:]
observations = Array{Float64}(undef, length(obslocs))

	coords, neighbors, areasoverlengths, _ = DPFEHM.regulargrid2d(mins, maxs, ns, 1.0)
	coords = vcat(coords, zeros(size(coords, 2))')#add a z component
	coords[3, :] = coords[2, :]#swap the y and z components
	coords[2, :] .= 0.0
	Qs = zeros(size(coords, 2))
	boundaryhead(x, y, z) = 0.0 - z
	dirichletnodes = Int[]
	dirichleths = zeros(size(coords, 2))
	for i = 1:size(coords, 2)
		if coords[3, i] == mins[2]
			push!(dirichletnodes, i)
			dirichleths[i] = boundaryhead(coords[:, i]...)
		elseif coords[3, i] == maxs[2] && coords[1, i] > maxs[1] * 1 / 4 && coords[1, i] < maxs[1] * 3 / 4
			Qs[i] = 1e-2
		end
	end

	lambda = 10.0#meters -- correlation length of log-conductivity
	sigma = 1.0#standard deviation of log-conductivity
	mu = 0.0#mean of log-permeability can be arbitrary because this is steady-state and there are no fluid sources

	#cov = GaussianRandomFields.CovarianceFunction(2, GaussianRandomFields.Matern(lambda, 1; σ=sigma))
	#anisotropy = [1 / lambda 0.1 / lambda;  / lambda 1 / lambda]
	anisotropy = [1 0.9; 0.9 1] / lambda
	cov = GaussianRandomFields.CovarianceFunction(2, GaussianRandomFields.AnisotropicExponential(anisotropy, σ=sigma))
	x_pts = range(mins[1], maxs[1]; length=ns[1])
	z_pts = range(mins[2], maxs[2]; length=ns[2])
	#@time grf = GaussianRandomFields.GaussianRandomField(cov, GaussianRandomFields.KarhunenLoeve(num_eigenvectors), x_pts, z_pts)
	@time grf = GaussianRandomFields.GaussianRandomField(cov, GaussianRandomFields.KarhunenLoeve(num_eigenvectors), x_pts, x_pts)
	@time logKs = GaussianRandomFields.sample(grf)
	alphas = fill(0.375, length(neighbors))
	Ns = fill(1.25, length(neighbors))
	parameterization = copy(grf.data.eigenfunc)
	sigmas = copy(grf.data.eigenval)

	#plot a realization
	fig, ax = PyPlot.subplots()
	ax.imshow(logKs, origin="lower")
	ax.title.set_text("Random Permeability Field")
	display(fig)
	println()
	PyPlot.close(fig)

	#plot the first few eigenvectors and some random ones
	fig, axs = PyPlot.subplots(2, 4, figsize=(16, 8))
	rp = Random.randperm(size(parameterization, 2))
	for (i, ax) in enumerate(axs)
		if i <= length(axs) / 2
			ax.imshow(reshape(parameterization[:, i], ns...), origin="lower")
			ax.title.set_text("Eigenvector $i")
		else
			ax.imshow(reshape(parameterization[:, rp[i]], ns...), origin="lower")
			ax.title.set_text("Eigenvector $(rp[i])")
		end
	end
	display(fig)
	println()
	PyPlot.close(fig)

	#logKs2Ks_neighbors(Ks) = map(neighbor->exp(0.5 * (Ks[neighbor[1]] + Ks[neighbor[2]])), neighbors)
	#logKs2Ks_neighbors(Ks) = exp.([0.5 * (Ks[i1] + Ks[i2]) for (i1, i2) in neighbors])
	logKs2Ks_neighbors(Ks) = exp.(0.5 * (Ks[map(p->p[1], neighbors)] .+ Ks[map(p->p[2], neighbors)]))#Zygote differentiates this efficiently but the definitions above are ineffecient with Zygote
	function solveforh(logKs, dirichleths)
		@assert length(logKs) == length(Qs)
		if maximum(logKs) - minimum(logKs) > 25
			return fill(NaN, length(Qs))#this is needed to prevent the solver from blowing up if the line search takes us somewhere crazy
		else
			Ks_neighbors = logKs2Ks_neighbors(logKs)
			psi0 = map(i->boundaryhead(coords[:, i]...), 1:size(coords, 2))
			return DPFEHM.richards_steadystate(psi0, Ks_neighbors, neighbors, areasoverlengths, dirichletnodes, dirichleths, coords, alphas, Ns, Qs)
		end
	end
	solveforheigs(x) = solveforh(x2logKs(x), dirichleths)
	x2logKs(x) = reshape(parameterization * (sigmas .* x), ns...)
	logKs_true = x2logKs(x_true)
	h_true = solveforheigs(x_true)

	obsnodes = getobsnodes(coords, obslocs)
	obssigma = 1e-3
	observations .= h_true[obsnodes]#set up the observations
	f(logKs) = sum((solveforh(logKs, dirichleths)[obsnodes] - observations) .^ 2 ./ obssigma ^ 2)
	print("forward solve time")
	@time h = solveforh(logKs, dirichleths)
	print("gradient time")
	@time zg = Zygote.gradient(f, logKs)[1]

	#plot the solution, the difference between the solution and the solution without recharge, and the logKs
	fig, axs = PyPlot.subplots(1, 4, figsize=(16, 4))
	ims = axs[1].imshow(reshape(h, ns...), origin="lower")
	axs[1].title.set_text("Head for Random Permeability")
	fig.colorbar(ims, ax=axs[1])
	ims = axs[2].imshow(DPFEHM.effective_saturation.(alphas[1], reshape(h, ns...), Ns[1]), origin="lower")
	axs[2].title.set_text("Effective Saturation\nfor Random Permeability")
	fig.colorbar(ims, ax=axs[2])
	ims = axs[3].imshow(logKs, origin="lower")
	axs[3].title.set_text("Random Permeability Field")
	fig.colorbar(ims, ax=axs[3])
	ims = axs[4].imshow(zg, origin="lower")
	axs[4].title.set_text("Gradient of Loss\nw.r.t. Permeability")
	fig.colorbar(ims, ax=axs[4])
	fig.tight_layout()
	display(fig)
	println()
	PyPlot.close(fig)


text:
function initial_location_of(path)
    # In the case where an argument tensor is permenantly assigned to DRAM, it will
    # have no `path`.
    #
    # Therefore, the absense of `path` implies an initial location of DRAM
    isempty(path) && return DRAM

    initial_location = first(path).location
    if initial_location == LOC_PMEM
        return PMEM
    elseif isdram(initial_location)
        return DRAM
    else
        error("$(initial_location)???")
    end
end

# Overload for a Frame
function configure!(f::nGraph.NFunction, frame::Frame)
    # Get initial schedules for the frame
    initial_schedule = get_schedule(frame)

    # Convert this into an appropriate format for the inner `configure!`
    schedule = Dict(
        t => (initial_location_of(path), getactions(tensor_graph, path))
        for (t, (tensor_graph, path)) in initial_schedule
    )

    # Add in the fixed tensors
    for t in fixed_tensors(frame)
        schedule[t] = (first(locations(t)), MoveAction[])
    end

    # TODO: Move this into the innermost `configure!`
    data = frame.profile_data
    for node in nodes(data)
        if nGraph.Lib.can_select_algo(nGraph.getpointer(unx(node)))
            # Only access this once we know that there is at least one node where the
            # algorithm can be decided.
            #
            # Otherwise, `frame.model[:algo_var]` will not be defined.
            algo_var = frame.model[:algo_var]
            count = 0
            local algo_enum
            for enum in enums(gettime(node))
                if approx_one(algo_var[node, enum])
                    count += 1
                    algo_enum = enum
                end
            end

            # Only one algorithm should be selected
            @assert count == 1
            nGraph.Lib.set_algo(
                nGraph.getpointer(unx(node)),
                convert(UInt, algo_enum),
                convert(UInt, Profiler.bytesat(gettime(node), algo_enum))
            )
        end
    end

    @info "Calling Inner Configure"
    return configure!(f, schedule, data)
end

# Get the path of the tensor traced through the graph
function get_schedule(F::Frame)
    data = F.profile_data
    model_graphs = F.model[:tensor_graphs]

    schedule = Dict{XTensor{XNode}, Tuple{MetaGraph, Vector{VertexMetadata}}}()

    for tensor in free_tensors(F)
        g = graph(descriptor(F, tensor))

        # Trace the route taken through the graph
        v = find_vertex((g, v) -> getmeta(g, v).location == LOC_SOURCE, g)

        path = [getmeta(g, v)]
        seen = Int[]
        while getmeta(g, v).location != LOC_SINK
            if isempty(outedges(g, v)) || in(v, seen)
                error("""
                $tensor
                $(getmeta(g, v))
                """)
            end

            push!(seen, v)
            for e in outedges(g, v)
                if approx_one(model_graphs[tensor, e])
                    v = dst(e)
                    break
                end
            end
            push!(path, getmeta(g, v))
        end
        # Drop the first source element and last sink element
        popfirst!(path)
        pop!(path)

        schedule[tensor] = (g, path)
    end

    return schedule
end

# Consume all of the PKEEP nodes.
function getkeeps(vertices::Vector{VertexMetadata}, index)
    keeps = XNode[]
    while checkbounds(Bool, vertices, index) && isdram(vertices[index].location)
        vertex = vertices[index]
        if vertex.isuser
            push!(keeps, vertices[index].op)
        end
        index += 1
    end
    return unique(keeps)
end

function isasync(tensor_graph, a::VertexMetadata, b::VertexMetadata)
    # Get the vertex number from the metadata - construct the edge
    src = a.vertex_number
    dst = b.vertex_number
    edge_metadata = getmeta(tensor_graph, edgetype(tensor_graph)(src, dst))
    return isasync(edge_metadata)
end

# Return `true` if there is an implied write to
write_to_pmem(a, b) = a == LOC_DRAM_PRE && ispmem(b)
read_from_pmem(a, b) = ispmem(a) && isdram(b)

function getactions(tensor_graph, vertices::Vector{VertexMetadata})
    actions = MoveAction[]

    # Same as for the initial location - vertices is empty, this implies that the tensor
    # is an argument and will thus have no actions.
    isempty(vertices) && return actions
    written_to_pmem = false

    for i in Iterators.drop(eachindex(vertices), 1)
        src = vertices[i-1]
        dst = vertices[i]
        a, b = src.location, dst.location

        # Determine whether this is an asynchronous move
        if isasync(tensor_graph, src, dst)
            concurrent = src.op
        else
            concurrent = nothing
        end

        if write_to_pmem(a, b)
            if written_to_pmem
                @show actions
                error()
            end
            # All downstream users are consumers
            consumers = unique(vertices[i].op for i in i:length(vertices) if vertices[i].isuser)

            # One solution to the ILP is to move data back and forth if it does not cause
            # any additional overhead.
            #
            # Here, we just filter out these movements by checking if the consumers of a
            # move is empty
            if !isempty(consumers)
                push!(actions, MoveAction(consumers, PMEM, true, concurrent))
                written_to_pmem = true
            end
        end

        if read_from_pmem(a, b)
            consumers = getkeeps(vertices, i)
            if !isempty(consumers)
                push!(actions, MoveAction(consumers, DRAM, false, concurrent))
            end
        end
    end

    # Need to filter out the first op from showing up in the actions because the first op
    # doesn't actually use the tensor - it produces it.
    producing_op = first(vertices).op

    for action in actions
        filter!(!isequal(producing_op), action.consumers)
    end
    return actions
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
function smooth_poly_max(x::Float64)
    if x<=zero(Float64)
        return(zero(Float64))
    elseif x>= one(Float64)
        return(x)
    else
        return( 2x^2 - x^3)
    end
end

function smooth_poly_max(x, ϵ, κ)
    mult_λ = (κ+1)*ϵ
    smooth_poly_max((x-ϵ)/mult_λ)*mult_λ + ϵ
end

function supersmooth(x::Float64)
   if x<=zero(Float64)
        return(zero(Float64))
   else
        return(exp(-1/x))
   end
end

function supersmooth_max(x)
    x*(supersmooth(x)/(supersmooth(x) + supersmooth(1-x)))
end

function supersmooth_max(x, ϵ, κ)
    mult_λ = (κ+1)*ϵ
    supersmooth_max((x-ϵ)/mult_λ)*mult_λ + ϵ
end

beyonce_palette = ["#424395" "#EB549A" "#5EC2DA" "#EBC915" "#018AC4"  "#550133"]


text:
using Base.Test
using QuantumOptics

mutable struct test_lazytensor <: Operator
    basis_l::Basis
    basis_r::Basis
    data::Matrix{Complex128}
    test_lazytensor(b1::Basis, b2::Basis, data) = length(b1) == size(data, 1) && length(b2) == size(data, 2) ? new(b1, b2, data) : throw(DimensionMismatch())
end

@testset "operators-lazytensor" begin

srand(0)

D(op1::Operator, op2::Operator) = abs(tracedistance_nh(full(op1), full(op2)))
D(x1::StateVector, x2::StateVector) = norm(x2-x1)

b1a = GenericBasis(2)
b1b = GenericBasis(3)
b2a = GenericBasis(1)
b2b = GenericBasis(4)
b3a = GenericBasis(6)
b3b = GenericBasis(5)

b_l = b1a⊗b2a⊗b3a
b_r = b1b⊗b2b⊗b3b

op1 = randoperator(b1a, b1b)
op2 = randoperator(b2a, b2b)
op3 = randoperator(b3a, b3b)

# Test creation
@test_throws AssertionError LazyTensor(b_l, b_r, [1], [randoperator(b1a)])
@test_throws AssertionError LazyTensor(b_l, b_r, [1, 2], [op1])
@test_throws AssertionError LazyTensor(b_l, b_r, [1, 2], [op1, sparse(randoperator(b_l, b_l))])
@test_throws AssertionError LazyTensor(b_l, b_r, [1, 2], [randoperator(b_r, b_r), sparse(op2)])

@test LazyTensor(b_l, b_r, [2, 1], [op2, op1]) == LazyTensor(b_l, b_r, [1, 2], [op1, op2])
x = randoperator(b2a)
@test LazyTensor(b_l, 2, x) == LazyTensor(b_l, b_l, [2], [x])

# Test copy
x = 2*LazyTensor(b_l, b_r, [1,2], [randoperator(b1a, b1b), sparse(randoperator(b2a, b2b))])
x_ = copy(x)
@test x == x_
@test !(x === x_)
x_.operators[1].data[1,1] = complex(10.)
@test x.operators[1].data[1,1] != x_.operators[1].data[1,1]
x_.factor = 3.
@test x_.factor != x.factor
x_.indices[2] = 100
@test x_.indices != x.indices

# Test full & sparse
I2 = identityoperator(b2a, b2b)
x = LazyTensor(b_l, b_r, [1, 3], [op1, sparse(op3)], 0.3)
@test 1e-12 > D(0.3*op1⊗full(I2)⊗op3, full(x))
@test 1e-12 > D(0.3*sparse(op1)⊗I2⊗sparse(op3), sparse(x))

# Test suboperators
@test operators_lazytensor.suboperator(x, 1) == op1
@test operators_lazytensor.suboperator(x, 3) == sparse(op3)
@test operators_lazytensor.suboperators(x, [1, 3]) == [op1, sparse(op3)]

# Arithmetic operations
# =====================
subop1 = randoperator(b1a, b1b)
subop2 = randoperator(b2a, b2b)
subop3 = randoperator(b3a, b3b)
I1 = full(identityoperator(b1a, b1b))
I2 = full(identityoperator(b2a, b2b))
I3 = full(identityoperator(b3a, b3b))
op1 = LazyTensor(b_l, b_r, [1, 3], [subop1, sparse(subop3)], 0.1)
op1_ = 0.1*subop1 ⊗ I2 ⊗ subop3
op2 = LazyTensor(b_l, b_r, [2, 3], [sparse(subop2), subop3], 0.7)
op2_ = 0.7*I1 ⊗ subop2 ⊗ subop3
op3 = 0.3*LazyTensor(b_l, b_r, 3, subop3)
op3_ = 0.3*I1 ⊗ I2 ⊗ subop3

x1 = Ket(b_r, rand(Complex128, length(b_r)))
x2 = Ket(b_r, rand(Complex128, length(b_r)))
xbra1 = Bra(b_l, rand(Complex128, length(b_l)))
xbra2 = Bra(b_l, rand(Complex128, length(b_l)))

# Addition
@test_throws ArgumentError op1 + op2
@test_throws ArgumentError op1 - op2
@test 1e-14 > D(-op1_, -op1)

# Test multiplication
@test_throws bases.IncompatibleBases op1*op2
@test 1e-11 > D(op1*(x1 + 0.3*x2), op1_*(x1 + 0.3*x2))
@test 1e-11 > D((xbra1 + 0.3*xbra2)*op1, (xbra1 + 0.3*xbra2)*op1_)
@test 1e-11 > D(op1*x1 + 0.3*op1*x2, op1_*x1 + 0.3*op1_*x2)
@test 1e-12 > D(dagger(x1)*dagger(0.3*op2), dagger(x1)*dagger(0.3*op2_))
@test 1e-12 > D(op1_*dagger(0.3*op2), op1_*dagger(0.3*op2_))
@test 1e-12 > D(dagger(0.3*op2)*op1_, dagger(0.3*op2_)*op1_)
@test 1e-12 > D(dagger(0.3*op2)*op1, dagger(0.3*op2_)*op1_)

# Test division
@test 1e-14 > D(op1/7, op1_/7)

# Test identityoperator
Idense = identityoperator(DenseOperator, b_r)
I = identityoperator(LazyTensor, b_r)
@test isa(I, LazyTensor)
@test full(I) == Idense
@test 1e-11 > D(I*x1, x1)
@test I == identityoperator(LazyTensor, b1b) ⊗ identityoperator(LazyTensor, b2b) ⊗ identityoperator(LazyTensor, b3b)

Idense = identityoperator(DenseOperator, b_l)
I = identityoperator(LazyTensor, b_l)
@test isa(I, LazyTensor)
@test full(I) == Idense
@test 1e-11 > D(xbra1*I, xbra1)
@test I == identityoperator(LazyTensor, b1a) ⊗ identityoperator(LazyTensor, b2a) ⊗ identityoperator(LazyTensor, b3a)

# Test trace and normalize
subop1 = randoperator(b1a)
I2 = full(identityoperator(b2a))
subop3 = randoperator(b3a)
op = LazyTensor(b_l, b_l, [1, 3], [subop1, sparse(subop3)], 0.1)
op_ = 0.1*subop1 ⊗ I2 ⊗ subop3

@test trace(op) ≈ trace(op_)
op_normalized = normalize(op)
@test trace(op_) ≈ trace(op)
@test 1 ≈ trace(op_normalized)
op_copy = deepcopy(op)
normalize!(op_copy)
@test trace(op) != trace(op_copy)
@test 1 ≈ trace(op_copy)

# Test partial trace
subop1 = randoperator(b1a)
I2 = full(identityoperator(b2a))
subop3 = randoperator(b3a)
op = LazyTensor(b_l, b_l, [1, 3], [subop1, sparse(subop3)], 0.1)
op_ = 0.1*subop1 ⊗ I2 ⊗ subop3

@test 1e-14 > D(ptrace(op_, 3), ptrace(op, 3))
@test 1e-14 > D(ptrace(op_, 2), ptrace(op, 2))
@test 1e-14 > D(ptrace(op_, 1), ptrace(op, 1))

@test 1e-14 > D(ptrace(op_, [2,3]), ptrace(op, [2,3]))
@test 1e-14 > D(ptrace(op_, [1,3]), ptrace(op, [1,3]))
@test 1e-14 > D(ptrace(op_, [1,2]), ptrace(op, [1,2]))

@test_throws ArgumentError ptrace(op, [1,2,3])

# Test expect
state = Ket(b_l, rand(Complex128, length(b_l)))
@test expect(op, state) ≈ expect(op_, state)

state = DenseOperator(b_l, b_l, rand(Complex128, length(b_l), length(b_l)))
@test expect(op, state) ≈ expect(op_, state)

# Permute systems
subop1 = randoperator(b1a, b1b)
subop2 = randoperator(b2a, b2b)
subop3 = randoperator(b3a, b3b)
I2 = full(identityoperator(b2a, b2b))
op = LazyTensor(b_l, b_r, [1, 3], [subop1, sparse(subop3)])*0.1
op_ = 0.1*subop1 ⊗ I2 ⊗ subop3

@test 1e-14 > D(permutesystems(op, [1, 3, 2]), permutesystems(op_, [1, 3, 2]))
@test 1e-14 > D(permutesystems(op, [2, 1, 3]), permutesystems(op_, [2, 1, 3]))
@test 1e-14 > D(permutesystems(op, [2, 3, 1]), permutesystems(op_, [2, 3, 1]))
@test 1e-14 > D(permutesystems(op, [3, 1, 2]), permutesystems(op_, [3, 1, 2]))
@test 1e-14 > D(permutesystems(op, [3, 2, 1]), permutesystems(op_, [3, 2, 1]))

# Test gemv
subop1 = randoperator(b1a, b1b)
subop2 = randoperator(b2a, b2b)
subop3 = randoperator(b3a, b3b)
I2 = full(identityoperator(b2a, b2b))
op = LazyTensor(b_l, b_r, [1, 3], [subop1, sparse(subop3)])*0.1
op_ = 0.1*subop1 ⊗ I2 ⊗ subop3

state = Ket(b_r, rand(Complex128, length(b_r)))
result_ = Ket(b_l, rand(Complex128, length(b_l)))
result = deepcopy(result_)
operators.gemv!(complex(1.), op, state, complex(0.), result)
@test 1e-13 > D(result, op_*state)

result = deepcopy(result_)
alpha = complex(1.5)
beta = complex(2.1)
operators.gemv!(alpha, op, state, beta, result)
@test 1e-13 > D(result, alpha*op_*state + beta*result_)

state = Bra(b_l, rand(Complex128, length(b_l)))
result_ = Bra(b_r, rand(Complex128, length(b_r)))
result = deepcopy(result_)
operators.gemv!(complex(1.), state, op, complex(0.), result)
@test 1e-13 > D(result, state*op_)

result = deepcopy(result_)
alpha = complex(1.5)
beta = complex(2.1)
operators.gemv!(alpha, state, op, beta, result)
@test 1e-13 > D(result, alpha*state*op_ + beta*result_)

# Test gemm
b_l2 = GenericBasis(17)
b_r2 = GenericBasis(13)
subop1 = randoperator(b1a, b1b)
subop2 = randoperator(b2a, b2b)
subop3 = randoperator(b3a, b3b)
I2 = full(identityoperator(b2a, b2b))
op = LazyTensor(b_l, b_r, [1, 3], [subop1, sparse(subop3)])*0.1
op_ = 0.1*subop1 ⊗ I2 ⊗ subop3

state = randoperator(b_r, b_r2)
result_ = randoperator(b_l, b_r2)
result = deepcopy(result_)
operators.gemm!(complex(1.), op, state, complex(0.), result)
@test 1e-12 > D(result, op_*state)

result = deepcopy(result_)
alpha = complex(1.5)
beta = complex(2.1)
operators.gemm!(alpha, op, state, beta, result)
@test 1e-12 > D(result, alpha*op_*state + beta*result_)

state = randoperator(b_l2, b_l)
result_ = randoperator(b_l2, b_r)
result = deepcopy(result_)
operators.gemm!(complex(1.), state, op, complex(0.), result)
@test 1e-12 > D(result, state*op_)

result = deepcopy(result_)
alpha = complex(1.5)
beta = complex(2.1)
operators.gemm!(alpha, state, op, beta, result)
@test 1e-12 > D(result, alpha*state*op_ + beta*result_)

# Test calling gemv with non-complex factors
state = Ket(b_r, rand(Complex128, length(b_r)))
result_ = Ket(b_l, rand(Complex128, length(b_l)))
result = deepcopy(result_)
operators.gemv!(1, op, state, 0, result)
@test 1e-13 > D(result, op_*state)

state = Bra(b_l, rand(Complex128, length(b_l)))
result_ = Bra(b_r, rand(Complex128, length(b_r)))
result = deepcopy(result_)
operators.gemv!(1, state, op, 0, result)
@test 1e-13 > D(result, state*op_)

# Test gemm errors
test_op = test_lazytensor(b1a, b1a, rand(2, 2))
test_lazy = LazyTensor(tensor(b1a, b1a), [1, 2], [test_op, test_op])
test_ket = Ket(tensor(b1a, b1a), rand(4))

@test_throws ArgumentError operators.gemv!(alpha, test_lazy, test_ket, beta, copy(test_ket))
@test_throws ArgumentError operators.gemv!(alpha, dagger(test_ket), test_lazy, beta, copy(dagger(test_ket)))

end # testset


text:
# ------------------------------------------------------------------------------------------
# ## Learn how to build a simple multi-layer-perceptron on the MNIST dataset
#
# MNIST from: https://github.com/FluxML/model-zoo/blob/master/mnist/mlp.jl
#
# Let's start by loading `Flux`, importing a few things from `Flux` explicitly, and bringing
# the `repeated` function into our scope.
# ------------------------------------------------------------------------------------------

using Flux, Flux.Data.MNIST
using Flux: onehotbatch, argmax, crossentropy, throttle
using Base.Iterators: repeated

# ------------------------------------------------------------------------------------------
# We can now store all the MNIST images in `imgs` and take a peak into this vector to see
# what the data looks like
# ------------------------------------------------------------------------------------------

imgs = MNIST.images()
imgs[3]

# ------------------------------------------------------------------------------------------
# Let's look at the type of an individual image.
# ------------------------------------------------------------------------------------------

typeof(imgs[3])

# ------------------------------------------------------------------------------------------
# #### Reorganizing our array of images
#
# We see this is a 2D array that stores `ColorTypes`. To work more easily with this data,
# let's convert all `ColorTypes` to floating point numbers.
# ------------------------------------------------------------------------------------------

fpt_imgs = float.(imgs)

# ------------------------------------------------------------------------------------------
# Now we can see what `imgs[3]` looks like as an array of floats, rather than as an array of
# colors!
# ------------------------------------------------------------------------------------------

fpt_imgs[3]

# ------------------------------------------------------------------------------------------
# **Let's stack the images to create one large 2D array, `X`, that stores the data for each
# image as a column.**
#
# To do this, we can **first** use `reshape` to unravel each image, creating a 1D array
# (`Vector`) of floats from a 2D array (`Matrix`) of floats.
# ------------------------------------------------------------------------------------------

unraveled_fpt_imgs = reshape.(fpt_imgs, :);
typeof(unraveled_fpt_imgs)

# ------------------------------------------------------------------------------------------
# (Note that `Vector` is an alias for a 1D `Array`.)
# ------------------------------------------------------------------------------------------

Vector

# ------------------------------------------------------------------------------------------
# This makes `unraveled_fpt_imgs` a `Vector` of `Vector`s where `imgs[3]` is now
# ------------------------------------------------------------------------------------------

unraveled_fpt_imgs[3]

# ------------------------------------------------------------------------------------------
# After using `reshape` to get a `Vector` of `Vector`s, we can use `hcat` to build a
# `Matrix`, `X`, from `unraveled_fpt_imgs` where the `Vector`s stored in
# `unraveled_fpt_imgs` will become the columns of `X`.
#
# Note that we're using the "splat" command below, `...`, which allows you to pass all the
# elements of an object to a function, rather than just passing the object itself.
# ------------------------------------------------------------------------------------------

X = hcat(unraveled_fpt_imgs...)

# ------------------------------------------------------------------------------------------
# #### How to go back to images from this 2D `Array`
#
# So now each column in X is an image reshaped to a vector of floating points. Let's pick
# one column and see what the digit is.
#
# Let's try to view the second image in the original array, `imgs`, by taking the second
# column of `X`
# ------------------------------------------------------------------------------------------

onefigure = X[:,2]

# ------------------------------------------------------------------------------------------
# We'll `reshape` this array to a 2D, 28x28 array,
# ------------------------------------------------------------------------------------------

t1 = reshape(onefigure,28,28)

# ------------------------------------------------------------------------------------------
# and finally use `colorview` from the `Images` package to view the handwritten digit.
# ------------------------------------------------------------------------------------------

using Images

colorview(Gray, t1)

# ------------------------------------------------------------------------------------------
# *Our data is in working order!*
#
# For our machine to learn the digit with which each image is associated, we'll need to
# train it using correct answers. Therefore we'll make use of the `labels` associated with
# these images from MNIST.
# ------------------------------------------------------------------------------------------

labels = MNIST.labels() # the true labels

# ------------------------------------------------------------------------------------------
# One-hot-encode the labels with `onehotbatch`
# ------------------------------------------------------------------------------------------

Y = onehotbatch(labels, 0:9)

# ------------------------------------------------------------------------------------------
# which gives a binary indicator vector for each figure
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Build the network
# ------------------------------------------------------------------------------------------

m = Chain(
  Dense(28^2, 32, relu),
  Dense(32, 10),
  softmax)

# ------------------------------------------------------------------------------------------
# Define the loss functions and accuracy
# ------------------------------------------------------------------------------------------

loss(x, y) = Flux.crossentropy(m(x), y)
accuracy(x, y) = mean(argmax(m(x)) .== argmax(y))

# ------------------------------------------------------------------------------------------
# Use `X` to create our training data and then declare our evaluation function:
# ------------------------------------------------------------------------------------------

dataset = repeated((X, Y), 200)
evalcb = () -> @show(loss(X, Y))
opt = ADAM(Flux.params(m))

# ------------------------------------------------------------------------------------------
# So far, we have defined our training data and our evaluation functions.
#
# Let's take a look at the function signature of Flux.train!
# ------------------------------------------------------------------------------------------

?Flux.train!

# ------------------------------------------------------------------------------------------
# **Now we can train our model and look at the accuracy thereafter.**
# ------------------------------------------------------------------------------------------

Flux.train!(loss, dataset, opt, cb = throttle(evalcb, 10))

accuracy(X, Y)

# ------------------------------------------------------------------------------------------
# Now that we've trained our model, let's create test data, `tX`,
# ------------------------------------------------------------------------------------------

tX = hcat(float.(reshape.(MNIST.images(:test), :))...)

# ------------------------------------------------------------------------------------------
# and run our model on one of the images from `tX`
# ------------------------------------------------------------------------------------------

test_image = m(tX[:,1])

indmax(test_image) - 1

# ------------------------------------------------------------------------------------------
# The largest element of `test_image` is the 8th element, so our model says that test_image
# is a "7".
#
# Now we can look at the original image.
# ------------------------------------------------------------------------------------------

using Images
t1 = reshape(tX[:,1],28,28)
colorview(Gray, t1)

# ------------------------------------------------------------------------------------------
# and there we have it!
# ------------------------------------------------------------------------------------------

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
import Term.Dendograms: link, Dendogram

@testset "DENDOGRAM" begin
    mydendo = Dendogram("awesome", "this", :is, "a", "dendogram!")
    otherdendo = Dendogram("head", "these", "are", "colorful", "leaves")
    smalldendo = Dendogram("head", [1, 2])

    @test string(mydendo) ==
          "                     \e[38;2;255;171;145mawesome\e[39m                      \n\e[1m\e[2m\e[38;2;144;202;249m     ┌────────────┬──────┴─────┬────────────┐     \e[22m\e[22m\e[39m\n\e[38;2;176;190;197m   this          is            a       dendogram! \e[39m"

    large = link(mydendo, otherdendo; title = "{red}superdendo{/red}")
    @test string(large) ==
          "                                             \e[38;2;255;171;145m\e[31msuperdendo\e[39m\e[38;2;255;171;145m\e[39m                                             \n\e[1m\e[2m\e[38;2;144;202;249m                         ┌────────────────────────┴───────────────────────┐                         \e[22m\e[22m\e[39m\n                     \e[38;2;255;171;145mawesome\e[39m                                             \e[38;2;255;171;145mhead\e[39m                       \n\e[1m\e[2m\e[38;2;144;202;249m     ┌────────────┬──────┴─────┬────────────┐     \e[22m\e[22m\e[39m\e[1m\e[2m\e[38;2;144;202;249m     ┌────────────┬──────┴─────┬────────────┐     \e[22m\e[22m\e[39m\n\e[38;2;176;190;197m   this          is            a       dendogram! \e[39m\e[38;2;176;190;197m   these         are       colorful      leaves   \e[39m"

    # nested = link(smalldendo, link(smalldendo, link(smalldendo, smalldendo; title="a level"); title="another level"); title="first level")
    # @test string(nested) == fromfile("./txtfiles/dendogram_nested.txt")

    D = link(mydendo, mydendo, mydendo)
    @test D.measure.w == 150
end


text:
# ---
# title: 16. 3Sum Closest
# id: problem16
# author: zhwang
# date: 2022-03-16
# difficulty: Medium
# categories: Array, Two Pointers
# link: <https://leetcode.com/problems/3sum-closest/description/>
# hidden: true
# ---
#
# Given an array `nums` of _n_ integers and an integer `target`, find three
# integers in `nums` such that the sum is closest to `target`. Return the sum of
# the three integers. You may assume that each input would have exactly one
# solution.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [-1,2,1,-4], target = 1
#     Output: 2
#     Explanation: The sum that is closest to the target is 2. (-1 + 2 + 1 = 2).
#
#
#
#
# **Constraints:**
#
#   * `3 <= nums.length <= 10^3`
#   * `-10^3 <= nums[i] <= 10^3`
#   * `-10^4 <= target <= 10^4`
#
#
## @lc code=start
using LeetCode

three_sum_closest(nums::Vector{Int}, target::Int) = three_sum_closest!(copy(nums), target)
function three_sum_closest!(nums::Vector{Int}, target::Int)::Int
    mindiff, n, res = Inf, length(sort!(nums)), 0
    for (i, first) in enumerate(nums)
        i > 1 && nums[i - 1] == first && continue
        left, right, newtarget = i + 1, n, target - first
        while left < right
            newsum = nums[left] + nums[right]
            newdiff = abs(newsum - newtarget)
            if newdiff < mindiff ## updates minimal difference
                newdiff == 0 && return target
                mindiff, res = newdiff, newsum + first
            elseif newsum > newtarget
                right -= 1
            else
                left += 1
            end
        end
    end
    res
end
## @lc code=end

## @lc test=start
@testset "16.3sum-closest.jl" begin
    @test three_sum_closest([-1,2,1,-4], 1) == 2
    @test three_sum_closest([0, 0, 0], 1) == 0
end
## @lc test=end

## @lc info=start
# link: [solution to 16]()
## @lc info=end


text:
using PyCall
function TextAlert(number,message)

py"""

def txt(num,mes):
    import smtplib
    from email.mime.text import MIMEText
    from email.mime.multipart import MIMEMultipart

    email = "*"
    pas = "*"

    sms_gateway = (str(num)+'@tmomail.net')
    smtp = "smtp.gmail.com"
    port = 587
    server = smtplib.SMTP(smtp,port)
    server.starttls()
    server.login(email,pas)

    #MIME module to structure message.
    msg = MIMEMultipart()
    msg['From'] = email
    msg['To'] = sms_gateway
    msg['Subject'] = "Julia Notification\n"
    body = (str(mes)+"\n")
    msg.attach(MIMEText(body, 'plain'))

    sms = msg.as_string()
    server.sendmail(email,sms_gateway,sms)
    server.quit()
"""
py"txt"(number,message)

end

max_stars_repo_path: Julia Monte Carlo Fitting/textalert.jl

max_issues_repo_path: Julia Monte Carlo Fitting/textalert.jl

max_forks_repo_path: Julia Monte Carlo Fitting/textalert.jl


text:
using Plots, Crayons.Box

using LinearAlgebra, Random, Optim, ForwardDiff

using Printf, ProgressMeter

using Distributions, Statistics

include("util/plot_ellipse.jl")
include("util/util.jl")
include("util/entropy.jl")
include("VI.jl")
include("VIdiag.jl")

gr()


text:
# AMDGPUnative run-time library
#
# This module defines method instances that will be compiled into a device-specific image
# and will be available to the AMDGPUnative compiler to call after Julia has generated code.

module Runtime

using ..AMDGPUnative
using LLVM
using LLVM.Interop

## representation of a runtime method instance

struct RuntimeMethodInstance
    def::Function

    return_type::Type
    types::Tuple
    name::Symbol

    # LLVM types cannot be cached, so we can't put them in the runtime method instance.
    # the actual types are constructed upon accessing them, based on a sentinel value:
    #  - nothing: construct the LLVM type based on its Julia counterparts
    #  - function: call this generator to get the type (when more control is needed)
    llvm_return_type::Union{Nothing, Function}
    llvm_types::Union{Nothing, Function}
    llvm_name::String
end

function Base.getproperty(rt::RuntimeMethodInstance, field::Symbol)
    value = getfield(rt, field)
    if field == :llvm_types
        if value == nothing
            LLVMType[convert.(LLVMType, typ) for typ in rt.types]
        else
            value()
        end
    elseif field == :llvm_return_type
        if value == nothing
            convert.(LLVMType, rt.return_type)
        else
            value()
        end
    else
        return value
    end
end

const methods = Dict{Symbol,RuntimeMethodInstance}()
get(name::Symbol) = methods[name]

# Register a Julia function `def` as a runtime library function identified by `name`. The
# function will be compiled upon first use for argument types `types` and should return
# `return_type`. Use `Runtime.get(name)` to get a reference to this method instance.
#
# The corresponding LLVM types `llvm_types` and `llvm_return_type` will be deduced from
# their Julia counterparts. To influence that conversion, pass a callable object instead;
# this object will be evaluated at run-time and the returned value will be used instead.
#
# When generating multiple runtime functions from a single definition, make sure to specify
# different values for `name`. The LLVM function name will be deduced from that name, but
# you can always specify `llvm_name` to influence that. Never use an LLVM name that starts
# with `julia_` or the function might clash with other compiled functions.
function compile(def, return_type, types, llvm_return_type=nothing, llvm_types=nothing;
                 name=typeof(def).name.mt.name, llvm_name="roc_$name")
    meth = RuntimeMethodInstance(def,
                                 return_type, types, name,
                                 llvm_return_type, llvm_types, llvm_name)
    if haskey(methods, name)
        error("Runtime function $name has already been registered!")
    end
    methods[name] = meth
    meth
end

## exception handling

function report_exception(ex)
    #= FIXME
    @rocprintf("""
        ERROR: a %s was thrown during kernel execution.
               Run Julia on debug level 2 for device stack traces.
        """, ex)
    =#
    return
end

compile(report_exception, Nothing, (Ptr{Cchar},))

function report_exception_name(ex)
    #= FIXME
    @rocprintf("""
        ERROR: a %s was thrown during kernel execution.
        Stacktrace:
        """, ex)
    =#
    return
end

function report_exception_frame(idx, func, file, line)
    # FIXME: @rocprintf(" [%i] %s at %s:%i\n", idx, func, file, line)
    return
end

compile(report_exception_frame, Nothing, (Cint, Ptr{Cchar}, Ptr{Cchar}, Cint))
compile(report_exception_name, Nothing, (Ptr{Cchar},))

# NOTE: no throw functions are provided here, but replaced by an LLVM pass instead
#       in order to provide some debug information without stack unwinding.

## GC

@enum AddressSpace begin
    Generic         = 1
    Tracked         = 100
    Derived         = 101
    CalleeRooted    = 102
    Loaded          = 103
end

# LLVM type of a tracked pointer
function T_prjlvalue()
    T_pjlvalue = convert(LLVMType, Any, true)
    LLVM.PointerType(eltype(T_pjlvalue), Tracked)
end

function gc_pool_alloc(sz::Csize_t)
    ptr = malloc(sz)
    return unsafe_pointer_to_objref(ptr)
end

compile(gc_pool_alloc, Any, (Csize_t,), T_prjlvalue)

## boxing and unboxing

const tag_type = UInt
const tag_size = sizeof(tag_type)

const gc_bits = 0x3 # FIXME

# get the type tag of a type at run-time
@generated function type_tag(::Val{type_name}) where type_name
    T_tag = convert(LLVMType, tag_type)
    T_ptag = LLVM.PointerType(T_tag)

    T_pjlvalue = convert(LLVMType, Any, true)

    # create function
    llvm_f, _ = create_function(T_tag)
    mod = LLVM.parent(llvm_f)

    # this isn't really a function, but we abuse it to get the JIT to resolve the address
    typ = LLVM.Function(mod, "jl_" * String(type_name) * "_type",
                        LLVM.FunctionType(T_pjlvalue))

    # generate IR
    Builder(JuliaContext()) do builder
        entry = BasicBlock(llvm_f, "entry", JuliaContext())
        position!(builder, entry)

        typ_var = bitcast!(builder, typ, T_ptag)

        tag = load!(builder, typ_var)

        ret!(builder, tag)
    end

    call_function(llvm_f, tag_type)
end

# we use `jl_value_ptr`, a Julia pseudo-intrinsic that can be used to box and unbox values

@generated function box(val, ::Val{type_name}) where type_name
    sz = sizeof(val)
    allocsz = sz + tag_size

    # type-tags are ephemeral, so look them up at run time
    #tag = unsafe_load(convert(Ptr{tag_type}, type_name))
    tag = :( type_tag(Val(type_name)) )

    quote
        Base.@_inline_meta

        ptr = malloc($(Csize_t(allocsz)))

        # store the type tag
        ptr = convert(Ptr{tag_type}, ptr)
        Core.Intrinsics.pointerset(ptr, $tag | $gc_bits, #=index=# 1, #=align=# $tag_size)

        # store the value
        ptr = convert(Ptr{$val}, ptr+tag_size)
        Core.Intrinsics.pointerset(ptr, val, #=index=# 1, #=align=# $sz)

        unsafe_pointer_to_objref(ptr)
    end
end

@inline function unbox(obj, ::Type{T}) where T
    ptr = ccall(:jl_value_ptr, Ptr{Cvoid}, (Any,), obj)

    # load the value
    ptr = convert(Ptr{T}, ptr)
    Core.Intrinsics.pointerref(ptr, #=index=# 1, #=align=# sizeof(T))
end

# generate functions functions that exist in the Julia runtime (see julia/src/datatype.c)
for (T, t) in [Int8   => :int8,  Int16  => :int16,  Int32  => :int32,  Int64  => :int64,
               UInt8  => :uint8, UInt16 => :uint16, UInt32 => :uint32, UInt64 => :uint64]
    box_fn   = Symbol("box_$t")
    unbox_fn = Symbol("unbox_$t")
    @eval begin
        $box_fn(val)   = box($T(val), Val($(QuoteNode(t))))
        $unbox_fn(obj) = unbox(obj, $T)

        compile($box_fn, Any, ($T,), T_prjlvalue; llvm_name=$"jl_$box_fn")
        compile($unbox_fn, $T, (Any,); llvm_name=$"jl_$unbox_fn")
    end
end

end


text:
using CorticalSpecrtralTemporalResponses
using AxisArrays
using DSP
const reasonable_response_maximum = 100

function bound(x,min,max)
  y = 1/(1+exp(-8((x - min)/(max - min) - 0.5)))
  miny = 1/(1+exp(-8(-0.5)))
  maxy = 1/(1+exp(-8(0.5)))
  min + (max-min)*(clamp(y,miny,maxy) - miny)/(maxy - miny)
end

function meanabs(A,n)
  result = reduce((x,y) -> x+abs(y),A,dims=n,init=zero(real(eltype(A))))
  result ./= size(A,n)
  dropdims(result,dims=n)
end

function findweights(condition,x,normalize_start)
  weights = if condition == :scales
    AxisArray(meanabs(x,axisdim(x,Axis{:freq})) .*
              ustrip.(uconvert.(cycoct,CorticalSpecrtralTemporalResponses.scales(x)))',
              AxisArrays.axes(x,Axis{:time}), AxisArrays.axes(x,Axis{:scale}))
  elseif condition ∈ [:freqs,:track]
    x
  end

  if normalize_start > 0s
    after_buildup = normalize_start .. last(times(weights))
    weights .-= minimum(weights[after_buildup])
    weights ./= maximum(weights[after_buildup])
  end
  weights
end

function remove_key_prefix!(prefix,dict)
  for key in keys(dict)
    if startswith(string(key),prefix)
      dict[Symbol(string(key)[length(prefix)+1:end])] = dict[key]
    end
  end
end

apply_bistable(x,args...;kwds...) = apply_bistable!(deepcopy(x),args...;kwds...)
function apply_bistable!(x,condition,params;
                         input_bound=(0.005,0.1),
                         lowpass=1.5,lowpass_order=3,
                         interactive=false,
                         intermediate_results=interactive,
                         normalize_start=-1.0s,
                         progressbar=interactive)

  if condition == :freqs
    remove_key_prefix!("f_",params)
  elseif condition == :scales
    remove_key_prefix!("s_",params)
  elseif condition == :track
    remove_key_prefix!("t_",params)
  end

  if iszero(params[:c_a]) && iszero(params[:c_m]) && iszero(params[:c_σ])
    return (result=x,)
  end

  noise_params = Dict(
    :τ_σ => params[:τ_σ], :c_σ => params[:c_σ]
  )
  adapt_params = Dict(
    :c_m => params[:c_m], :τ_m => params[:τ_m],
    :c_a => params[:c_a], :τ_a => params[:τ_a],
    :c_x => params[:c_x], :τ_x => params[:τ_x]
  )

  weights = findweights(condition,x,asseconds(normalize_start))
  weights .= bound.(weights,input_bound...)
  input_weights = copy(weights)

  wn = drift(weights;noise_params...,progressbar=progressbar)
  wna,a,m = adaptmi(wn;W_m=weighting(x,condition,params),
                    shape_y = x -> max(0,x),progressbar=progressbar,
                    adapt_params...)

  low = digitalfilter(Lowpass(lowpass;fs=ustrip(uconvert(Hz,1/Δt(wna)))),
                      Butterworth(lowpass_order))
  wna_low = filt!(similar(wna),low,wna)

  # shouldn't the below instead be replaced with:
  x .*= wna_low./weights

  if intermediate_results
    (result=x,inweights=input_weights,outweights=wna_low,adapt=a,inhibit=m)
  else
    (result=x,)
  end

end

max_stars_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_issues_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_forks_repo_name: JuliaTagBot/AuditoryBistabilityLE


text:
export null, isnull, all_ones, UndefValue, PointerNull

null(typ::LLVMType) = Value(API.LLVMConstNull(ref(typ)))

all_ones(typ::LLVMType) = Value(API.LLVMConstAllOnes(ref(typ)))

isnull(val::Value) = convert(Core.Bool, API.LLVMIsNull(ref(val)))

@checked struct UndefValue <: User
    ref::reftype(User)
end
identify(::Type{Value}, ::Val{API.LLVMUndefValueValueKind}) = UndefValue

UndefValue(typ::LLVMType) = UndefValue(API.LLVMGetUndef(ref(typ)))

abstract type Constant <: User end

# forward declarations
@checked struct Module
    ref::API.LLVMModuleRef
end
abstract type Instruction <: User end

@checked struct PointerNull <: Constant
    ref::reftype(User)
end
identify(::Type{Value}, ::Val{API.LLVMConstantPointerNullValueKind}) = PointerNull

PointerNull(typ::PointerType) = PointerNull(API.LLVMConstPointerNull(ref(typ)))

## scalar

export ConstantInt, ConstantFP

@checked struct ConstantInt <: Constant
    ref::reftype(Constant)
end
identify(::Type{Value}, ::Val{API.LLVMConstantIntValueKind}) = ConstantInt

# NOTE: fixed set for dispatch, also because we can't rely on sizeof(T)==width(T)
const WideInteger = Union{Int64, UInt64}
ConstantInt(typ::IntegerType, val::WideInteger, signed=false) =
    ConstantInt(API.LLVMConstInt(ref(typ), reinterpret(Culonglong, val),
                convert(Bool, signed)))
const SmallInteger = Union{Int8, Int16, Int32, UInt8, UInt16, UInt32}
ConstantInt(typ::IntegerType, val::SmallInteger, signed=false) =
    ConstantInt(typ, convert(Int64, val), signed)

function ConstantInt(typ::IntegerType, val::Integer, signed=false)
    valbits = ceil(Int, log2(abs(val))) + 1
    numwords = ceil(Int, valbits / 64)
    words = Vector{Culonglong}(undef, numwords)
    for i in 1:numwords
        words[i] = (val >> 64(i-1)) % Culonglong
    end
    return ConstantInt(API.LLVMConstIntOfArbitraryPrecision(ref(typ), numwords, words))
end

# NOTE: fixed set where sizeof(T) does match the numerical width
const SizeableInteger = Union{Int8, Int16, Int32, Int64, Int128, UInt8, UInt16, UInt32, UInt64, UInt128}
function ConstantInt(val::T, ctx::Context=GlobalContext()) where T<:SizeableInteger
    typ = IntType(sizeof(T)*8, ctx)
    return ConstantInt(typ, val, T<:Signed)
end

Base.convert(::Type{T}, val::ConstantInt) where {T<:Unsigned} =
    convert(T, API.LLVMConstIntGetZExtValue(ref(val)))

Base.convert(::Type{T}, val::ConstantInt) where {T<:Signed} =
    convert(T, API.LLVMConstIntGetSExtValue(ref(val)))

@checked struct ConstantFP <: Constant
    ref::reftype(Constant)
end
identify(::Type{Value}, ::Val{API.LLVMConstantFPValueKind}) = ConstantFP

ConstantFP(typ::LLVMDouble, val::Real) =
    ConstantFP(API.LLVMConstReal(ref(typ), Cdouble(val)))

Base.convert(::Type{Float64}, val::ConstantFP) =
    API.LLVMConstRealGetDouble(ref(val), Ref{API.LLVMBool}())

## aggregate

export ConstantAggregateZero

@checked struct ConstantAggregateZero <: Constant
    ref::reftype(Constant)
end
identify(::Type{Value}, ::Val{API.LLVMConstantAggregateZeroValueKind}) = ConstantAggregateZero

# there currently seems to be no function in the LLVM-C interface which returns a
# ConstantAggregateZero value directly, but values can occur through calls to LLVMConstNull

## constant expressions

export ConstantExpr, ConstantAggregate, ConstantArray, ConstantStruct, ConstantVector, InlineAsm

@checked struct ConstantExpr <: Constant
    ref::reftype(Constant)
end
identify(::Type{Value}, ::Val{API.LLVMConstantExprValueKind}) = ConstantExpr

abstract type ConstantAggregate <: Constant end

@checked struct ConstantArray <: ConstantAggregate
    ref::reftype(Constant)
end
identify(::Type{Value}, ::Val{API.LLVMConstantArrayValueKind}) = ConstantArray

@checked struct ConstantStruct <: ConstantAggregate
    ref::reftype(Constant)
end
identify(::Type{Value}, ::Val{API.LLVMConstantStructValueKind}) = ConstantStruct

@checked struct ConstantVector <: ConstantAggregate
    ref::reftype(Constant)
end
identify(::Type{Value}, ::Val{API.LLVMConstantVectorValueKind}) = ConstantVector

@checked struct InlineAsm <: Constant
    ref::reftype(Constant)
end
identify(::Type{Value}, ::Val{API.LLVMInlineAsmValueKind}) = InlineAsm

InlineAsm(typ::FunctionType, asm::String, constraints::String,
          side_effects::Core.Bool, align_stack::Core.Bool=false) =
    InlineAsm(API.LLVMConstInlineAsm(ref(typ), asm, constraints,
                                     convert(Bool, side_effects),
                                     convert(Bool, align_stack)))

## global values

abstract type GlobalValue <: Constant end

export GlobalValue,
       isdeclaration,
       linkage, linkage!,
       section, section!,
       visibility, visibility!,
       dllstorage, dllstorage!,
       unnamed_addr, unnamed_addr!,
       alignment, alignment!

parent(val::GlobalValue) = Module(API.LLVMGetGlobalParent(ref(val)))

isdeclaration(val::GlobalValue) = convert(Core.Bool, API.LLVMIsDeclaration(ref(val)))

linkage(val::GlobalValue) = API.LLVMGetLinkage(ref(val))
linkage!(val::GlobalValue, linkage::API.LLVMLinkage) =
    API.LLVMSetLinkage(ref(val), linkage)

function section(val::GlobalValue)
  #=
  The following started to fail on LLVM 4.0:
    Context() do ctx
      LLVM.Module("SomeModule", ctx) do mod
        st = LLVM.StructType("SomeType", ctx)
        ft = LLVM.FunctionType(st, [st])
        fn = LLVM.Function(mod, "SomeFunction", ft)
        section(fn) == ""
      end
      end
  =#
  section_ptr = API.LLVMGetSection(ref(val))
  return section_ptr != C_NULL ? unsafe_string(section_ptr) : ""
end
section!(val::GlobalValue, sec::String) = API.LLVMSetSection(ref(val), sec)

visibility(val::GlobalValue) = API.LLVMGetVisibility(ref(val))
visibility!(val::GlobalValue, viz::API.LLVMVisibility) =
    API.LLVMSetVisibility(ref(val), viz)

dllstorage(val::GlobalValue) = API.LLVMGetDLLStorageClass(ref(val))
dllstorage!(val::GlobalValue, storage::API.LLVMDLLStorageClass) =
    API.LLVMSetDLLStorageClass(ref(val), storage)

unnamed_addr(val::GlobalValue) = convert(Core.Bool, API.LLVMHasUnnamedAddr(ref(val)))
unnamed_addr!(val::GlobalValue, flag::Core.Bool) =
    API.LLVMSetUnnamedAddr(ref(val), convert(Bool, flag))

const AlignedValue = Union{GlobalValue,Instruction}   # load, store, alloca
alignment(val::AlignedValue) = API.LLVMGetAlignment(ref(val))
alignment!(val::AlignedValue, bytes::Integer) = API.LLVMSetAlignment(ref(val), bytes)

## global variables

abstract type GlobalObject <: GlobalValue end

export GlobalVariable, unsafe_delete!,
       initializer, initializer!,
       isthreadlocal, threadlocal!,
       threadlocalmode, threadlocalmode!,
       isconstant, constant!,
       isextinit, extinit!

@checked struct GlobalVariable <: GlobalObject
    ref::reftype(GlobalObject)
end
identify(::Type{Value}, ::Val{API.LLVMGlobalVariableValueKind}) = GlobalVariable

GlobalVariable(mod::Module, typ::LLVMType, name::String) =
    GlobalVariable(API.LLVMAddGlobal(ref(mod), ref(typ), name))

GlobalVariable(mod::Module, typ::LLVMType, name::String, addrspace::Integer) =
    GlobalVariable(API.LLVMAddGlobalInAddressSpace(ref(mod), ref(typ),
                                                   name, addrspace))

unsafe_delete!(::Module, gv::GlobalVariable) = API.LLVMDeleteGlobal(ref(gv))

initializer(gv::GlobalVariable) =
  Value(API.LLVMGetInitializer(ref(gv)))
initializer!(gv::GlobalVariable, val::Constant) =
  API.LLVMSetInitializer(ref(gv), ref(val))

isthreadlocal(gv::GlobalVariable) = convert(Core.Bool, API.LLVMIsThreadLocal(ref(gv)))
threadlocal!(gv::GlobalVariable, bool) =
  API.LLVMSetThreadLocal(ref(gv), convert(Bool, bool))

isconstant(gv::GlobalVariable) = convert(Core.Bool, API.LLVMIsGlobalConstant(ref(gv)))
constant!(gv::GlobalVariable, bool) =
  API.LLVMSetGlobalConstant(ref(gv), convert(Bool, bool))

threadlocalmode(gv::GlobalVariable) = API.LLVMGetThreadLocalMode(ref(gv))
threadlocalmode!(gv::GlobalVariable, mode) =
  API.LLVMSetThreadLocalMode(ref(gv), mode)

isextinit(gv::GlobalVariable) =
  convert(Core.Bool, API.LLVMIsExternallyInitialized(ref(gv)))
extinit!(gv::GlobalVariable, bool) =
  API.LLVMSetExternallyInitialized(ref(gv), convert(Bool, bool))


text:
using Test

using Plots, LatPhysPlottingPlots
using LatPhysUnitcellLibrary, LatPhysLatticeConstruction

# 2D
uc = getUnitcellHoneycomb()
lt = getLatticePeriodic(uc, 5)
plot(lt)

# 3D
uc = getUnitcell_10_3_b()
lt = getLatticePeriodic(uc, 5)
plot(lt)


text:
# https://juliagraphics.github.io/Gtk.jl/latest/manual/textwidgets/#Label-1

using Gtk

# A GtkLabel is the most basic text widget that has already been used behind
# the scene in any previous example involving a GtkButton. A GtkLabel is
# constructed by calling
#
# label = GtkLabel("My text")
# The text of a label can be changed using
#
# GAccessor.text(label,"My other text")

win = GtkWindow("Text Widgets", 400, 200)
label = GtkLabel("My text")
push!(win, label)
showall(win)

GAccessor.text(label, "My other text")

# Furthermore, a label has limited support for adding formatted text. This is
# done using the markup function:
#
# GAccessor.markup(label,"""My bold text\n
#                           GTK+ website""")
# The syntax for this markup text is borrowed from html and explained here:
# https://developer.gnome.org/pygtk/stable/pango-markup-language.html

GAccessor.markup(label,
"""My bold text\n
GTK websitee""")

#
# A label can be made selectable using
#
# GAccessor.selectable(label,true)
# This can be used if the user should be allowed to copy the text.
#

GAccessor.selectable(label, true)

# The justification of a label can be changed in the following way:
#
# GAccessor.justify(label,Gtk.GConstants.GtkJustification.RIGHT)
# Possible values of the enum GtkJustification are LEFT,RIGHT,CENTER, and FILL.

GAccessor.justify(label, Gtk.GConstants.GtkJustification.RIGHT)

#
# Automatic line wrapping can be enabled using
#
# GAccessor.text(label,repeat("Very long text! ",20))
# GAccessor.line_wrap(label,true)
# Note that this will only happen, if the size of the widget is limited using
# layout constraints.

GAccessor.text(label, repeat("Very long text! ", 20))
GAccessor.line_wrap(label, true)


text:
const INFINITY = -1
const INVALID_NODE_INDEX = -1

abstract type AbstractNode end

struct VariableNode <: AbstractNode
    index::Int
end

struct ConstraintNode <: AbstractNode
    index::Int
end

struct ObjectiveNode <: AbstractNode
    index::Int
end

abstract type AbstractEdge end

struct Edge <: AbstractEdge
    bridge_index::Int
    added_variables::Vector{VariableNode}
    added_constraints::Vector{ConstraintNode}
end

struct ObjectiveEdge <: AbstractEdge
    bridge_index::Int
    added_variables::Vector{VariableNode}
    added_constraints::Vector{ConstraintNode}
    added_objective::ObjectiveNode
end

mutable struct Graph
    variable_edges::Vector{Vector{Edge}}
    variable_constraint_node::Vector{ConstraintNode}
    variable_constraint_cost::Vector{Int}
    # variable node index -> Number of bridges that need to be used
    variable_dist::Vector{Int}
    # variable node index -> Index of bridge to be used
    variable_best::Vector{Int}
    variable_last_correct::Int
    constraint_edges::Vector{Vector{Edge}}
    # constraint node index -> Number of bridges that need to be used
    constraint_dist::Vector{Int}
    # constraint node index -> Index of bridge to be used
    constraint_best::Vector{Int}
    constraint_last_correct::Int
    objective_edges::Vector{Vector{ObjectiveEdge}}
    # objective node index -> Number of bridges that need to be used
    objective_dist::Vector{Int}
    # objective node index -> Index of bridge to be used
    objective_best::Vector{Int}
    objective_last_correct::Int
end

function Graph()
    return Graph(
        Vector{Edge}[],
        ConstraintNode[],
        Int[],
        Int[],
        Int[],
        0,
        Vector{Edge}[],
        Int[],
        Int[],
        0,
        Vector{ObjectiveEdge}[],
        Int[],
        Int[],
        0,
    )
end

function Base.show(io::IO, graph::Graph)
    print(io, "Bridge graph with ")
    print(io, length(graph.variable_best), " variable nodes, ")
    print(io, length(graph.constraint_best), " constraint nodes and ")
    return print(io, length(graph.objective_best), " objective nodes.")
end

function variable_nodes(graph::Graph)
    return MOIU.LazyMap{VariableNode}(
        i -> VariableNode(i),
        eachindex(graph.variable_best),
    )
end

function constraint_nodes(graph::Graph)
    return MOIU.LazyMap{ConstraintNode}(
        i -> ConstraintNode(i),
        eachindex(graph.constraint_best),
    )
end

function objective_nodes(graph::Graph)
    return MOIU.LazyMap{ObjectiveNode}(
        i -> ObjectiveNode(i),
        eachindex(graph.objective_best),
    )
end

# After `add_bridge(b, BT)`, some constrained variables `(S,)` in
# `keys(b.variable_best)` or constraints `(F, S)` in `keys(b.constraint_best)`
# or `(F,)` in `keys(b.objective_best)` may be bridged
# with less bridges than `b.variable_dist[(S,)]`,
# `b.constraint_dist[(F, S)]` or `b.objective_dist[(F,)]` using `BT`.
# We could either recompute the distance from every node or clear the
# dictionary so that the distance is computed lazily at the next `supports_constraint`
# call. We prefer clearing the dictionaries so as this is called for each
# bridge added and recomputing the distance for each bridge would be a wase
# if several bridges are added consecutively.
function Base.empty!(graph::Graph)
    empty!(graph.variable_edges)
    empty!(graph.variable_constraint_node)
    empty!(graph.variable_constraint_cost)
    empty!(graph.variable_dist)
    empty!(graph.variable_best)
    graph.variable_last_correct = 0
    empty!(graph.constraint_edges)
    empty!(graph.constraint_dist)
    empty!(graph.constraint_best)
    graph.constraint_last_correct = 0
    empty!(graph.objective_edges)
    empty!(graph.objective_dist)
    empty!(graph.objective_best)
    graph.objective_last_correct = 0
    return
end

function add_edge(graph::Graph, node::VariableNode, edge::Edge)
    push!(graph.variable_edges[node.index], edge)
    return
end

function add_edge(graph::Graph, node::ConstraintNode, edge::Edge)
    push!(graph.constraint_edges[node.index], edge)
    return
end

function add_edge(graph::Graph, node::ObjectiveNode, edge::ObjectiveEdge)
    push!(graph.objective_edges[node.index], edge)
    return
end

function add_variable_node(graph::Graph)
    push!(graph.variable_edges, Edge[])
    # Use an invalid index so that the code errors instead return something
    # incorrect in case `set_variable_constraint_node` is not called.
    push!(graph.variable_constraint_node, ConstraintNode(INVALID_NODE_INDEX))
    push!(graph.variable_constraint_cost, 0)
    push!(graph.variable_dist, INFINITY)
    push!(graph.variable_best, 0)
    return VariableNode(length(graph.variable_best))
end

function set_variable_constraint_node(
    graph::Graph,
    variable_node::VariableNode,
    constraint_node::ConstraintNode,
    cost::Int,
)
    graph.variable_constraint_node[variable_node.index] = constraint_node
    graph.variable_constraint_cost[variable_node.index] = cost
    return
end

function add_constraint_node(graph::Graph)
    push!(graph.constraint_edges, Edge[])
    push!(graph.constraint_dist, INFINITY)
    push!(graph.constraint_best, 0)
    return ConstraintNode(length(graph.constraint_best))
end

function add_objective_node(graph::Graph)
    push!(graph.objective_edges, ObjectiveEdge[])
    push!(graph.objective_dist, INFINITY)
    push!(graph.objective_best, 0)
    return ObjectiveNode(length(graph.objective_best))
end

function bridging_cost(graph::Graph, node::AbstractNode)
    bellman_ford!(graph)
    dist = _dist(graph, node)
    return dist == INFINITY ? Inf : float(dist)
end

function bridge_index(graph::Graph, node::VariableNode)
    bellman_ford!(graph)
    return graph.variable_best[node.index]
end

function bridge_index(graph::Graph, node::ConstraintNode)
    bellman_ford!(graph)
    return graph.constraint_best[node.index]
end

function bridge_index(graph::Graph, node::ObjectiveNode)
    bellman_ford!(graph)
    return graph.objective_best[node.index]
end

"""
    is_variable_edge_best(graph::Graph, node::VariableNode)

Return a `Bool` indicating whether the value of `_dist(graph, node)` can be
achieved with a variable bridge.
"""
function is_variable_edge_best(graph::Graph, node::VariableNode)
    bellman_ford!(graph)
    return graph.variable_dist[node.index] == _dist(graph, node)
end

# Update `b.variable_dist`, `b.constraint_dist` `b.variable_best` and
# `b.constraint_best` for constrained variable types in `variables` and
# constraint types in `constraints`.
function updated_dist(graph::Graph, current::Int, edges::Vector{<:AbstractEdge})
    bridge_index = 0
    for edge in edges
        if supports_added_no_update(graph, edge)
            # `added_dist != INFINITY` because `supports_added_no_update`
            # returned `true`.
            dist = 1 + added_dist(graph, edge)
            if current == INFINITY || dist < current
                current = dist
                bridge_index = edge.bridge_index
            end
        end
    end
    return current, bridge_index
end

function bellman_ford!(graph::Graph)
    # Bellman-Ford algorithm
    changed = true # Has b.constraint_dist changed in the last iteration ?
    while changed
        changed = false
        for i in (graph.variable_last_correct+1):length(graph.variable_best)
            dist, best = updated_dist(
                graph,
                graph.variable_dist[i],
                graph.variable_edges[i],
            )
            if !iszero(best)
                graph.variable_dist[i] = dist
                graph.variable_best[i] = best
                changed = true
            end
        end
        for i in (graph.constraint_last_correct+1):length(graph.constraint_best)
            dist, best = updated_dist(
                graph,
                graph.constraint_dist[i],
                graph.constraint_edges[i],
            )
            if !iszero(best)
                graph.constraint_dist[i] = dist
                graph.constraint_best[i] = best
                changed = true
            end
        end
        for i in (graph.objective_last_correct+1):length(graph.objective_best)
            dist, best = updated_dist(
                graph,
                graph.objective_dist[i],
                graph.objective_edges[i],
            )
            if !iszero(best)
                graph.objective_dist[i] = dist
                graph.objective_best[i] = best
                changed = true
            end
        end
    end
    graph.variable_last_correct = length(graph.variable_best)
    graph.constraint_last_correct = length(graph.constraint_best)
    graph.objective_last_correct = length(graph.objective_best)
    return
end

function _dist(graph::Graph, node::VariableNode)
    if iszero(node.index)
        return 0
    end
    constraint_node = graph.variable_constraint_node[node.index]
    # If free variables are bridged but the functionize bridge was not added,
    # constraint_node is `ConstraintNode(INVALID_NODE_INDEX)`.
    dc =
        constraint_node.index == INVALID_NODE_INDEX ? INFINITY :
        _dist(graph, constraint_node)
    dv = graph.variable_dist[node.index]
    if dc == INFINITY
        return dv
    else
        dc += graph.variable_constraint_cost[node.index]
        if dv == INFINITY
            return dc
        else
            return min(dv, dc)
        end
    end
end

function _dist(graph::Graph, node::ConstraintNode)
    return iszero(node.index) ? 0 : graph.constraint_dist[node.index]
end

function _dist(graph::Graph, node::ObjectiveNode)
    return iszero(node.index) ? 0 : graph.objective_dist[node.index]
end

function _sum_dist(graph::Graph, nodes::Vector{<:AbstractNode})
    return mapreduce(node -> _dist(graph, node), +, nodes, init = 0)
end

function added_dist(graph::Graph, edge::Edge)
    return _sum_dist(graph, edge.added_variables) +
           _sum_dist(graph, edge.added_constraints)
end

function added_dist(graph::Graph, edge::ObjectiveEdge)
    return _sum_dist(graph, edge.added_variables) +
           _sum_dist(graph, edge.added_constraints) +
           _dist(graph, edge.added_objective)
end

function supports_no_update(graph::Graph, node::AbstractNode)
    return iszero(node.index) || _dist(graph, node) != INFINITY
end

function _supports_all_no_update(graph::Graph, nodes::Vector{<:AbstractNode})
    return all(node -> supports_no_update(graph, node), nodes)
end

function supports_added_no_update(graph::Graph, edge::Edge)
    return _supports_all_no_update(graph, edge.added_variables) &&
           _supports_all_no_update(graph, edge.added_constraints)
end

function supports_added_no_update(graph::Graph, edge::ObjectiveEdge)
    return _supports_all_no_update(graph, edge.added_variables) &&
           _supports_all_no_update(graph, edge.added_constraints) &&
           supports_no_update(graph, edge.added_objective)
end