text:
function raw_ref(model::BilevelModel, idx::Int)
    if haskey(model.ctr_upper, idx)
        return model.ctr_upper[idx]
    elseif haskey(model.ctr_lower, idx)
        return model.ctr_lower[idx]
    else
        error("Index $(idx) does not belong to BilevelModel")
    end
end
raw_ref(cref::BilevelConstraintRef) = raw_ref(cref.model, cref.index)

function BilevelConstraintRef(model, idx)
    raw = raw_ref(model, idx)
    return JuMP.ConstraintRef(model, idx, raw.shape)
end
JuMP.constraint_type(::AbstractBilevelModel) = BilevelConstraintRef
level(cref::BilevelConstraintRef) = cref.model.ctr_level[cref.index]
function JuMP.add_constraint(::BilevelModel, ::JuMP.AbstractConstraint, ::String="")
    error(
        "Can't add constraint directly to the bilevel model `m`, "*
        "attach the constraint to the upper or lower model "*
        "with @constraint(Upper(m), ...) or @constraint(Lower(m), ...)")
end
function JuMP.constraint_object(con_ref::ConstraintRef{BilevelModel, Int})
    raw = raw_ref(con_ref)
    return JuMP.constraint_object(raw)
end
# JuMP.add_constraint(m::UpperModel, c::JuMP.VectorConstraint, name::String="") =
#     error("no vec ctr")
function JuMP.add_constraint(m::InnerBilevelModel, c::Union{JuMP.ScalarConstraint{F,S},JuMP.VectorConstraint{F,S}}, name::String="") where {F,S}
    blm = bilevel_model(m)
    blm.nextconidx += 1
    cref = JuMP.ConstraintRef(blm, blm.nextconidx, JuMP.shape(c))
    func = JuMP.jump_function(c)
    level_func = replace_variables(func, bilevel_model(m), mylevel_var_list(m), level(m))
    level_c = JuMP.build_constraint(error, level_func, c.set)
    level_cref = JuMP.add_constraint(mylevel_model(m), level_c, name)
    blm.ctr_level[cref.index] = level(m)
    mylevel_ctr_list(m)[cref.index] = level_cref
    blm.constraints[cref.index] = c
    blm.ctr_info[cref.index] = empty_info(c)
    JuMP.set_name(cref, name)
    cref
end

JuMP.is_valid(m::BilevelModel, cref::BilevelConstraintRef) = cref.index in keys(m.constraints)
JuMP.is_valid(m::InnerBilevelModel, cref::BilevelConstraintRef) =
    JuMP.is_valid(bilevel_model(m), cref) && level(cref) == level(m)
function JuMP.constraint_object(cref::BilevelConstraintRef, F::Type, S::Type)
    c = cref.model.constraints[cref.index]
    # `TypeError` should be thrown is `F` and `S` are not correct
    # This is needed for the tests in `constraints.jl`
    c.func::F
    c.set::S
    c
end
function empty_info(c::JuMP.ScalarConstraint{F,S}) where {F,S}
    return ConstraintInfo{Float64}()
end
function empty_info(c::JuMP.VectorConstraint{F,S}) where {F,S}
    return ConstraintInfo{Vector{Float64}}(MOI.dimension(c.set))
end

function JuMP.set_dual_start_value(cref::BilevelConstraintRef, value::T) where T<:Number
    cref.model.ctr_info[cref.index].start = value
end
function JuMP.set_dual_start_value(cref::BilevelConstraintRef, value::T) where T<:Vector{S} where S
    array = cref.model.ctr_info[cref.index].start
    @assert length(array) == length(value)
    copyto!(array, value)
end
function JuMP.dual_start_value(cref::BilevelConstraintRef)
    cref.model.ctr_info[cref.index].start
end
function set_dual_upper_bound(cref::BilevelConstraintRef, value::T) where T<:Number
    cref.model.ctr_info[cref.index].upper = value
end
function set_dual_upper_bound(cref::BilevelConstraintRef, value::T) where T<:Vector{S} where S
    array = cref.model.ctr_info[cref.index].upper
    @assert length(array) == length(value)
    copyto!(array, value)
end
function get_dual_upper_bound(cref::BilevelConstraintRef)
    cref.model.ctr_info[cref.index].upper
end
function set_dual_lower_bound(cref::BilevelConstraintRef, value::T) where T<:Number
    cref.model.ctr_info[cref.index].lower = value
end
function set_dual_lower_bound(cref::BilevelConstraintRef, value::T) where T<:Vector{S} where S
    array = cref.model.ctr_info[cref.index].lower
    @assert length(array) == length(value)
    copyto!(array, value)
end
function get_dual_lower_bound(cref::BilevelConstraintRef)
    cref.model.ctr_info[cref.index].lower
end

function set_primal_upper_bound_hint(vref::BilevelVariableRef, value::T) where T<:Number
    vref.model.var_info[vref.idx].upper = value
end
function get_primal_upper_bound_hint(vref::BilevelVariableRef)
    vref.model.var_info[vref.idx].upper
end
function set_primal_lower_bound_hint(vref::BilevelVariableRef, value::T) where T<:Number
    vref.model.var_info[vref.idx].lower = value
end
function get_primal_lower_bound_hint(vref::BilevelVariableRef)
    vref.model.var_info[vref.idx].lower
end

function JuMP.value(cref::BilevelConstraintRef; result::Int = 1)
    if level(cref) == BilevelJuMP.LOWER_ONLY
        # Constraint index on the lower model
        con_lower_idx = cref.model.ctr_lower[cref.index].index
        # Single bilevel model constraint associated with the lower level constraint
        con_sblm_idx = cref.model.lower_to_sblm[con_lower_idx]
    else
        # Constraint index on the lower model
        con_upper_idx = cref.model.ctr_upper[cref.index].index
        # Single bilevel model constraint associated with the lower level constraint
        con_sblm_idx = cref.model.upper_to_sblm[con_upper_idx]
    end
    # Solver constraint associated with the single bilevel model constraint
    con_solver_idx = cref.model.sblm_to_solver[con_sblm_idx]
    return MOI.get(cref.model.solver, MOI.ConstraintPrimal(result), con_solver_idx)
end
# variables again (duals)
# code for using dual variables associated with lower level constraints
# in the upper level

function JuMP.num_constraints(model::BilevelModel)
    return length(model.ctr_level)
end
function JuMP.num_constraints(model::LowerModel)
    return length(model.m.ctr_lower)
end
function JuMP.num_constraints(model::UpperModel)
    return length(model.m.ctr_upper)
end
function JuMP.num_constraints(model::InnerBilevelModel, f, s)
    return JuMP.num_constraints(mylevel_model(model), f, s)
end
function JuMP.num_constraints(model::BilevelModel, f, s)
    return JuMP.num_constraints(Upper(model), f, s) +
        JuMP.num_constraints(Lower(model), f, s)
end

struct DualOf
    ci::BilevelConstraintRef
end
struct DualVariableInfo
    info::JuMP.VariableInfo
    ci::BilevelConstraintRef
end
function JuMP.build_variable(
    _error::Function,
    info::JuMP.VariableInfo,
    dual_of::DualOf;
    extra_kw_args...,
)

    if level(dual_of.ci) != LOWER_ONLY
        error("Variables can only be tied to LOWER level constraints, got $(dual_of.ci.level) level")
    end
    for (kwarg, _) in extra_kw_args
        _error("Unrecognized keyword argument $kwarg")
    end

    if info.has_lb
        set_dual_lower_bound(dual_of.ci, info.lower_bound)
        # info.has_lb = false
        # info.lower_bound = NaN
    end
    if info.has_ub
        set_dual_upper_bound(dual_of.ci, info.upper_bound)
        # info.has_ub = false
        # info.upper_bound = NaN
    end
    info.has_fix   && _error("Dual variable does not support fixing")
    if info.has_start
        JuMP.set_dual_start_value(dual_of.ci, info.start)
        # info.has_start = false
        # info.start = NaN
    end
    info.binary    && _error("Dual variable cannot be binary")
    info.integer   && _error("Dual variable cannot be integer")

    info = JuMP.VariableInfo(false, NaN, false, NaN,
                             false, NaN, false, NaN,
                             false, false)

    return DualVariableInfo(
        info,
        dual_of.ci
    )
end
function JuMP.add_variable(inner::UpperModel, dual_info::DualVariableInfo, name::String="")
    # TODO vector version
    m = bilevel_model(inner)
    m.nextvaridx += 1
    vref = BilevelVariableRef(m, m.nextvaridx, DUAL_OF_LOWER)
    v_upper = JuMP.add_variable(m.upper, JuMP.ScalarVariable(dual_info.info), name)
    m.var_upper[vref.idx] = v_upper
    m.var_level[vref.idx] = DUAL_OF_LOWER
    m.upper_var_to_lower_ctr_link[v_upper] = m.ctr_lower[dual_info.ci.index] # TODO improve this
    m.variables[vref.idx] = JuMP.ScalarVariable(dual_info.info)
    JuMP.set_name(vref, name)
    m.var_info[vref.idx] = empty_info(dual_info)
    vref
end
function empty_info(::DualVariableInfo)
    return VariableInfo()
end

function get_constrain_ref(vref::BilevelVariableRef)
    model = vref.model
    ctr_ref = model.upper_var_to_lower_ctr_link[model.var_upper[vref.idx]]
    idx = -1
    for (ind, ref) in model.ctr_lower
        if ref == ctr_ref
            idx = ind
        end
    end
    @assert idx != -1
    return BilevelConstraintRef(model, idx)
end

function JuMP.dual(cref::BilevelConstraintRef)
    # Right now this code assumes there is no possibility for vectorized constraints
    if level(cref) == BilevelJuMP.LOWER_ONLY
        # Constraint index on the lower model
        con_lower_ref = cref.model.ctr_lower[cref.index]
        con_lower_idx = con_lower_ref.index
        # Dual variable associated with constraint index
        model_var_idxs = cref.model.lower_primal_dual_map.primal_con_dual_var[con_lower_idx]
        # Single bilevel model variable associated with the dual variable
        sblm_var_idxs = MOI.VariableIndex[]
        for vi in model_var_idxs
            push!(sblm_var_idxs, cref.model.lower_dual_to_sblm[vi])
        end
        # Solver variable associated withe the sblm model
        solver_var_idxs = MOI.VariableIndex[]
        for vi in sblm_var_idxs
            push!(solver_var_idxs, cref.model.sblm_to_solver[vi])
        end
        pre_duals = MOI.get(cref.model.solver, MOI.VariablePrimal(), solver_var_idxs)
        return JuMP.reshape_vector(
            pre_duals,
            JuMP.dual_shape(con_lower_ref.shape)
            )
    elseif level(cref) == BilevelJuMP.UPPER_ONLY
        m = cref.model
        con_upper_ref = cref.model.ctr_upper[cref.index]
        solver_ctr_idx = m.sblm_to_solver[m.upper_to_sblm[JuMP.index(con_upper_ref)]]
        pre_duals = MOI.get(cref.model.solver, MOI.ConstraintDual(), solver_ctr_idx)
        return JuMP.reshape_vector(
            pre_duals,
            JuMP.dual_shape(con_upper_ref.shape)
            )
    else
        error("Dual solutions of upper level constraints are not available. Either the solution method does nto porvide duals or or the solver failed to get one.")
    end
end

function JuMP.normalized_rhs(cref::BilevelConstraintRef)
    return JuMP.normalized_rhs(raw_ref(cref))
end

function JuMP.set_normalized_rhs(cref::BilevelConstraintRef, val)
    return JuMP.set_normalized_rhs(raw_ref(cref), val)
end

function JuMP.add_to_function_constant(cref::BilevelConstraintRef, val)
    return JuMP.add_to_function_constant(raw_ref(cref), val)
end

function JuMP.normalized_coefficient(cref::BilevelConstraintRef, var::BilevelVariableRef)
    cidx = cref.index
    model = cref.model
    level = model.ctr_level[cidx]
    vidx = var.idx
    level_var = if level == UPPER_ONLY
        model.var_upper[vidx]
    else
        model.var_lower[vidx]
    end
    return JuMP.normalized_coefficient(raw_ref(cref), level_var)
end

function JuMP.set_normalized_coefficient(
    cref::BilevelConstraintRef, var::BilevelVariableRef, val)
    cidx = cref.index
    model = cref.model
    level = model.ctr_level[cidx]
    vidx = var.idx
    level_var = if level == UPPER_ONLY
        model.var_upper[vidx]
    else
        model.var_lower[vidx]
    end
    return JuMP.set_normalized_coefficient(raw_ref(cref), level_var, val)
end

function JuMP.list_of_constraint_types(
    model::InnerBilevelModel,
)::Vector{Tuple{DataType,DataType}}
    return JuMP.list_of_constraint_types(mylevel_model(model))
end
function JuMP.list_of_constraint_types(
    model::BilevelModel,
)::Vector{Tuple{DataType,DataType}}
    return unique!(vcat(
        JuMP.list_of_constraint_types(Upper(model)),
        JuMP.list_of_constraint_types(Lower(model)),
    ))
end

function JuMP.all_constraints(model::BilevelModel, f, s)
    return unique!(vcat(
        JuMP.all_constraints(Upper(model), f, s),
        JuMP.all_constraints(Lower(model), f, s),
    ))
end

function JuMP.all_constraints(model::InnerBilevelModel, f, s)
    build_reverse_ctr_map!(model)
    m = mylevel_model(model)
    list = JuMP.all_constraints(m, f, s)
    get_reverse_ctr_map.(model, list)
end
function build_reverse_ctr_map!(um::UpperModel)
    m = bilevel_model(um)
    m.ctr_upper_rev = Dict{JuMP.ConstraintRef, JuMP.ConstraintRef}()
    for (idx, ref) in m.ctr_upper
        m.ctr_upper_rev[ref] = BilevelConstraintRef(m, idx)
    end
end
function build_reverse_ctr_map!(lm::LowerModel)
    m = bilevel_model(lm)
    m.ctr_lower_rev = Dict{JuMP.ConstraintRef, JuMP.ConstraintRef}()
    for (idx, ref) in m.ctr_lower
        m.ctr_lower_rev[ref] = BilevelConstraintRef(m, idx)
    end
    return nothing
end
get_reverse_ctr_map(m::UpperModel, idx) = m.m.ctr_upper_rev[idx]
get_reverse_ctr_map(m::LowerModel, idx) = m.m.ctr_lower_rev[idx]

function JuMP.delete(mod::BilevelModel, cref::BilevelConstraintRef)
    model = cref.model
    @assert model === mod
    idx = cref.index
    delete!(model.constraints, idx)
    delete!(model.connames, idx)
    model.need_rebuild_names_ctr = true
    delete!(model.ctr_level, idx)
    if haskey(model.ctr_upper, idx)
        c_up = model.ctr_upper[idx]
        delete!(model.ctr_upper, idx)
        JuMP.delete(model.upper, c_up)
    end
    if haskey(model.ctr_lower, idx)
        c_lo = model.ctr_lower[idx]
        delete!(model.ctr_lower, idx)
        JuMP.delete(model.lower, c_lo)
    end
    delete!(model.ctr_info, idx)
    model.ctr_upper_rev = nothing
    model.ctr_lower_rev = nothing
    return nothing
end


text:
## Basic heuristics ##
export NullHeuristic, GoalCountHeuristic, ManhattanHeuristic

"Null heuristic that always returns zero."
struct NullHeuristic <: Heuristic end

Base.hash(::NullHeuristic, h::UInt) = hash(NullHeuristic, h)

compute(h::NullHeuristic, domain::Domain, state::State, spec::Specification) = 0

"Heuristic that counts the number of goals un/satisfied."
struct GoalCountHeuristic <: Heuristic
    dir::Symbol # :forward or :backward
    GoalCountHeuristic() = new(:forward)
    GoalCountHeuristic(dir) = new(dir)
end

Base.hash(::GoalCountHeuristic, h::UInt) = hash(GoalCountHeuristic, h)

function compute(h::GoalCountHeuristic,
                 domain::Domain, state::State, spec::Specification)
    goals = get_goal_terms(spec)
    count = sum([!state[domain, g] for g in goals])
    return h.dir == :backward ? length(goals) - count : count
end

"Computes Manhattan distance to the goal for the specified numeric fluents."
mutable struct ManhattanHeuristic <: Heuristic
    fluents::Vector{Term}
    goal_state::State
    pre_key::UInt64 # Key to check if information needs to be precomputed again
    ManhattanHeuristic(fluents) = new(fluents)
end

Base.hash(heuristic::ManhattanHeuristic, h::UInt) =
    hash(heuristic.fluents, hash(ManhattanHeuristic, h))

function precompute!(h::ManhattanHeuristic,
                     domain::Domain, state::State, spec::Specification)
    h.goal_state = State(get_goal_terms(spec))
    h.pre_key = objectid(spec)
    return h
end

function is_precomputed(h::ManhattanHeuristic,
                        domain::Domain, state::State, spec::Specification)
    return (isdefined(h, :goal_state) && objectid(spec) == h.pre_key)
end

function compute(h::ManhattanHeuristic,
                 domain::Domain, state::State, spec::Specification)
    # Precompute if necessary
    if !is_precomputed(h, domain, state, spec)
        precompute!(h, domain, state, spec) end
    # Compute L1 distance between fluents in the current and goal state
    goal_vals = [h.goal_state[domain, f] for f in h.fluents]
    curr_vals = [state[domain, f] for f in h.fluents]
    dist = sum(abs.(goal_vals - curr_vals))
    return dist
end

ext: jl


text:
"""
	printLog(simulator::Quantum2dFractalSimulator; printlog="none")

# option `printtlog`:

* "none" : default - print nothing
* "coef"
* "norm"
* "label" : print the labels for coeficients
* "mag"
* "magmute" : calculate magnetization for every renorm steps without printing

"""
function printLog(simulator::Quantum2dFractalSimulator; printlog="none")
	if printlog=="coef"
		printCoefficients(simulator)
	elseif printlog=="norm"
		printNormalizationFactor(simulator)
	elseif printlog=="label"
		printCoefficientsLabel(simulator)
	elseif printlog=="mag"
		magnetization = getExpectationValue(simulator)
		println(magnetization)
	elseif printlog=="magmute"
		magnetization = getExpectationValue(simulator)
	end
end

function printCoefficientsLabel(simulator::Quantum2dFractalSimulator)
	print("count", "\t")
	for i = 1:getLengthOfNormList(simulator)
		print(getListOfNormFactors(simulator)[i],"\t")
	end
	println()
end

function printCoefficients(simulator::Quantum2dFractalSimulator)
	iteration = getWholeiteration(simulator) +2 -getCount(simulator)
	lengCoef = getLengthOfNormList(simulator)
	print(iteration, "\t")
	for i = 1:lengCoef
		@printf "%.5e\t" simulator.coefficients[iteration,i]
	end
	println()
end

function printNormalizationFactor(simulator::Quantum2dFractalSimulator)
	iteration = getCount(simulator)
	lengCoef = getLengthOfNormList(simulator)
	print(iteration, "\t")
	for i = 1:lengCoef
		@printf "%.5e\t" simulator.normalizationFactor[iteration,i]
	end
	println()
end


text:
module OptMaps

using Random
using Base.Iterators
using ColorSchemes
using ColorTypes
using Images
using Luxor
using CSV
using DataFrames
using MultivariateStats

export OptMap, create_map!, draw_solution!, save_all, save_map, save_base_img
export MAP_SIDE

const MAP_SIDE = 800

mutable struct OptMap
    variable_names::Array{String,1}
    variable_domains::Array{Tuple{R,R},1} where R <: Number
    constraint_names::Array{String,1}
    objective_f::Function
    constraints_f
    sample_points_budget::Int
    verbose::Bool
    min_2d::AbstractArray{Float64,1}
    max_2d::AbstractArray{Float64,1}
    X_nd::AbstractArray{Float64,2}
    C::AbstractArray{Float64,2}
    Z::AbstractArray{Float64,2}
    M::Array{RGB{Normed{UInt8,8}},2}
    S::Array{RGB{Normed{UInt8,8}},2}
    pca #::MultivariateStats.PCA{Float64}
end

function OptMap(variable_names::Array{String,1},
                variable_domains::Array{Tuple{R,R},1} where R <: Number,
                objective_f::Function;
                constraint_names::Array{String,1}=["none"],
                constraints_f=[(x) -> true],
                sample_points_budget=5000000,
                verbose=false)::OptMap

    O = OptMap( variable_names,
                variable_domains,
                constraint_names,
                objective_f,
                constraints_f,
                sample_points_budget,
                verbose,
                zeros(1),
                zeros(1),
                zeros(1,1),
                zeros(1,1),
                zeros(1,1),
                zeros(RGB{Normed{UInt8, 8}}, MAP_SIDE, MAP_SIDE),
                zeros(RGB{Normed{UInt8, 8}}, MAP_SIDE, MAP_SIDE),
                nothing)

    return O
end

function create_map!(O::OptMap, constraint_brightness::Float64=0.7; colormap=:viridis)
    cmap = colorschemes[colormap]

    fit!(O)

    n_constraints = length(O.constraints_f)

    O.verbose && println("Creating OptMap")

    Z = zeros(Float64, MAP_SIDE, MAP_SIDE)
    A = ones(Float64, MAP_SIDE, MAP_SIDE)
    C = ones(Float64, MAP_SIDE^2, n_constraints)

    @inbounds for i in 1:MAP_SIDE
        @inbounds for j in 1:MAP_SIDE
            k = ((i - 1) * MAP_SIDE + j)

            @inbounds for c in 1:n_constraints
                C[k,c] = O.constraints_f[c](O.X_nd[:,k])
            end

            is_feasible_point = *(C[k,:]...) # product of all constraints
            Z[i,j] = O.objective_f(O.X_nd[:,k])
            A[i,j] = Bool(is_feasible_point) ? 1.0 : constraint_brightness
        end
    end

    M = get(cmap, Z, (minimum(Z), maximum(Z)))
    M = RGB.(M .* A)
    M = convert(Array{RGB{N0f8}}, M)

    O.Z = Z
    O.M = M
    O.C = C
end

function draw_solution!(O::OptMap, path_to_solution::AbstractArray{Float64,2}; point_radius::Int=50, line_width::Float64=3.0)
    O.verbose && println("Drawing OptMap")

    path_2d_tmp = MultivariateStats.transform(O.pca, path_to_solution)
    path_2d = zeros(Float64, 2, size(path_to_solution, 2))

    @inbounds for i in 1:size(path_to_solution, 2)
        path_2d[:,i] = coords_to_pixels(path_2d_tmp[:,i], O.min_2d, O.max_2d)
    end

    tmp_map_name = tempname() * ".png"
    save(tmp_map_name, O.M)
    img = readpng(tmp_map_name)

    w = img.width
    h = img.height

    Drawing(w, h, tmp_map_name)
    placeimage(img, 0, 0, 1.0)
    setline(line_width)

    sethue("white")
    circle(path_2d[1,1], path_2d[2,1], point_radius, :stroke)
    sethue("red")
    circle(path_2d[1,end], path_2d[2,end], point_radius, :stroke)

    sethue(1.0, 0.0, 1.0)
    @inbounds for i in 2:size(path_2d, 2)
        p1 = Point((path_2d[:,i-1])...)
        p2 = Point((path_2d[:,i])...)

        circle(p1, line_width*2, :fill)
        circle(p2, line_width*2, :fill)
        line(p1, p2, :stroke)
    end

    finish()

    O.S = load(tmp_map_name)
    rm(tmp_map_name)
end

function save_base_img(O::OptMap, file_name::String)
    O.verbose && println("Saving image to $file_name")
    save(file_name, O.M)
end

function save_map(O::OptMap, file_prefix::String, path::String=".")
    map_name = "$(path)/$(file_prefix)_optmap.png"

    O.verbose && println("Saving map to $map_name")
    save(map_name, O.S)
end

function save_all(O::OptMap, file_prefix::String, path::String=".")
    save_map(O, file_prefix, path)

    X_csv = "$(path)/$(file_prefix)_optmap_X.csv"
    C_csv = "$(path)/$(file_prefix)_optmap_C.csv"
    Z_csv = "$(path)/$(file_prefix)_optmap_Z.csv"

    O.verbose && println("Saving X_nd to $X_csv")
    CSV.write(X_csv, DataFrame(O.X_nd'), header=O.variable_names)

    O.verbose && println("Saving C to $C_csv")
    CSV.write(C_csv, DataFrame(O.C), header=O.constraint_names)

    O.verbose && println("Saving Z to $Z_csv")
    CSV.write(Z_csv, DataFrame(O.Z), writeheader=false)
end

function fit!(O::OptMap)
    O.verbose && println("Sampling nD points")
    X = sample_nd_points(O.variable_domains, O.sample_points_budget)

    O.verbose && println("Projecting nD -> 2D ($(size(X)))")
    X_2d, p = project_2d(X)

    O.verbose && println("Unprojecting 2D point grid")
    X_nd = inverse_project(X_2d, p)

    O.pca = p
    O.X_nd = X_nd

    O.min_2d = [minimum(X_2d[1,:]), minimum(X_2d[2,:])]
    O.max_2d = [maximum(X_2d[1,:]), maximum(X_2d[2,:])]
end

function sample_nd_points(  variables_domains::Array{Tuple{R,R},1} where R <: Number,
                            point_budget::Int)::AbstractArray{Float64,2}
    n_vars = length(variables_domains)

    max_points = convert(Int, floor(point_budget^(1 / n_vars)))

    variable_ranges = []

    @inbounds for vd in variables_domains
        if typeof(vd[1]) == Float64
            points = max_points
        else
            points = vd[2] - vd[1] + 1
            points = min(points, max_points)
        end

        var_range = range(vd[1], vd[2], length=points)
        push!(variable_ranges, var_range)
    end

    tmp = vec(collect(product(variable_ranges...)))
    X = zeros(n_vars, length(tmp))

    @inbounds for j in 1:n_vars
        for i in 1:length(tmp)
            X[j,i] = tmp[i][j]
        end
    end

    return X
end

function project_2d(X::AbstractArray{Float64,2})
    p = fit(PCA, X; maxoutdim=2, method=:svd)
    X_2d = MultivariateStats.transform(p, X)

    return X_2d, p
end

function inverse_project(X_2d::AbstractArray{Float64,2}, p)::AbstractArray{Float64,2}
    tmp_grid = vec(collect(product(
                            range(minimum(X_2d[1,:]), maximum(X_2d[1,:]), length=MAP_SIDE),
                            range(minimum(X_2d[2,:]), maximum(X_2d[2,:]), length=MAP_SIDE))))
    X_2d_grid = zeros(Float64, 2, MAP_SIDE^2)

    @inbounds for i in 1:size(tmp_grid)[1]
        X_2d_grid[1,i] = tmp_grid[i][1]
        X_2d_grid[2,i] = tmp_grid[i][2]
    end

    X_nd = reconstruct(p, X_2d_grid)
    return X_nd
end

function coords_to_pixels(point::AbstractArray{Float64,1}, min_2d::Array{Float64,1}, max_2d::Array{Float64,1}; mingrid=1, maxgrid=MAP_SIDE)
    p = (point - min_2d) ./ (max_2d - min_2d)
    p_trunc = floor.(p .* (maxgrid - mingrid))

    return p_trunc
end

end

max_stars_repo_name: nmaterise/optmap

max_issues_repo_name: nmaterise/optmap

max_forks_repo_name: nmaterise/optmap


text:
"""
    sv = selectElement(v, sy, dw=[1,1,...,1], tri=0)

Select elements of a multi-element DAC to minimize the selection error,
subject to the constraint that the nominal DAC output is `v`, i.e. `v` = (`sv`)ᵀ(`dw`).
Assume that the preferred usage order is that given by `dw`.

Input
- `v`:    DSM output in -M:2:M
- `sy`:   M×1 desired usage vector
- `dw`:   M×1 vector of element weights
- `tri`:  Flag indicating the elements are tri-level, i.e. `sv` = 0 is allowed. NOT SUPPORTED YET
Output
- `sv`:   M×1 selection vector. ±1 if `tri` = 0, {0,±1} if `tri` = 1.
"""
function selectElement(v, sy, dw=[], tri=0)
    M = length(sy)
    if isempty(dw) || all(dw .== 1)
        dw_is_one = true
        sdw = M
    else
        dw_is_one = false
        sdw = sum(dw)
    end

    if abs(v) > sdw
        error("|v| is too large.");
    elseif v == sdw
        return ones(M)
    elseif v == -sdw
        return -ones(M)
    end

    sv = -ones(Int, M)
    possibilities = sortperm(-sy)  # Determine the element priority

    # Speed up selection in the usual case where all element weights are one.
    # Suggested 1997-03-05 by J. A. Cherry.
    if dw_is_one
        sv[possibilities[1:(M+v)÷2]] .= 1
        return sv
    end

    # Go through sv possibilities one by one, until one which meets the
    # v = sv'*dw constraint is found.
    i = 1                    # Selection level
    pointer = ones(Int, M)   # Array of pointers to selected elements
    selected = zeros(Int, M) # Selected elements
    dw0 = [0; dw]
    while true
        while pointer[i] > M
            # backtrack
            i = i - 1
            if i == 0
                break   # failure!
            end
            pointer[i] = pointer[i] + 1;
            selected[i+1:end] .= 0
        end
        selected[i] = possibilities[pointer[i]]
        dv = 2 * sum(dw0[selected.+1]) - sdw
        if dv == v
            break       #success!
        elseif dv < v
            # Proceed to the next level of selection
            pointer[i+1] = pointer[i] + 1
            i = i + 1
        else
            # Try the next element at the current level.
            pointer[i] = pointer[i] + 1
        end
    end
    selected = selected[selected .!= 0]
    sv[selected] .= 1

    return sv
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
module BigCleaner

include("Pipeline.jl")

export Pipeline
export add_step!

end


text:
struct GitBlob <: GitObject
    data::Vector{UInt8}
end

function Base.sizeof(blob::GitBlob)
    length(blob.data)
end

function oid(blob::GitBlob)
    ctx = SHA.SHA1_CTX()
    SHA.update!(ctx, Vector{UInt8}("blob $(sizeof(blob))\0"))
    SHA.update!(ctx, blob.data)
    return SHA1(SHA.digest!(ctx))
end

max_stars_repo_name: simonbyrne/GitX

max_issues_repo_name: simonbyrne/GitX

max_forks_repo_name: simonbyrne/GitX


text:
using MultivariateStats
using Base.Test

srand(34568)

## data

d = 3
n = 5

X = rand(d, n)
mv = vec(mean(X, 2))
C = cov(X, 2)
C0 = copy(C)

emax = maximum(eigvals(C))
rc = 0.01
Cr = C + (emax * rc) * eye(d)

## cov_whitening

cf = cholfact(C, :U)
W = cov_whitening(cf)
@test istriu(W)
@test W'C * W ≈ eye(d)

cf = cholfact(C, :L)
W = cov_whitening(cf)
@test istriu(W)
@test W'C * W ≈ eye(d)

W = cov_whitening(C)
@test C == C0
@test istriu(W)
@test W'C * W ≈ eye(d)

W = cov_whitening(C, 0.0)
@test C == C0
@test istriu(W)
@test W'C * W ≈ eye(d)

W = cov_whitening(C, rc)
@test istriu(W)
@test W'Cr * W ≈ eye(d)

# Whitening from data

f = fit(Whitening, X)
W = f.W
@test isa(f, Whitening{Float64})
@test mean(f) === f.mean
@test istriu(W)
@test W'C * W ≈ eye(d)
@test transform(f, X) ≈ W' * (X .- f.mean)

f = fit(Whitening, X; regcoef=rc)
W = f.W
@test W'Cr * W ≈ eye(d)

f = fit(Whitening, X; mean=mv)
W = f.W
@test W'C * W ≈ eye(d)

f = fit(Whitening, X; mean=0)
Cx = (X * X') / (n - 1)
W = f.W
@test W'Cx * W ≈ eye(d)

# invsqrtm

R = invsqrtm(C)
@test C == C0
@test R ≈ inv(sqrtm(C))


text:
####################################################
# Script to train a patchy PCA/ICA network with 1 hidden layer (SGD)
# patch connectivity is controled via "n_patches" & "n_hidden_per_patch" variables

###################################################
# Include libraries and core

using Pkg; Pkg.activate("./../core/"); Pkg.instantiate()
include("./../core/src/ratenets/core.jl")
include("./../core/src/ratenets/helpers_patchy.jl")

###################################################
# Define task and network

# This keyword decides whether PCA or ICA is used for feature extraction:
featurefunction = getPCAscores
# OR: getICAscores

rectifyfeatures = true # Otherwise patchy PCA/ICA is not more powerful than global PCA/ICA (since vanilla PCA/ICA is linear)
dataset = "MNIST" # dataset to be used: "MNIST" or "CIFAR10"
colorimage = dataset == "CIFAR10"
patch_size = 10 # (linear) patch size of receptive fields in first layer p ∈ [1,28] or [1,32] for MNIST/CIFAR10 resp.
n_patches = 500 # number of feature extraction patches
n_hidden_per_patch = 10 # hidden neurons/features per patch. Overall number of hidden neurons is: n_patches * n_hidden_per_patch
n_inits = 3 # number of initializations per net (for averaging)
iterations = 10^6 # iterations per initialization
learningrates = 1e-3*[1; 1] # learningrates per layer
nonlinearity = [relu!] # nonlinearities
nonlinearity_diff = [relu_diff!]

###################################################
# Load and preprocess data

smallimgs, labels, smallimgstest, labelstest = import_data(dataset)
subtractmean!(smallimgs) # substract pixel-wise mean
subtractmean!(smallimgstest)

# Extract PCA/ICA features per patch
smallimgs, smallimgstest = getpatchyhiddenreps(smallimgs, smallimgstest, n_patches, n_hidden_per_patch;
    featurefunction = featurefunction, patchsize = patch_size, colorimage = colorimage)

# Rectify features if applicable
rectifyfeatures && for imgs in [smallimgs,smallimgstest]
     clamp!(imgs,0,Inf64)
end

##############################################################################
# Train network

print(string("train classifier on (", n_patches * n_hidden_per_patch,") unsupervised features (without hidden layer) ...\n"))
for j in 1:n_inits
  net = Network([size(smallimgs)[1], 10])
  learn!(net, getsmallimg, getlabel, iterations, learningrates, 0.,
                     nonlinearity = nonlinearity, nonlinearity_diff = nonlinearity_diff)
  error_train = geterrors(net, smallimgs, labels, nonlinearity = nonlinearity)
  error_test = geterrors(net, smallimgstest, labelstest, nonlinearity = nonlinearity)
  print(string(100 * error_train, " % on training set\n"))
  print(string(100 * error_test, " % on test set\n"))
end


text:
module Problox

using PyCall
using Conda
using MacroTools
using MacroTools: @capture

# Install ProbLog.
if PyCall.conda
	Conda.add("pip")
	pip = joinpath(Conda.BINDIR, "pip")
	run(`$pip install problog`)
else
	try
        pyimport("problog")
        pyimport("problog.program")
        pyimport("problog.logic")
        pyimport("problog.tasks")
        pyimport("problog.tasks.sample")
	catch ee
		typeof(ee) <: PyCall.PyError || rethrow(ee)
		warn("""
Python Dependencies not installed
Please either:
 - Rebuild PyCall to use Conda, by running in the julia REPL:
    - `ENV[PYTHON]=""; Pkg.build("PyCall"); Pkg.build("Problox")`
 - Or install the depencences, eg by running pip
	- `pip install problog`
	"""
		)
	end
end

problog = pyimport("problog")
pp = pyimport("problog.program")
pl = pyimport("problog.logic")
pt = pyimport("problog.tasks")
sample = pyimport("problog.tasks.sample")

# User-defined terms.
function term(str::String)
    return pl.Term(str)
end

SimpleProgram() = pp.SimpleProgram()
Var(str::String) =  pl.Var(str)
variable(s::String) = Var(String(s))
variable(s::Symbol) = Var(String(s))
Constant(str::String) =  pl.Constant(str)
AnnotatedDisjunction(p::Array{PyObject}, q::PyObject) = pl.AnnotatedDisjunction(p, q)
AnnotatedDisjunction(p::NTuple{N, PyObject}, q::PyObject) where N = pl.AnnotatedDisjunction([p...], q)

evaluate(p::PyObject) = problog.get_evaluatable().create_from(p).evaluate()

include("dsl.jl")

# DSL.
export @problox

# Lower-level API.
export SimpleProgram, Var, Constant, AnnotatedDisjunction, term

# Convenience.
query = term("query")
export query, evaluate, variable

end # module


text:
import Base.show

using Colors
using CSV
using LinearAlgebra
using Plots

"""
    getpoints(path::String, [usingnames::Bool])

Read the x and y coordinates from a list in a CSV file at `path`.

Also read the names if `usingnames`.

# Example
```julia-repl
julia> getpoints("points.txt")
```
"""
function getpoints(path::String, usingnames::Bool = false)
    xs = Int64[]
    ys = Int64[]
    names = String[]
    for row ∈ CSV.Rows(path)
        push!(xs,parse(Int64,row.x))
        push!(ys,parse(Int64,row.y))
        if usingnames
            push!(names,row.name)
        end
    end
    return usingnames ? (names, xs, ys) : (xs, ys)
end

"""
    bruteforcevoronoi(xs::Vector{Int64}, ys::Vector{Int64})

Draw a Voronoi diagram from lists of coordinates `xs` and `ys` by individually colouring every square of a 1000x1000 grid.

# Example
```julia-repl
julia> bruteforcevoronoi(rand(1:1000, 10), rand(1:1000, 10))
```
"""
function bruteforcevoronoi(xs::Vector{Int64}, ys::Vector{Int64})
    m = length(xs)
    colours = distinguishable_colors(m, lchoices = 50:5:100)

    pointcolours = fill(RGB(0,0,0), (1000,1000))
    for y ∈ 1:1000
        for x ∈ 1:1000
            pointcolours[x,y] = colours[findmin(hypot.(xs .- x, ys .- y))[2]]
        end
    end

    diagram = scatter(
        [(x,y) for x ∈ 1:1000, y ∈ 1:1000][:];
        c = pointcolours[:],
        markersize = 1,
        markerstrokewidth = 0,
        legend = false,
        showaxis = false,
        ticks = false,
        size = (1000, 1000)
    )
    return scatter!(
        diagram,
        xs,
        ys;
        c = :black,
        markersize = 2,
        markerstrokewidth = 0
    )
end

"""
    montecarlovoronoi(n::Int64, xs::Vector{Int64}, ys::Vector{Int64})

Sample `n` times from a 1000x1000 grid to approximate a Voronoi diagram from lists of coordinates `xs` and `ys`.

# Example
```julia-repl
julia> montecarlovoronoi(100000, rand(1:1000, 10), rand(1:1000, 10))
```
"""
function montecarlovoronoi(n::Int64, xs::Vector{Int64}, ys::Vector{Int64})
    m = length(xs)
    colours = distinguishable_colors(m, lchoices = 50:5:100)

    randxs = Int64[]
    randys = Int64[]
    randcolours = RGB[]

    for _ ∈ 1:n
        x, y = rand(1:1000), rand(1:1000)
        push!(randxs, x)
        push!(randys, y)
        push!(randcolours, colours[findmin(hypot.(xs .- x, ys .- y))[2]])
    end

    diagram = scatter(
        randxs,
        randys;
        c = randcolours,
        markersize = 1,
        markerstrokewidth = 0,
        legend = false,
        showaxis = false,
        ticks = false,
        size = (1000, 1000)
    )
    return scatter!(
        diagram,
        xs,
        ys;
        c = :black,
        markersize = 2,
        markerstrokewidth = 0,
    )
end

"""
    VoronoiPoint(x::Real, y::Real, neighbours::Vector{VoronoiPoint})
    VoronoiPoint(x::Real, y::Real)

Construct a `VoronoiPoint` at the coordinates `x`, `y`, representing a point in a Voronoi diagram.

If no neighbours are specified, default to an empty `Vector{VoronoiPoint}`.

# Examples
```julia-repl
julia> p = VoronoiPoint(3,4)
(3,4)

julia> q = VoronoiPoint(4,5,[p])
(4,5)

julia> p.neighbours
VoronoiPoint[]

julia> q.neighbours
1-element Vector{VoronoiPoint}:
 (3,4)
```
"""
mutable struct VoronoiPoint
    x::Real
    y::Real
    neighbours::Vector{VoronoiPoint}
end
VoronoiPoint(x::T, y::T) where T <: Real = VoronoiPoint(x, y, VoronoiPoint[])

show(io::IO, p::VoronoiPoint) = print(io, "(", p.x, ",", p.y, ")")

"""
    pos(p::VoronoiPoint)

Return the `x` and `y` coordinates of the `VoronoiPoint` `p` as a `Vector`.

# Example
```julia-repl
julia> pos(VoronoiPoint(3,4))
2-element Vector{Int64}:
 3
 4
```
"""
pos(p::VoronoiPoint) = [p.x, p.y]

"""
    perpbisector(p::VoronoiPoint, q::VoronoiPoint)

Return the midpoint of `VoronoiPoint`s `p` and `q`, and direction vector along the perpendicular bisector between of the same length as `p - q` and direction such that `p` is to the right and `q` is to the left.

# Example
```julia-repl
julia> perpbisector(VoronoiPoint(0,1), VoronoiPoint(2,1))
([1.0, 1.0], [0, -2])
```
"""
function perpbisector(p::VoronoiPoint, q::VoronoiPoint)
    return ( (pos(p) + pos(q))/2 , [[0,1] [-1,0]] * (pos(p) - pos(q)) )
end

"""
    circumcentre(p::VoronoiPoint, q::VoronoiPoint, r::VoronoiPoint)

Return the circumcentre of the three `VoronoiPoint`s `p`, `q`, and `r`.

If they are collinear, return `[Inf, Inf]`.

# Example
```julia-repl
julia> circumcentre(VoronoiPoint(0,0), VoronoiPoint(1,0), VoronoiPoint(0,1))
2-element Vector{Float64}:
 0.5
 0.5
```
"""
function circumcentre(p::VoronoiPoint, q::VoronoiPoint, r::VoronoiPoint)
    (x₁, v₁) = perpbisector(p,q)
    (x₂, v₂) = perpbisector(p,r)
    det([v₁ v₂]) == 0 && return [Inf, Inf]
    M = [v₁ v₂]^-1
    t₁ = (M * (x₂ - x₁))[1]
    return x₁ + t₁*v₁
end

"""
    findnextpoint(p::VoronoiPoint, q::VoronoiPoint, z::Vector{T}, reversedirection::Bool) where T <: Real

Find the `r` from `q.neighbours` where the circumcentre of `p`, `q`, and `r` is closest to `z` in the direction determined by `perpbisector`.

If `reversedirection == true`, find `r` based on the opposite direction.

# Examples
```julia-repl
julia> p = VoronoiPoint(0,0)
(0,0)

julia> q = VoronoiPoint(0,1,[p, VoronoiPoint(3,4), VoronoiPoint(-3,4)])
(0,1)

julia> findnextpoint(p, q, [0, 0.5], false)
(3,4)

julia> findnextpoint(p, q, [0, 0.5], true)
(-3,4)
```
"""
function findnextpoint(p::VoronoiPoint, q::VoronoiPoint, z::Vector{T}, reversedirection::Bool) where T <: Real
    directionvector = perpbisector(p,q)[2] .* (-1)^reversedirection

    mindisplacement = Inf
    nextpoint = VoronoiPoint(Inf,Inf)
    for r ∈ q.neighbours
        r ∈ [p,q] && continue
        displacement = dot(directionvector, circumcentre(p,q,r) - z)
        displacement > 0 && displacement < mindisplacement && (mindisplacement = displacement; nextpoint = r)
    end
    return nextpoint
end

"""
    addpoint(x::Real, y::Real, allpoints::Vector{VoronoiPoint})

Add a new `VoronoiPoint` at coordinates `x` and `y`, including calculating all neighbours from `allpoints`.

Also, update neighbour lists of other points in `allpoints` to agree with the new addition.

# Example
```julia-repl
julia> p = VoronoiPoint(1,1)
(1,1)

julia> q = VoronoiPoint(2,2)
(2,2)

julia> p.neighbours = [VoronoiPoint(Inf,Inf), q, VoronoiPoint(Inf,Inf)]
3-element Vector{VoronoiPoint}:
 (Inf,Inf)
 (2,2)
 (Inf,Inf)

julia> q.neighbours = [VoronoiPoint(Inf,Inf), p, VoronoiPoint(Inf,Inf)]
3-element Vector{VoronoiPoint}:
 (Inf,Inf)
 (1,1)
 (Inf,Inf)

julia> r = addpoint(1, 2, [p,q])
(1,2)

julia> p.neighbours
4-element Vector{VoronoiPoint}:
 (Inf,Inf)
 (2,2)
 (1,2)
 (Inf,Inf)

julia> q.neighbours
4-element Vector{VoronoiPoint}:
 (Inf,Inf)
 (1,2)
 (1,1)
 (Inf,Inf)

julia> r.neighbours
4-element Vector{VoronoiPoint}:
 (Inf,Inf)
 (1,1)
 (2,2)
 (Inf,Inf)
```
"""
function addpoint(x::Real, y::Real, allpoints::Vector{VoronoiPoint})
    p = VoronoiPoint(x, y)
    nearestpoint = allpoints[findmin([hypot(x-q.x, y-q.y) for q ∈ allpoints])[2]]

    z = perpbisector(p, nearestpoint)[1]
    q = nearestpoint
    reversedirection = false
    while true
        push!(p.neighbours, q)
        r = findnextpoint(p, q, z, reversedirection)

        # Edge goes to infinity
        if pos(r) == [Inf,Inf]
            push!(p.neighbours, r)
            if reversedirection
                break
            else
                reversedirection = true
                q = nearestpoint
                z = perpbisector(p, nearestpoint)[1]
                pop!(reverse!(p.neighbours))
                continue
            end

        # Back to start
        elseif pos(r) == pos(nearestpoint)
            reverse!(p.neighbours)
            break

        else
            z = circumcentre(p,q,r)
            q = r
        end
    end

    n = length(p.neighbours)
    for i ∈ 1:n
        q = p.neighbours[i]
        pos(q) == [Inf, Inf] && continue
        j = findlast(r -> pos(r) == pos(p.neighbours[mod(i-1,1:n)]), q.neighbours)
        k = findfirst(r -> pos(r) == pos(p.neighbours[mod(i+1,1:n)]), q.neighbours)
        if j < k
            q.neighbours = vcat(q.neighbours[j:k], p)
        else
            q.neighbours = vcat(q.neighbours[1:k], p, q.neighbours[j:end])
        end
    end

    return p
end

"""
    voronoipoints(xs::Vector{T}, ys::Vector{T}) where T <: Real

Calculate the lists of neighbouring regions of a Voronoi diagram with points at coordinates given by `xs` and `ys`.

Return a list of `VoronoiPoint`s containing this information, with `Rational` parameters.

# Example
```julia-repl
voronoipoints([1,2,3],[2,1,3])
julia> points = voronoipoints([1,2,3],[2,1,3])
3-element Vector{VoronoiPoint}:
 (1//1,2//1)
 (2//1,1//1)
 (3//1,3//1)

julia> points[1].neighbours
4-element Vector{VoronoiPoint}:
 (Inf,Inf)
 (2//1,1//1)
 (3//1,3//1)
 (Inf,Inf)

julia> points[2].neighbours
4-element Vector{VoronoiPoint}:
 (Inf,Inf)
 (3//1,3//1)
 (1//1,2//1)
 (Inf,Inf)

julia> points[3].neighbours
4-element Vector{VoronoiPoint}:
 (Inf,Inf)
 (1//1,2//1)
 (2//1,1//1)
 (Inf,Inf)
```
"""
function voronoipoints(xs::Vector{T}, ys::Vector{T}) where T <: Real
    # Rational numbers!
    xs, ys = Rational{BigInt}.(xs), Rational{BigInt}.(ys)
    points = [VoronoiPoint(xs[1],ys[1],[VoronoiPoint(Inf,Inf)])]
    for i ∈ 2:length(xs)
        push!(points,addpoint(xs[i],ys[i],points))
    end
    return points
end

"""
    voronoi(xs::Vector{T}, ys::Vector{T}; [names::Vector{String}]) where T <: Real

Draw a Voronoi diagram from points at coordinates `xs` and `ys` in the region (x,y) ∈ [0,1000]².

Points given `names` in a legend if specified.

# Example
```julia-repl
julia> voronoi(rand(1:1000, 10), rand(1:1000, 10))
```
"""
function voronoi(xs::Vector{T}, ys::Vector{T}; names::Vector{String} = String[]) where T <: Real
    points = voronoipoints(xs,ys)
    m = length(xs)
    colours = distinguishable_colors(m, lchoices = 50:5:100)
    usingnames = !isempty(names)

    polygons = Shape[]
    for p ∈ points
        n = length(p.neighbours)
        vertices = Tuple{Rational{Int64}, Rational{Int64}}[]

        if pos(p.neighbours[1]) == [Inf,Inf]
            startindex, endindex = 2, n-2

            x₁, v₁ = perpbisector(p.neighbours[n-1],p)
            x₁, v₁ = Rational{Int64}.(x₁), Rational{Int64}.(v₁)
            # Values of t₁ where the line x₁ + t₁v₁ intersects the borders
            t₁vals = [-x₁[2]/v₁[2], (1000-x₁[1])/v₁[1], (1000-x₁[2])/v₁[2], -x₁[1]/v₁[1]]
            t₁ = minimum(t₁vals[t₁vals .> 0])
            push!(vertices,Tuple(x₁ + t₁*v₁))

            x₂, v₂ = perpbisector(p,p.neighbours[2])
            x₂, v₂ = Rational{Int64}.(x₂), Rational{Int64}.(v₂)
            t₂vals = [-x₂[2]/v₂[2], (1000-x₂[1])/v₂[1], (1000-x₂[2])/v₂[2], -x₂[1]/v₂[1]]
            t₂ = minimum(t₂vals[t₂vals .> 0])

            t₁index = findfirst(==(t₁),t₁vals)
            t₂index = findfirst(==(t₂),t₂vals)
            indexdifference = mod(t₂index - t₁index,4)
            corners = [(0,0), (1000,0), (1000,1000), (0,1000)]
            for j ∈ 1:indexdifference
                push!(vertices,corners[mod(t₁index+j, 1:4)])
            end

            push!(vertices,Tuple(x₂ + t₂*v₂))

        else
            startindex, endindex = 1, n
        end

        for j ∈ startindex:endindex
            q,r = p.neighbours[j], p.neighbours[mod(j+1,1:n)]
            push!(vertices,Tuple(circumcentre(p,q,r)))
        end

        push!(polygons, Shape(vertices))
    end

    diagram = plot(
        polygons,
        c = colours',
        label = (usingnames ? hcat(names...) : false),
        showaxis = false,
        ticks = false,
        size = (1000, 1000),
        xlims = (0, 1000),
        ylims = (0, 1000)
    )
    return scatter!(diagram, xs, ys, c = :black, markersize = 2, label = false)
end

"""
    voronoidistribution(points::Vector{VoronoiPoint})

Count the number of edges of each of the polygons from a Voronoi diagram, as represented by `points`.

`points` should usually be the output of the function `voronoipoints`.

Display the counts in a bar chart.

# Example
```julia-repl
voronoidistribution(voronoipoints(rand(1:1000, 10), rand(1:1000, 10)))
```
"""
function voronoidistribution(points::Vector{VoronoiPoint})
    neighbourcount = broadcast(p -> pos(p.neighbours[1]) == [Inf, Inf] ? 0 : length(p.neighbours), points)

    m = maximum(neighbourcount)
    total = count(!=(0), neighbourcount)
    return bar(
        1:m,
        [count(==(i), neighbourcount)//total for i ∈ 1:m],
        legend = false,
        xticks = 1:m,
        yticks = 0:0.1:1
    )
end

"""
    randomlattice(shape::Symbol, n::Int64, σ::Rational{Int64})

Create a lattice of `n^2` points, with each perturbed randomly both vertically and horizontally by a normal distribution of standard deviation `σ`.

If `shape == :square`, the lattice is a square lattice.

If `shape == :hexagonal`, the lattice is an approximately equilateral hexagonal lattice.

Return a bar chart of the distribution of edges of the polygons in the Voronoi diagram generated by these points.

# Examples
```julia-repl
julia> randomlattice(:square, 10, 1)

julia> randomlattice(:hexagonal, 20, 1//10)
```
"""
function randomlattice(shape::Symbol, n::Int64, σ::Rational{Int64})
	if shape == :square
		coords = [(x,y) for x ∈ 1:n, y ∈ 1:n]
	elseif shape == :hexagonal
		coords = hcat([(x,(6//7)*y) for x ∈ 1:n, y ∈ 1:2:n],
			[(x+1//2,(6//7)*y) for x ∈ 1:n, y ∈ 2:2:n])
	else
		error("Invalid shape")
	end

	xs,ys = [xy[1] for xy ∈ coords[:]], [xy[2] for xy ∈ coords[:]]
	randxs = xs .+ σ * round.(Int64, n * randn(n^2))//n
	randys = ys .+ σ * round.(Int64, n * randn(n^2))//n
	return voronoidistribution(voronoipoints(randxs,randys))
end


text:
using FdeSolver
using Test
using Statistics

@testset "FdeSolver.jl" begin

tSpan = [0, 10]     # [intial time, final time]
x0 = [1 1]             # intial value
β = 2            # order of the derivative
par = [16.0, 4.0] # [spring constant for a mass on a spring, inertial mass]
h = 0.01

function F(t, x, par)

      K = par[1]
      m = par[2]

      - K ./ m .* x

end

t, Yapp = FDEsolver(F, tSpan, x0, β, par, h = h)

a = x0[1] .* map(cos, sqrt(par[1] / par[2]) .* t) .+ x0[2] ./ sqrt(par[1] / par[2]) .* map(sin, sqrt(par[1] / par[2]) .* t)

@test @isdefined(t)
@test @isdefined(Yapp)
@test abs((mean(Yapp) - mean(a))) < 0.05

end


text:
"""
**SeisRemove**

*Delete a seis file (deletes the text file, binary data file, and binary header file if there is one)*

**IN**
* in: input filename

**OUT**

*Credits: AS, 2015*

"""
function SeisRemove(filename::String)

	filename_data = ParseDataName(filename)
	filename_headers = ParseHeaderName(filename)
	rm(filename);
	rm(filename_data);
	if filename_headers != "NULL"
		rm(filename_headers);
	end

end

function SeisRemove(filename::Array{String,1})

	for j = 1 : length(filename)
		SeisRemove(filename[j])
	end

end


text:
module Sherlogs

    export Sherlog64,Sherlog32,Sherlog16,
            DrWatson64,DrWatson32,DrWatson16,
            LogBook,reset_logbook,reset_logbooks,
            get_logbook,get_logbooks,
            get_stacktrace,get_stacktraces, get_st_index,
            entropy

    import Base: +, -, *, /, ^, promote_rule, round, rem
    import StatsBase: entropy
    import Random
    import Random: rand, randn, randexp

    include("logbook.jl")
    include("SherlogN.jl")
    include("DrWatsonN.jl")
    include("stacktraces.jl")

end


text:
function guess_infinity(diagram, infinity)
    if isnothing(infinity)
        if hasproperty(diagram, :threshold)
            thresh = threshold(diagram)
        else
            thresh = nothing
        end
        if isnothing(thresh) || ismissing(thresh) || isinf(thresh)
            return 1.1 * maximum(Iterators.filter(isfinite, Iterators.flatten(diagram)))
        else
            return thresh
        end
    else
        return infinity
    end
end

function AbstractPlotting.convert_arguments(
    ::Type{<:Scatter}, intervals::AbstractVector{<:PersistenceInterval}
)
    return ([Point2f0(birth(int), death(int)) for int in intervals],)
end

@recipe(DiagramBackground) do scene
    return Theme(; persistence=false, gapwidth=0.1)
end

function AbstractPlotting.plot!(p::DiagramBackground)
    t_lo = p[1]
    t_hi = p[2]
    inf_val = p[3]
    persistence = p[:persistence]
    gap = @lift ($t_hi - $t_lo) * $(p[:gapwidth])

    persistence_line = lift(t_lo, t_hi, gap, persistence) do t_lo, t_hi, gap, persistence
        if persistence
            return [Point2f0(t_lo - gap, 0), Point2f0(t_hi + gap, 0)]
        else
            return [Point2f0(t_lo - gap, t_lo - gap), Point2f0(t_hi + gap, t_hi + gap)]
        end
    end
    lines!(p, persistence_line)

    inf_line = @lift [Point2f0($t_lo - $gap, $inf_val), Point2f0($t_hi + $gap, $inf_val)]
    return lines!(p, inf_line; linestyle=:dot, color=:gray)
end

@recipe(DiagramPlot, diagram) do scene
    return Theme(; palette=DEFAULT_PALETTE, color=1, infinity=nothing, persistence=false)
end

function AbstractPlotting.plot!(p::DiagramPlot)
    diag = p[1]
    infinity = @lift guess_infinity($diag, $(p[:infinity]))
    points = lift(diag, infinity, p[:persistence]) do diag, infinity, persistence
        pts = convert_arguments(Scatter, diag)[1]
        for i in eachindex(pts)
            if persistence
                pts[i] -= Point2f0(0, pts[i][1])
            end
            if isinf(pts[i][2])
                pts[i] = Point2f0(pts[i][1], infinity)
            end
        end
        pts
    end
    return scatter!(p, points; color=get_color(p, p[:color]))
end

function AbstractPlotting.default_theme(
    scene::SceneLike, ::Type{<:Plot(PersistenceDiagram)}
)
    return Theme(; color=:gray, infinity=nothing, persistence=false)
end

function _default_colors(diags)
    if all(d -> hasproperty(d, :dim), diags) && allunique(dim.(diags))
        return dim.(diags) .+ 1
    else
        return 1:length(diags)
    end
end

function plot_diagram!(
    scene,
    diags;
    infinity=nothing,
    persistence=false,
    palette=DEFAULT_PALETTE,
    time=Observable(nothing),
    gapwidth=0.1,
    colors=_default_colors(diags),
)
    if !(time isa Observable)
        time = Observable(time)
    end
    t_lo, t_hi, inf_val = PersistenceDiagrams.limits(diags, infinity)
    width = t_hi - t_lo
    gap = width * gapwidth

    diagrambackground!(scene, t_lo, t_hi, inf_val; persistence, gapwidth)

    for (i, diag) in enumerate(diags)
        !isempty(diag) && diagramplot!(
            scene, diag; infinity=infinity, persistence=persistence, color=colors[i]
        )
    end
    xlims!(scene, t_lo - gap, t_hi + gap)
    ylims!(scene, t_lo - gap, t_hi + gap)
    if scene isa LAxis
        scene.xlabel = "birth"
        scene.ylabel = persistence ? "persistence" : "death"
    else
        xlabel!(scene, "birth")
        ylabel!(scene, persistence ? "persistence" : "death")
    end

    time_line = lift(time) do t
        if isnothing(t)
            [Point2f0(0, 0)]
        else
            [Point2f0(t_lo - gap, t), Point2f0(t, t), Point2f0(t, t_hi + gap)]
        end
    end
    lines!(scene, time_line; linestyle=:dash)

    return scene
end
function plot_diagram!(scene, diag::PersistenceDiagram; kwargs...)
    return plot_diagram!(scene, (diag,); kwargs...)
end
plot_diagram(diags; kwargs...) = plot_diagram!(Scene(), diags; kwargs...)

for T in (
    AbstractVector{<:PersistenceDiagram},
    NTuple{<:Any,PersistenceDiagram},
    PersistenceDiagram,
)
    @eval AbstractPlotting.plot(diags::$T; kwargs...) = plot_diagram(diags; kwargs...)
end

@recipe(Bars) do scene
    Theme(; color=:black, ystart=1, linewidth=3, infinity=nothing)
end
function AbstractPlotting.plot!(p::Bars)
    infinity = lift(guess_infinity, p[1], p[:infinity])
    bar_points = lift(p[1], p[:ystart], infinity) do points, ystart, infinity
        res = Point2f0[]
        for (i, p) in enumerate(points)
            y = i - 1 + ystart
            x1, x2 = p
            x2 = isinf(x2) ? infinity : x2
            append!(res, (Point2f0(x1, y), Point2f0(x2, y)))
        end
        if isempty(res)
            [Point2f0(0, 0)]
        else
            res
        end
    end
    return linesegments!(p, bar_points; color=p[:color], linewidth=p[:linewidth])
end

function plot_barcode!(
    scene,
    diags;
    infinity=Observable(nothing),
    palette=DEFAULT_PALETTE,
    linewidth=3,
    time=Observable(nothing),
    colors=_default_colors(diags),
)
    cscheme = PlotUtils.get_colorscheme(palette)
    lims = @lift PersistenceDiagrams.limits(diags, $infinity)
    t_lo, t_hi, inf_val = to_value(lims)
    infinity = @lift $lims[3]
    width = t_hi - t_lo
    gap = width * 0.1
    n_bars = sum(length, diags)
    ygap = n_bars * 0.1

    # Inf line
    lines!(
        scene,
        [Point2f0(inf_val, 1 - ygap), Point2f0(inf_val, n_bars + ygap)];
        linestyle=:dot,
        color=:gray,
    )

    ystart = 1
    for (i, diag) in enumerate(diags)
        if !isempty(diag)
            color = cscheme[colors[i]]
            bars!(
                scene,
                diag;
                linewidth=linewidth,
                color=color,
                ystart=ystart,
                infinity=infinity,
            )
            ystart += length(diag)
        end
    end
    xlims!(scene, t_lo, t_hi)
    ylims!(scene, 0, n_bars)
    if scene isa LAxis
        scene.xlabel = "t"
        scene.ylabel = "i"
    else
        xlabel!(scene, "t")
        ylabel!(scene, "i")
    end

    time_line = lift(time) do t
        if isnothing(t)
            [Point2f0(0, 0)]
        else
            [Point2f0(t, 1 - ygap), Point2f0(t, n_bars + ygap)]
        end
    end
    lines!(scene, time_line; linestyle=:dash)

    return scene
end
function plot_barcode!(scene, diag::PersistenceDiagram; kwargs...)
    return plot_barcode!(scene, (diag,); kwargs...)
end
plot_barcode(diags; kwargs...) = plot_barcode!(Scene(), diags; kwargs...)


text:
const c = Center()
const f = Face()

function build_condition(Topo, side, dim)
    if Topo == :Bounded
        return :(($side < 1) | ($side > grid.$dim))
    elseif Topo == :LeftConnected
        return :(($side > grid.$dim))
    else
        return :(($side < 1))
    end
end

Topos = [:Bounded, :LeftConnected, :RightConnected]

#####
##### Exterior node and peripheral node
#####

"""
    inactive_cell(i, j, k, grid)

Return `true` when the tracer cell at `i, j, k` is lies outside the "active domain" of
the grid in `Bounded` directions. Otherwise, return `false`.
"""
@inline inactive_cell(i, j, k, grid) = false

for Topo in Topos

    xcondition = build_condition(Topo, :i, :Nx)
    ycondition = build_condition(Topo, :j, :Ny)
    zcondition = build_condition(Topo, :k, :Nz)

    @eval begin
        XBoundedGrid = AbstractGrid{<:Any, <:$Topo}
        YBoundedGrid = AbstractGrid{<:Any, <:Any, <:$Topo}
        ZBoundedGrid = AbstractGrid{<:Any, <:Any, <:Any, <:$Topo}

        @inline inactive_cell(i, j, k, grid::XBoundedGrid) = ifelse($xcondition, true, false)
        @inline inactive_cell(i, j, k, grid::YBoundedGrid) = ifelse($ycondition, true, false)
        @inline inactive_cell(i, j, k, grid::ZBoundedGrid) = ifelse($zcondition, true, false)
    end
    for OtherTopo in Topos
        xycondition = :( $xcondition | $(build_condition(OtherTopo, :j, :Ny)))
        xzcondition = :( $xcondition | $(build_condition(OtherTopo, :k, :Nz)))
        yzcondition = :( $ycondition | $(build_condition(OtherTopo, :k, :Nz)))

        @eval begin
            XYBoundedGrid = AbstractGrid{<:Any, <:$Topo, <:$OtherTopo}
            XZBoundedGrid = AbstractGrid{<:Any, <:$Topo, <:Any, <:$OtherTopo}
            YZBoundedGrid = AbstractGrid{<:Any, <:Any, <:$Topo, <:$OtherTopo}

            @inline inactive_cell(i, j, k, grid::XYBoundedGrid) = ifelse($xycondition, true, false)
            @inline inactive_cell(i, j, k, grid::XZBoundedGrid) = ifelse($xzcondition, true, false)
            @inline inactive_cell(i, j, k, grid::YZBoundedGrid) = ifelse($yzcondition, true, false)
        end
        for LastTopo in Topos
            xyzcondition = :( $xycondition | $(build_condition(LastTopo, :k, :Nz)))

            @eval begin
                XYZBoundedGrid = AbstractGrid{<:Any, <:$Topo, <:$OtherTopo, <:$LastTopo}

                @inline inactive_cell(i, j, k, grid::XYZBoundedGrid) = ifelse($xyzcondition, true, false)
            end
        end
    end
end

"""
    inactive_node(LX, LY, LZ, i, j, k, grid)

Return `true` when the location `(LX, LY, LZ)` is "inactive" and thus not directly
associated with an "active" cell.

For `Face` locations, this means the node is surrounded by `inactive_cell`s:
the interfaces of "active" cells are _not_ `inactive_node`.

For `Center` locations, this means the direction is `Bounded` and that the
cell or interface centered on the location is completely outside the active
region of the grid.
"""
@inline inactive_node(LX, LY, LZ, i, j, k, grid) = inactive_cell(i, j, k, grid)

@inline inactive_node(::Face, LY, LZ, i, j, k, grid) = inactive_cell(i, j, k, grid) & inactive_cell(i-1, j, k, grid)
@inline inactive_node(LX, ::Face, LZ, i, j, k, grid) = inactive_cell(i, j, k, grid) & inactive_cell(i, j-1, k, grid)
@inline inactive_node(LX, LY, ::Face, i, j, k, grid) = inactive_cell(i, j, k, grid) & inactive_cell(i, j, k-1, grid)

@inline inactive_node(::Face, ::Face, LZ, i, j, k, grid) = inactive_node(c, f, c, i, j, k, grid) & inactive_node(c, f, c, i-1, j, k, grid)
@inline inactive_node(::Face, LY, ::Face, i, j, k, grid) = inactive_node(c, c, f, i, j, k, grid) & inactive_node(c, c, f, i-1, j, k, grid)
@inline inactive_node(LX, ::Face, ::Face, i, j, k, grid) = inactive_node(c, f, c, i, j, k, grid) & inactive_node(c, f, c, i, j, k-1, grid)

@inline inactive_node(::Face, ::Face, ::Face, i, j, k, grid) = inactive_node(c, f, f, i, j, k, grid) & inactive_node(c, f, f, i-1, j, k, grid)

"""
    peripheral_node(LX, LY, LZ, i, j, k, grid)

Return `true` when the location `(LX, LY, LZ)`, is _either_ inactive or
lies on the boundary between inactive and active cells in a `Bounded` direction.
"""
@inline peripheral_node(LX, LY, LZ, i, j, k, grid) = inactive_cell(i, j, k, grid)

@inline peripheral_node(::Face, LY, LZ, i, j, k, grid) = inactive_cell(i, j, k, grid) | inactive_cell(i-1, j, k, grid)
@inline peripheral_node(LX, ::Face, LZ, i, j, k, grid) = inactive_cell(i, j, k, grid) | inactive_cell(i, j-1, k, grid)
@inline peripheral_node(LX, LY, ::Face, i, j, k, grid) = inactive_cell(i, j, k, grid) | inactive_cell(i, j, k-1, grid)

@inline peripheral_node(::Face, ::Face, LZ, i, j, k, grid) = peripheral_node(c, f, c, i, j, k, grid) | peripheral_node(c, f, c, i-1, j, k, grid)
@inline peripheral_node(::Face, LY, ::Face, i, j, k, grid) = peripheral_node(c, c, f, i, j, k, grid) | peripheral_node(c, c, f, i-1, j, k, grid)
@inline peripheral_node(LX, ::Face, ::Face, i, j, k, grid) = peripheral_node(c, f, c, i, j, k, grid) | peripheral_node(c, f, c, i, j, k-1, grid)

@inline peripheral_node(::Face, ::Face, ::Face, i, j, k, grid) = peripheral_node(c, f, f, i, j, k, grid) | peripheral_node(c, f, f, i-1, j, k, grid)

"""
    boundary_node(LX, LY, LZ, i, j, k, grid)

Return `true` when the location `(LX, LY, LZ)` lies on a boundary.
"""
@inline boundary_node(LX, LY, LZ, i, j, k, grid) = peripheral_node(LX, LY, LZ, i, j, k, grid) & !inactive_node(LX, LY, LZ, i, j, k, grid)


text:
using DashBootstrapComponents

card_content = [
    dbc_cardheader("Card header"),
    dbc_cardbody([
        html_h5("Card title", className = "card-title"),
        html_p("This is some card content that we'll reuse", className = "card-text"),
    ]),
];

row_1 = dbc_row(
    [
        dbc_col(dbc_card(card_content, color = "primary", outline = true)),
        dbc_col(dbc_card(card_content, color = "secondary", outline = true)),
        dbc_col(dbc_card(card_content, color = "info", outline = true)),
    ],
    className = "mb-4",
);

row_2 = dbc_row(
    [
        dbc_col(dbc_card(card_content, color = "success", outline = true)),
        dbc_col(dbc_card(card_content, color = "warning", outline = true)),
        dbc_col(dbc_card(card_content, color = "danger", outline = true)),
    ],
    className = "mb-4",
);

row_3 = dbc_row([
    dbc_col(dbc_card(card_content, color = "light", outline = true)),
    dbc_col(dbc_card(card_content, color = "dark", outline = true)),
]);

cards = html_div([row_1, row_2, row_3]);

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
SQ._collect{R<:Relation}(_src::Tuple{R}, q::SQ.Node{SQ.GroupBy}) =
    Relations._collect(first(_src), q)

function _collect{S<:Attributed}(_src::Relation{S}, q::SQ.Node{SQ.GroupBy})
    # Make a copy, don't mutate
    src = copy(_src)
    # Map from predicates (e.g. `A > .5`) to aliasing column names (e.g. `pred_1`)
    predicates = Dict{Symbol, Expr}()
    # pre_groupby!(src, q, predicates)
    new_src, groupbys = pre_groupby!(src, q, predicates)
    return groupby(new_src, groupbys, predicates)
end

function pre_groupby!{S<:Attributed}(src::Relation{S}, q, predicates)#::Tuple{Relation, Vector{Symbol}}
    _groupbys = copy(q.args)
    groupbys = Vector{Symbol}(length(_groupbys))
    dos = Vector{SQ.Select}()
    i = 1 # count number of predicates encountered
    for (d, (j, groupby)) in zip(q.dos, enumerate(_groupbys))
        is_predicate, f, ai = SQ.parts(d)
        if is_predicate
            alias = Symbol("pred_$i")
            predicates[alias] = groupby
            groupbys[j] = alias
            # include the result of applying the predicate as a new column
            push!(dos, SQ.Select(alias, f, ai))
            i += 1
        else
            if isa(groupby, Expr) && groupby.head == :.
                groupbys[j] = groupby.args[2].args[1]
            elseif isa(groupby, Symbol)
                groupbys[j] = groupby
            end
        end
    end

    pred_cols = _collect(src, SQ.Node{SQ.Select}((SQ.DataNode(src),), [], dos))
    new_src = together(src, pred_cols)

    return new_src, groupbys
end

function groupby{S<:Attributed}(src::Relation{S}, groupbys, predicates)::SQ.Grouped
    metadata = Dict{Symbol, Any}()
    metadata[:predicates] = predicates
    # obtain the field names of the groupbys (either the names of the selected
    # columns or the predicate aliases given to groupby predicates)
    group_indices, group_levels = build_group_data(src, groupbys)
    return SQ.Grouped(src, group_indices, groupbys, group_levels, metadata)
end

function build_group_data(src, groupbys)
    p = project(src, groupbys...)
    group_indices = Dict{typeof(p).parameters[2], Vector{Int}}()
    group_levels = grow_indices!(group_indices, p)
    return group_indices, group_levels
end

@noinline function grow_indices!{S,T}(group_indices, src::Relation{S,T})::Vector{T}
    group_levels = Vector{T}()
    for (i, row) in enumerate(src) # the row is the group level
        if haskey(group_indices, row)
            push!(group_indices[row], i)
        else
            group_indices[row] = [i]
            push!(group_levels, row)
        end
    end
    return group_levels
end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

#
# package typedb.client.test.behaviour.concept.type.entitytype;
#
# public class EntityTypeSteps {}


text:
const width = height = 1000.0
const trunklength = 400.0
const scalefactor = 0.6
const startingangle = 1.5 * pi
const deltaangle = 0.2 * pi

function tree(fh, x, y, len, theta)
   if len >= 1.0
       x2 = x + len * cos(theta)
       y2 = y + len * sin(theta)
       write(fh, "<line x1='$x' y1='$y' x2='$x2' y2='$y2' style='stroke:rgb(0,0,0);stroke-width:1'/>\n")
       tree(fh, x2, y2, len * scalefactor, theta + deltaangle)
       tree(fh, x2, y2, len * scalefactor, theta - deltaangle)
    end
end

outsvg = open("tree.svg", "w")
write(outsvg,
    """<?xml version='1.0' encoding='utf-8' standalone='no'?>
    <!DOCTYPE svg PUBLIC '-//W3C//DTD SVG 1.1//EN'
    'http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd'>
    <svg width='100%%' height='100%%' version='1.1'
    xmlns='http://www.w3.org/2000/svg'>\n""")

tree(outsvg, 0.5 * width, height, trunklength, startingangle)

write(outsvg, "</svg>\n") # view file tree.svg in browser

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


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

"""
    BlockCrossValidation(sides; loss=Dict())

Cross-validation with blocks of given `sides`. Optionally,
specify `loss` function from `LossFunctions.jl` for some
of the variables. If only one side is provided, then blocks
become cubes.

## References

* Roberts et al. 2017. [Cross-validation strategies for data with
  temporal, spatial, hierarchical, or phylogenetic structure]
  (https://onlinelibrary.wiley.com/doi/10.1111/ecog.02881)
* Pohjankukka et al. 2017. [Estimating the prediction performance
  of spatial models via spatial k-fold cross-validation]
  (https://www.tandfonline.com/doi/full/10.1080/13658816.2017.1346255)
"""
struct BlockCrossValidation{S} <: AbstractErrorEstimator
  sides::S
  loss::Dict{Symbol,SupervisedLoss}
end

BlockCrossValidation(sides; loss=Dict()) =
  BlockCrossValidation{typeof(sides)}(sides, loss)

function Base.error(solver::AbstractLearningSolver,
                    problem::LearningProblem,
                    eestimator::BlockCrossValidation)
  sdata = sourcedata(problem)
  ovars = outputvars(task(problem))
  loss  = eestimator.loss
  sides = eestimator.sides
  for var in ovars
    if var ∉ keys(loss)
      loss[var] = defaultloss(sdata[1,var])
    end
  end

  bsides    = length(sides) > 1 ? sides : ntuple(i->sides, ndims(sdata))
  blocks    = partition(sdata, BlockPartitioner(bsides))
  bsubsets  = subsets(blocks)
  neighbors = blocks[:neighbors]
  allblocks = 1:length(blocks)

  solutions = pmap(allblocks) do b
    # source and target blocks
    sblocks = setdiff(allblocks, vcat(neighbors[b], b))
    tblocks = [b]

    # source and target indices
    sinds = reduce(vcat, bsubsets[sblocks])
    tinds = reduce(vcat, bsubsets[tblocks])

    # setup and solve learning sub-problem
    subproblem = LearningProblem(view(sdata, sinds),
                                 view(sdata, tinds),
                                 task(problem))
    solve(subproblem, solver)
  end

  result = pmap(ovars) do var
    losses = map(allblocks) do b
      dview = view(sdata, bsubsets[b])
      y = dview[var]
      ŷ = solutions[b][var]
      value(loss[var], y, ŷ, AggMode.Mean())
    end
    var => mean(losses)
  end

  Dict(result)
end


text:
using Simulacrum
using HyperGraphs, Symbolics

# cloning one hypergraph
chx = ChemicalHyperGraph{Num}()
@test (clone(chx, 1) == chx) && !(clone(chx, 1) === chx)

# cloning n hypergraphs
@variables t X(t) Y(t) k
che = ChemicalHyperEdge([X], [Y], k)
n = rand(1:100)
cloned_chxs = clone(ChemicalHyperGraph(che), n)
@test length(cloned_chxs) == n

# rename_weights, where all rates are symbolic variables
chxs = clone(ChemicalHyperGraph(che), n, rename_weights = true)
weight.([hyperedges(chx)[1] for chx in chxs])

# rename_weights, where all rates are integers (nothing should happen)
chxs = clone(ChemicalHyperGraph(ChemicalHyperEdge([X], [Y], 2)), n, rename_weights = true)
@test all(weight.([hyperedges(chx)[1] for chx in chxs]) .== 2)

# rename_weights, where some weights are symbolic variables and some weights are integers
chxs = clone(ChemicalHyperGraph([ChemicalHyperEdge([X], [Y], 2), ChemicalHyperEdge([X], [Y], k)]), n, rename_weights = true)
@test all(weight.([hyperedges(chx)[1] for chx in chxs]) .== 2)
weight.([hyperedges(chx)[2] for chx in chxs])


text:
const ENTITIES_URL = "https://s3.us-east-2.amazonaws.com/dgl.ai/dataset/"
const ENTITIES_DATASETS = [:aifb, :am, :mutag, :bgs]

entities_init() = register(DataDep(
    "Entities",
    """
    The relational entities networks "AIFB", "MUTAG", "BGS" and "AM" from
    the `"Modeling Relational Data with Graph Convolutional Networks"
    <https://arxiv.org/abs/1703.06103>`_ paper.
    Training and test splits are given by node indices.
    """,
    [joinpath(ENTITIES_URL, "$(d).tgz") for d in ENTITIES_DATASETS],
    "e58bcfddd240d9bbc830bcae74e9854f1f778a96d3072930395f47d3c8e6f342";
    post_fetch_method=preprocess_entities,
))

function preprocess_entities(local_path)
    for dataset in ENTITIES_DATASETS
        tgz_file = @datadep_str "Entities/$(dataset).tgz"
        dir = joinpath(dirname(tgz_file), "$(dataset)")
        unzip_tgz(tgz_file, dir)

        # nt_file = joinpath(dir, "$(dataset)_stripped.nt.gz")
        # train_file = joinpath(dir, "trainingSet.tsv")
        # test_file = joinpath(dir, "testSet.tsv")

        # py"""
        # import rdflib as rdf
        # import gzip
        # g = rdf.Graph()
        # with gzip.open($nt_file, 'rb') as file:
        #     g.parse(file, format='nt')
        # """

        # train = CSV.read(train_file, delim='\t')
        # test = CSV.read(test_file, delim='\t')
        # uri = HTTP.URI("http://data.bgs.ac.uk/id/Geochronology/Division/CN")
    end
end

struct Entities <:Dataset
    dataset::Symbol

    function Entities(dataset::Symbol)
        @assert dataset in ENTITIES_DATASETS "`dataset` should be one of :aifb, :am, :mutag, :bgs."
        new(dataset)
    end
end

# function traindata(ent::Entities)
#     file = @datadep_str "Entities/$(ent.dataset).train.jld2"
#     @load file graph train_X train_y
#     graph, train_X, train_y
# end

# function testdata(ent::Entities)
#     file = @datadep_str "Entities/$(ent.dataset).test.jld2"
#     @load file graph test_X test_y
#     graph, test_X, test_y
# end


text:
#############Test Metropolis-Hastings samplers####################

fv1 = [0.1,0.2]
fv2 = [0.3,0.4]
fr1 = samples(:base,2,1,Float64,Float64) ; copy!(fr1.values,fv1)
fr2 = samples(:base,2,2,Float64,Float64) ; copy!(fr2.values,3,fv2,1,2)

for args in [((:mh,:normal,0.1,eye(2)),4,"Scalable Metropolis-Hastings","Metropolis"),
             ((:mh,:normal,0.2,ones(2)),2,"Scalable Metropolis-Hastings","Metropolis"),
             ((:mh,:lognormal,0.1,2),1,"Scalable Metropolis-Hastings","Metropolis-Hastings"),
             ((:mh,:bactrian,0.1,0.5*ones(2),:normal,0.95),3,"Scalable Metropolis-Hastings","Metropolis")]
    s = sampler(args[1]...)
    t = samplerstate(s,args[2],Float64,Float64)
    @test numparas(t.from) == 2
    @test numsamples(t.from) == 1
    @test numparas(t.proposals) == 2
    @test numsamples(t.proposals) == args[2]
    @test numparas(t) == 2
    @test numsamples(t) == args[2]
    @test from(t) == t.from
    @test proposals(t) == t.proposals
    @test acceptance(t) == zeros(args[2])
    @test GeneralizedMetropolisHastings.samplername(s) == args[3]
    @test GeneralizedMetropolisHastings.samplerstatename(t) == args[4]
    setfrom!(t,fr1)
    @test from(t).values == fv1
    setfrom!(t,fr2,1)
    @test from(t).values == zeros(2)
    setfrom!(t,fr2,2)
    @test from(t).values == fv2
end

function testprepare(t,f1,f2,x1,x2,s)
    @test t.from.values == x1
    @test_approx_eq f1(t.density.distribution) x2
    @test_approx_eq f2(t.density.distribution) s
end

function testsym(t,a)
    f = t.from.loglikelihood[1] + t.from.logprior[1]
    @test a == t.proposals.loglikelihood + t.proposals.logprior - f
end

function testasym(t,a)
    f = t.from.loglikelihood[1] + t.from.logprior[1]
    for i=1:length(a)
        d21 = deepcopy(t.density.distribution)
        d21 = recenter(d21,t.proposals.values[:,i])
        k12 = Distributions.logpdf(t.density.distribution,t.proposals.values[:,i])
        k21 = Distributions.logpdf(d21,t.from.values)
        @test_approx_eq a[i] (t.proposals.loglikelihood[i] +t.proposals.logprior[i] -f + k12 - k21)
    end
end

function testsmmala(t,a)
    f = t.from.loglikelihood[1] + t.from.logprior[1]
    for i=1:length(a)
        d21 = deepcopy(t.density.distribution)
        d21 = update(d21,GeneralizedMetropolisHastings._meancov(t,i)...)
        k12 = Distributions.logpdf(t.density.distribution,t.proposals.values[:,i])
        k21 = Distributions.logpdf(d21,t.from.values)
        @test_approx_eq a[i] (t.proposals.loglikelihood[i] +t.proposals.logprior[i] -f + k12 - k21)
    end
end

#test prepare! functions
for p in [(:normal,mean,cov,identity,testsym),
          (:lognormal,Distributions.location,Distributions.scale,log,testasym)]
    s = sampler(:mh,p[1],0.1,eye(2))
    tind = samplerstate(s,1,Float64,Float64)
    taux = samplerstate(s,3,Float64,Float64)
    copy!(tind.proposals.values,fv1)
    copy!(taux.proposals.values,3,2fv1,1,2)
    #test the prepare! functions
    prepare!(tind)
    testprepare(tind,p[2],p[3],zeros(2),p[4](zeros(2)),0.01eye(2))
    prepare!(tind,true)
    testprepare(tind,p[2],p[3],fv1,p[4](fv1),0.01eye(2))
    prepareauxiliary!(tind,taux)
    testprepare(taux,p[2],p[3],fv1,p[4](fv1),0.01eye(2))
    prepareindicator!(tind)
    testprepare(tind,p[2],p[3],fv1,p[4](fv1),0.01eye(2))
    prepareindicator!(tind,taux,2)
    testprepare(tind,p[2],p[3],2fv1,p[4](2fv1),0.01eye(2))
    prepare!(tind,true)
    testprepare(tind,p[2],p[3],fv1,p[4](fv1),0.01eye(2))
    #test the propose! function
    d = distribution(p[1],fv1,eye(2)) ; d = rescale(d,0.1)
    srand(46732) ; propose!(tind)
    srand(46732) ; r1 = rand!(d,similar(tind.proposals.values))
    srand(46733) ; propose!(taux)
    srand(46733) ; r2 = rand!(d,similar(taux.proposals.values))
    @test tind.proposals.values == r1
    @test taux.proposals.values == r2
    #test the acceptance function
    copy!(tind.from.loglikelihood,[-1.0])
    copy!(tind.from.logprior,[-5.0])
    copy!(tind.proposals.loglikelihood,[-0.5])
    copy!(tind.proposals.logprior,[-4.0])
    copy!(taux.from.loglikelihood,[-1.0])
    copy!(taux.from.logprior,[-5.0])
    copy!(taux.proposals.loglikelihood,[-0.1,-0.2,-0.3])
    copy!(taux.proposals.logprior,[-3.0,-3.0,-3.0])
    a1 = acceptance!(tind)
    a2 = acceptance!(taux)
    @test length(a1) == 1
    @test length(a2) == 3
    p[5](tind,a1)
    p[5](taux,a2)
    tune!(tind,2.0)
    tune!(taux,0.5)
    @test_approx_eq p[3](tind.density.distribution) 0.04*eye(2)
    @test_approx_eq p[3](taux.density.distribution) 0.0025*eye(2)
    println("Test show() function for ",p[1])
    show(s)
    show(tind)
    show(taux)
    println("======================================")
end

#############Test Adaptive Metropolis samplers####################

ru1 = GeneralizedMetropolisHastings._runningstate(Val{:normal},2,0.1)
v = [1.0,2.0]
@test ru1.iteration == 0
@test ru1.runningmean == zeros(2)
@test ru1.runningcov == zeros(2,2)
@test ru1.scale == 0.1
@test size(ru1.tempdiff) == (2,)
@test size(ru1.densitycov) == (2,2)

@test GeneralizedMetropolisHastings._covupdatefactor(1,Float64) == 0.0
@test GeneralizedMetropolisHastings._covupdatefactor(2,Float64) == 0.0
@test GeneralizedMetropolisHastings._covupdatefactor(3,Float64) == 0.5

@test_approx_eq GeneralizedMetropolisHastings._updatetempdiff!(zeros(2),v/3,v) 2v/3
@test_approx_eq GeneralizedMetropolisHastings._updaterunningcov!(eye(2),2v/3,3) eye(2)/2+2v/3*2v'/3/3
@test_approx_eq GeneralizedMetropolisHastings._updaterunningmean!(v/3,v,3) (2*v/3+v)/3
@test_approx_eq GeneralizedMetropolisHastings._updatedensitycov!(zeros(2,2),ones(2,2),0.1) ones(2,2)+0.1*0.1*eye(2)

niterations = 10
rv = rand!(zeros(2,niterations))
for i=1:niterations
    GeneralizedMetropolisHastings.updaterunning!(ru1,rv[:,i])
end

#test showing that the covariance is indeed the correct sample covariance
@test ru1.iteration == niterations
@test_approx_eq ru1.runningmean mean(rv,2)
@test_approx_eq ru1.runningcov cov(rv,vardim=2)

as1 = sampler(:adaptive,0.1,2)
atind = samplerstate(as1,1,Float64,Float64)
ataux = samplerstate(as1,3,Float64,Float64,true)

@test numparas(as1) == 2
@test GeneralizedMetropolisHastings.samplername(as1) == "Adaptive Normal Metropolis"
@test GeneralizedMetropolisHastings.samplerstatename(atind) == "Adaptive Normal Metropolis"
@test isnull(atind.runningstate) == false
@test isnull(ataux.runningstate) == true

copy!(atind.from.values,v)
copy!(ataux.from.values,v)
GeneralizedMetropolisHastings.updaterunning!(atind)
GeneralizedMetropolisHastings.updaterunning!(ataux)
@test mean(atind.density.distribution) == zeros(2)
@test_approx_eq cov(atind.density.distribution) get(atind.runningstate).densitycov
@test mean(ataux.density.distribution) == zeros(2)
@test_approx_eq cov(ataux.density.distribution) 0.01*eye(2)

copy!(atind.proposals.values,2.5v)
copy!(ataux.proposals.values,3,3.5v,1,2)
#test the prepare! functions
prepare!(atind)
@test get(atind.runningstate).iteration == 2
testprepare(atind,mean,cov,v,zeros(2),0.01eye(2))
prepare!(atind,true)
@test get(atind.runningstate).iteration == 3
testprepare(atind,mean,cov,2.5v,zeros(2),get(atind.runningstate).densitycov)
prepareauxiliary!(atind,ataux)
testprepare(ataux,mean,cov,2.5v,zeros(2),get(atind.runningstate).densitycov)
prepareindicator!(atind)
@test get(atind.runningstate).iteration == 4
testprepare(atind,mean,cov,2.5v,zeros(2),get(atind.runningstate).densitycov)
prepareindicator!(atind,ataux,2)
@test get(atind.runningstate).iteration == 5
testprepare(atind,mean,cov,3.5v,zeros(2),get(atind.runningstate).densitycov)
#test the propose! function
dind = distribution(:normal,atind.from.values,cov(atind.density.distribution))
daux = distribution(:normal,ataux.from.values,cov(ataux.density.distribution))
srand(46732) ; propose!(atind)
srand(46732) ; r1 = rand!(dind,similar(atind.proposals.values))
srand(46733) ; propose!(ataux)
srand(46733) ; r2 = rand!(daux,similar(ataux.proposals.values))
@test atind.proposals.values == r1
@test ataux.proposals.values == r2
#test the acceptance function
copy!(atind.from.loglikelihood,[-1.0])
copy!(atind.from.logprior,[-5.0])
copy!(atind.proposals.loglikelihood,[-0.5])
copy!(atind.proposals.logprior,[-4.0])
copy!(ataux.from.loglikelihood,[-1.0])
copy!(ataux.from.logprior,[-5.0])
copy!(ataux.proposals.loglikelihood,[-0.1,-0.2,-0.3])
copy!(ataux.proposals.logprior,[-3.0,-3.0,-3.0])
a1 = acceptance!(atind)
a2 = acceptance!(ataux)
@test length(a1) == 1
@test length(a2) == 3
testsym(atind,a1)
testsym(ataux,a2)
tune!(atind,2.0)
tune!(ataux,0.5)
@test get(atind.runningstate).scale == 0.2
println("Test show() function for ",:adaptive)
show(as1)
show(atind)
show(ataux)
println("======================================")

########################
### Tests of SmMALA
for p in [(false,:full,:standard,:full),
          (true,:full,:tr,:full),
          (true,:tangent,:tr,:tangent)]
    smpol = GeneralizedMetropolisHastings.SmMALAPolicy(p[1],p[2])
    @test traitvalue(smpol.smmala) == p[3]
    @test traitvalue(smpol.tensor) == p[4]
end

@test_approx_eq GeneralizedMetropolisHastings._adddiag!(eye(2),0.5) 1.5*eye(2)
@test_approx_eq GeneralizedMetropolisHastings._addfrom!(ones(2),0.5ones(2)) 1.5*ones(2)

from1 = 0.5ones(2)
grad1 = [-0.5,0.3]
ten1 = [1.4 1.2;1.2 1.4]
m1,c1 = GeneralizedMetropolisHastings._meancov(Val{:standard},from1,grad1,ten1,0.1)
m2,c2 = GeneralizedMetropolisHastings._meancov(Val{:tr},from1,grad1,ten1,0.1)
@test_approx_eq m1 from1+ten1\grad1*0.1*0.1/2
@test_approx_eq m2 from1+(ten1+eye(2)/0.1)\grad1
@test_approx_eq c1 inv(ten1)*0.01
@test_approx_eq c2 inv(ten1+eye(2)/0.1)

for p in [(false,"Simplified mMALA",:standard,:full),
          (true,"Trust Region Simplified mMALA",:tr,:full),
          (10,"Trust Region Simplified mMALA with Tangent Tensor",:tr,:tangent)]
    d = density(:normal,zeros(2),eye(2))
    sm = sampler(:smmala,0.1,2,p[1])
    smind = samplerstate(sm,1,Float64,Float64)
    smaux = samplerstate(sm,3,Float64,Float64)

    @test numparas(sm) == 2
    @test GeneralizedMetropolisHastings.samplername(sm) == p[2]
    @test GeneralizedMetropolisHastings.samplerstatename(smind) == p[2]

    @test traitvalue(smind.policy.smmala) == p[3]
    @test traitvalue(smind.policy.tensor) == p[4]

    copy!(smind.from.values,from1)
    copy!(smind.from.gradloglikelihood,grad1)
    copy!(smind.from.tensorloglikelihood,ten1)
    copy!(smind.proposals.values,2from1)
    copy!(smind.proposals.gradloglikelihood,grad1)
    copy!(smind.proposals.tensorloglikelihood,ten1)
    copy!(smaux.from.values,3from1)
    copy!(smaux.from.gradloglikelihood,grad1)
    copy!(smaux.from.tensorloglikelihood,ten1)
    copy!(smaux.proposals.values,[from1 2from1 3from1])
    copy!(smaux.proposals.gradloglikelihood,[grad1 0.5grad1 0.3grad1])
    copy!(smaux.proposals.tensorloglikelihood,cat(3,3ten1,1.5ten1,ten1))

    @test GeneralizedMetropolisHastings._meancov(smind) == GeneralizedMetropolisHastings._meancov(Val{p[3]},from1,grad1,ten1,0.1)
    @test GeneralizedMetropolisHastings._meancov(smind,1) == GeneralizedMetropolisHastings._meancov(Val{p[3]},2from1,grad1,ten1,0.1)
    @test GeneralizedMetropolisHastings._meancov(smaux) == GeneralizedMetropolisHastings._meancov(Val{p[3]},3from1,grad1,ten1,0.1)
    @test GeneralizedMetropolisHastings._meancov(smaux,2) == GeneralizedMetropolisHastings._meancov(Val{p[3]},2from1,0.5grad1,1.5ten1,0.1)

    GeneralizedMetropolisHastings.update!(d,GeneralizedMetropolisHastings._meancov(Val{p[3]},from1,grad1,ten1,0.1)...)
    GeneralizedMetropolisHastings._updatedensity!(smind)
    @test_approx_eq mean(smind.density.distribution) mean(d.distribution)
    @test_approx_eq cov(smind.density.distribution) cov(d.distribution)
    GeneralizedMetropolisHastings.update!(d,GeneralizedMetropolisHastings._meancov(Val{p[3]},2from1,0.5grad1,1.5ten1,0.1)...)
    GeneralizedMetropolisHastings._updatedensity!(smaux,2)
    @test_approx_eq mean(smaux.density.distribution) mean(d.distribution)
    @test_approx_eq cov(smaux.density.distribution) cov(d.distribution)

    #test the prepare! functions
    prepare!(smind)
    testprepare(smind,mean,cov,from1,GeneralizedMetropolisHastings._meancov(Val{p[3]},from1,grad1,ten1,0.1)...)
    prepare!(smind,true)
    testprepare(smind,mean,cov,2from1,GeneralizedMetropolisHastings._meancov(Val{p[3]},2from1,grad1,ten1,0.1)...)
    prepareauxiliary!(smind,smaux)
    testprepare(smaux,mean,cov,2from1,GeneralizedMetropolisHastings._meancov(Val{p[3]},2from1,grad1,ten1,0.1)...)
    prepareindicator!(smind)
    testprepare(smind,mean,cov,2from1,GeneralizedMetropolisHastings._meancov(Val{p[3]},2from1,grad1,ten1,0.1)...)
    prepareindicator!(smind,smaux,2)
    testprepare(smind,mean,cov,2from1,GeneralizedMetropolisHastings._meancov(Val{p[3]},2from1,0.5grad1,1.5ten1,0.1)...)
    prepare!(smind,true)
    testprepare(smind,mean,cov,2from1,GeneralizedMetropolisHastings._meancov(Val{p[3]},2from1,grad1,ten1,0.1)...)

    #test the propose! function
    dind = distribution(:normal,GeneralizedMetropolisHastings._meancov(Val{p[3]},2from1,grad1,ten1,0.1)...)
    daux = distribution(:normal,GeneralizedMetropolisHastings._meancov(Val{p[3]},2from1,grad1,ten1,0.1)...)
    srand(46732) ; propose!(smind)
    srand(46732) ; r1 = rand!(dind,similar(smind.proposals.values))
    srand(46733) ; propose!(smaux)
    srand(46733) ; r2 = rand!(daux,similar(smaux.proposals.values))
    @test smind.proposals.values == r1
    @test smaux.proposals.values == r2

    #test the acceptance function
    copy!(smind.from.loglikelihood,[-1.0])
    copy!(smind.from.logprior,[-5.0])
    copy!(smind.proposals.loglikelihood,[-0.5])
    copy!(smind.proposals.logprior,[-4.0])
    copy!(smaux.from.loglikelihood,[-1.0])
    copy!(smaux.from.logprior,[-5.0])
    copy!(smaux.proposals.loglikelihood,[-0.1,-0.2,-0.3])
    copy!(smaux.proposals.logprior,[-3.0,-3.0,-3.0])
    a1 = acceptance!(smind)
    a2 = acceptance!(smaux)
    @test length(a1) == 1
    @test length(a2) == 3
    testsmmala(smind,a1)
    testsmmala(smaux,a2)

    println("Test show() function for smmala $(p[3]) $(p[4])")
    show(sm)
    show(smind)
    show(smaux)
    println("======================================")
end


text:
# Small (currently handmaintained) list of abbreviations for units in commonly-found orders
# of magnitude. There's a delicate balance to be struck between polluting the namespace
# and being useful. We might in the future give up and just generate all the orders of
# magnitude for everything.
module ShortUnits
    export µg, mg, g, kg, µm, mm, cm, m, km, mJ, J, MJ, HJ, fF, pF, nF, µF, mF,
        F, pΩ, nΩ, µΩ, mΩ, Ω, kΩ, MΩ, mV, V, MV, mHz, Hz, kHz, MHz, GHz, THz, fA, pA,
        nA, µA, mA, A, C, N, mol, ns, ms, µs, s, S, nW, µW, mW, W, kW, MW, GW, eV,
        K, mPa, Pa, kPa, MPa, GPa, nm, torr, mtorr, atm, rad, deg, sr,
        MeV, GeV, TeV

    using SIUnits

    const µg    = Micro * Gram
    const mg    = Milli * Gram
    const g     = Gram
    const kg    = KiloGram

    const nm    = Nano * Meter
    const µm    = Micro * Meter
    const mm    = Milli * Meter
    const cm    = CentiMeter
    const m     = Meter
    const km    = Kilo*Meter

    const mJ    = Milli*Joule
    const J     = Joule
    const MJ    = Mega*Joule
    const GJ    = Giga*Joule

    const fF    = Femto*Farad
    const pF    = Pico*Farad
    const nF    = Nano*Farad
    const µF    = Micro*Farad
    const mF    = Milli*Farad
    const F     = Farad

    const pΩ    = Pico*Ohm
    const nΩ    = Nano*Ohm
    const µΩ    = Micro*Ohm
    const mΩ    = Milli*Ohm
    const Ω     = Ohm
    const kΩ    = Kilo*Ohm
    const MΩ    = Mega*Ohm

    const mV    = Milli*Volt
    const V     = Volt
    const MV    = Mega*Volt

    const ns    = Nano*Second
    const ms    = Milli*Second
    const µs    = Micro*Second
    const s     = Second

    const S     = Siemens

    const mHz   = Milli*Hertz
    const Hz    = Hertz
    const kHz   = Kilo*Hertz
    const MHz   = Mega*Hertz
    const GHz   = Giga*Hertz
    const THz   = Tera*Hertz

    const fA    = Femto*Ampere
    const pA    = Pico*Ampere
    const nA    = Nano*Ampere
    const µA    = Micro*Ampere
    const mA    = Milli*Ampere
    const A     = Ampere

    const mPa   = Milli*Pascal
    const Pa    = Pascal
    const kPa   = Kilo*Pascal
    const MPa   = Mega*Pascal
    const GPa   = Giga*Pascal

    const C     = Coulomb

    const N     = Newton

    const nW    = Nano*Watt
    const µW    = Micro*Watt
    const mW    = Milli*Watt
    const W     = Watt
    const MW    = Mega*Watt
    const kW    = Kilo*Watt
    const GW    = Giga*Watt

    const mol   = Mole

    const K     = Kelvin

    const eV    = ElectronVolt
    const MeV   = Mega*ElectronVolt
    const GeV   = Giga*ElectronVolt
    const TeV   = Tera*ElectronVolt

    const torr  = Torr
    const mtorr = Milli*Torr

    const atm   = Atmosphere

    const rad   = Radian
    const deg   = Degree
    const sr    = Steradian

    # Support the other unicode codepoint for $\mu$
    # (namely the one we get on the REPL)
    export μg, μm, μF, μA, μW
    const μg = µg
    const μm = µm
    const μF = µF
    const μA = µA
    const μW = µW

end


text:
# Preventing polarization from forming as a function of number of attributes.
# Figure 1
using DrWatson
quickactivate(@__DIR__)

using PolarizationFramework

ns = [5, 9]#[9,15,25]
reps = [100, 100]
reps_dict = Dict(zip(ns, reps))
# two parts of simulations
# gs = [1:2:21...]
gs = [25, 29, 33, 37, 41, 45, 49, 55, 61, 67, 73, 81, 89, 97]
threshold = 0.5;
vs = [4, @onlyif("attr_types" == "OA", 1000)] #includes CA

attr_types = ["BA", "UA", "OA", "UPA"]
# attr_types = ["BA", "OA", "UPA"]
# attr_types = ["OA"]

gammas = [0.5, 1.5, 4]

all_params = @strdict(ns, gs, threshold, vs, attr_types)
dicts = dict_list(all_params)
[d["reps"] = reps_dict[d["ns"]] for d in dicts] #adding reps
# [d["gammas"] = gammas for d in dicts] #adding gammas

for params in dicts
    n, g, attr_type, v, rep = let
        @unpack ns, threshold, reps, gs, attr_types, vs = params
        ns, gs, attr_types, vs, reps
    end

    println("Started n=$n and g=$g and attr_type=", attr_type, " and v=$v.")

    if attr_type == "UA"
        attr = UnorderedAttributes(g, threshold, v)
    elseif attr_type == "BA"
        attr = BinaryAttributes(g)
    elseif attr_type == "OA"
        attr = OrderedAttributes(g, threshold, v)
    elseif attr_type == "UPA"
        attr = UnorderedPositiveAttributes(g, threshold, v)
    else
        throw(attr_type)
    end

    r = using_heider_attr(
        n,
        attr,
        gammas,
        rep,
        3000.0,
        "Heider7!";
        disp_each = 0,
        disp_more_every = 600,
        save_each = 600,
        files_folder = ["data", "sims"],
        filename_prefix = "NumerFig1",
    )
end

max_stars_repo_name: pjgorski/PolarizationFramework

max_issues_repo_name: pjgorski/PolarizationFramework

max_forks_repo_name: pjgorski/PolarizationFramework


text:
module DataBagsExamples

using DataBags

export
    BagEx1,
    BagEx2,
    BagEx3

#------------------------------------------------------------------------------
# EXAMPLE 1

# Define a concrete sub-type of `DataBags.AbstractDataBag`.
struct BagEx1{K,V,D<:AbstractDict{K,V}} <: DataBags.AbstractDataBag{K,V,D}
    data::D # object used to store key-value pairs
end

# Override `DataBags.contents` to yield the dictionary that stores the data.
DataBags.contents(A::BagEx1) = Base.getfield(A, :data)

#------------------------------------------------------------------------------
# EXAMPLE 2

# Define a concrete sub-type of `DataBags.AbstractDataBag`.
struct BagEx2{K,V,D<:AbstractDict{K,V}} <: DataBags.AbstractDataBag{K,V,D}
    # Object used to store key-value pairs.
    data::D

    # Explicitely define inner constructor to avoid outer constructor
    # automatically created by Julia.
    BagEx2{K,V,D}(data::D) where {K,V,D<:AbstractDict{K,V}} = new{K,V,D}(data)
end

# Override `DataBags.contents` to yield the dictionary that stores the data.
DataBags.contents(A::BagEx2) = Base.getfield(A, :data)

# Override `DataBags.wrap` to create an instance of `BagEx2` that stores
# its data in a given dictionary.
DataBags.wrap(::Type{BagEx2}, data::D) where {K,V,D<:AbstractDict{K,V}} =
    BagEx2{K,V,D}(data)

# Outer constructor.
BagEx2(args...; kdws...) =
    wrap(BagEx2, contents(Dict{Any,Any}, args...; kdws...))

#------------------------------------------------------------------------------
# EXAMPLE 3

struct BagEx3 <: DataBags.AbstractDataBag{Symbol,Any,Dict{Symbol,Any}}
    data::Dict{Symbol,Any}
    BagEx3(args...; kwds...) =
        new(DataBags.contents(Dict{Symbol,Any}, args...; kwds...))
end
DataBags.contents(A::BagEx3) = Base.getfield(A, :data)

end


text:
################################################################################
#                      TEST: Kinematics using Quaternion
################################################################################

for i = 1:samples
    # Create a random quaternion.
    qba0 = Quaternion(randn(), randn(), randn(), randn())
    qba0 = qba0/norm(qba0)

    # Create a random velocity vector.
    wba_a = randn(3)

    # Propagate the initial DCM using the sampled velocity vector.
    Δ   = 0.0001
    qba = copy(qba0)
    for k = 1:1000
        dqba = dquat(qba,vect(qba\wba_a*qba))

        qba  = qba + dqba*Δ
    end

    # In the end, the vector aligned with w must not change.
    v1 = vect(qba0\wba_a*qba0)
    v2 = vect( qba\wba_a*qba)

    @test norm(v2-v1) < 1e-10
end


text:
using PGFPlots
using QHull
using JLD2
using LazySets
include("models/problems.jl")
include("OverApprox/src/OA_relational_util.jl")
include("MIP/src/overt_to_mip.jl")
include("MIP/src/mip_utils.jl")
include("models/single_pendulum/single_pend.jl")

controller_type = "small" # change to "small" or "big" as needed
if controller_type == "small"
    JLD2.@load "examples/jmlr/data/single_pend_-0.2167/single_pendulum_reachability_small_controller_data_1step.jld2"
    one_step_state_sets = symbolic_state_sets
    JLD2.@load "examples/jmlr/data/single_pend_-0.2167/single_pendulum_reachability_small_controller_data.jld2"
elseif controller_type == "big"
    JLD2.@load "examples/jmlr/data/single_pend_-0.2167/single_pendulum_reachability_big_controller_data_1step.jld2"
    one_step_state_sets = symbolic_state_sets
    JLD2.@load "examples/jmlr/data/single_pend_-0.2167/single_pendulum_reachability_big_controller_data.jld2"
end

# monte_carlo_simulate
mc_state_sets, xvec, x0, mc_meas_sets, yvec, y0 = monte_carlo_simulate(query, input_set, n_sim=1000000)

# clean up sets
init_state, reachable_state_sets = clean_up_sets(concrete_state_sets, symbolic_state_sets, concretization_intervals)

define_color("concrete_color", 0x139EAB)
define_color("symbolic_color", 0x9BFF85)
define_color("mc_color", 0x38caff)
# input_set_style = "dashed, blue, mark=none"
concrete_line_style = "solid, concrete_color, mark=none"
symbolic_line_style = "solid, symbolic_color, very thick, mark=none"
mc_line_style = "solid, mc_color, ultra thick, mark=none"
conc_style = "solid, concrete_color, mark=none, fill=concrete_color"
sym_style = "solid, symbolic_color, mark=none, fill=symbolic_color"
mc_style = "solid, mc_color, mark=none, fill=mc_color"
conc_style_trans = conc_style*", fill opacity=0.5"
sym_style_trans = sym_style*", fill opacity=0.5"
mc_style_trans = mc_style*", fill opacity=0.5"

###################################################
# Plot 1: 1D plots
###################################################

# plot intervals using lines in PGF plots
fig = PGFPlots.Axis(style="width=10cm, height=10cm", ylabel="\$x_1\$", xlabel="timesteps", title="Single Pendulum Reachable Sets, $controller_type Controller")

dims = [1]
conc_subsets = get_interval_subsets(one_step_state_sets, dims)
sym_subsets = get_interval_subsets(reachable_state_sets, dims)
mc_subsets = get_interval_subsets(mc_state_sets, dims)
deleteat!(mc_subsets, 1) # pop init set off start
# plot init set
push!(fig, PGFPlots.Plots.Linear([0.,0.], [low(input_set)[dims]..., high(input_set)[dims]...], style=concrete_line_style, legendentry="Concrete Sets"))
# plot init set
push!(fig, PGFPlots.Plots.Linear([0.,0.], [low(input_set)[dims]..., high(input_set)[dims]...], style=symbolic_line_style, legendentry="OVERT Hybrid Symbolic Sets"))
# plot init set
push!(fig, PGFPlots.Plots.Linear([0.,0.], [low(input_set)[dims]..., high(input_set)[dims]...], style=mc_line_style, legendentry="Monte Carlo Simulations"))
# avoid set
push!(fig, PGFPlots.Plots.Linear([-1., 26], [-0.2617, -0.2617], style="red, dashed", mark="none", legendentry="Unsafe Set"))
# plot the rest of the sets
for t in 1:query.ntime
    # if t == query.ntime # last timestep
    #     # include legend entry
    #     timestep = Float64[t,t]
    #     push!(fig, PGFPlots.Plots.Linear(timestep, conc_subsets[t], style=concrete_line_style, legendentry="Concrete Sets"))
    #     push!(fig, PGFPlots.Plots.Linear(timestep,sym_subsets[t], style=symbolic_line_style, markSize=1, legendentry="OVERT Hybrid Symbolic Sets"))
    #     push!(fig, PGFPlots.Plots.Linear(timestep,mc_subsets[t], style=mc_line_style, markSize=1, legendentry="Monte Carlo Simulations"))
    # else
        timestep = Float64[t,t]
        push!(fig, PGFPlots.Plots.Linear(timestep, conc_subsets[t], style=concrete_line_style))
        push!(fig, PGFPlots.Plots.Linear(timestep,sym_subsets[t], style=symbolic_line_style, markSize=1))
        push!(fig, PGFPlots.Plots.Linear(timestep,mc_subsets[t], style=mc_line_style, markSize=1))
    # end
end

fig.legendStyle = "at={(1.05,1.0)}, anchor=north west"

PGFPlots.save("examples/jmlr/plots/single_pendulum_$(controller_type).tex", fig)
PGFPlots.save("examples/jmlr/plots/single_pendulum_$(controller_type).pdf", fig)

###################################################
# Plot 2: 2D plots
###################################################
dims = [1,2]
fig = PGFPlots.Axis(style="width=10cm, height=10cm, axis equal image", xlabel="\$x_1\$", ylabel="\$x_2\$", title="Single Pendulum Reachable Sets, $controller_type Controller")

# plot init set in both concrete and hybrid colors
inputx, inputy = get_rectangle(input_set, dims)
push!(fig, PGFPlots.Plots.Linear(inputx, inputy, style=conc_style, legendentry="Concrete Sets"))
push!(fig, PGFPlots.Plots.Linear(inputx, inputy, style=sym_style, legendentry="OVERT Hybrid Symbolic Sets"))
push!(fig, PGFPlots.Plots.Linear(inputx, inputy, style=mc_style, legendentry="Monte Carlo Simulations"))
# avoid set
push!(fig, PGFPlots.Plots.Linear([-0.2617, -0.2617], [-0.9, 0.5], style="red, dashed", mark="none", legendentry="Unsafe Set"))

for t in 1:query.ntime
    if t == query.ntime
        push!(fig, PGFPlots.Plots.Linear( get_rectangle(one_step_state_sets[t], dims)..., style=conc_style))
        push!(fig, PGFPlots.Plots.Linear( get_rectangle(reachable_state_sets[t], dims)..., style=sym_style))
        # push!(fig, PGFPlots.Plots.Linear( get_rectangle(mc_state_sets[t+1], dims)..., style=mc_style_transparent))
        # plot mc_sim points
        points = xvec[:, t, dims]
        #plot hull
        border_idx = chull(points).vertices
        x = points[border_idx, 1]
        y = points[border_idx, 2]
        push!(fig, PGFPlots.Plots.Linear([x..., x[1]], [y..., y[1]], style=mc_style))
    else
        push!(fig, PGFPlots.Plots.Linear( get_rectangle(one_step_state_sets[t], dims)..., style=conc_style_trans))
        push!(fig, PGFPlots.Plots.Linear( get_rectangle(reachable_state_sets[t], dims)..., style=sym_style_trans))
        # push!(fig, PGFPlots.Plots.Linear( get_rectangle(mc_state_sets[t+1], dims)..., style=mc_style_transparent))
        # plot mc simulations hull
        points = xvec[:, t, dims]
        border_idx = chull(points).vertices
        x = points[border_idx, 1]
        y = points[border_idx, 2]
        push!(fig, PGFPlots.Plots.Linear([x..., x[1]], [y..., y[1]], style=mc_style))
    end
end

# plot avoid set again because it gets covered up
push!(fig, PGFPlots.Plots.Linear([-0.2617, -0.2617], [-0.9, 0.25], style="red, dashed", mark="none"))

fig.legendStyle =  "at={(1.05,1.0)}, anchor=north west"

PGFPlots.save("examples/jmlr/plots/single_pendulum/single_pendulum_$(controller_type)_x12.tex", fig)
PGFPlots.save("examples/jmlr/plots/single_pendulum/single_pendulum_$(controller_type)_x12.pdf", fig)