text:
# --------------------------------------------------------------------------- #
## General
# --------------------------------------------------------------------------- #

function to_string(model::Model, block::T) where {T <: AbstractBlock}
    error("Not implemented yet")
end

function create_string_vector(n::Int)
    return Vector{String}(undef, n)
end

struct StringBuildingBlocks
    eq::String
    ws::String
    sep::String

    function StringBuildingBlocks(;
        eq = " = ",
        ws = " ",
        sep = "_")

        return new(eq, ws, sep)
    end
end

const stringbuildingblocks = StringBuildingBlocks()

function find_block_name(model::Model, block::T) where {T <: AbstractBlock}
    block_index = model.map_node_to_index[block]
    block_name  = model.names[block_index]
    return block_name
end

function std_unary_outputpin_name()
    output_pin = "y"
    return output_pin
end

function std_unary_inputpin_name()
    input_pin = "x"
    return input_pin
end

function std_numbered_inputpin_name(n::Int)
    input_pin = "x" * "$(n)"
    return input_pin
end

function std_numbered_outputpin_name(n::Int)
    output_pin = "y" * "$(n)"
    return output_pin
end

function build_equation_assignment(lhs::String, rhs::String)
    str = lhs * stringbuildingblocks.eq * rhs
    return str
end

# --------------------------------------------------------------------------- #
## AbstractUnaryTransformer
# --------------------------------------------------------------------------- #
#=
Identity{T}  # OK
Gain{T} # OK
    g::T
    x::InputPin{T}
    y::OutputPin{T}

Sin{T} # OK
    x::InputPin{T}
    y::OutputPin{T}

UserDefinedUnaryFunction{T} # OK
    fcn
    x::InputPin{T}
    y::OutputPin{T}
=#

function to_string(model::Model, block::Identity)

    block_name = find_block_name(model, block)

    input_pin  = std_unary_inputpin_name()
    output_pin = std_unary_outputpin_name()

    eq_parameter_str = nothing

    lhs = block_name * stringbuildingblocks.sep * output_pin
    rhs = block_name * stringbuildingblocks.sep * input_pin

    eq_str    = create_string_vector(1)
    eq_str[1] = build_equation_assignment(lhs, rhs)

    return (eq_parameter_str, eq_str)
end

function to_string(model::Model, block::Gain)

    block_name = find_block_name(model, block)

    input_pin  = std_unary_inputpin_name()
    output_pin = std_unary_outputpin_name()

    eq_parameter_str_lhs_1 = block_name * stringbuildingblocks.sep * "g"
    eq_parameter_str_rhs_1 = "$(block.g)"

    eq_parameter_str    = create_string_vector(1)
    eq_parameter_str[1] = build_equation_assignment(eq_parameter_str_lhs_1, eq_parameter_str_rhs_1)

    lhs = block_name * stringbuildingblocks.sep * output_pin
    rhs = block_name * stringbuildingblocks.sep * input_pin * " * " * eq_parameter_str_lhs_1

    eq_str    = create_string_vector(1)
    eq_str[1] = build_equation_assignment(lhs, rhs)

    return (eq_parameter_str, eq_str)
end

function to_string(model::Model, block::Sin)

    block_name = find_block_name(model, block)

    input_pin  = std_unary_inputpin_name()
    output_pin = std_unary_outputpin_name()

    eq_parameter_str = nothing

    lhs = block_name * stringbuildingblocks.sep * output_pin
    rhs = "sin(" * block_name * stringbuildingblocks.sep * input_pin * ")"

    eq_str    = create_string_vector(1)
    eq_str[1] = build_equation_assignment(lhs, rhs)

    return (eq_parameter_str, eq_str)

end

function to_string(model::Model, block::UserDefinedUnaryFunction)

    block_name = find_block_name(model, block)
    block_index = model.map_node_to_index[block]

    eq_parameter_str = nothing

    input_pin  = std_unary_inputpin_name()
    output_pin = std_unary_outputpin_name()

    lhs = block_name * stringbuildingblocks.sep * output_pin

    arg = block_name * stringbuildingblocks.sep * input_pin
    rhs = "model.nodes[$(block_index)].fcn($(arg))"

    eq_str    = create_string_vector(1)
    eq_str[1] = build_equation_assignment(lhs, rhs)

    return (eq_parameter_str, eq_str)
end

# --------------------------------------------------------------------------- #
## AbstractBinaryTransformer
# --------------------------------------------------------------------------- #
#=
Plus
    x1::InputPin{T}
    x2::InputPin{T}
    y::OutputPin{T}

Minus

Multiply
Divide
=#

function to_string_binary_template(model::Model, block::T, operand) where {T <: AbstractBinaryTransformer}

    block_name = find_block_name(model, block)

    input_pin_1  = std_numbered_inputpin_name(1)
    input_pin_2  = std_numbered_inputpin_name(2)
    output_pin   = std_unary_outputpin_name()

    eq_parameter_str = nothing

    lhs = block_name * stringbuildingblocks.sep * output_pin

    x1 = block_name * stringbuildingblocks.sep * input_pin_1
    x2 = block_name * stringbuildingblocks.sep * input_pin_2
    rhs = x1 * " " * operand * " " * x2

    eq_str    = create_string_vector(1)
    eq_str[1] = build_equation_assignment(lhs, rhs)

    return (eq_parameter_str, eq_str)

end

function to_string(model::Model, block::Plus)
    return to_string_binary_template(model, block, "+")
end

function to_string(model::Model, block::Minus)
    return to_string_binary_template(model, block, "-")
end

function to_string(model::Model, block::Multiply)
    return to_string_binary_template(model, block, "*")
end

function to_string(model::Model, block::Divide)
    return to_string_binary_template(model, block, "/")
end

# --------------------------------------------------------------------------- #
## AbstractSource
# --------------------------------------------------------------------------- #
#=
Constant{T}  # OK
    value::T
    y::OutputPin{T}

SinGenerator{T} # OK
    amplitude::T
    frequency::T
    phase::T
    y::OutputPin{T}

Step{T} # OK
=#

function to_string(model::Model, block::Constant)

    block_name = find_block_name(model, block)

    output_pin = std_unary_outputpin_name()

    eq_parameter_str = nothing

    lhs = block_name * stringbuildingblocks.sep * output_pin
    rhs = "$(block.value)"

    eq_str    = create_string_vector(1)
    eq_str[1] = build_equation_assignment(lhs, rhs)

    return (eq_parameter_str, eq_str)
end

function to_string(model::Model, block::SinGenerator)
    block_name = find_block_name(model, block)

    output_pin = std_unary_outputpin_name()

    eq_parameter_str_lhs_1 = block_name * stringbuildingblocks.sep * "amplitude"
    eq_parameter_str_lhs_2 = block_name * stringbuildingblocks.sep * "frequency"
    eq_parameter_str_lhs_3 = block_name * stringbuildingblocks.sep * "phase"

    eq_parameter_str_rhs_1 = "$(block.amplitude)"
    eq_parameter_str_rhs_2 = "$(block.frequency)"
    eq_parameter_str_rhs_3 = "$(block.phase)"

    eq_parameter_str    = create_string_vector(3)
    eq_parameter_str[1] = build_equation_assignment(eq_parameter_str_lhs_1, eq_parameter_str_rhs_1)
    eq_parameter_str[2] = build_equation_assignment(eq_parameter_str_lhs_2, eq_parameter_str_rhs_2)
    eq_parameter_str[3] = build_equation_assignment(eq_parameter_str_lhs_3, eq_parameter_str_rhs_3)

    lhs = block_name * stringbuildingblocks.sep * output_pin
    rhs = "$(eq_parameter_str_lhs_1) * sin($(eq_parameter_str_lhs_2) * time + $(eq_parameter_str_lhs_3))"

    eq_str    = create_string_vector(1)
    eq_str[1] = build_equation_assignment(lhs, rhs)

    return (eq_parameter_str, eq_str)
end

function to_string(model::Model, block::Step)

    block_name = find_block_name(model, block)

    output_pin = std_unary_outputpin_name()

    eq_parameter_str_lhs_1 = block_name * stringbuildingblocks.sep * "amplitude"
    eq_parameter_str_lhs_2 = block_name * stringbuildingblocks.sep * "time_delay"
    eq_parameter_str_rhs_1 = "$(block.amplitude)"
    eq_parameter_str_rhs_2 = "$(block.time_delay)"

    eq_parameter_str    = create_string_vector(2)
    eq_parameter_str[1] = build_equation_assignment(eq_parameter_str_lhs_1, eq_parameter_str_rhs_1)
    eq_parameter_str[2] = build_equation_assignment(eq_parameter_str_lhs_2, eq_parameter_str_rhs_2)

    lhs = block_name * stringbuildingblocks.sep * output_pin
    rhs = "time >= $(eq_parameter_str_lhs_2) ? $(eq_parameter_str_lhs_1) : zero(typeof(time))"

    eq_str    = create_string_vector(1)
    eq_str[1] = build_equation_assignment(lhs, rhs)

    return (eq_parameter_str, eq_str)
end

function to_string(model::Model, block::Curve)

    block_name  = find_block_name(model, block)
    block_index = model.map_node_to_index[block]

    eq_parameter_str = nothing

    output_pin = std_unary_outputpin_name()

    lhs = block_name * stringbuildingblocks.sep * output_pin

    arg = "time"
    rhs = "eval(model.nodes[$(block_index)], $(arg))"

    eq_str    = create_string_vector(1)
    eq_str[1] = build_equation_assignment(lhs, rhs)

    return (eq_parameter_str, eq_str)
end

# --------------------------------------------------------------------------- #
## AbstractIntegrator
# --------------------------------------------------------------------------- #
#=
Integrator # OK
=#

function to_string(model::Model, block::Integrator)

    block_name          = find_block_name(model, block)
    integrator_indicies = find_integrators(model)
    block_index         = model.map_node_to_index[block]
    integrator_index    = findfirst(integrator_indicies .== block_index)

    output_pin = std_unary_outputpin_name()

    eq_parameter_str    = nothing

    lhs = block_name * stringbuildingblocks.sep * output_pin
    rhs = "q[" * "$(integrator_index)" * "]"

    eq_str    = create_string_vector(1)
    eq_str[1] = build_equation_assignment(lhs, rhs)

    return (eq_parameter_str, eq_str)

end

# --------------------------------------------------------------------------- #
## Model
# --------------------------------------------------------------------------- #

function model_to_string(
    model::Model,
    forest,
    integrator_indices,
    tree_node_vecvec,
    evalulation_indices_vecvec)

    # for each node, get the string represenation
    # sort the parameter equation at the front, the assignments at the end
    # or better into two separate lists
    # if nothing == no parameter in that block -> skip ode_solver
    #
    # for each equation == for each element in the forest
    # from the bottom up propagate each connections in a form of:
    #
    # receiver_block_name_inputpinname = sender_block_name_outpinname

    # ode function structure:
    # section 1: node parameter assignments
    # section 2: for each equation:
    # from bottom nodes (source or integrator) its inner represenation
    # assignment to its outputpin
    # assignement to its receiver pins
    #

    # create the strings for each node
    N_nodes = length(model.nodes)

    equation_param_vec      = Vector{Any}(undef, N_nodes)
    equation_assignment_vec = Vector{Any}(undef, N_nodes)

    for ii = 1:N_nodes
        (eq_parameter_str, eq_str)  = to_string(model, model.nodes[ii])
        equation_param_vec[ii]      = eq_parameter_str
        equation_assignment_vec[ii] = eq_str
    end

    # get rid of the nothing assignments
    for ii = N_nodes:-1:1
        if isnothing(equation_param_vec[ii])
            popat!(equation_param_vec, ii)
        end
    end

    # flatten the parameters
    L_eq_params = length(equation_param_vec)
    Ls          = map(length, equation_param_vec)
    N_param_eqs = sum(Ls)
    flat_equation_param_vec = create_string_vector(N_param_eqs)

    kk = 0
    for ii = 1:length(equation_param_vec)
        for jj = 1:Ls[ii]
            kk += 1
            flat_equation_param_vec[kk] = equation_param_vec[ii][jj]
        end
    end

    # for all trees:
    # bottom up create the connection strings
    n_trees        = length(forest)
    # connection_set = Set{String}()
    has_been_added = falses(N_nodes)

    eq_vecvec = datastructs.Stack{Any}(N_nodes) # preallocate with N_nodes

    for tt = 1:n_trees
        eval_indices = evalulation_indices_vecvec[tt]
        L = length(eval_indices)

        for ii = 1:L
            node_idx = eval_indices[ii]
            node     = model.nodes[node_idx]

            # add the node equation
            if has_been_added[node_idx]
                continue
            else
                has_been_added[node_idx] = true
                Leq = length(equation_assignment_vec[node_idx])
                for bb = 1:Leq
                    push!(eq_vecvec, equation_assignment_vec[node_idx][bb])
                end

                # add the outgoing connection equations
                rhs_block_name = model.names[node_idx]

                symbol_list = get_output_pin_names(node)
                for sym in symbol_list
                    out_pin = getfield(node, sym)

                    # handle outgoing signals
                    signal = out_pin.left_child_signal

                    while true
                        cont_obj       = signal.end_pin.containing_object
                        cont_obj_index = model.map_node_to_index[cont_obj]
                        lhs_name       = model.names[cont_obj_index]

                        # how to get the name of the input pin?
                        cont_obj_input_pin_symbol_list = get_input_pin_names(cont_obj)

                        # find the input pin
                        inputpin_name = ""
                        for aa = 1:length(cont_obj_input_pin_symbol_list)

                            if getfield(cont_obj, cont_obj_input_pin_symbol_list[aa]) == signal.end_pin
                                inputpin_name = string(cont_obj_input_pin_symbol_list[aa])
                                break
                            end
                        end
                        if isempty(inputpin_name)
                            error("Input pin could not be found")
                        end

                        conn = lhs_name * "_" * inputpin_name * " = " * rhs_block_name * "_" * string(sym)
                        push!(eq_vecvec, conn)

                        if has_right_sibling(signal)
                            signal = signal.right_sibling_signal
                        else
                            break
                        end
                    end # while outgoing signals
                end # for output pins
            end # if has been added
        end # for nodes in one tree

        # handle last assignemt to der_q[idx] = rhs
        node_idx             = eval_indices[end]
        integrator_index     = findfirst(node_idx .== integrator_indices)
        integrator_node_name = model.names[node_idx]

        rhs_str = integrator_node_name * "_x"
        der_eq = "der_q[" * "$(integrator_index)" * "] = " * rhs_str
        push!(eq_vecvec, der_eq)
    end # for trees

    # flatten eq_vecvec
    # eq_vecvec_data = eq_vecvec.data[1:eq_vecvec.ptr]
    # L_eq_vecvec = length(eq_vecvec_data)
    # Ls          = map(length, eq_vecvec_data)
    # N_param_eqs = sum(Ls)
    # flat_equation_assignments = create_string_vector(N_param_eqs)

    # kk = 0
    # for ii = 1:L_eq_vecvec
    #     for jj = 1:Ls[ii]
    #         kk += 1
    #         flat_equation_assignments[kk] = eq_vecvec_data[ii][jj]
    #     end
    # end

    flat_equation_assignments = eq_vecvec.data[1:eq_vecvec.ptr]

    # create function signature
    newline   = "\n"
    str_c     = ", "
    str_der_q = "der_q"
    str_time  = "time"
    str_q     = "q"
    str_model = "model"
    fcn_name  = "function generated_inner_ode_fcn"
    fcn_signature = fcn_name * "(" * str_der_q * str_c * str_time * str_c * str_q * str_c * str_model * ")"
    fcn_end = "return der_q\nend"

    # create function body
    fcn_body = [flat_equation_param_vec; flat_equation_assignments]
    fcn_body = map(x -> x * "\n", fcn_body)

    fcn_str = [fcn_signature; newline; newline; fcn_body; newline; newline; fcn_end]
    fcn_total_str = join(fcn_str)

    return fcn_total_str # (equation_param_vec, equation_assignment_vec, fcn_str)
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
#Defining a generator does not make a lot of sense, since for a chosen board_size, there is only one possible instance.
#We define it to se how much our agent can learn on a single instance.

struct NQueensGenerator <: AbstractModelGenerator
    board_size::Int
    #density::Real
    function NQueensGenerator(board_size)
        @assert board_size > 1
        new(board_size)
    end
end

"""
    fill_with_generator!(cpmodel::CPModel, gen::GraphColoringGenerator)
Fill a CPModel with the variables and constraints generated. We fill it directly instead of
creating temporary files for efficiency purpose.

This generator create graps for the NQueens problem.

"""
function fill_with_generator!(cpmodel::CPModel, gen::NQueensGenerator; seed=nothing)
    cpmodel.limit.numberOfSolutions = 1
    if !isnothing(seed)
        Random.seed!(seed)
    end

    #density = gen.density
    board_size = gen.board_size

    #nb_edges = floor(Int64, density * nb_nodes)

    rows = Vector{SeaPearl.AbstractIntVar}(undef, board_size)
    for i = 1:board_size
        rows[i] = SeaPearl.IntVar(1, board_size, "row_"*string(i), cpmodel.trailer)
        SeaPearl.addVariable!(cpmodel, rows[i]; branchable=true)
    end

    rows_plus = Vector{SeaPearl.AbstractIntVar}(undef, board_size)
    for i = 1:board_size
        rows_plus[i] = SeaPearl.IntVarViewOffset(rows[i], i, rows[i].id*"+"*string(i))
        #SeaPearl.addVariable!(model, rows_plus[i]; branchable=false)
    end

    rows_minus = Vector{SeaPearl.AbstractIntVar}(undef, board_size)
    for i = 1:board_size
        rows_minus[i] = SeaPearl.IntVarViewOffset(rows[i], -i, rows[i].id*"-"*string(i))
        #SeaPearl.addVariable!(model, rows_minus[i]; branchable=false)
    end

    SeaPearl.addConstraint!(cpmodel, SeaPearl.AllDifferent(rows, cpmodel.trailer))
    SeaPearl.addConstraint!(cpmodel, SeaPearl.AllDifferent(rows_plus, cpmodel.trailer))
    SeaPearl.addConstraint!(cpmodel, SeaPearl.AllDifferent(rows_minus, cpmodel.trailer))
    return nothing
    #return model
end

ext: jl


text:
using Documenter, MathOptInterface

makedocs(
    sitename = "MathOptInterface",
    # See https://github.com/JuliaDocs/Documenter.jl/issues/868
    html_prettyurls = get(ENV, "CI", nothing) == "true",
    # See https://github.com/JuliaOpt/JuMP.jl/issues/1576
    strict = true,
    pages = [
        "Introduction" => "index.md",
        "Manual" => "apimanual.md",
        "Reference" => "apireference.md",
    ]
)

deploydocs(
    repo   = "github.com/JuliaOpt/MathOptInterface.jl.git",
)


text:
# Aliases: I = 1, V = 5, X = 10, L = 50, C = 100, D = 500, M = 1000.
function to_roman(number::Integer)
    number < 1 && error("Invalid number")
    alias = [
        ["I", "II", "III", "IV", "V", "VI", "VII", "VIII", "IX"],
        ["X", "XX", "XXX", "XL", "L", "LX", "LXX", "LXXX", "XC"],
        ["C", "CC", "CCC", "CD", "D", "DC", "DCC", "DCCC", "CM"],
        ["M", "MM", "MMM"],
    ]
    i = 1
    result = []
    while number > 0
        j = rem(number, 10)
        if j > 0
            push!(result, alias[i][j])
        end
        i += 1
        number = div(number, 10)
    end
    join(reverse(result))
end


text:
@testset "$(rpad("Stochastic equations",80))" begin

    ################################################################################
    # Test SDE: Stochastic Differential Equation
    ################################################################################

    x₁ₛ = [0.5  0.0 -0.5; 0.0  0.5  0.0]

    function sde_v(λ, t, q, v)
        v[1] = λ*q[1]
        v[2] = λ*q[2]
    end

    function sde_u(μ, t, q, u)
        u[1] = μ*q[1]
        u[2] = μ*q[2]
    end

    λ  = 2.
    μ  = 1.

    sde_v_params = (t, q, v) -> sde_v(λ, t, q, v)
    sde_u_params = (t, q, v) -> sde_u(μ, t, q, v)

    sde  = SDE(1, 1, sde_v_params, sde_u_params, t₀, x₀)
    sde1 = SDE(1, 1, sde_v_params, sde_u_params, x₀)
    sde2 = SDE(1, 3, sde_v_params, sde_u_params, x₀)
    sde3 = SDE(1, 1, sde_v_params, sde_u_params, x₁ₛ)

    @test sde == sde1
    @test sde != sde2
    @test sde != sde3

    @test hash(sde) == hash(sde1)
    @test hash(sde) != hash(sde2)
    @test hash(sde) != hash(sde3)

    @test sde1.d == 2
    @test sde2.d == 2
    @test sde3.d == 2

    @test sde1.ni == 1
    @test sde2.ni == 1
    @test sde3.ni == 3

    @test sde1.ns == 1
    @test sde2.ns == 3
    @test sde3.ns == 1

    @test sde == similar(sde, t₀, x₀, 1)
    @test sde == similar(sde, t₀, x₀)
    @test sde == similar(sde, x₀, 1)
    @test sde == similar(sde, x₀)

    ################################################################################
    # Test PSDE: Partitioned Stochastic Differential Equation
    ################################################################################

    noise_intensity = 0.1

    q₁ₛ = [0.5  0.0 -0.5]
    p₁ₛ = [0.0  0.5  0.0]

    function psde_v(t, q, p, v)
        v[1] = p[1]
    end

    function psde_f(t, q, p, f)
        f[1] = - q[1]
    end

    function psde_B(t, q, p, B)
        B[1,1] = noise_intensity * p[1]
    end

    function psde_G(t, q, p, G)
        G[1,1] = - noise_intensity * q[1]
    end

    psde  = PSDE(1, 1, psde_v, psde_f, psde_B, psde_G, t₀, q₀, p₀)
    psde1 = PSDE(1, 1, psde_v, psde_f, psde_B, psde_G, q₀, p₀)
    psde2 = PSDE(1, 3, psde_v, psde_f, psde_B, psde_G, q₀, p₀)
    psde3 = PSDE(1, 1, psde_v, psde_f, psde_B, psde_G, q₁ₛ, p₁ₛ)

    @test psde == psde1
    @test psde != psde2
    @test psde != psde3

    @test hash(psde) == hash(psde1)
    @test hash(psde) != hash(psde2)
    @test hash(psde) != hash(psde3)

    @test psde1.d == 1
    @test psde2.d == 1
    @test psde3.d == 1

    @test psde1.ni == 1
    @test psde2.ni == 1
    @test psde3.ni == 3

    @test psde1.ns == 1
    @test psde2.ns == 3
    @test psde3.ns == 1

    @test psde == similar(psde, t₀, q₀, p₀, 1)
    @test psde == similar(psde, q₀, p₀, 1)

    @test psde != similar(psde, t₀, q₁ₛ, p₁ₛ)
    @test psde != similar(psde, q₁ₛ, p₁ₛ)

    ################################################################################
    # Test SPSDE: Split Partitioned Stochastic Differential Equation
    ################################################################################

    noise_intensity = 0.1

    q₁ₛ = [0.5  0.0 -0.5]
    p₁ₛ = [0.0  0.5  0.0]

    function spsde_v(t, q, p, v)
        v[1] = p[1]
    end

    function spsde_f1(t, q, p, f)
        f[1] = - q[1]
    end

    function spsde_f2(t, q, p, f)
        f[1] = 0
    end

    function spsde_B(t, q, p, B)
        B[1,1] = noise_intensity * p[1]
    end

    function spsde_G1(t, q, p, G)
        G[1,1] = - noise_intensity * q[1]
    end

    function spsde_G2(t, q, p, G)
        G[1,1] = 0
    end

    spsde  = SPSDE(1, 1, spsde_v, spsde_f1, spsde_f2, spsde_B, spsde_G1, spsde_G2, t₀, q₀, p₀)
    spsde1 = SPSDE(1, 1, spsde_v, spsde_f1, spsde_f2, spsde_B, spsde_G1, spsde_G2, q₀, p₀)
    spsde2 = SPSDE(1, 3, spsde_v, spsde_f1, spsde_f2, spsde_B, spsde_G1, spsde_G2, q₀, p₀)
    spsde3 = SPSDE(1, 1, spsde_v, spsde_f1, spsde_f2, spsde_B, spsde_G1, spsde_G2, q₁ₛ, p₁ₛ)

    @test spsde == spsde1
    @test spsde != spsde2
    @test spsde != spsde3

    @test hash(spsde) == hash(spsde1)
    @test hash(spsde) != hash(spsde2)
    @test hash(spsde) != hash(spsde3)

    @test spsde1.d == 1
    @test spsde2.d == 1
    @test spsde3.d == 1

    @test spsde1.ni == 1
    @test spsde2.ni == 1
    @test spsde3.ni == 3

    @test spsde1.ns == 1
    @test spsde2.ns == 3
    @test spsde3.ns == 1

    @test spsde == similar(spsde, t₀, q₀, p₀, 1)
    @test spsde == similar(spsde, q₀, p₀, 1)

    @test spsde != similar(spsde, t₀, q₁ₛ, p₁ₛ)
    @test spsde != similar(spsde, q₁ₛ, p₁ₛ)

end


text:
module DataReader
  export readSequences, FastaRecord, readMatrix, ScoreMatrix

  import Base.==

  immutable FastaRecord
    description :: ASCIIString
    sequence :: ASCIIString

  end

  ==(a::FastaRecord, b::FastaRecord) = a.description == b.description && a.sequence == b.sequence

  immutable ScoreMatrix
    keys :: Dict{Char, Int}
    hsh  :: Array{Array{Float64, 1}, 1}
  end
  ScoreMatrix() = ScoreMatrix(Dict{Char, Int}(), [])

  function readSequences(input_file_name :: AbstractString)
    records = FastaRecord[]
    temp_desc = ""
    temp_str = ""
    input_file = open(input_file_name,"r")
    while !eof(input_file)
      s = rstrip(readline(input_file), ['\r','\n'])
      if s[1] == '>'
        if length(temp_str) > 0
          push!(records, FastaRecord(temp_desc, temp_str))
          temp_str = ""
        end
        temp_desc = s[2:end]
      else
        temp_str = string(temp_str, s)
      end
    end

    if (temp_str != "")
      push!(records, FastaRecord(temp_desc, temp_str))
    end
    records
  end

  #helper function to avoid problems with readdlm from standard lib (empty lines skipping, \r\n, etc.)
  function my_readdlm(input_file :: AbstractString)
    removeEol = str-> strip(str, ['\r', '\n'])
    notCommentOrEmptyLine = x -> !startswith(x, "#") && !isempty(x)
    getNumbersAndChars = n -> ismatch(r"[^A-Z*]", n) ? float(n) : n[1]

    linesProcessor = lines -> map( x-> map(getNumbersAndChars, split(x)),
                              filter(notCommentOrEmptyLine, map(removeEol, lines)))

    f = open(input_file)
    lines = linesProcessor(readlines(f))
    close(f)
    lines
  end

  function readMatrix(input_file :: AbstractString)
    result = my_readdlm(input_file)
    #result = readdlm(input_file)
    # thus dirty fix needed to avoid error on empty lines in input file. in julia 0.3 readdlm doesn't skip them
    k = [ getindex(getindex(result, 1), inner_key)[1] => inner_key for inner_key = 1:24]
    h = [
        [
            getindex(getindex(result, outer_key), inner_key + 1) :: Float64
            for inner_key = 1 : outer_key - 1
        ]
        for outer_key = 2 : 25
    ]
    ScoreMatrix(k, h)
  end

end


text:
using DataFrames, Statistics, StatsBase

function analyze_results_stats(doc)
    df = readtable(doc, header=true, separator=' ')
    categorical!(df, :indicator)
    categorical!(df, :criteria)

    df_ind_dist = filter(r -> any(occursin.(["dist"], r.indicator)), df)
    df_ind_num = filter(r -> any(occursin.(["number"], r.indicator)), df)
    df_crit_dist = filter(r -> any(occursin.(["dist"], r.criteria)), df)
    df_crit_num = filter(r -> any(occursin.(["number"], r.criteria)), df)

    df_ind_crit_dist = filter(r -> any(occursin.(["dist"], r.indicator)), df_crit_dist)
    df_ind_crit_num = filter(r -> any(occursin.(["number"], r.indicator)), df_crit_num)

    results = Array{Array{Any, 1},1}()
    stats = summarystats(df_ind_dist[:makespan])
    t_ind_dist = [stats.mean, stats.min, stats.q25, stats.median, stats.q75, stats.max, "ind_dist"]
    push!(results, t_ind_dist)

    stats = summarystats(df_ind_num[:makespan])
    t_ind_num = [stats.mean, stats.min, stats.q25, stats.median, stats.q75, stats.max, "ind_num"]
    push!(results, t_ind_num)

    stats = summarystats(df_crit_dist[:makespan])
    t_crit_dist = [stats.mean, stats.min, stats.q25, stats.median, stats.q75, stats.max, "crit_dist"]
    push!(results, t_crit_dist)

    stats = summarystats(df_crit_num[:makespan])
    t_crit_num = [stats.mean, stats.min, stats.q25, stats.median, stats.q75, stats.max, "crit_num"]
    push!(results, t_crit_num)

    stats = summarystats(df_ind_crit_dist[:makespan])
    t_ind_crit_dist = [stats.mean, stats.min, stats.q25, stats.median, stats.q75, stats.max, "ind_crit_dist"]
    push!(results, t_ind_crit_dist)

    stats = summarystats(df_ind_crit_num[:makespan])
    t_ind_crit_num = [stats.mean, stats.min, stats.q25, stats.median, stats.q75, stats.max, "ind_crit_num"]
    push!(results, t_ind_crit_num)

    return(results)
end

function init_results_stats_file(output::String)

    open(output, "w") do o
        println(o, "Mean Minimum 1stQuartile Median 3rdQuartile Maximum Type")
    end
end

function write_results_stats_file(output::String, results::Array{Array{Any,1},1})

    open(output, "a+") do o
        for t in results
            for s in t
                print(o, string(s)*" ")
            end
            print(o, "\n")
        end
    end
end

#RESULTS MAIN PROGRAM
for file in reverse(readdir("./results/"))
    println("-------------")
    println(file)
    if file != ".ipynb_checkpoints" && file != "60C_stats.txt" && file != "240C_stats.txt" && file != "500C_stats.txt" && file != "1000C_stats.txt"
    # if file == "60C_25Type_UniformDense_2QC.txt"
        doc = "./results/"*file
        results = analyze_results_stats(doc)

        if file[1:2] == "60"
            out_doc = "60C_stats.txt"
        elseif file[1:3] == "240"
            out_doc = "240C_stats.txt"
        elseif file[1:3] == "500"
            out_doc = "500C_stats.txt"
        else
            out_doc = "1000C_stats.txt"
        end

        if !(out_doc in readdir("./results/"))
            init_results_stats_file("./results/"*out_doc)
        end

        write_results_stats_file("./results/"*out_doc, results)

    end
end

max_stars_repo_licenses: Zlib MIT

max_issues_repo_licenses: Zlib MIT

max_forks_repo_licenses: Zlib MIT


text:
# This file is part of Kpax3. License is MIT.

"""
Logposterior minus the log(normalizing constant)
"""
function logposterior(R::Vector{Int},
                      C::Matrix{UInt8},
                      priorR::PriorRowPartition,
                      priorC::PriorColPartition,
                      k::Int)
  # TODO: write the function
  nothing
end


text:
using AdventOfCodeSolutions
using Test
using OffsetArrays

function input(puzzle::Puzzle{2019, 2, n}) where n
    io = openInput(puzzle)
    inputs = split(read(io, String), ',', keepempty=false)
    return parse.(Int, inputs)
end

function runIntCode(memory)
    location = 0

    while (true)
        opcode = memory[location]

        if (opcode == 1)
            (left, right, output) = memory[location + 1:location + 3]
            memory[output] = sum(memory[[left, right]])
        end
        if (opcode == 2)
            (left, right, output) = memory[location + 1:location + 3]
            memory[output] = prod(memory[[left, right]])
        end

        if (opcode == 99)
            break
        end
        location += 4
    end

    return memory
end

function runWithInput(memory, noun, verb)
    memory[1] = noun
    memory[2] = verb
    runIntCode(memory)
end

function solve(::Puzzle{2019, 2, 1}, input)
    runWithInput(OffsetVector(input, -1), 12, 2)[0]
end

function solve(::Puzzle{2019, 2, 2}, input)
    (noun, verb) = Iterators.filter(
        ((noun, verb),) -> runWithInput(OffsetVector(copy(input), -1), noun, verb)[0] == 19690720,
        Iterators.product(0:99, 0:99)
    ) |>first

    return 100 * noun + verb
end

@testset "Day 2" begin
    @test runIntCode(OffsetVector([1,0,0,0,99], -1)) == OffsetVector([2,0,0,0,99], -1)
    @test runIntCode(OffsetVector([2,3,0,3,99], -1)) == OffsetVector([2,3,0,6,99], -1)
    @test runIntCode(OffsetVector([2,4,4,5,99,0], -1)) == OffsetVector([2,4,4,5,99,9801], -1)
    @test runIntCode(OffsetVector([1,1,1,4,99,5,6,0,99], -1)) == OffsetVector([30,1,1,4,2,5,6,0,99], -1)

    @test (solve(Puzzle(2019, 2, 1), input(Puzzle(2019, 2, 1))) == 5866714)
    @test (solve(Puzzle(2019, 2, 2), input(Puzzle(2019, 2, 2))) == 5208)
end


text:
# Test Hamiltonian increment
h = ham(1e-2, 0*q('z'))
@test h(0.0, q('x')) ≈ q('x')
h = ham(1e-2, t->0*q('z'))
@test h(0.0, q('x')) ≈ q('x')
h = ham(π/2, q('x'))
@test h(0.0, ground(q)) ≈ q(1,1)
h = ham(1e-2, 0*q('z'), ket=true)
@test h(0.0, xplus) ≈ xplus
h = ham(π/2, q('x'), ket=true)
@test h(0.0, groundvec(q)) ≈ -im*q(1)

# Test superoperator Hamiltonian increment
h = sham(1e-2, 0*q('z'))
@test h(0.0, superket(q('x'))) ≈ superket(q('x'))
h = sham(1e-2, t->0*q('z'))
@test h(0.0, superket(q('x'))) ≈ superket(q('x'))
h = sham(π/2, q('x'))
@test h(0.0, superket(ground(q))) ≈ superket(q(1,1))

# Test RK4 Hamiltonian increment
h = ham_rk4(1e-2, 0*q('z'))
@test h(0.0, q('x')) ≈ q('x')
h = ham_rk4(1e-2, t->0*q('z'))
@test h(0.0, q('x')) ≈ q('x')
h = ham_rk4(1e-2, 0*q('z'), ket=true)
@test h(0.0, xplus) ≈ xplus

# Test jump/no-jump Lindblad increment
l = lind(1e-6, 0*q('z'), clist=[q('d')])
@test l(0.0, q(1,1))[2,2] ≈ exp(-1e-6)
l = lind(1e-6, clist=[q('d')])
@test l(0.0, q(1,1))[2,2] ≈ exp(-1e-6)

# Test RK4 Lindblad increment
l = lind_rk4(1e-6, 0*q('z'), clist=[q('d')])
@test l(0.0, q(1,1))[2,2] ≈ exp(-1e-6)
l = lind_rk4(1e-6, clist=[q('d')])
@test l(0.0, q(1,1))[2,2] ≈ exp(-1e-6)

# Test Superoperator Lindblad increment
l = slind(1e-6, 0*q('z'), clist=[q('d')])
@test unsuperket(l(0.0, superket(q(1,1))))[2,2] ≈ exp(-1e-6)
l = slind(1e-6, t->0*q('z'), clist=[q('d')])
@test unsuperket(l(0.0, superket(q(1,1))))[2,2] ≈ exp(-1e-6)

# Test trajectory code
t = trajectory(ham(0.0, q('z')), ground(q), (0.0, 1.0), ρ->ρ[1,1], dt=1e-3, points=100, verbose=false)
@test t[1].t == range(0.0,stop=1.0,length=100)
@test t[1][end] ≈ QComp(1)
@test t[1](1.0) ≈ QComp(1)
@test size(t[1]) == size(t[1].v)
t = trajectory(ham(1e-4, (2π/4).*q('y')), ground(q), (0.0, 1.0), ρ->ρ[2,2], dt=1e-4, points=1000, verbose=false)
@test t[1].t == range(0.0,stop=1.0,length=1000)
@test t[1][end] ≈ QComp(1) atol=1e-4
@test t[1](1.0) ≈ QComp(1) atol=1e-4

# Test trajectory point truncation
t = trajectory(ham(0.0, q('z')), ground(q), (0.0, 1.0), ρ->ρ[1,1], dt=1e-2, points=1000, verbose=false)
@test t[1].t == range(0.0,stop=1.0,length=101)
@test t[1][end] ≈ QComp(1)

# Test stochastic trajectory code
t = trajectory(meas(1e-3, q('z'), clist=[q('d')], mclist=[(q('z'), 2.0, 1.0)]), ground(q), (0.0, 1.0), ρ->ρ[1,1], dt=1e-3, points=100, verbose=false)
@test t[1].t == range(0.0,stop=1.0,length=100)
@test length(t) == 2
t = trajectory(meas(1e-3, q('z'), flist=[t->q('d')], mflist=[(t->q('z'), 2.0, 1.0)]), ground(q), (0.0, 1.0), ρ->ρ[1,1], dt=1e-3, points=100, verbose=false)
@test t[1].t == range(0.0,stop=1.0,length=100)
@test length(t) == 2

# Test stochastic ensemble code
e = Ensemble(t[1])
@test e(1).t == range(0.0,stop=1.0,length=100)
@test size(e(1).v) == size(t[1].v)
@test size(e) == size(e.a)
@test length(e) == length(e.a)
@test ndims(e) == ndims(e.a)
e = Ensemble(t[1].t, 1, eltype(t[1].v))
#e[1] = t[1].v
#@test e[1][1] == t[1].v[1]
e = Ensemble(t[1].t, t[1].v)
#@test e[1][1] == t[1].v[1]
e = ensemble(2, meas(1e-3, q('z'), mclist=[(q('z'), 2.0, 1.0)], clist=[q('d')]), ground(q), (0.0, 1.0), ρ->real(ρ[1,1]), dt=1e-3, points=100, verbose=false)
#@test e[1].n == 2
@test length(mean(e[1])) == length(e[1](1))
@test length(std(e[1])) == length(e[1](1))
@test length(median(e[1])) == length(e[1](1))


text:
using SafeTestsets
using ABM4OSN

@safetestset "ABM4OSN Tests" begin
    @safetestset "Network Tests" begin
        include(joinpath("test", "network_tests.jl"))
    end
    @safetestset "Property Function Tests" begin
        include(joinpath("test", "property_function_tests.jl"))
    end
    @safetestset "Post Tests" begin
        include(joinpath("test", "post_tests.jl"))
    end
end

# using ABM4OSN
# using Test
# using JLD
# using JLD2
# using CSV
# using LightGraphs

# batch = Config[]
# for ac in 50:50:200
#     push!(batch,Config(network=cfg_net(agent_count=ac)))
# end

# batch

# simulate_batch(batch, batch_desc="Agentcount")

# convert_results()
# result[1]

# result = load("results/_tmpstate.jld2")[""]

# using Statistics
# using LightGraphs
# mean(agent.desired_input_count for agent in result[3][2])

# std(indegree(result[3][1]))

# ABM4OSN.generate_desired_input_count(100,25)

# for i in 1:1000
#     push!(randenn, randn)
# end

# using Plots
# plot(randenn)
# histogram(randenn)

# sim2 = load("sim500.jld2")
# net_evolution = deepcopy(sim2["sim500"][1][3])

# CSV.write("dataframe.csv", sim2["sim500"][1][1])

# bla = "bla"

# using Plots
# plot(outdegrees, alpha=0.2)
# outdegrees
# arrtest = []
# last(arrtest)

# network = ABM4OSN.create_network(100,10)
# agents = ABM4OSN.create_agents(network)
# ABM4OSN.update_network!((network,agents),Config())
# agents
# agents[12].active = false
# agents
# active_agents = [agent.id for agent in agents if !agent.active]
# find(agents, active_agents)
# testset = indexin([agent for agent in agents if !agent.active],agents)

# testset = [indexin([1.2, 1.3], [1,1.1,1.2,1.3])][1]
# testset2 = Int64[]
# for elem in testset[1]
#     push!(testset2, elem)
# end

# testset2

# for elem in testset[1]
# indexin(1.2, [1,1.1,1.2,1.3])
# indexin

# using LightGraphs
# rem_vertices!(network, testset2)

# @time sim2 = simulate(Config(network = cfg_net(agent_count = 100), simulation=cfg_sim(n_iter=200)));

# sim[2]

# mean([weight for weight in sim[1][2].Weight if weight>=0])
# sim[1][1][:, [1,2,5,6,7]]

# sim
# save("sim500.jld2", "sim500", sim)

# testing = load("sim500.jld2")
# testing["sim500"]

# length(sim[3][2])

# @testset "ABM4OSN.jl" begin
#     # Write your own tests here.
#     @test typeof(Config()) == Config
#     @test Config(network=cfg_net(agent_count=100)).network.agent_count == 100
#     @test length(simulate()[3][2]) == 100
# end

# maximum([1,2,3])


text:
export logenergy, ℓ

"""
`logenergy(ω)`

Unnormalized joint log density
"""
function logenergy(ω::AbstractΩ)
  reduce(ω.data; init = 0.0) do logpdf_, (id, (dist, val))
    logpdf_ + logpdf(dist, val)
  end
end

const ℓ = logenergy

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
using GLPlot, GLAbstraction, ModernGL

window = createdisplay(w=1000,h=1000,eyeposition=Vec3(1.,1.,1.), lookat=Vec3(0.,0.,0.));
color  = Vec4(1,0,0,1)
n=100
h=1./n
r=h:h:1.
t=(-1:h:1+h)*π
x=float32(r*cos(t)')
y=float32(r*sin(t)')

f(x,y)  = exp(-10x.^2-20y.^2)  # arbitrary function of f
z       = Float32[float32(f(x[k,j],y[k,j])) for k=1:size(x,1),j=1:size(x,2)]
obj     = glplot(z, xrange=x, yrange=y, color="xyz.z>0 ? vec4(.1,.1,0.5+3*xyz.z, 1.0) : vec4(.1,.1-3*xyz.z,0.5,  1.0);")

renderloop(window)


text:
struct DynamicOracle{T} <: AbstractOracle{T}
    system::T
    oracle_function
    arc
end

"""
    DynamicOracle(system, oracle_function = dynamic_oracle; arc = untyped)

Create a dynamic oracle for dependency parsing (see [Goldberg & Nivre, 2012](https://aclweb.org/anthology/C/C12/C12-1059.pdf)).

Dynamic oracles allow multiple gold transitions from each parser configuration.
"""
DynamicOracle(system, oracle_function = dynamic_oracle; arc = untyped) =
    DynamicOracle(system, oracle_function, arc)

system(oracle::DynamicOracle) = oracle.system

(oracle::DynamicOracle)(tree::DependencyTree; kwargs...) =
    TreeOracle(oracle, tree; kwargs...)

function oracle_state(gold::TreeOracle{<:DynamicOracle}, cfg)
    arc, oracle = gold.oracle.arc, gold.oracle.oracle_function
    A = possible_transitions(cfg, gold.tree, arc)
    isoptimal = t -> oracle(t, cfg, gold.tree)
    G = filter(isoptimal, A)
    return OracleState(cfg, A, G)
end

function gold_transitions(oracle::DynamicOracle, cfg, gold::DependencyTree)
    isoptimal(t) = oracle.oracle_function(t, cfg, gold)
    A = possible_transitions(cfg, gold, oracle.arc)
    return filter(isoptimal, A)
end


text:
__precompile__()
module DiffWrappers

using ContinuousTransformations: TransformLogLikelihood
using DiffResults: GradientResult, DiffResult, ImmutableDiffResult
using DocStringExtensions: SIGNATURES
using ForwardDiff: gradient!, GradientConfig, Chunk
using Parameters: @unpack
import Base: length

export ForwardGradientWrapper

# forward

struct ForwardGradientWrapper{Tf, Tc, Tg}
    f::Tf
    config::Tc
    gr::Tg
end

"""
    $SIGNATURES

Create a wrapper for ``f: ℝⁿ→ℝ`` that returns the value `f(z)` and the gradient
``∇f(z)`` when called with vector `z`, using `ForwardDiff`. `x` is used only for
ascertaining the size for the buffer, its value is ignored.

Additional arguments are passed to `ForwardDiff.GradientConfig`.

Values returned by the wrapper don't share structure.
"""
function ForwardGradientWrapper(f, x; chunk = Chunk(x))
    config = GradientConfig(f, x, chunk)
    gr = GradientResult(x)
    ForwardGradientWrapper(f, config, gr)
end

ForwardGradientWrapper(t::TransformLogLikelihood; args...) =
    ForwardGradientWrapper(t, zeros(length(t)); args...)

"""
    $SIGNATURES

Helper function to copy shared structure if necessary.
"""
_ensure_unshared(gr::ImmutableDiffResult) = gr

# FIXME: revisit once https://github.com/JuliaDiff/DiffResults.jl/issues/4
# is fixed
_ensure_unshared(gr::DiffResult) = deepcopy(gr)

function (fgw::ForwardGradientWrapper)(x)
    @unpack f, config, gr = fgw
    _ensure_unshared(gradient!(gr, f, x, config))
end

length(fgw::ForwardGradientWrapper) = length(fgw.gr.derivs[1])

end # module


text:
using BenchmarkTools

data11 = """
8548335644
6576521782
1223677762
1284713113
6125654778
6435726842
5664175556
1445736556
2248473568
6451473526
"""

function day11()
    part = [0, 0]
    day11lines = filter(!isempty, strip.(split(data11, "\n")))
    nrows, ncols = length(day11lines), length(first(day11lines))
    mat = zeros(Int, nrows, ncols)
    for i in 1:nrows
        mat[i, :] .= parse.(Int, collect(day11lines[i]))
    end

    for step in 1:1000
        mat .+= 1
        while any(mat .> 9)
            for c in CartesianIndices(mat)
                if mat[c] > 9
                    for k in max(1, c[1] - 1):min(c[1] + 1, nrows),
                        l in max(c[2] - 1, 1):min(c[2] + 1, ncols)

                        if mat[k, l] > 0
                            mat[k, l] += 1
                        end
                    end
                    mat[c] = 0
                    step < 101 && (part[1] += 1)
                end
            end
        end
        if all(==(first(mat)), mat)
            part[2] == 0 && (part[2] = step)
            step >= 100 && break
        end
    end
    return part
end

part = day11()
println("Part 1:", part[1])
println("Part 2:", part[2])

@btime day11()

#=
Part 1:1601
Part 2:368
  7.539 ms (86240 allocations: 2.71 MiB)
=#


text:
using OrdinaryDiffEq, Zygote
using DiffEqSensitivity
using Random

function dt!(du, u, p, t)
  x, y = u
  α, β, δ, γ = p
  du[1] = dx = α*x - β*x*y
  du[2] = dy = -δ*y + γ*x*y
end

n_par=3
Random.seed!(2)
u0=rand(2,n_par)
u0[:,1] = [1.0,1.0]
tspan = (0.0, 10.0)
p = [2.2, 1.0, 2.0, 0.4]
prob_ode = ODEProblem(dt!, u0[:,1], tspan)

function test_loss(p1,prob)

	function prob_func(prob, i, repeat)
		@show i
		remake(prob,u0=u0[:,i])
	end

	#define ensemble problem
	ensembleprob = EnsembleProblem(prob,prob_func = prob_func)

	u = Array(solve(ensembleprob, Tsit5(), EnsembleThreads(), trajectories=n_par,
	p=p,
	sensealg = ForwardDiffSensitivity(),
	saveat = 0.1, dt=0.001))[:,end,:]
	loss=sum(u)
	return loss
end

@time gs = Zygote.gradient(p) do p
	test_loss(p,prob_ode)
end


text:
#promote_rule(::Type{S}, ::Type{T}) where {S<:AbstractIrrational,T<:SNumber} =
#    promote_type(promote_type(S, real(T)), T)


text:
using DrWatson
quickactivate(@__DIR__)
@show projectdir()
##
using CSV
using DataFrames
using Dates
using Plots
using BSON
using Underscores
import Downloads
using ColorSchemes
using RData
using CovidSurvey
using StatsBase
colors = ColorSchemes.tableau_10.colors
# using Makie
# using GLMakie
# GLMakie.activate!()
plotlyjs()
##

function parsetime(sraw)
    svec = split(sraw, "-")
    if length(svec) == 3
        d,m,y = svec
        sclean = join([d,m,"20$y"], "-")
        return Date(sclean, "d-u-y")
    else
        d = sraw[1:2]
        m = sraw[3:5]
        y = sraw[6:9]
        sclean = join([d,m,y], "-")
        return Date(sclean, "d-u-y")
    end
end

function readcontacts( fname )
    @_ fname |>
        CSV.read(__, DataFrame) |>
        DataFrames.select(__,
            :dataRound_glb => :batch,
            :edate => ByRow(x->parsetime(x)) => :date,
            "country_string" => :region,
            "qNew_3_2_resp" => :family,
            "qNew_3_3_resp" => :colleagues,
            "qNew_3_4_resp" => :friends,
            "qNew_3_5_resp" => :strangers,
        ) |>
    sort(__, :date)
end

function plot_rt(fname)
    p = plot(title="Reproduction Number")
    data = RData.load(fname)
    dates      = data["time"]
    values     = @_ data["draws"] |>
                    transpose |>
                    Array |>
                    __[:, 2001 : end]
    recipe = Regional.RtPlottingRecipe1( nothing, nothing, nothing, (; dates, values), "Rt")
    return plot!(p, recipe)
end

function firefox(p; fname="tmp.html", fdir=normpath(homedir(), ".tmp"))
    # mkpath(dir)
    filepath = normpath(fdir, fname)
    savefig(p, filepath)
    run(`firefox $(filepath)`, wait=false)
    return p
end

##
let
    ctype = "strangers"
    rawcontacts = readcontacts( projectdir("data/contacts/int/international_contacts.csv") )
    regions = unique(rawcontacts.region)
    for region in regions
        # plot raw data
        plts = Vector{Plots.Plot}()
        raw = filter(:region => ==(region), rawcontacts)
        p = scatter(raw.date, raw.strangers, mα=0.1, mc=colors[2], msc=:match; title="raw contacts to $ctype in $region", lab=nothing)
        push!(plts, p)
        # plotweeks!(p1, raw.date)

        # plot preprocessed data
        batched  = CSV.read(projectdir("data/contacts/int/mean_contact_rates_region=$(region)_aggregation=batch.csv"), DataFrame)
        smoothed = CSV.read(projectdir("data/contacts/int/smoothed_contact_rates_region=$(region).csv"), DataFrame)
        suffix = region != "Denmark" ? "" : "_low_resolution"

        p = plot(title="pre-processed contacts to $ctype in $region")
        scatter!(batched.date, batched[:,ctype], mα=1, mc=colors[1], msc=:match; lab=nothing)
        plot!(smoothed.date, smoothed[:,ctype*suffix]; label=false, ribbon=3smoothed[:,"$(ctype)_std"suffix], hover=smoothed.date, lw=4)
        # plotweeks!(p2, batched.date)
        push!(plts, p)

        # plot case counts
        epidemic = CSV.read(projectdir("data/epidata/epidata_region=$(region).csv"), DataFrame)
        is = findfirst(>(22), epidemic.cases_per_100000)
        ds = epidemic.date[is]
        de1 = epidemic.date[is:end][findfirst(<(10), epidemic.cases_per_100000[is:end])]
        de2 = epidemic.date[findfirst(>(0.05), epidemic.fraction_vaccinated)]
        de  = min(de1, de2)
        p = plot(title="observed cases")
        vspan!([ds, de], label="active epidemic", α=0.2, fc=:midnightblue)
        plot!(epidemic.date, epidemic.deaths, fc=colors[1], lc=colors[1], label="smoothed cases / 100_000", lw=4)#, hover=epidemic.date)
        #plot!(twinx(), epidemic.date, epidemic.cases, fc=colors[1], lc=colors[1], label="cases", lw=0.8, α=0.5)
        push!(plts, p)

        # plot predictors
        mobility = CSV.read(projectdir("data/mobility/mobility_country=$(region).csv"), DataFrame)
        p = plot(title="predictors")
        vspan!([ds, de], label="active epidemic", α=0.2, fc=:midnightblue)
        plot!(p, smoothed.date, zscore(smoothed[:,ctype]); label="contact to $ctype", hover=smoothed.date, lw=4)
        plot!(p, mobility.date, zscore(mobility.apple); label="apple", lw=1, hover=mobility.date)
        plot!(p, mobility.date, zscore(mobility.google); label="google", lw=1, hover=mobility.date)
        push!(plts, p)

        fname = projectdir("reports/stanfit/2021-10-05/sampled_rt_region=$(region).rds")
        if isfile(fname)
            p = plot_rt(fname)
            vspan!([ds, de], label="active epidemic", α=0.2, fc=:midnightblue)
            push!(plts, p)
        end

        p = plot(title="fraction of vaccinated")
        vspan!([ds, de], label="active epidemic", α=0.2, fc=:midnightblue)
        plot!(epidemic.date, epidemic.fraction_vaccinated * 100, fc=colors[1], lc=colors[1], label="vaccinated", lw=4)#, hover=epidemic.date)
        push!(plts, p)

        p = plot(plts..., layout=(length(plts),1), size=(1000, length(plts)*250), sharex=true, link=:x)
        firefox(p, fname="$region.html")
    end
end

## =========================================================================
#         tests
pr = let
    data = RData.load(projectdir("reports/stanfit/2021-10-05/sampled_predictions.rds"))
    dates      = data["time"]
    values     = @_ data["draws"] |>
                transpose |>
                Array |>
                __[:, 2001 : end]
    (; dates, values )
end
o = ( dates = epidemic.date, values = epidemic.deaths )

p = plot()
plot!(p, o.dates, o.values, α=0.5, lc=:match, lab="observed", c=:midnightblue, lw=4, ylab="cases")
plot_confidence_timeseries!(p, pr.dates, pr.values; label = "predicted (95%CI)", c=:midnightblue) #Dict(hover=>strdates)
plot_confidence_timeseries!(p, e.dates, e.values; label = "predicted (50%CI)") #Dict(hover=>strdates)
plot!(epidemic.date, epidemic.deaths, fc=colors[1], lc=colors[1], label="deaths", lw=4)#, hover=epidemic.date)
##
# plot case counts
let
    epidemic = CSV.read(projectdir("data/epidata/epidata_region=$(region).csv"), DataFrame)
    is = findfirst(>(22), epidemic.cases_per_100000)
    ds = epidemic.date[is]
    de1 = epidemic.date[is:end][findfirst(<(10), epidemic.cases_per_100000[is:end])]
    de2 = epidemic.date[findfirst(>(0.05), epidemic.fraction_vaccinated)]
    de  = min(de1, de2)
    p1 = plot(title="observed deaths")
    vspan!([ds, de], label="active epidemic", α=0.2, fc=:midnightblue)
    plot!(epidemic.date, epidemic.deaths, fc=colors[1], lc=colors[1], label="deaths", lw=4)#, hover=epidemic.date)
    #plot!(twinx(), epidemic.date, epidemic.cases, fc=colors[1], lc=colors[1], label="cases", lw=0.8, α=0.5)

    p2 = plot(recipe)

    # plot predictors
    mobility = CSV.read(projectdir("data/mobility/mobility_country=$(region).csv"), DataFrame)
    p3 = plot(title="predictors")
    vspan!([ds, de], label="active epidemic", α=0.2, fc=:midnightblue)
    plot!(p3, smoothed.date, zscore(smoothed[:,ctype]); label="contact to $ctype", hover=smoothed.date, lw=4)
    plot!(p3, mobility.date, zscore(mobility.apple); label="apple", lw=1, hover=mobility.date)
    plot!(p3, mobility.date, zscore(mobility.google); label="google", lw=1, hover=mobility.date)

    p4 = plot(title="fraction of vaccinated")
    vspan!([ds, de], label="active epidemic", α=0.2, fc=:midnightblue)
    plot!(epidemic.date, epidemic.fraction_vaccinated * 100, fc=colors[1], lc=colors[1], label="vaccinated", lw=4)#, hover=epidemic.date)

    p = plot(p1, p2, p3, p4, layout=(4,1), size=(1000, 4*250), sharex=true, link=:x)
end

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	load_fuels_data(setup::Dict, path::AbstractString, inputs_fuel::Dict)

Function for reading input parameters related to fuel costs and CO$_2$ content of fuels
"""
function load_fuels_data(setup::Dict, path::AbstractString, inputs_fuel::Dict)

	# Fuel related inputs - read in different files depending on if time domain reduction is activated or not
	#data_directory = chop(replace(path, pwd() => ""), head = 1, tail = 0)
	data_directory = joinpath(path, setup["TimeDomainReductionFolder"])
	if setup["TimeDomainReduction"] == 1  && isfile(joinpath(data_directory,"Load_data.csv")) && isfile(joinpath(data_directory,"Generators_variability.csv")) && isfile(joinpath(data_directory,"Fuels_data.csv")) # Use Time Domain Reduced data for GenX
		fuels_in = DataFrame(CSV.File(joinpath(data_directory,"Fuels_data.csv"), header=true), copycols=true)
	else  # Run without Time Domain Reduction OR Getting original input data for Time Domain Reduction
		fuels_in = DataFrame(CSV.File(joinpath(path,"Fuels_data.csv"), header=true), copycols=true)
	end

	# Fuel costs .&  CO2 emissions rate for each fuel type (stored in dictionary objects)
	fuels = names(fuels_in)[2:end] # fuel type indexes
	costs = Matrix(fuels_in[2:end,2:end])
	# New addition for variable fuel price
	CO2_content = fuels_in[1,2:end] # tons CO2/MMBtu
	fuel_costs = Dict{AbstractString,Array{Float64}}()
	fuel_CO2 = Dict{AbstractString,Float64}()
	for i = 1:length(fuels)
		if setup["ParameterScale"] == 1
			fuel_costs[string(fuels[i])] = costs[:,i] / ModelScalingFactor
			fuel_CO2[string(fuels[i])] = CO2_content[i] / ModelScalingFactor # kton/MMBTU
		else
			fuel_costs[string(fuels[i])] = costs[:,i]
			fuel_CO2[string(fuels[i])] = CO2_content[i] # ton/MMBTU
		end
	end

	inputs_fuel["fuels"] = fuels
	inputs_fuel["fuel_costs"] = fuel_costs
	inputs_fuel["fuel_CO2"] = fuel_CO2

	println("Fuels_data.csv Successfully Read!")

	return inputs_fuel, fuel_costs, fuel_CO2
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
function preproccsing!(df::DataFrames.DataFrame)
    df[!, :dX] .= 0.0
    df[!, :dY] .= 0.0
    df.dX[2:end] .= diff(df.corrected_x)
    df.dY[2:end] .= diff(df.corrected_y)
    df[df.FRAME2.==1, [:dX, :dY]] .= NaN
    df.dR2 = abs2.(df.dX) + abs2.(df.dY)
    df.dR = sqrt.(df.dR2)
    start_point = findall(x -> x == 2, df.FRAME2)
    end_point = start_point[2:end] .- 2
    terminus = size(df, 1)
    append!(end_point, terminus)
    track_length = Integer[]
    track_length = end_point .- start_point .+ 1
    track_num = maximum(df.TrackID)
    max_length = Int64(maximum(df.FRAME2))

    return track_length, track_num, max_length, start_point
end

function chi(
    df;
    idlabel = :TrackID,
    datalabel = :POSITION_X,
    δ::Int = 1
)
    dx = Float64[]
    @inbounds for n = sort(collect(Set(df[!, idlabel])))
        m = extract(df, n, idlabel, [datalabel])
        @simd for t = 1:size(m, 1)-δ
            append!(dx, m[t+δ] - m[t])
        end
    end
    dx
end

function displacement(
    df;
    idlabel = :TrackID,
    xlabel = :POSITION_X,
    ylabel = :POSITION_Y,
    δ::Int = 1
)
    dr = Float64[]
    @inbounds for n = sort(collect(Set(df[!, idlabel])))
        m = extract(df, n, idlabel, [xlabel, ylabel])
        @simd for t = 1:size(m, 1)-δ
            append!(dr, displacement(m, t, δ))
        end
    end
    dr
end

function xy2matrix(df, id, x, y, frame)
    N = maximum(df[!, id])
    xy = Array{Union{Nothing,Float64}}(nothing, maximum(df[!, frame]), 2, N + 1 - minimum(df[!, id]))
    @inbounds for (n, i) = enumerate(minimum(df[!, id]):maximum(df[!, id]))
    # @inbounds for n = minimum(df[!, id]):N
        m = extract(df, i, id, [x, y])
        @simd for t = 1:size(m, 1)
            xy[t, 1, n] = m[t, 1]
            xy[t, 2, n] = m[t, 2]
        end
    end
    xy
end

nomapround(b, d) = (x -> round.(x, digits = d)).(b)

"""
    group_files(dir) -> GroupedDataFrame{DataFrame}

Make grouped dataframe that exist in dir.
# TODO: raise an error if no csv files exist in the directory.
"""
function group_files(dir::AbstractString)
    cd(dir)
    genotypes = readdir()
    filter!(isdir, genotypes)
    df = []
    for genotype in genotypes
        files = glob("./$genotype/*.csv")
        for file in files
            tmp_df = CSV.read(file)
            tmp_df[:, :Genotype] .= genotype
            push!(df, tmp_df)
        end
    end
    df = vcat(df...)
    grouped_df = groupby(df, [:Genotype])
    return grouped_df
end

function label_mean_displacement!(
    gd::GroupedDataFrame;
    displacement = :mean_dis,
    max_dis::Int = 1,
)
    @argcheck displacement in [:mean_dis, :corrected_mean_dis]
    @argcheck max_dis > 0
    bin_array = range(0, stop = max_dis, length = max_dis * 20)
    transform!(
        gd,
        :mean_dis => (x -> int(cut(x, bin_array, extend = true), type = Int)) => :bin,
        ungroup = false,
    )
    transform!(
        gd,
        :corrected_mean_dis =>
            (x -> int(cut(x, bin_array, extend = true), type = Int)) => :corrected_bin,
        ungroup = false,
    )
    return collect(bin_array)
end

function time_average(df::DataFrame, f::Function, id::Symbol, datalabel)
    result = []
    N = maximum(df[:, id])
    @inbounds for n = 1:N
        data = Matrix(df[df[!, id].==n, datalabel])
        T = size(data, 1)
        for δ = 1:T-1
            c = 0.0
            @simd for t = 1:T-δ
                c += f(data, t, δ)
            end
            c /= (T - δ)
            push!(result, [n, c, δ])
        end
    end
    result
end

"""
    add_noise!(m) -> Matrix
Update m by adding a random uniform noise(± 1e-4), if m[t, (1 or 2)] ≈ ,[t-1, (1 or 2)].
This is needed when diff(coordinates) == 0 cause an error.
```
"""
function add_noise!(m)
    @inbounds @simd for t = 2:size(m, 1)
        if m[t, 1] ≈ m[t-1, 1]
            m[t, 1] += rand(Uniform(-1e-4, 1e-4))
        end
        if m[t, 2] ≈ m[t-1, 2]
            m[t, 2] += rand(Uniform(-1e-4, 1e-4))
        end
    end
end

function extract(df::DataFrame, n::Int, id::Symbol, x::Symbol)
    Vector(@where(df, cols(id) .== n)[!, x])
    # Vector{Float64}(df[df[:, id] .== n, x])
end

function extract(df::DataFrame, n::Int, id::Symbol, label::Vector{Symbol})
    Matrix(@where(df, cols(id) .== n)[!, label])
end

# https://discourse.julialang.org/t/a-good-way-to-filter-an-array-with-keeping-its-dimension/42603
function remove_nothing(x::AbstractArray)
    x[.!any.(isnothing, eachrow(x)), :]
end

# function frameaverage(data; L = 2)
#     new_data = zero(size(data))
#     @inbounds for t = 1:L-1
#         new_data[t, ]
#     @inbounds for t = 1:size(data, 1)

# end

# function displacement(r::AbstractMatrix, δ::Int)
#     sqrt((r[1+δ, 2] - r[1, 2])^2 + (r[1+δ, 1] - r[1, 1])^2)
# end

# function displacement(r::AbstractMatrix, t::Int, δ::Int)
#     sqrt((r[t+δ, 2] - r[t, 2])^2 + (r[t+δ, 1] - r[t, 1])^2)
# end

# function squared_displacement(r::AbstractMatrix, t::Int, δ::Int)
#     displacement(r, t, δ)^2
# end

max_stars_repo_name: SosUts/spt

max_issues_repo_name: SosUts/spt

max_forks_repo_name: SosUts/spt


text:
module GeoArrays

using ArchGDAL
using GDAL
using CoordinateTransformations
using StaticArrays

include("geoarray.jl")
include("geoutils.jl")
include("io.jl")
include("utils.jl")
include("interpolate.jl")
include("crs.jl")

export GeoArray

export coords
export centercoords
export indices

export compose!

export interpolate!

export epsg!

end


text:
using Flux
using Functors

struct ConvBlock
    conv
end

@functor ConvBlock
function ConvBlock(
    kernel_size::Int,
    in_channels::Int,
    out_channels::Int,
    act::Bool = true,
    down::Bool = true;
    kwargs...
)
return ConvBlock(
    Chain(
        if down
            Conv((kernel_size, kernel_size), in_channels => out_channels; kwargs...)
        else
            ConvTranspose((kernel_size, kernel_size), in_channels => out_channels; kwargs...)
        end,
        InstanceNorm(out_channels),
        if act
           x -> relu.(x)
        else
            identity
        end)
)  |> gpu
end

function (net::ConvBlock)(x)
    return net.conv(x)
end

struct ResidualBlock
    block
end

@functor ResidualBlock
function ResidualBlock(
    in_channels::Int
) return ResidualBlock(
    Chain(
        ConvBlock(3, in_channels, in_channels,true,true; pad=1),
        ConvBlock(3, in_channels, in_channels,false,true;pad=1)
    )
)   |> gpu
end

function(net::ResidualBlock)(x)
    return x + net.block(x)
end

struct Block
    conv
end

@functor Block
function Block(
    in_channels::Int,
    out_channels::Int;
    stride::Int
) return Block(
    Chain(
        Conv((4,4), in_channels=>out_channels; stride=stride, pad=1),
        InstanceNorm(out_channels),
        x -> leakyrelu.(x, 0.2)
    )
) |> gpu
end

function(net::Block)(x)
    return net.conv(x)
end


text:
IntOrString = Union{Int, AbstractString}

1::IntOrString
"Hello"::IntOrString
# error
# 1.0::IntOrString

Option = Union{T, Nothing}

max_stars_repo_name: mkanenobu/trashbox

max_issues_repo_name: mkanenobu/trashbox

max_forks_repo_name: mkanenobu/trashbox


text:
struct Dielectric <: Material
    refractive_index::AbstractFloat
end

function schlick(cosine::AbstractFloat, refractive_index::AbstractFloat)
    r0::AbstractFloat = ((1.0 - refractive_index) / (1.0 + refractive_index))^2
    r0 + (1.0 - r0) * (1.0 - cosine)^5
end

function scatter(dielectric::Dielectric, ray::Ray, hp::HitPoint)::Tuple{Bool, Union{Nothing, AbstractVec}, Union{Nothing, Ray}}
    attenuation::AbstractVec = Vec(1.0, 1.0, 1.0)
    local scattered::Ray

    succeeded::Bool, refracted::Union{Nothing, AbstractVec} = refract(ray.direction, hp.normal, 1.0 / dielectric.refractive_index)
    if succeeded
        cosine::AbstractFloat = -dot(ray.direction, hp.normal) / norm(ray.direction)
        reflect_prob::AbstractFloat = schlick(cosine, refractive_index)
        if rand() >= reflect_prob
            scattered = Ray(hp.point, refracted)
            return true, attenuation, scattered
        end
    end

    reflected::AbstractVec = reflect(ray.direction, hp.normal)
    scattered = Ray(hp.point, reflected)
    true, attenuation, scattered
end

max_stars_repo_name: goooxu/JuliaRayTracing

max_issues_repo_name: goooxu/JuliaRayTracing

max_forks_repo_name: goooxu/JuliaRayTracing


text:
# using Flux
using Pkg

Pkg.activate(".")

using JuliaRL
using JuliaRL.LinearRL
# using JuliaRL.Environments
using Random
using ProgressMeter

import JuliaRL.step!, JuliaRL.start!

# TileCoder = JuliaRL.FeatureCreators.TileCoder

abstract type AbstractQPolicy <: AbstractPolicy end

mutable struct EpsilonGreedyQPolicy <: AbstractQPolicy
    ϵ::Float64
    actions::AbstractArray
end

function get_action(policy::AbstractQPolicy, values::AbstractArray; rng=Random.GLOBAL_RNG)
    action = findmax(values)[2]
    if rand(rng) < policy.ϵ
        action = rand(rng, policy.actions)
    end
    return action
end

mutable struct TileCoderAgent <: AbstractAgent
    Q::SparseQFunction
    opt::Optimizer
    fc::FeatureCreators.TileCoder
    # iht::TileCoder.IHT
    # tilings::Integer
    # tiles::Integer
    actions::AbstractArray
    action::Integer
    γ::Float64
    policy::AbstractQPolicy
    ϕ_t::Array{Int64, 1}
    ϕ_tp1::Array{Int64, 1}
    function TileCoderAgent(opt::Optimizer, size_env_state::Integer, num_actions::Integer, tilings::Integer, tiles::Integer, γ::Float64, policy::AbstractQPolicy)

        num_features_per_action = (tilings*(tiles+1)^size_env_state)
        num_features = num_features_per_action*num_actions
        fc = FeatureCreators.TileCoder(tilings, tiles, size_env_state; wrap=false, wrapwidths=0.0)
        Q = SparseQFunction(
            num_features,
            num_features_per_action,
            num_actions)
        # opt = WatkinsQ(α)
        actions = 1:3
        return new(Q, opt, fc, actions, -1, γ, policy, zeros(Int64, tilings), zeros(Int64, tilings))
    end
end

get_action(agent::TileCoderAgent, ϕ; rng=Random.GLOBAL_RNG) = get_action(agent.policy, [agent.Q(ϕ, a) for a = agent.actions]; rng=rng)

function start!(agent::TileCoderAgent, env_s_tp1; rng=Random.GLOBAL_RNG, kwargs...)
    # agent.ϕ_t .= TileCoder.tiles!(
    #     agent.iht,
    #     agent.tilings,
    #     env_s_tp1.*agent.tiles)
    agent.ϕ_t .= FeatureCreators.create_features(agent.fc, env_s_tp1)

    agent.action = get_action(agent, agent.ϕ_t; rng=rng)
    # if rand(rng) < agent.ϵ
    #     agent.action = rand(1:3)
    # end
    return agent.action
end

function step!(agent::TileCoderAgent, env_s_tp1, r, terminal; rng=Random.GLOBAL_RNG, kwargs...)
    # agent.ϕ_tp1 .= TileCoder.tiles!(
    #     agent.iht,
    #     agent.tilings,
    #     env_s_tp1.*agent.tiles)

    agent.ϕ_tp1 .= FeatureCreators.create_features(agent.fc, env_s_tp1)
    # println(agent.ϕ_t)

    LinearRL.update!(
        agent.Q,
        agent.opt,
        agent.ϕ_t,
        agent.ϕ_tp1,
        r,
        agent.γ,
        1.0,
        terminal,
        agent.action)

    agent.ϕ_t .= agent.ϕ_tp1
    # println(agent.ϕ_t)
    agent.action = get_action(agent, agent.ϕ_t; rng=rng)
    # println(agent.ϕ_t)
    # println([agent.Q(agent.ϕ_t, a) for a = agent.actions])
    # println(agent.Q)
    # if rand() < agent.ϵ
    #     agent.action = rand(rng, 1:3)
    # end

    return agent.action

end

function run(α=0.5/8, ϵ=0.1, γ=1.0, tilings=32, tiles=2; n_episodes=1000, seed=1002349)

    ϵ = 0.1
    α = 0.5/tilings
    # weights = zeros((tilings*(tiles+1)^2)*3)
    # ag = TileCoderAgent(SparseQFunction((tilings*(tiles+1)^2)*3, ((tilings*(tiles+1)^2)), 3),
    #                     TileCoder.IHT(tilings*(tiles+1)^2),
    #                     1:3)
    opt = WatkinsQ(α)
    agent = TileCoderAgent(opt, 2, 3, tilings, tiles, γ, EpsilonGreedyQPolicy(ϵ, 1:3))

    rng = Random.MersenneTwister(seed)

    env = MountainCar(rng)
    cumulative_reward_array = zeros(Int64, n_episodes)
    # @showprogress 0.1 "Episode: " for episode = 1:n_episodes
    for episode = 1:n_episodes

        terminal = false
        num_steps = 0
        cumulative_reward = 0
        _, state = start!(env; rng=rng)
        action = start!(agent, state)
        while !terminal

            _, state_prime, reward, terminal = step!(env, action-1)

            action = step!(agent, state_prime, reward, terminal)

            num_steps += 1
            cumulative_reward += reward

        end
        cumulative_reward_array[episode] = cumulative_reward
    end
    return cumulative_reward_array
end

max_stars_repo_name: mkschleg/JuliaRL

max_issues_repo_name: mkschleg/JuliaRL

max_forks_repo_name: mkschleg/JuliaRL


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualNetworkRuleListResult <: SwaggerModel
    value::Any # spec type: Union{ Nothing, Vector{VirtualNetworkRule} } # spec name: value
    nextLink::Any # spec type: Union{ Nothing, String } # spec name: nextLink

    function VirtualNetworkRuleListResult(;value=nothing, nextLink=nothing)
        o = new()
        validate_property(VirtualNetworkRuleListResult, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        validate_property(VirtualNetworkRuleListResult, Symbol("nextLink"), nextLink)
        setfield!(o, Symbol("nextLink"), nextLink)
        o
    end
end # type VirtualNetworkRuleListResult

const _property_map_VirtualNetworkRuleListResult = Dict{Symbol,Symbol}(Symbol("value")=>Symbol("value"), Symbol("nextLink")=>Symbol("nextLink"))
const _property_types_VirtualNetworkRuleListResult = Dict{Symbol,String}(Symbol("value")=>"Vector{VirtualNetworkRule}", Symbol("nextLink")=>"String")
Base.propertynames(::Type{ VirtualNetworkRuleListResult }) = collect(keys(_property_map_VirtualNetworkRuleListResult))
Swagger.property_type(::Type{ VirtualNetworkRuleListResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualNetworkRuleListResult[name]))}
Swagger.field_name(::Type{ VirtualNetworkRuleListResult }, property_name::Symbol) =  _property_map_VirtualNetworkRuleListResult[property_name]

function check_required(o::VirtualNetworkRuleListResult)
    true
end

function validate_property(::Type{ VirtualNetworkRuleListResult }, name::Symbol, val)
end


text:
export Flat

struct Flat
end

Distributions.logpdf(::Flat, x) = 0.0

struct HalfFlat
end

Distributions.logpdf(::HalfFlat, x) = 0.0


text:
module HTSLib

import BioGenerics:
    BioGenerics,
    seqname,
    sequence,
    metadata,
    leftposition,
    rightposition,
    isoverlapping,
    metainfoval,
    metainfotag,
    header,
    IO.tryread!,
    IO.stream

# BioGenerics reexport
export seqname,
    sequence,
    metadata,
    leftposition,
    rightposition,
    isoverlapping,
    metainfoval,
    metainfotag,
    header,
    tryread!,
    stream

export sequence!,
    setsequence!,
    seqlength,
    mappingquality,
    refid,
    refname,
    quality,
    quality!,
    setquality!,
    queryname,
    next_refid,
    next_refpos

export htslib
export eachrecord

using Printf: @printf
using Setfield: @set

import HttpIO: HttpIO, PoorGCloudAuth, GCSClient, GCSFileIO, HttpFileIO

import TranscodingStreams

include("htslib/htslib.jl")
include("recordflag.jl")
include("bamrecord.jl")
include("htsopen.jl")
include("iteration.jl")

"""
The version of htslib.
"""
const HTSLIB_VERSION = VersionNumber(unsafe_string(htslib.hts_version()))

end # module


text:
# This file includes utilities

export Q, dbtoval, valtodb, energy

"""
    $SIGNATURES

Computes the energy of the discrete time signal `s` sampled with `ts` (defaults to 1) seconds.
"""
energy(s::AbstractVector, ts=1.) = sum(abs.(s).^2) * ts

"""
    $SIGNATURES

Converts `γ` from `dB` value to its real value.

# Example
```julia
julia> dbtoval(2)  # γ = 2 dB
1.5848931924611136
```
"""
dbtoval(γ) = 10^(γ / 10)

"""
    $SIGNATURES

Converts `val` tı dB scale.
"""
valtodb(val) = 10 * log10(val)

"""
    $SIGNATURES

Q-function defined as.
```math
    Q(x) = \\int_{x}^{\\infty} exp(-\\dfrac{x^2}{2}) dx
```
"""
Q(x) = 1 / 2 * erfc(x / sqrt(2))


text:
export
    rows, cols, diags, antidiags

rows(m::AbstractArray{T,2}) where {T} = [vec(m[i,:]) for i in 1:size(m, 1)]
cols(m::AbstractArray{T,2}) where {T} = [vec(m[:,j]) for j in 1:size(m, 2)]
diags(m::AbstractArray{T,2}) where {T} = [diag(m, k) for k in (1 - size(m, 1)):(size(m, 2) - 1)]
antidiags(m::AbstractArray{T,2}) where {T} = rotl90(m) |> diags


text:
#===================================================================================================
Multiobjective differential evolution (DE) with DEMO: a hybridization of DE and NSGA-II

Robič, Tea, and Bogdan Filipič. "Differential evolution for multiobjective optimization."
International conference on evolutionary multi-criterion optimization. Springer, Berlin, Heidelberg,
2005.
===================================================================================================#

"""
    demo!(de::StdMoDE, evaluator::Function, ngen; show_progress::Bool=true)::Vector{Vector{Int}}

DEMO algorithm for multiobjective differential evolution. `evaluator` is the objective function that
should return a vector of multiple values.

As recommended in the paper, the DEMO/parent variant is implemented. The mutation operator is
DE/rand/1 and the recombination is conducted with binomial crossover.

Robič, Tea, and Bogdan Filipič. "Differential evolution for multiobjective optimization."
International conference on evolutionary multi-criterion optimization. Springer, Berlin, Heidelberg,
2005.
"""
function demo!(de::StdMoDE, evaluator::Function, ngen; show_progress::Bool=true)::Vector{Vector{Int}}
    n = nindividuals(de)
    d = ndims(de)
    ss = senses(de)
    # initial evaluation
    de.fitness .= hcat(evaluator.(eachcol(de.pop))...)
    # in DEMO, we will have an intermediate population of size [n, 2n]
    im_pop = zeros(eltype(de.pop), d, 2*n)
    im_fitness = zeros(eltype(de.fitness), nobjectives(de), 2*n)
    im_indices = collect(1:2*n)
    U = de.pop[:, 1] # copy to avoid allocation
    V = de.pop[:, 1]
    F = de.options.F
    Cr = de.options.Cr
    if show_progress
        pg = Progress(ngen, 0.1)
        ProgressMeter.ijulia_behavior(:clear)
    end
    for g = 1:ngen
        im_pop[:, 1:n] .= de.pop
        im_fitness[:, 1:n] .= de.fitness
        im_indices .= 1: 2*n
        cursor = n # the actual number of individuals in `im_pop`
        for i = 1:n  # produce n offspring
            mutate_rand_1!(V, i, im_pop, F, @view im_indices[1:cursor]) # donor
            X = view(im_pop, :, i)  # target (base)
            bounce_back!(V, X, de.bounds)
            crossover_binomial!(U, X, V, Cr)  # trial
            fit = evaluator(U)
            parent_fit = @view im_fitness[:, i]
            if dominate(fit, parent_fit, ss) # better, replace the parent
                im_pop[:, i] .= U
                im_fitness[:, i] .= fit
            elseif !dominate(parent_fit, fit, ss)  # equally good, add to the population
                cursor += 1
                im_pop[:, cursor] .= U
                im_fitness[:, cursor] .= fit
            end
        end
        # now we have a population of size `cursor` ∈ [n, 2n]
        # follow NSGA-II to reduce the population to size n
        if cursor > n
            fronts = nondominated_sort(@view(im_fitness[:, 1:cursor]), ss)
            i_ind = 0
            i_front = 1
            # copy individuals in preferential fronts to de
            while i_ind + length(fronts[i_front]) <= n
                front = fronts[i_front]
                l = length(front)
                de.pop[:, i_ind+1:i_ind+l] .= @view im_pop[:, front]
                de.fitness[:, i_ind+1:i_ind+l] .= @view im_fitness[:, front]
                i_ind += l
                i_front += 1
            end
            # if there are still vacancies
            # we need only a portion of the next front accroding to crowding distance
            if i_ind < n
                front = fronts[i_front]
                l = length(front)
                im_indices[1:l] .= 1:l
                cds_i_front = assign_crowding_distance(view(im_fitness, :, front), [@view im_indices[1:l]])
                # arg-sort of the front in decending order w.r.t crowding distance
                ix = @view im_indices[1:l] # just a buffer
                sortperm!(ix, cds_i_front; rev=true)
                # pick the ones with larger crowding distance
                for i = 1:n-i_ind
                    picked = front[ix[i]]
                    de.pop[:, i_ind+i] .= @view im_pop[:, picked]
                    de.fitness[:, i_ind+i] .= @view im_fitness[:, picked]
                end
            end
        else # we have coincidentally n individuals
            front = Vector{Int}[]
            de.pop .= @view im_pop[:, 1:n]
            de.fitness .= @view im_fitness[:, 1:n]
        end

        if show_progress
            ProgressMeter.next!(pg; showvalues=[(:gen, "$g/$ngen")])
        end
    end
    # finally, perform nondominated sorting again to obtain fronts
    return nondominated_sort(de.fitness, ss)
end


text:
# will be executed on workspace process

module PlutoFormatter

using Base
using Markdown
import Markdown: html, htmlinline, LaTeX, withtag, htmlesc
import Distributed
import Base64

# We add a method for the Markdown -> HTML conversion that takes a LaTeX chunk from the Markdown tree and adds our custom span
function htmlinline(io::IO, x::LaTeX)
    withtag(io, :span, :class => "tex") do
        print(io, '$')
        htmlesc(io, x.formula)
        print(io, '$')
    end
end

# this one for block equations: (double $$)
function html(io::IO, x::LaTeX)
    withtag(io, :p, :class => "tex") do
        print(io, '$', '$')
        htmlesc(io, x.formula)
        print(io, '$', '$')
    end
end

"The `IOContext` used for converting arbitrary objects to pretty strings."
iocontext = IOContext(stdout, :color => false, :compact => true, :limit => true, :displaysize => (18, 120))

"""Format `val` using the richest possible output, return formatted string and used MIME type.

Currently, the MIME type is one of `text/html` or `text/plain`, the former being richest."""
function format_output(val::Any)::Tuple{String, MIME}
    # in order of coolness
    # text/plain always matches
    mime = let
        mimes = [MIME("text/html"), MIME("image/svg+xml"), MIME("image/png"), MIME("image/jpg"), MIME("text/plain")]
        first(filter(m->Base.invokelatest(showable, m, val), mimes))
    end

    if val === nothing
        "", mime
    else
        try
            result = Base.invokelatest(repr, mime, val; context = iocontext)
            # MIME rewrites for types other than text/plain or text/html
            if mime ∈ [MIME("image/png"), MIME("image/jpg")]
                result = "<img src=\"data:$(mime);base64,$(Base64.base64encode(result))\" />"
                mime = MIME("text/html")
            end
            result, mime
        catch ex
            "Failed to show value: \n" * sprint(showerror, ex, stacktrace(backtrace())), MIME("text/plain")
        end
    end
end

function format_output(ex::Exception, bt::Array{Any, 1})::Tuple{String, MIME}
    sprint(showerror, ex, bt), MIME("text/plain")
end

function format_output(ex::Exception)::Tuple{String, MIME}
    sprint(showerror, ex), MIME("text/plain")
end

function format_output(val::CapturedException)::Tuple{String, MIME}

    ## We hide the part of the stacktrace that belongs to Pluto's evalling of user code.

    bt = try
        new_bt = val.processed_bt
        # If this is a ModuleWorkspace, then that's everything starting from the last `eval`.
        # For a ProcessWorkspace, it's everything starting from the 2nd to last `eval`.
        howdeep = Distributed.myid() == 1 ? 1 : 2

        for _ in 1:howdeep
            until = findfirst(b -> b[1].func == :eval, reverse(new_bt))
            new_bt = until === nothing ? new_bt : new_bt[1:(length(new_bt) - until)]
        end

        # We don't include the deepest item of the stacktrace, since it is always
        # `top-level scope at none:0`
        new_bt[1:end-1]
    catch ex
        val.processed_bt
    end

    format_output(val.ex, bt)
end

end


text:
function e010(N)
    return sum_of_primes_below(N)
end

function solution010()
    N = 2000000
    answer = e010(N)
    println("The answer to problem 010 is $answer")
end

function sum_of_primes_below(n::Integer)
    s = 0
    for x in Channel(c -> get_primes_below(c,n))
        s += x
    end
    return s
end

max_stars_repo_name: mithodin/ttpe

max_issues_repo_name: mithodin/ttpe

max_forks_repo_name: mithodin/ttpe


text:
"""
    RungeKuttaConvolutionQuadrature{T,N,NN}

T: the value type of the basis function.
N: the number of stages.
NN: N*N.

Performs a convolution quadrature on a laplaceKernel to represent an operator
in time domain using an implicit Runge-Kutta method.

laplaceKernel: function of the Laplace variable s that returns an IntegralOperator.
A, b: Coefficient matrix and vectors from the Butcher tableau.
Δt: time step.
zTransformedTermCount: Number of terms in the inverse Z-transform.
contourRadius: radius of circle used as integration contour for the inverse Z-transform.
"""
struct RungeKuttaConvolutionQuadrature{LK, T, N, NN}
	laplaceKernel         :: LK # function of s that returns an IntegralOperator
	A                     :: SArray{Tuple{N,N},T,2,NN}
	b                     :: SVector{N,T}
	Δt                    :: T
	zTransformedTermCount :: Int
	contourRadius         :: T
end
scalartype(rkcq::RungeKuttaConvolutionQuadrature{LK, T, N, NN}) where {LK, T, N, NN} = Complex{T};

# M = H*diagm(D)*invH
struct DiagonalizedMatrix{T,N,NN}
	H    :: SArray{Tuple{N,N},Complex{T},2,NN}
	invH :: SArray{Tuple{N,N},Complex{T},2,NN}
	D    :: SVector{N,Complex{T}}
end

# M = H*diagm(D)*invH
function diagonalizedmatrix(M :: SArray{Tuple{N,N},Complex{T},2,NN}) where {T,N,NN}
	ef = eigfact(Array{Complex{T},2}(M));

	efValues = SVector{N,Complex{T}}(ef.values)  :: SVector{N,Complex{T}};
	efVectors = SArray{Tuple{N,N},Complex{T},2,NN}(ef.vectors) :: SArray{Tuple{N,N},Complex{T},2,NN};
	return DiagonalizedMatrix(efVectors, inv(efVectors), efValues);
end

function assemble(rkcq :: RungeKuttaConvolutionQuadrature,
                  testfns :: StagedTimeStep,
                  trialfns :: StagedTimeStep)

	laplaceKernel = rkcq.laplaceKernel;
	A = rkcq.A;
	b = rkcq.b;
	Δt = rkcq.Δt;
	Q = rkcq.zTransformedTermCount;
	rho = rkcq.contourRadius;
	p = length(b); # stage count

	test_spatial_basis  = testfns.spatialBasis;
	trial_spatial_basis = trialfns.spatialBasis;

	# Compute the Z transformed sequence.
	# Assume that the operator applied on the conjugate of s is the same as the
	# conjugate of the operator applied on s,
	# so that only half of the values are computed
	Qmax = Q>>1+1;
	M = numfunctions(test_spatial_basis);
	N = numfunctions(trial_spatial_basis);
	Tz = promote_type(scalartype(rkcq), scalartype(testfns), scalartype(trialfns));
	Zz = Vector{Array{Tz,2}}(undef,Qmax);
	blocksEigenvalues = Vector{Array{Tz,2}}(undef,p);
	tmpDiag = Vector{Tz}(undef,p);
	for q = 0:Qmax-1
		# Build a temporary matrix for each eigenvalue
		s = laplace_to_z(rho, q, Q, Δt, A, b);
		sFactorized = diagonalizedmatrix(s);
		for (i,sD) in enumerate(sFactorized.D)
			blocksEigenvalues[i] = assemble(laplaceKernel(sD), test_spatial_basis, trial_spatial_basis);
		end

		# Compute the Z transformed matrix by block
		Zz[q+1] = zeros(Tz, M*p, N*p);
		for m = 1:M
			for n = 1:N
				for i = 1:p
					tmpDiag[i] = blocksEigenvalues[i][m,n];
				end
				D = SVector{p,Tz}(tmpDiag);
				Zz[q+1][(m-1)*p+(1:p),(n-1)*p+(1:p)] = sFactorized.H * diagm(D) * sFactorized.invH;
			end
		end
	end

	# return the inverse Z transform
	kmax = Q;
	T = real(Tz);
	Z = zeros(T, M*p, N*p, kmax)
	for q = 0:kmax-1
		Z[:,:,q+1] = real_inverse_z_transform(q, rho, Q, Zz);
	end
	return Z

end


text:
using SIIPExamples

using PowerSystems
using PowerSimulations
const PSI = PowerSimulations
using PowerSystemCaseBuilder

using Dates
using DataFrames

using Cbc #solver

sys = build_system(PSITestSystems, "modified_RTS_GMLC_DA_sys")

template_uc = OperationsProblemTemplate()

print_tree(PSI.AbstractDeviceFormulation)

set_device_model!(template_uc, Line, StaticBranch)
set_device_model!(template_uc, Transformer2W, StaticBranch)
set_device_model!(template_uc, TapTransformer, StaticBranch)

set_device_model!(template_uc, ThermalStandard, ThermalStandardUnitCommitment)
set_device_model!(template_uc, RenewableDispatch, RenewableFullDispatch)
set_device_model!(template_uc, PowerLoad, StaticPowerLoad)
set_device_model!(template_uc, HydroDispatch, FixedOutput)
set_device_model!(template_uc, HydroEnergyReservoir, HydroDispatchRunOfRiver)
set_device_model!(template_uc, RenewableFix, FixedOutput)

set_service_model!(template_uc, VariableReserve{ReserveUp}, RangeReserve)
set_service_model!(template_uc, VariableReserve{ReserveDown}, RangeReserve)

set_transmission_model!(template_uc, CopperPlatePowerModel)

solver = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 1, "ratioGap" => 0.5)

op_problem = OperationsProblem(template_uc, sys; optimizer = solver, horizon = 24)

build!(op_problem, output_dir = mktempdir())

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
function addCloudVert!(fgl::FactorGraph, exvert::Graphs.ExVertex)
  cv = CloudGraphs.exVertex2CloudVertex(exvert);
  CloudGraphs.add_vertex!(fgl.cg, cv);
  fgl.cgIDs[exvert.index] = cv.neo4jNodeId
  IncrementalInference.addGraphsVert!(fgl, exvert)
end

# Return Graphs.ExVertex type containing data according to id
function getExVertFromCloud(fgl::FactorGraph, fgid::Int64; bigdata::Bool=false)
  neoID = fgl.cgIDs[fgid]
  cvr = CloudGraphs.get_vertex(fgl.cg, neoID, false)
  CloudGraphs.cloudVertex2ExVertex(cvr)
end

function getExVertFromCloud(fgl::FactorGraph, lbl::ASCIIString; bigdata::Bool=false)
  getExVertFromCloud(fgl, fgl.IDs[lbl],bigdata=bigdata)
end

function updateFullCloudVertData!(fgl::FactorGraph,
    nv::Graphs.ExVertex)

  # TODO -- this get_vertex seems excessive, but we need the CloudVertex
  neoID =fgl.cgIDs[nv.index]
  # println("updateFullCloudVertData! -- trying to get $(neoID)")
  vert = CloudGraphs.get_vertex(fgl.cg, neoID, false)
  vert.packed = nv.attributes["data"]
  # TODO -- ignoring other properties
  fgl.g.vertices[nv.index].attributes["data"] = nv.attributes["data"]
  CloudGraphs.update_vertex!(fgl.cg, vert)
end

function makeAddCloudEdge!(fgl::FactorGraph, v1::Graphs.ExVertex, v2::Graphs.ExVertex)
  cv1 = CloudGraphs.get_vertex(fgl.cg, fgl.cgIDs[v1.index], false)
  cv2 = CloudGraphs.get_vertex(fgl.cg, fgl.cgIDs[v2.index], false)
  ce = CloudGraphs.CloudEdge(cv1, cv2, "DEPENDENCE");
  retrel = CloudGraphs.add_edge!(fgl.cg, ce);

  # TODO -- keep this edge id in function node data, must refactor
  push!(v2.attributes["data"].edgeIDs, retrel.id) # TODO -- not good way to do this
  updateFullCloudVertData!(fgl, v2)

  IncrementalInference.makeAddEdge!(fgl, v1, v2, saveedgeID=false)
  retrel.id
end

# return list of neighbors as Graphs.ExVertex type
function getCloudOutNeighbors(fgl::FactorGraph, vert::Graphs.ExVertex)
  # @show vert.index
  # @show fgl.cgIDs
  cgid = fgl.cgIDs[vert.index]
  cv = CloudGraphs.get_vertex(fgl.cg, cgid, false)
  neighs = CloudGraphs.get_neighbors(fgl.cg, cv)
  neExV = Graphs.ExVertex[]
  for n in neighs
    push!(neExV,  CloudGraphs.cloudVertex2ExVertex(n))
  end
  return neExV
end

function getEdgeFromCloud(fgl::FactorGraph, id::Int64)
  println("getting id=$(id)")
  CloudGraphs.get_edge(fgl.cg, id)
end

function deleteCloudVertex!(fgl::FactorGraph, vert::Graphs.ExVertex)
  neoID = fgl.cgIDs[vert.index]
  cvr = CloudGraphs.get_vertex(fgl.cg, neoID, false)
  CloudGraphs.delete_vertex!(fgl.cg, cvr)
end

function deleteCloudEdge!(fgl::FactorGraph, edge::CloudEdge)
  CloudGraphs.delete_edge!(fgl.cg, edge)
end

function setCloudDataLayerAPI!()
  # cgapi = DataLayerAPI(addCloudVert!,            # addvertex
  #                      dlapi.getvertex,          # getvertex
  #                      makeAddCloudEdge!,        # makeaddedge
  #                      graphsGetEdge,           # getedge
  #                      dlapi.outneighbors,       # outneighbors
  #                      +, +, +, + )
  dlapi.addvertex! = addCloudVert!
  dlapi.getvertex = getExVertFromCloud
  dlapi.makeaddedge! = makeAddCloudEdge!
  dlapi.getedge = getEdgeFromCloud
  dlapi.updatevertex! = updateFullCloudVertData!
  dlapi.outneighbors = getCloudOutNeighbors
  dlapi.deletevertex! = deleteCloudVertex!
  dlapi.deleteedge! = deleteCloudEdge!
  dlapi.cgEnabled = true

  println("Changed internal API calls to use CloudGraphs in appropriate places.")
  nothing
end

  #


text:
"""
    log_embeddings(logger::TBLogger, name::AbstractString, mat::AbstractMatrix; metadata, metadata_header, img_labels, step=step(logger))

Log embedding data to tensorboard and visualize in 3-D or 2-D with PCA, t-SNE or UMAP.
- mat: 2-D `Matrix` of data, with rows representing the samples and columns representing the features
- metadata: Array of labels for each sample. Each element across 1st dimenstion will be converted to `string`
- metadata_header: 1-D Array. Useful when samples have multiple labels. size should be same as size of each row in `metadata`
- img_labels: TBImages object representing image labels for each sample.
    - each number of images (N) must be equal to number of samples in `mat`
    - each image must be a square (H == W)
    - the value √N * W must be less than or equal to 8192 because of tensorboard restrictions.
"""
function log_embeddings(logger::TBLogger, name::AbstractString, mat::AbstractMatrix; metadata=nothing, metadata_header=nothing, img_labels=nothing, step=nothing)
    ndims(mat) == 2 || throw(DimensionMismatch("Embedding matrix must be 2-Dimensional. ("*string(ndims(mat))*" ≠ 2)"))
    matrix_path = joinpath(logger.logdir, repr(step), name)
    mkpath(matrix_path)
    if metadata != nothing
        size(metadata, 1) == size(mat, 1) || throw(ErrorException("#labels must be equal to #samples. ("*string(length(metadata))*" ≠ "*string(size(mat, 1))*")"))
        if metadata_header == nothing
            metadata = [string(x) for x in metadata]
        else
            length(metadata_header) == size(metadata, 2) || throw(ErrorException("length of header must be equal to the number of columns in metadata. ("*string(length(metadata_header))*" ≠ "*string(size(metadata, 2))*")"))
            metadata = [join(metadata_header, '\t'); [join(metadata[i, :], '\t') for i in 1:size(metadata, 1)]]
        end
        write_metadata(metadata, matrix_path)
    end
    if img_labels != nothing
        img_labels = convert_to_NCHW(content(img_labels), img_labels.format)
        size(img_labels, 1) == size(mat, 1) || throw(ErrorException("#label images must be equal to #samples. ("*string(size(img_labels, 1))*" ≠ "*string(size(mat, 1))*")"))
        size(img_labels, 3) == size(img_labels, 4) || throw(ErrorException("#label images must be square"))
        write_sprite(img_labels, matrix_path)
    end
    write_matrix(mat, matrix_path)
    write_pbtext(name, logger.logdir, matrix_path, metadata, img_labels, step)
end

function write_matrix(mat::AbstractMatrix, matrix_path)
    matrix_path = joinpath(matrix_path, "tensor.tsv")
    mat = convert(Array{Float64,2}, mat)
    open(matrix_path, "w") do file
        for i in 1:size(mat, 1)
            write(file, join(mat[i, :], '\t'))
            write(file, '\n')
        end
    end
end

function write_metadata(metadata::AbstractArray, matrix_path)
    matrix_path = joinpath(matrix_path, "metadata.tsv")
    open(matrix_path, "w") do file
        for x in metadata
            write(file, string(x)*'\n')
        end
    end
end

function write_sprite(img_labels::AbstractArray, matrix_path)
    n, _, _, w = size(img_labels)
    sqrt(n)*w <= 8192 || throw(ErrorException("the value √N * W must be less than or equal to 8192 because of tensorboard restrictions"))
    total_pixels = size(img_labels, 1)*size(img_labels, 3)*size(img_labels, 4)
    pixels_one_side = sqrt(total_pixels)
    number_of_images_per_row = Integer(ceil(pixels_one_side/size(img_labels, 4)))
    arranged_img_CHW = make_grid_of_images(img_labels, number_of_images_per_row)
    sprite_size = size(arranged_img_CHW, 3)
    arranged_augment_square_CHW = zeros((3, sprite_size, sprite_size))
    arranged_augment_square_CHW[:, 1:size(arranged_img_CHW, 2), :] = arranged_img_CHW
    sprite_path = joinpath(matrix_path, "sprite.png")
    open(sprite_path; write=true) do io
        save(Stream{format"PNG"}(io), colorview(RGB, arranged_augment_square_CHW))
    end
end

function make_grid_of_images(img_labels::AbstractArray, ncols::Integer)
    if size(img_labels, 2) == 1
        img_labels = cat(img_labels, img_labels, img_labels, dims = 2)
    end
    n, c, h, w = size(img_labels)
    ncols = min(ncols, n)
    nrows = Integer(ceil(n/ncols))
    grid = zeros(c, h*nrows, w*ncols)
    i = 0
    for y in 1:nrows
        for x in 1:ncols
            if i >= n
                break
            end
            grid[:, (y-1)*h+1:y*h, (x-1)*w+1:x*w] = img_labels[i+1, :, :, :]
            i = i+1
        end
    end
    grid
end

function write_pbtext(name::AbstractString, path, matrix_path, metadata, img_labels, step)
    metadata_path = joinpath(matrix_path, "metadata.tsv")
    img_labels_path = joinpath(matrix_path, "sprite.png")
    matrix_path = joinpath(matrix_path, "tensor.tsv")
    path = joinpath(path, "projector_config.pbtxt")
    isfile(path) || write(path, "") # workaround https://github.com/JuliaCloud/AWSS3.jl/issues/173
    open(path, "a") do file
        write(file, "embeddings {\n")
        write(file, "tensor_name: \""*name*":"*repr(step)*"\"\n")
        write(file, string("tensor_path: \"", matrix_path, "\"\n"))
        if metadata != nothing
            write(file, string("metadata_path: \"", metadata_path, "\"\n"))
        end
        if img_labels != nothing
            write(file, "sprite {\n")
            write(file, string("image_path: \"", img_labels_path, "\"\n"))
            write(file, "single_image_dim: "*string(size(img_labels, 4))*"\n")
            write(file, "single_image_dim: "*string(size(img_labels, 3))*"\n")
            write(file, "}\n")
        end
        write(file, "}\n")
    end
    return
end


text:
# example of reading properties from chemsep .xml database with PhysProps

using PhysProps
using CSV # for writing out the data frame

fileName = String("data\\chemsep1.xml")

dbRoot = getDBRoot(fileName)

compound = PhysProps.parse_compound(firstelement(dbRoot)) #returns a Dict


text:
function create_space(config::Config)
    base_space = LightGraphs.barabasi_albert(config.n_nodes, config.n0, config.k)
    for new_node in (config.n_nodes + 1):1:(config.n_nodes + config.agentcount)
        sampled_connection_node = StatsBase.sample(1:config.n_nodes)
        add_vertex!(base_space)
        add_edge!(base_space, new_node, sampled_connection_node)
    end
    return Agents.GraphSpace(base_space)
end

function initialize_model(config::Config, space::Agents.GraphSpace)
    properties = get_model_properties(config)
    model = Agents.AgentBasedModel(Citizen, space; properties=properties, scheduler=random_activation)
    add_public_spaces!(model, config.n_public_spaces)
    populate_space!(model, config.initial_infectives)
    initialize_statecounts!(model)
    return model
end

function get_model_properties(config::Config)
    return Dict(
        :agentcount => config.agentcount,
        :transmission_prob => config.transmission_prob,
        :test_prob => config.test_prob,
        :quarantine_duration => config.quarantine_duration,
        :behavior_reduction_factor => config.behavior_reduction_factor,
        :cost_matrix => (player1 = [1 3; 2 2], player2 = [1 2; 3 2]),
        :behavior_activated => config.behavior_activated,
        :authority_recommendation => 0.0,
        :tick_authority_recommendation => config.tick_authority_recommendation,
        :tick_authority_policy => config.tick_authority_policy,
        :homes => Int64[],
        :public_spaces => Int64[],
        :behavior_policy => false,
        :clock => 0,
        :S_count => config.agentcount,
        :E_count => 0,
        :Iu_count => 0,
        :Id_count => 0,
        :I_count => 0,
        :R_count => 0,
        :B_count => 0,
        :Id_cumulative => 0,
        :I_cumulative => 0,
        :new_cases_real => 0,
        :new_cases_detected => 0,
    )
end

function add_public_spaces!(model::Agents.ABM, n_public_spaces::Int64)
    degree_dict = Dict()
    for v in LightGraphs.vertices(model.space.graph)
        degree_dict[v] = LightGraphs.degree(model.space.graph, v)
    end
    vertices_by_degree = sort(collect(degree_dict), by=last, rev=true)
    for i in 1:n_public_spaces
        push!(model.public_spaces, popfirst!(vertices_by_degree)[1])
    end
    return model
end

function populate_space!(model::Agents.ABM, initial_infectives::Int64)
    highest_home = deepcopy(LightGraphs.nv(model.space.graph))
    lowest_home = highest_home - model.agentcount + 1
    positions_initial_infectives = StatsBase.sample(lowest_home:highest_home, initial_infectives)
    for (id, pos) in enumerate(collect(highest_home:-1:lowest_home))
        if !(pos in positions_initial_infectives)
            Agents.add_agent_pos!(Citizen(id, pos, :S), model)
            push!(model.homes, pos)
        else
            Agents.add_agent_pos!(Citizen(id, pos, :E), model)
            push!(model.homes, pos)
        end
    end
    return model
end

function initialize_statecounts!(model::Agents.ABM)
    model.S_count = sum([model.agents[i].state == :S for i in 1:model.agentcount])
    model.E_count = sum([model.agents[i].state == :E for i in 1:model.agentcount])
    model.I_count = deepcopy(model.E_count)
    model.I_cumulative = deepcopy(model.I_count)
    model.new_cases_real = deepcopy(model.I_count)
    return model
end


text:
module TestErrorTrace

using Test
using Try
using Try.Internal: ErrorTrace

function module_context(f)
    old = Try.enable_errortrace(true)
    try
        Base.invokelatest(f)
    finally
        Try.enable_errortrace(old)
    end
end

@noinline f1(x) = x ? Ok(nothing) : Err(ErrorException("nope"))
@noinline f2(x) = f1(x)
@noinline f3(x) = f2(x)

function test_errtrace()
    err = f1(false)
    @test err isa Err
    errmsg = sprint(show, "text/plain", err)
    @test occursin("f1(x::Bool)", errmsg)

    exc = try
        Try.unwrap(err)
        nothing
    catch x
        x
    end
    @test exc isa ErrorTrace
    excmsg = sprint(showerror, exc)
    @test occursin("f1(x::Bool)", excmsg)
end

end  # module


text:
"""
"""
function population_to_spectrum(u, m, x, p)
    A, ω, Γ = p

end


text:
@testset "For.jl" begin
    filename = "data\\test.for"
    df = DDR2import.For.read(filename)
    @test df[1,:RANK] == "Base"
    @test df[2,:Dep] == "ALANDISLANDS"
    @test df[3,:Des] == "ESTONIA"
    @test df[4,Symbol("2017")] == 2.536986301
    @test df[5,Symbol("2024")] == 0.013151227
end